水文分析與計算——頻率分析及洪水流量過程

[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