水文分析与计算——频率分析及洪水流量过程

[cpp] 
//频率分析.h 
void PinLvFenXi() 
{ 
    using namespace std; 
    const int NYear = EndYear - 1862 + 1,//调查考证期(1862-2000) 
        L = 1,//实测期特大洪水个数 
        a = 2,//考证期及实测期特大洪水个数 
        YearSuperYearW[a] = {1903, 1975},//特大洪水年份 
        SuperYearW[a] = {21000, 12300};//特大洪水值 
    int order;//降序的序号 
    double AverageMaxYearW[J] = {0}, 
        HundredYearW[J] = {0},//不同年极值的百年一遇设计值 
        K,//模比系数 
        Cv = 0;//变差系数近似无偏估计值 
        ofstream outfile; 
        outfile.open("outfile_YearW_P.txt"); 
 
 
    outfile<<setw(10)<<"年份"<<setw(10)<<"最大1日" 
            <<setw(10)<<"降序序号"<<setw(10)<<"经验频率"<<endl; 
    outfile<<"调查考证期："<<endl; 
    for(int M = 0; M < a; M++)//最大1日特大洪水经验频率 
        outfile<<setw(10)<<YearSuperYearW[M]<<setw(10)<<SuperYearW[M] 
            <<setw(10)<<M+1<<setw(10)<<(double)(M+1)/(NYear+1)<<endl; 
    outfile<<"实测期："<<endl; 
    for(int i = 0; i < Y; i++) 
    {//最大1日实测洪水经验频率 
        for(int M = 0; M < a; M++) 
            if(StartYear+i == YearSuperYearW[M]) 
            {                
                outfile<<setw(10)<<YearSuperYearW[M]<<setw(10)<<SuperYearW[M] 
                    <<setw(10)<<"已计入特大洪水"<<endl; 
                i++; 
                M = 0;//再检验i++ 
            } 
        if(i == Y) break; 
        AverageMaxYearW[0] += MaxYearW[i][0]; 
        outfile<<setw(10)<<StartYear+i<<setw(10)<<MaxYearW[i][0]; 
        order = 1;//从1开始降序排序的(忽略L个特大洪水) 
        for(int j = 0; j < Y; j++) 
            if(MaxYearW[i][0]<MaxYearW[j][0]) order++; 
        outfile<<setw(10)<<order - L<<setw(10)<<(double)(order - L)/(Y+1)<<endl; 
    } 
    AverageMaxYearW[0] *= ((NYear - a)/(Y - L)); 
    for(int M = 0; M < a; M++) 
        AverageMaxYearW[0] += SuperYearW[M]; 
    AverageMaxYearW[0] /= NYear;//得到最大1日洪水的多年平均值 
    for(int i = 0; i < Y; i++) 
    {//计算最大1日洪水的变差系数 
        for(int M = 0; M < a; M++) 
            if(StartYear+i == YearSuperYearW[M]) 
            {                
                i++; 
                M = 0;//再检验i++ 
            } 
        if(i == Y) break; 
        Cv += pow(MaxYearW[i][0] - AverageMaxYearW[0], 2); 
    } 
    Cv *= ((NYear - a)/(Y - L)); 
    for(int M = 0; M < a; M++) 
        Cv += pow(SuperYearW[M] - AverageMaxYearW[0], 2); 
    Cv = pow(Cv/(NYear - 1), 0.5)/AverageMaxYearW[0]; 
 
 
    outfile<<endl<<endl<<setw(10)<<"年份"; 
    for(int k = 1; k < J; k++)//其他年极值洪水经验频率 
        outfile<<setw(7)<<"最大"<<NumJ[k]<<"日" 
            <<setw(10)<<"降序序号"<<setw(10)<<"经验频率"; 
    outfile<<endl; 
    for(int i = 0; i < Y; i++) 
    {//算经验频率 
        outfile<<setw(10)<<StartYear+i; 
        for(int k = 1; k < J; k++) 
        { 
            outfile<<setw(10)<<MaxYearW[i][k]; 
            order = 1;//从1开始降序排序的 
            for(int j = 0; j < Y; j++) 
                if(MaxYearW[i][k]<MaxYearW[j][k]) order++; 
            outfile<<setw(10)<<order<<setw(10)<<(double)order/(Y+1); 
        } 
        outfile<<endl; 
    } 
    outfile.close(); 
 
    for(int k = 1; k < J; k++) 
    {//算其他年极值的多年平均值 
        for(int i = 0; i < Y; i++) 
            AverageMaxYearW[k] += MaxYearW[i][k]; 
        AverageMaxYearW[k] /= Y; 
    } 
    cout<<"年极值系列频率分析——理论频率曲线参数初值如下："<<endl; 
    for(int k = 0; k < J; k++) 
    {//计算其他年极值的变差系数 
        if(k != 0) 
        { 
            Cv = 0; 
            for(int i = 0; i < Y; i++) 
            { 
                K = MaxYearW[i][k]/AverageMaxYearW[k]; 
                Cv += pow(K - 1, 2); 
            } 
            Cv = pow(Cv/(Y - 1), 0.5); 
        } 
        cout<<"最大"<<NumJ[k]<<"日"<<"多年平均值=  "<<AverageMaxYearW[k]<<endl 
            <<"变差系数近似无偏估计值Cv = "<<Cv<<endl<<endl; 
    } 
 
 
 
    //推求设计洪水过程线 
    double ClassicalMonthQ[31],//典型洪水月流量 
        ForcastQ[31],//设计洪水过程 
        MaxClassicalYearW[J],//典型洪水的极值洪量 
        RatioW[J];//同频率倍比 
    int date[J];//记录典型年极值时段，判断是否为包含关系 
    ifstream infile; 
    infile.open("infile_ClassicalWprocess.txt"); 
    for(int i = 0; i < 31; i++) 
        infile>>ClassicalMonthQ[i]; 
    infile.close(); 
    outfile.open("outfile_ForecastWprocess.txt"); 
 
    for(int k = 0; k<J; k++) 
    {//求典型洪水的极值洪量 
         MaxClassicalYearW[k] = 0;//初始化 
         for(int j = 0; j < 31 - NumJ[k] + 1; j++) 
         { 
             temp_W = 0;//初始化 
             for(int n = 0; n < NumJ[k]; n++) 
                 temp_W += ClassicalMonthQ[j + n]; 
             if(temp_W >MaxClassicalYearW[k]) 
             { 
                MaxClassicalYearW[k] = temp_W;//(m3/s*月) 
                date[k] = j; 
             } 
         } 
         cout<<"最大"<<NumJ[k]<<"日的时段：1975年8月" 
             <<date[k]+1<<"日至"<<date[k]+NumJ[k]<<"日"<<endl; 
    } 
    cout<<"由以上时段包含关系判断典型洪水选取是否合适！" 
        <<"合适请输入1，否则请关闭！"<<endl; 
    cin>>order; 
    if(order == 1) 
    { 
        cout<<"输入由理论频率曲线得到的百年一遇设计值："<<endl; 
        for(int k = 0; k < J; k++) 
        {//求同频率放大倍比 
            cout<<"最大"<<NumJ[k]<<"日  " 
                <<"典型洪水："<<MaxClassicalYearW[k] 
                <<"     设计洪水："; 
            cin>>HundredYearW[k]; 
            if( k == 0) RatioW[k] = HundredYearW[k]/MaxClassicalYearW[k]; 
            else RatioW[k] = (HundredYearW[k] - HundredYearW[k - 1])/(MaxClassicalYearW[k] - MaxClassicalYearW[k - 1]); 
        } 
        for(int k = 0; k < J; k++) 
        {//求设计洪水过程 
            for(int i = date[k]; i < date[k]+NumJ[k]; i++) 
            { 
                if(k == 0)  ForcastQ[i] = ClassicalMonthQ[i]*RatioW[k]; 
                else if(i < date[k-1]|| i>=date[k-1]+NumJ[k-1]) 
                    ForcastQ[i] = ClassicalMonthQ[i]*RatioW[k]; 
            } 
        } 
        outfile<<setw(10)<<"典型洪水"<<setw(10)<<"设计洪水"<<endl; 
        for(int i = date[J - 1]; i < date[J - 1] + NumJ[J - 1]; i++) 
            outfile<<setw(10)<<ClassicalMonthQ[i]<<setw(10)<<ForcastQ[i]<<endl; 
        outfile.close(); 
        cout<<"频率分析已结束，请关闭！"<<endl; 
        cin>>Cv;//控制台暂停 
    } 
} 
作者：Superwen_go