/*
  M--庫容曲線插值節(jié)點數(shù)
  M--泄流曲線節(jié)點數(shù)
  K1--計算時段總數(shù)
  Z0--起調(diào)水位
  S[0]--起調(diào)水位相應(yīng)庫容
  S[K]--時段K的水庫庫容
  O[K]--時段K下泄流量
  Q[K]--對段K入庫流量（設(shè)計洪水過程）
  Z，V--庫容曲線的水位，庫容數(shù)組
  V，q--泄流曲線的庫容，流量數(shù)組
  QE--防洪控制點的安全泄流量。
  Y（K）--閘門全開時的泄流能力
  VF--防洪庫容
  ZF--防洪高水位
*/
#include<fstream.h>
#include<math.h>
#define K1 45
#define M  25
#define N  35
void main()
{
	int i,j,k,l;
	double Z0,V0,QE,DV,SMAX,ZF,VF,A,C,D;
	double  S[K1],O[K1],Q[K1],Z[K1],Y[K1];
	double Z1[M],V[M],q[N],Z2[N];
	ifstream ZV("D:\\水位庫容關(guān)系.txt");
	for(i=0;i<M;i++)
		ZV>>Z1[i]>>V[i];
	ZV.close;
	ifstream Zq("D:\\泄流能力曲線.txt");
	for(i=0;i<N;i++)
		Zq>>Z2[i]>>q[i];
	Zq.close;
	ifstream Qin("D:\\防洪標(biāo)準(zhǔn)洪水.txt");
	for(i=0;i<K1;i++)
		Qin>>Q[i];
	Qin.close;
	Z0=98.0;
	QE=18000;
	for(i=0;i<M-1;i++)
		if(Z0>=Z1[i]&&Z0<Z1[i+1])
			V0=V[i]+(Z0-Z1[i])*(V[i+1]-V[i])/(Z1[i+1]-Z1[i]);
	for(j=0;j<N-1;j++)
	    if(Z0>=Z2[j]&&Z0<Z2[j+1])
			Y[0]=q[j]+(Z0-Z2[j])*(q[j+1]-q[j])/(Z2[j+1]-Z2[j]);
    S[0]=V0;
	Z[0]=Z0;
	O[0]=Q[0];
		if(Q[0]>Y[0])
			O[0]=Y[0];
		else
			O[0]=Q[0];
		for(k=0;k<K1-1;k++)
		{
			if(Q[k+1]<=QE)
			{
				if(Q[k+1]<=Y[k])
				{
					O[k+1]=Q[k+1];
					Y[k+1]=Y[k];
					S[k+1]=S[k];
				}
			    else
				{
					    
						O[k+1]=O[k];
			            l=0;
						while(l<200)
						{
							S[k+1]=(Q[k]+Q[k+1]-O[k]-O[k+1])*6*3600/200000000+S[k];
							C=S[k+1];
							{
								for(i=0;i<M-1;i++)
									if(C>=V[i]&&C<V[i+1])
							           A=Z1[i]+(Z1[i+1]-Z1[i])*(C-V[i])/(V[i+1]-V[i]);
								if(C>V[M-1])
						            A=Z1[M-1]+(C-V[M-1])*(Z1[M-1]-Z1[M-2])/(V[M-1]-V[M-2]);
					            for(j=0;j<N-1;j++)
						            if(A>=Z2[j]&&A<Z2[j+1])
					                   D=q[j]+(q[j+1]-q[j])*(A-Z2[j])/(Z2[j+1]-Z2[j]);
									if(A>Z2[N-1])   
									   D=q[N-1];
							}
							if(fabs(O[k+1]-D)>0.0001)
							{
								O[k+1]=(O[k+1]+D)/2;
				            	l++;
							}
				            else
							{
				               	O[k+1]=D;
					            break;
							}					
						}	
						Y[k+1]=D;
				}
			}
			if(Q[k+1]>QE) 
			{
				
				DV=(Q[k]+Q[k+1]-O[k]-QE)*3600*6/200000000;
				C=S[k]+DV;
				for(i=0;i<M-1;i++)
					if(C>=V[i]&&C<V[i+1])
						A=Z1[i]+(C-V[i])*(Z1[i+1]-Z1[i])/(V[i+1]-V[i]);
                    if(C>V[M-1])
					    A=Z1[M-1]+(C-V[M-1])*(Z1[M-1]-Z1[M-2])/(V[M-1]-V[M-2]);
				for(j=0;j<N-1;j++)
					if(A>=Z2[j]&&A<Z2[j+1])
						D=q[j]+(A-Z2[j])*(q[j+1]-q[j])/(Z2[j+1]-Z2[j]);
					if(A>Z2[N-1])   
						D=q[N-1];
				if(D>QE)
				{
					O[k+1]=QE;
					Y[k+1]=D;
					S[k+1]=C;
				}
				else
				{
					    double O1;
						O[k+1]=O[k];
			            l=0;
						while(l<200)
						{
							S[k+1]=(Q[k]+Q[k+1]-O[k]-O[k+1])*6*3600/200000000+S[k];
							C=S[k+1];
							{
								for(i=0;i<M-1;i++)
									if(C>=V[i]&&C<V[i+1])
							           A=Z1[i]+(Z1[i+1]-Z1[i])*(C-V[i])/(V[i+1]-V[i]);
								if(C>V[M-1])
						            A=Z1[M-1]+(C-V[M-1])*(Z1[M-1]-Z1[M-2])/(V[M-1]-V[M-2]);
					            for(j=0;j<N-1;j++)
						            if(A>=Z2[j]&&A<Z2[j+1])
					                   D=q[j]+(q[j+1]-q[j])*(A-Z2[j])/(Z2[j+1]-Z2[j]);
								if(A>Z2[N-1])   
									D=q[N-1];
							    if(D>QE) 
									O1=QE;
								else 
									O1=D;
							}
							if(fabs(O[k+1]-O1)>0.0001)
							{
								O[k+1]=(O[k+1]+O1)/2;
					            l++;								
							}
				            else
					           	break;
						}
						Y[k+1]=D;					    
				}
			}
		}
		SMAX=0.0;
	    for(k=0;k<K1;k++)
		{
			if(SMAX<S[k])
			{
				SMAX=S[k];
				A=k;
			}
		}
		for(k=A;k<K1-1;k++)
		{
			if(S[k]>V0)
			{
		       double O1;
	           O[k+1]=O[k];
               l=0;
		       while(l<200)
			   {
			         S[k+1]=(Q[k]+Q[k+1]-O[k]-O[k+1])*6*3600/200000000+S[k];
			         C=S[k+1];
					 {
				        for(i=0;i<M-1;i++)
				        if(C>=V[i]&&C<V[i+1])
		                   A=Z1[i]+(Z1[i+1]-Z1[i])*(C-V[i])/(V[i+1]-V[i]);
				        if(C>V[M-1])
				           A=Z1[M-1]+(C-V[M-1])*(Z1[M-1]-Z1[M-2])/(V[M-1]-V[M-2]);
				        for(j=0;j<N-1;j++)
				        if(A>=Z2[j]&&A<Z2[j+1])
				           D=q[j]+(q[j+1]-q[j])*(A-Z2[j])/(Z2[j+1]-Z2[j]);
				        if(A>Z2[N-1])   
				           D=q[N-1];
				        if(D>QE) 
				           O1=QE;
				        else 
				           O1=D;
					 }
			         if(fabs(O[k+1]-O1)>0.0001)
					 {
				        O[k+1]=(O[k+1]+O1)/2;
			        	l++;								
					 }
			         else
					 {
				        O[k+1]=(Q[k]+O1)*0.4;
			            break;
					 }
			   }
		       Y[k+1]=D;
			}
            else
			{
				O[k+1]=Q[k+1];
				S[k+1]=S[k];
				Y[k+1]=Y[k];
			}
		}
        for(k=0;k<K1;k++)
		{
			for(i=0;i<M-1;i++)
				if(S[k]>=V[i]&&S[k]<V[i+1])
					Z[k]=Z1[i]+(Z1[i+1]-Z1[i])*(S[k]-V[i])/(V[i+1]-V[i]);
			if(S[k]>V[M-1])
				Z[k]=Z1[M-1]+(Z1[M-1]-Z1[M-2])*(S[k]-V[M-1])/(V[M-1]-V[M-2]);
		}
        for(i=0;i<M-1;i++)
		if(SMAX>=V[i]&&SMAX<V[i+1])
			ZF=Z1[i]+(Z1[i+1]-Z1[i])*(SMAX-V[i])/(V[i+1]-V[i]);
		VF=SMAX-S[0];
        ofstream file1("D:\\防洪標(biāo)準(zhǔn)調(diào)洪結(jié)果.txt");
		file1<<"來水過程"<<"\t"<<"下泄能力"<<"\t"<<"下泄過程"<<"\t"<<"水庫蓄水量"<<"\t"<<"水庫水位"<<endl;
			for(k=0;k<K1;k++)
				file1<<Q[k]<<"\t"<<Y[k]<<"\t"<<O[k]<<"\t"<<S[k]<<"\t"<<Z[k]<<endl;
        file1<<"防洪高水位水位"<<ZF<<endl;
		file1<<"防洪庫容"<<VF<<endl;
			file1.close();
}