/*
�@�@�@‚PŽŸŒ³‰Í�°•Ï“®ŒvŽZ


�@�@�@‰‰�K–â‘è‚P‚X
*/

#include<stdio.h>
#include<math.h>
#include<stdlib.h>


double CALWF(double S,double G,double D);
double CALH(double H[],double Z[],double B[],double Q,double SN,double H0,double G,double DX[],double EPS,double IE[],int NJ);
double CALIWA(double S,double G,double D);
double CALQSU(double S,double G,double D,int NJ,double US[],double TS[],double TSC,double WF,double QSU[]);
double CALC(int NJ,double X[],double B[],double QSU[],double WF,double C[],double Q,double H[],int IDT1,int IDT2);
void POUT(double Z[],double ZZ[],double H[],double TS[],double US[],double QSU[],double C[],int NJ,int TIM,int PTIM[],int ITM);
double CALDZ(int NJ,double RAMDA,int IDT2,double WF,double Z[],double QSU[],double C[]);
/*********************************************************************************************/
double H[60],IE[60];
double WF,USC,TSC;

/*********************************************************************************************/
int main(void){
	double X[60],DX[60],Z[60],B[60],US[60],TS[60],ZZ[60],QSU[60],C[60];
	//double IE[60]/*,HH[60]*/;
	int PTIM[6],OT[6];
	double Q=500,H0=12.0,SN=0.02,D=0.0001,S=1.65;
	double G=9.8,RAMDA=0.4,EPS=0.001;
	int NJ=51,IDT1=50,IDT2=100,ETIM=86400;

	OT[1]=0,OT[2]=1,OT[3]=4,OT[4]=8,OT[5]=12,OT[6]=24;


	FILE *fp;
	if((fp=fopen("out.dat","w"))==NULL){
		printf("�o—̓tƒ@ƒCƒ‹‚ðƒI�[ƒvƒ“‚Å‚«‚Ü‚¹‚ñ");
		exit(1);
	}

	int TIM=0,ITM=1,i,j;
	for(i=2;i<6;i++)
		PTIM[i]=OT[i]*3600;
	for(i=1;i<6;i++)
		printf("PTIM = %6d",PTIM[i]);
/********************************************************************************************/
	for(j=1;j<NJ;j++){
		X[j]=200*(j-1);
		B[j]=100;
		Z[j]=1/700*(j-1)*200;
		ZZ[j]=Z[j];
	}

	for(j=2;j<NJ;j++)
		DX[j]=X[j]-X[j-1];
	
/*********************************************************************************************/

	CALWF(S,G,D);
//----------------------------------------------------------------------------------------------
L1000:
//----------------------------------------------------------------------------------------------	
	CALH(H,Z,B,Q,SN,H0,G,DX,EPS,IE,NJ);
	
	for(j=1;j<NJ;j++){
		TS[j]=H[j]*IE[j]/(S*D);
		US[j]=sqrt(G*H[j]*IE[j]);

		printf("TS =%10.5lf,US =%10.5lf\n",TS[j],US[j]);
	}

	CALIWA(S,G,D);
	printf("OK1\n");
	CALQSU(S,G,D,NJ,US,TS,TSC,WF,QSU);
	printf("OK2\n");
	CALC(NJ,X,B,QSU,WF,C,Q,H,IDT1,IDT2);
	printf("OK3\n");
	POUT(Z,ZZ,H,TS,US,QSU,C,NJ,TIM,PTIM,ITM);
	printf("OK4\n");
	CALDZ(NJ,RAMDA,IDT2,WF,Z,QSU,C);

//
/*********************************************************************************************/
//
	TIM=TIM+IDT2;
	if(TIM<ETIM)goto L1000;

	fclose(fp);
	return 0;
}

/***********************************************************************************************/
/*-----WF‚ÌŒvŽZ-----*/
double CALWF(double S,double G,double D){
	double DD,DDD,A;

	DD=D*100;
	A=36*0.01*0.01/(S*G*100);
	DDD=DD*DD*DD;
	WF=(sqrt(2/3+A/DDD)-sqrt(A/DDD))*sqrt(S*G*100*DD)/100;

	return WF;
}
/*-----H‚ÌŒvŽZ-----*/
double CALH(double H[],double Z[],double B[],double Q,double SN,double H0,double G,double DX[],double EPS,double IE[],int NJ){
	//double H[60],Z[60],B[60],DX[60],IE[60];
	double QQ,SNN,HH,BB,H3,FD,FU,FH,DFDH1,DFDH2,DFDH;
	int j,jj;
	QQ=Q*Q;
	SNN=SN*SN;
	H[1]=H0-Z[1];

	for(j=1;j<NJ;j++){
		jj=j-1;
		if(j==1)goto L20;
		HH=H[jj]*H[jj];
		BB=B[jj]*B[jj];
		H3=pow(H[jj],10/3);
		FD=H[jj]+Z[jj]+QQ/(2*G*BB*HH)+SNN*QQ*DX[j]/(2*BB*H3);
		H[j]=H[jj];
		L30: 
		HH=H[j]*H[j];
		BB=B[j]*B[j];
		H3=pow(H[j],10/3);
		FU=H[j]+Z[j]+QQ/(2*G*BB*HH)-SNN*QQ*DX[j]/(2*BB*H3);
		FH=FU-FD;
		if(fabs(FH)>=EPS){
			DFDH1=QQ/(G*BB*H[j]*H[j]*H[j]);
			DFDH2=(5/3)*SNN*QQ*DX[j]/(BB*pow(H[j],13/3));
			DFDH=1-DFDH1+DFDH2;
			H[j]=H[j]-FH/DFDH;
			goto L30;
		}
	L20: 
	IE[j]=SNN*QQ/(B[j]*B[j]*pow(H[j],10/3));
	}
	
	return H[60];
}

/*-----USC,TSC‚ÌŒvŽZ-----*/
double CALIWA(double S,double G,double D){
/*** IWAGAKI METHOD ***/
double DD,USC2,USC,TSC;
DD=D*100;
if(DD>=0.303)USC2=80.9*DD;
if((0.118<=DD)&&(DD<0.303)) USC2=134.6*pow(DD,31/22);
if((0.0565<=DD)&&(DD<0.118)) USC2=55*DD;
if((0.0065<=DD)&&(DD<0.0565)) USC2=8.41*pow(DD,11/32);
if(DD<0.0065) USC2=226*DD;
USC=sqrt(USC2)/100;
TSC=USC*USC/(S*G*D);
return USC;
}
/*-----QSU‚ÌŒvŽZ-----*/
double CALQSU(double S,double G,double D,int NJ,double US[],double TS[],double TSC,double WF,double QSU[]){
//double TSC;
double A1=0.436184,A2=-0.120168,A3=0.937298;
double AS=0.14,BS=0.143,PK=0.008;
double RW=1,RS=2.65,ET=0.5;
double PI=3.141592;
double SSS=sqrt(S*G*D),AD,ADD,T,D1,D2,OMEGA;
int j;
for(j=1;j<NJ;j++){
	AD=BS/TS[j]-1/ET;
	ADD=fabs(AD);
	T=1/(1+0.33267*sqrt(2)*ADD);
	D1=1/(2*sqrt(PI))*exp(-ADD*ADD);
	D2=1/(2*sqrt(PI))*exp(-ADD*ADD)*(A1*T+A2*T*T+A3*T*T*T);
	if(AD<0.0)D2=1-D2;
	OMEGA=TS[j]/BS*D1/D2+TS[j]/BS/ET-1;
	QSU[j]=PK*(AS*RW/RS*OMEGA*SSS/sqrt(TS[j])-WF);
	if(QSU[j]<0)QSU[j]=0;
	if(US[j]/WF<1.08)QSU[j]=0;
	if(TS[j]<TSC) QSU[j]=0;
}
return QSU[60];
}
/*-----C‚ÌŒvŽZ-----*/
double CALC(int NJ,double X[],double B[],double QSU[],double WF,double C[],double Q,double H[],int IDT1,int IDT2){
	double CA[60],CC[60];
	double DXB;
	int TIM,j;
	CC[NJ]=QSU[NJ]/WF;
	for(j=1;j<NJ;j++){
		C[j]=CC[j];
	}
	for(TIM=IDT1;TIM<IDT2;TIM+=IDT1){
		for(j=1;j<NJ;j++){
			CA[j]=CC[j];
		}	
		for(j=NJ-1;j>1;j--){
			DXB=(X[j+1]-X[j])*B[j];
			CC[j]=CA[j]+IDT1/H[j]*(QSU[j]-WF*CA[j]-Q*(CA[j]-CA[j+1])/DXB);
		}

	}
	return C[60];
};
/*-----PRINT OUT-----*/
void POUT(double Z[],double ZZ[],double H[],double TS[],double US[],double QSU[],double C[],int NJ,int TIM,int PTIM[],int ITM){
double HH[60],DZ;
int j;
FILE *fp;
printf("TIME = %5d,ITM = %5d\n",TIM,ITM);
if(TIM<PTIM[ITM])goto L900;	
ITM=ITM+1;
fprintf(fp,"TIME= %5.0lf HR\n",TIM/3600);
fprintf(fp,"   NO   Z (m)      DZ (m)    HH (m)    H (m)       T*    U* (cm/s) QSU (cm/s)    C\n");
for(j=1;j<NJ;j++){
	HH[j]=H[j]+Z[j];
	DZ=Z[j]-ZZ[j];
	fprintf(fp,"%5d,%10.5lf,%10.5lf,%10.5lf,%10.5lf,%10.5lf,%10.5lf,%10.5lf,%10.5lf\n",j,Z[j],DZ,HH[j],H[j],TS[j],US[j]*100,QSU[j]*100,C[j]);
}
fprintf(fp,"\n");
L900:
return;
}
/*-----DZ‚ÌŒvŽZ-----*/
double CALDZ(int NJ,double RAMDA,int IDT2,double WF,double Z[],double QSU[],double C[]){
	int j;
	double DZ;
	for(j=1;j<NJ;j++){
		DZ=1/(1-RAMDA)*(WF*C[j]-QSU[j])*IDT2;
		Z[j]=Z[j]+DZ;
	}
	return Z[j];
}