cpp1.cpp

水文预报程序

#include<stdio.h>
#include<math.h>
#define NUMBER 485
#define F      686
#define WM     150
#define WUM    30
#define WLM    80
#define WDM    40
#define B      1.5
#define SM     9
#define EX     1.9
#define KG     0.109
#define KKG    0.964
#define KSS    0.156
#define KKSS   0.975
 float S[NUMBER],AU[NUMBER],FR[NUMBER],RS[NUMBER],RSS[NUMBER],RG[NUMBER];
 float SMM,WMM,b;
 struct shu
  {
   int month;
   int day;
   int begin;
   int finish;
   float p1,p2,p3,Q;
  }data[NUMBER];
 int ii[10];
 float E1[12]={0.71,0.68,0.93,1.57,2.29,2.65,3.41,2.76,2.44,1.75,1.21,0.87};
 float E2[12]={0.60,0.31,0.55,0.97,1.49,1.71,2.34,2.45,1.03,0.77,0.60,0.36};
 float  e[12]={1.0,1.0,1.0,1.36,1.09,0.90,1.68,1.47,0.78,1.0,1.0,1.0};
 float EP[NUMBER],P[NUMBER],PE[NUMBER],EU[NUMBER],EL[NUMBER],ED[NUMBER],E[NUMBER],WU[NUMBER],WL[NUMBER],WD[NUMBER],W[NUMBER],a[NUMBER],R[NUMBER];
 void main()

/*产流计算*/
 {
  float FE=0.8;
  float C=0.15;
  int i;
  FILE *fp;
  fp=fopen("Data.txt","r");
  for(i=0;i<NUMBER;i++)
  fscanf(fp,"%d    %d    %d    %d    %f    %f    %f     %f\n",&data[i].month,&data[i].day,&data[i].begin,&data[i].finish,&data[i].p1,&data[i].p2,&data[i].p3,&data[i].Q);
  fclose(fp);
 for(i=0;i<NUMBER;i++)
  {
   if((P[i]=(data[i].p1+data[i].p2+data[i].p3)/3)>=3)
   EP[i]=E2[data[i].month-1]/4.0*e[data[i].month-1];
   else
   EP[i]=E1[data[i].month-1]/4.0*e[data[i].month-1];
  }
  WU[0]=WUM*FE;
  WL[0]=WLM*FE;
  WD[0]=WDM*FE;
  W[0]=WM*FE;
  R[0]=0;
  /*FILE *fp1;
  fp1=fopen("Data1.txt","w");
  for(i=0;i<NUMBER;i++)
  fprintf(fp,"%f\n",EP[i]);*/
 for(i=1;i<NUMBER;i++)
  {
   if(WU[i-1]+P[i]>=EP[i])
     {EU[i]=EP[i];EL[i]=0;ED[i]=0;}
   if(WU[i-1]+P[i]<EP[i]&&WL[i-1]>=C*WLM)
     {EU[i]=WU[i-1]+P[i];
      EL[i]=(EP[i]-EU[i])*WL[i-1]/WLM;
      ED[i]=0;}
   if(WU[i-1]+P[i]<EP[i]&&C*(EP[i]-EU[i])<=WL[i-1]&&WL[i-1]<C*WLM)
     {EU[i]=WU[i-1]+P[i];
      EL[i]=C*(EP[i]-ED[i]);
      ED[i]=0;}
   if(WU[i-1]+P[i]<EP[i]&&WL[i-1]<C*(EP[i]-EU[i]))
     {EU[i]=WU[i-1]+P[i];
      EL[i]=WL[i-1];
      ED[i]=C*(EP[i]-EU[i])-EL[i];}
   E[i]=EU[i]+EL[i]+ED[i];
   PE[i]=P[i]-E[i];
   WMM=WM*(1+B);
   a[i]=WMM*(1-pow((1-W[i-1]/WM),1.0/(1+B)));
   if(PE[i]>=0)
     {
      if(a[i]+PE[i]<=WMM)
      R[i]=PE[i]+W[i-1]-WM+WM*pow((1-(PE[i]+a[i])/WMM),(1+B));
      else if(a[i]+PE[i]>WMM)
      R[i]=PE[i]+W[i-1]-WM;
     }
   else
      R[i]=0;
   WU[i]=WU[i-1]+PE[i]-R[i];
   if(WU[i]>=WUM)
     {WL[i]=WL[i-1]+WU[i]-WUM;WU[i]=WUM;}
   else
      WL[i]=WL[i-1];
   if(WL[i]>=WLM)
     {WD[i]=WD[i-1]+WL[i]-WLM;WL[i]=WLM;}
   else
      WD[i]=WD[i-1];
   if(WD[i]>=WDM)
      WD[i]=WDM;
   else
      WD[i]=WD[i-1];
   if(WU[i]<0)
      WU[i]=0;
   if(WL[i]<0)
      WL[i]=0;
   if(WD[i]<0)
      WD[i]=0;
   W[i]=WU[i]+WL[i]+WD[i];
 }
   FILE *fp1;
   fp1=fopen("result1.txt","w");
   fprintf(fp1,"          产流计算表          \n");
   fprintf(fp1," P  EP   EU   EL   ED   E   PE   WU   WL   WD   W   R   a\n");
   for(i=0;i<NUMBER;i++)
   fprintf(fp1,"%5.2f  %5.2f  %5.2f  %5.2f  %5.2f  %5.2f  %5.2f  %5.2f  %5.2f  %5.2f  %5.2f  %5.2f  %5.2f\n",
	   P[i],EP[i],EU[i],EL[i],ED[i],E[i],PE[i],WU[i],WL[i],WD[i],W[i],R[i],a[i]);



   SMM=SM*(1+EX);
 S[0]=0.8*SM;
 for(i=0;i<NUMBER;i++)
 {
  AU[i]=SMM*(1-pow((1-S[i]/SM),1.0/(1+EX)));
  if(PE[i]<=0)
    {
     FR[i]=1-pow((1-W[i]/WM),B/(1+B));
     RS[i]=0;
     RSS[i]=S[i]*KSS*FR[i];
     RG[i]=S[i]*KG*FR[i];
     S[i+1]=(1-KSS-KG)*S[i];
    }
  else
    {
     FR[i]=R[i]/PE[i];
     if(PE[i]+AU[i]<SMM)
       {
	b=SM*pow((1-(PE[i]+AU[i])/SMM),(1+EX));
	RS[i]=(PE[i]-SM+S[i]+b)*FR[i];
	RSS[i]=(SM-b)*KSS*FR[i];
	RG[i]=(SM-b)*KG*FR[i];
	S[i+1]=(1-KSS-KG)*(SM-b);
       }
      else
       {
	RS[i]=(PE[i]-SM+S[i])*FR[i];
	RSS[i]=SM*KSS*FR[i];
	RG[i]=SM*KG*FR[i];
	S[i+1]=(1-KSS-KG)*SM;
       }
    }
 }
  FILE *fp2;
  fp2=fopen("result2.txt","w");
  fprintf(fp2,"  RS     RSS     RG\n");
  fprintf(fp2,"                   \n");
  for(i=0;i<NUMBER;i++)
 {
  fprintf(fp2,"%5.2f     %5.2f     %5.2f\n",RS[i],RSS[i],RG[i]);
 }


  int j;
 int D;
 D=24/6;
 float UH[14]={0,12.4,102.4,51.9,39.7,33.5,24.4,18.8,13.7,9.4,6.3,3.4,1.7,0};
 float Q1[498],Q2[498],Q3[498],Qc[498];
 for(i=0;i<498;i++)
 Q1[i]=0;
 for(i=0;i<NUMBER;i++)
 for(j=0;j<14;j++)
 Q1[i+j]=Q1[i+j]+UH[j]*RS[i]/10;
 Q2[0]=KKG/(KKG+KKSS)*11.6;
 Q3[0]=KKSS/(KKG+KKSS)*11.6;
 for(i=1;i<498;i++)
    {
     Q2[i]=Q2[i-1]*pow(KKG,1.0/D)+RG[i-1]*(1-pow(KKG,1.0/D))*F/3.6/6;
     Q3[i]=Q3[i-1]*pow(KKSS,1.0/D)+RSS[i-1]*(1-pow(KKSS,1.0/D))*F/3.6/6;
    }
 FILE *fp3;
 fp3=fopen("result3.txt","w");
 for(i=0;i<498;i++)
    {
     Qc[i]=Q1[i]+Q2[i]+Q3[i];
     fprintf(fp3,"%.2f\n",Qc[i]);
    }

j=0;
 for(i=0;i<NUMBER;i++)
   {
    if(data[i].month==4 && data[i].day==11 && data[i].begin==2)
    ii[j++]=i;
    if(data[i].month==4 && data[i].day==18 && data[i].begin==20 && data[i].finish==2)
    ii[j++]=i;
    if(data[i].month==4 && data[i].day==19 && data[i].begin==2 && data[i].finish==8)
    ii[j++]=i;
    if(data[i].month==4 && data[i].day==27 && data[i].begin==14 && data[i].finish==20 )
    ii[j++]=i;
    if(data[i].month==4 && data[i].day==28 && data[i].begin==2 && data[i].finish==8)
    ii[j++]=i;
    if(data[i].month==5 && data[i].day==7 && data[i].begin==8 && data[i].finish==14)
    ii[j++]=i;
    if(data[i].month==5 && data[i].day==7 && data[i].begin==14 && data[i].finish==20)
    ii[j++]=i;
    if(data[i].month==5 && data[i].day==25 && data[i].begin==14 && data[i].finish==20)
    ii[j++]=i;
    if(data[i].month==6 && data[i].day==28 && data[i].begin==8 && data[i].finish==14)
    ii[j++]=i;
    if(data[i].month==7 && data[i].day==9  && data[i].begin==2 && data[i].finish==8)
    ii[j++]=i;
   }
 FILE *fp4;
 fp4=fopen("result4.txt","w");
 for(i=0;i<5;i++)
 fprintf(fp4,"第%d次洪水起止时间（6h一个时段）是%d----%d\n",i+1,ii[2*i]+1,ii[2*i+1]+1);



 int   k,k1[5],k2[5],dt[5];
 int temp1=0,temp2=0;

 float max1[5]={0};
 float max2[5]={0};
 float dq[5];
 
 for(i=0;i<5;i++)
 for(k=ii[2*i];k<ii[2*i+1]+1;k++)
 { 
    if(max1[i]<data[k].Q)
      {
       k1[i]=k;max1[i]=data[k].Q ;  
      }
    if(max2[i]<Qc[k])
      {
       k2[i]=k;max2[i]=Qc[k] ;
      }
 } 
 
 
 

 for(i=0;i<5;i++)
    {
     dt[i]=abs(k1[i]-k2[i]);
     dq[i]=abs(max1[i]-max2[i]);
    }
  FILE *fp5;
  fp5=fopen("result5.txt","w");
  
  for(i=0;i<5;i++)
     {
      fprintf(fp5,"%d,%d\n",k1[i],k2[i]);
      fprintf(fp5,"实测第%d次洪峰流量Q为%.2f,出现时间是(%d month,%d day,%d time)\n",i+1,data[k1[i]].Q,data[k1[i]].month,data[k1[i]].day,data[k1[i]].begin);
      fprintf(fp5,"计算的%d次洪峰流量Q为%.2f,出现时间是(%d month,%d day,%d time)\n",i+1,data[k2[i]].Q,data[k2[i]].month,data[k2[i]].day,data[k2[i]].begin);
      fprintf(fp5,"第%d次实测洪峰流量与计算洪峰流量之差为%.2f\n",i+1,dq[i]);
      fprintf(fp5,"第%d次实测洪峰时间与计算洪峰时间之差为%d\n",i+1,dt[i]);
     }

 float Rc[5],R0[5],dR[5];
 
  for(i=0;i<5;i++)
   {
	 Rc[i]=0,R0[i]=0;
	 for(k=ii[2*i];k<ii[2*i+1]+1;k++)
	 {
	   Rc[i]+=R[k];
	   R0[i]=R0[i]+data[k].Q;
	   }
	   R0[i]=(R0[i]-(data[ii[2*i]].Q+data[ii[2*i+1]].Q-data[ii[2*i]].Q))*6*3.6/F;
	   dR[i]=abs(Rc[i]-R0[i]);}
  FILE *fp6;
  fp6=fopen("result6.txt","w");
    for(i=0;i<5;i++)
	{ fprintf(fp6,"%.2f,%.2f\n",Rc[i],R0[i]);
	   fprintf(fp6,"第%d次实测径流深R0与计算径流深Rc的差值为dR[%d]=%.2f\n",i+1,i+1,dR[i]);}
  
 

    int a1=0,a2=0,a3=0,a4=0;
   float m1,m2,m3,m4;
   float s1=0,s2,s3=0;
   float dy;
   for(i=0;i<5;i++)
      {
       if((dt[i]*1.0)/(k1[i]-ii[i])<0.3)
       a1++;
      }
   m1=a1/5.0;
   for(i=0;i<5;i++)
      {
       if((dq[i]*1.0)/data[k1[i]].Q<0.2)
       a2++;
      }
   m2=a2/5.0;
   for(i=0;i<5;i++)
   {
	if(abs(Rc[i]-R0[i])/(R0[i]*1.0)<0.2&&abs(Rc[i]-R0[i])<=20)
	   a3++;
   }
   m3=a3/5.0;
   for(i=0;i<NUMBER;i++)
   {
    s1+=(Qc[i]-data[i].Q)*(Qc[i]-data[i].Q);
    s2+=data[i].Q/485;
    s3+=(Qc[i]-s2)*(Qc[i]-s2);
   }
   dy=1-s1/s3;
   for(i=0;i<NUMBER;i++)
   {
    if(abs(Qc[i]-data[i].Q)/data[i].Q<0.2)
    a4++;
   }
   m4=a4/485.0;
   FILE *fp7;
   fp7=fopen("result7.txt","w");
   fprintf(fp7,"洪峰出现时间预报合格率为%.2f\n",m1);
   fprintf(fp7,"洪峰流量预报合格率为%.2f\n",m2);
   fprintf(fp7,"计算流量确定性系数DC为%.2f\n",dy);
   fprintf(fp7,"径流深预报合格率为%.2f\n",m3);
   fprintf(fp7,"流量过程预报合格率为%.2f\n",m4);

}