//#pragma once
#include <stdio.h>
// TODO: プログラムに必要な追加ヘッダーをここで参照してください。
//#define _USE_MATH_DEFINES
#include <math.h>
#include <stdlib.h>
#define  q 10.0			//単位幅流量
#define  i0 0.01		//河床勾配
#define  dx 200			//計算区間
#define  n  0.050		//Manningの粗度係数

#define eps 1.0e-5		//許容誤差

#define D  0.05         //礫径(50mm)
#define R  1.65         //土砂の水中比重
#define RAMDA 0.4       //空隙率
#define TSC 0.05        //T*c
#define dt 10           //delta t (sec)
#define G 9.8           //重力加速度 

double HC();
double hn();

double function1(double h);
double function2(double h);
double function3(double h);

double function4(double h);
double function5(double TS);

int main(int argc, char *argv[])
{
	FILE *fp;
	if((fp=fopen("kasyohendo50-6.txt","w"))==NULL){
		printf("kasyohendo50-6.txtをオープンできません。\n");
		exit(1);
	}

/********************************************************/
	double h0,hc;
	h0=hn();	hc=HC();
	printf("等流水深,限界水深\n");
	printf("h0 = %f, hc = %f\n",h0,hc);
/*******************************************************/
	double h1,h2,z,Z[60];
	double UHEN;
//
	double H[60],TS[60],QB[60],DZ[60],Fr[60];
	int i;
	int TIME=0;
	for( i=0 ; i<42  ; i++){ 
		Z[i]=i*i0*dx;
		printf("%6.3f\n",Z[i]);        //初期河床
	}
	
	double ZZ[60];
	FILE *fi;
	if((fi=fopen("kasyohendo50-5.csv","r"))==NULL){
	  printf("ファイルをオープンできません。\n");
	  exit(1);
	}
	for(i=0;i<40;i++)
	  fscanf(fi,"%lf",&ZZ[i]);            //データの読み込み

	fclose(fi);

	for(i=0;i<40;i++){
	  printf("ZZ[%d] = %6.3f\n",i,ZZ[i]);
	  Z[i]=ZZ[i];                            //前の数値を代入
	}
	
	
	/*******************************************************/
L100:
	if(TIME>3600) goto end;          //終了の判定
	printf("TIME = %d sec\n",TIME);


	h1=50.0;									//ダム地点の水位
	H[0]=h1;
	printf("%6.3f\n",h1);
	//fprintf(fp,"縦断距離,河床高,水深,水位,速度水頭\n");
for(int j=0;j<41;j++){			//for先頭
	h2=h1;			//第一近似
Start:
	double Esokudo2=function1(h2);
	//printf("%f\n",Esokudo2);
	double Eloss2=function2(h2);
	//printf("%f\n",Eloss2);
		z=Z[j+1]-Z[j];
		//z=i0*dx;
		//printf("%f\n",z);
//
		double Esokudo1=function1(h1);
		//printf("%f\n",Esokudo1);
		double Eloss1=function2(h1);
		//printf("%f\n",Eloss1);

	UHEN=Esokudo2+h2-Eloss2-(Esokudo1+h1+Eloss1-z);

	if(fabs(UHEN)<eps){
		h1=h2;
		goto Next;
	}
	
	else if(fabs(UHEN)>=eps){
		double func=function3(h1);
		h2-=UHEN/func;
		goto Start;
	}

Next:
	printf("j= %d,z= %6.3f,h= %6.3f,Esokudo= %6.3f\n",j,i0*dx*j,h1,Esokudo1);
	//fprintf(fp," %d, %f, %f, %f, %f\n",(j+1)*dx,i0*dx*(j+1),h1,h1+i0*dx*(j+1),Esokudo1);
	Fr[0]=q/(H[0]*sqrt(G*H[0]));
	Fr[j+1]=q/(h1*sqrt(G*h1));
	printf("%8.5f\n",Fr[j+1]);

	if(Fr[j+1]<=1.0)
		H[j+1]=h1;             //常流 
	else
		H[j+1]=hc;              //射流の場合は限界水深
	TS[j]=function4(H[j]);
	if(TS[j]>TSC)
		QB[j]=function5(TS[j]);
	else
		QB[j]=0;
	printf("%8.5e  ,  %8.5e\n",TS[j],QB[j]);
	//i++;
}												//for終わり
/**************************************************************/
	int k;
	DZ[0]=0,DZ[40]=0;
	for(k=1;k<40;k++){
		DZ[k]=-dt/(1-RAMDA)*(QB[k]-QB[k-1])/dx;
		printf("DZ[%d] = %8.5f\n",k,DZ[k]);
	}
	for(k=0;k<41;k++)
		Z[k]+=DZ[k];

	if(TIME==0 || TIME==980|| TIME==3600 || TIME==43200 || TIME==86400 ||TIME ==172800 || TIME==345600){
		fprintf(fp,"TIME = %d sec\n",TIME);
		fprintf(fp," NO , KP,  Z  ,  H  , DZ(mm) ,  T*  ,  QB  ,Fr\n");
		for(k=0;k<41;k++){
		  fprintf(fp,"%d,  %d,  %6.3f,  %6.3f,  %6.3f,  %6.3f,  %6.3f,  %6.3f\n",k,k*dx,Z[k],H[k],DZ[k]*1000,TS[k],QB[k],Fr[k] );
		}
	}
	TIME+=dt;
	goto L100;


/***********************************************************/
end:
	fprintf(fp,"\n");
	fclose(fp);

	return 0;
}
/**********************************************************/
double HC(){
	return pow((1.1*q*q/G),(double)1/3);	//限界水深
}
double hn(){
	return pow(n*n*q*q/i0,(double)3/10);
}
double function1(double h){
	return q*q/(2*G*h*h);
}
double function2(double h){
	return n*n*q*q*dx/(2*pow(h,(double)10/3));
}
double function3(double h){
	return -q*q/(9.8*h*h*h)+1+10/3*n*n*q*q*dx/2*1/pow(h,13/3);
}
double function4(double h){
	return n*n*q*q/(R*D*pow(h,7/3));
}
double function5(double TS){
	return 8.0*sqrt(R*G*pow(D,3))*pow(TS-TSC,3/2);
}