xaj.txt

新安江模型是分散型性的水文流域模型

      6 :K1:=K6;
      7 :K1:=K7;
      8 :K1:=K8;
      else K1:=K;
    end;
    xmcl(K1,WM,B,IMP,P,EM,WU0+WL0+WD0,R);
    zsf(K1,WM,WUM,WLM,C,P,EM,WU0,WL0,WD0,R,WU0[i+1],WL0[i+1],WD0[i+1],E);
    syhf(K1,IMP,WM,B,KSS,KG,SM,EX,P,EM,WU0+WL0+WD0,R,S0,SDCount,RS,RSS,RG,S1);
    S0 :=S1;
  end; //end if
  S :=S1; WU :=WU0[N+1]; WL :=WL0[N+1]; WD :=WD0[N+1]; FR :=FR;
  //Memo1.Lines.Clear;
  for i :=1 to N do IU :=RS+RSS;//地面径流和壤中流之和
  //-------流域汇流计算
  m :=30;
  if SDCount=8 then Hour :=3 else if SDCount=1 then Hour:=24;
  SetLength(UH,m+1);
  UnitLine(NN,NK,Hour,UH); //计算单位线
  Confluence(IU,UH,QRS);//地面径流汇流计算
  QJL :=25;  //基流
  //地下汇流
  QRG[0] :=QRG0;
  for i :=1 to N do QRG :=QRG[i-1]*KKG+RG*(1-KKG)*U0;
  for i :=1 to N do Q :=QRG+QRS*U0+QJL;
  QRG0 :=Q[High(Q)]; //最后时刻的流量
end;

procedure Tfrm_XinAnJiangModel.XinAnJiang1(rq: array of TDateTime; K, K4,
  K5, K6, K7, K8, WM, B, IMP, WUM, WLM, C, SM, EX, KSS, KG, KKG, NN,
  NK: Double; N: Integer; P, EM: array of double; SDCount: Integer;
  U0: Double; var QRG0, FR: Double; var S,WU,WL,WD,E,R,RS,RSS,RG,Q: array of double);
var
  Year,Month,Day :Word;
  i,m ,hour:Integer;
  //RS,RSS,RG :array of double;
  S1,K1,QJL :Double;
  IU,QRS,QRG :array of double; //IU=RS+RSS--净雨
  UH :array of double;
begin
  //SetLength(RS,High(P)+1);    SetLength(RSS,High(P)+1);   SetLength(RG,High(P)+1);
  for i :=Low(P)+1 to High(P) do begin
    DecodeDate(rq,Year,Month,Day);
    case Month of
      4 :K1:=K4;
      5 :K1 :=K5;
      6 :K1:=K6;
      7 :K1:=K7;
      8 :K1:=K8;
      else K1:=K;
    end;
    xmcl(K1,WM,B,IMP,P,EM,WU[i-1]+WL[i-1]+WD[i-1],R);
    zsf(K1,WM,WUM,WLM,C,P,EM,WU[i-1],WL[i-1],WD[i-1],R,WU,WL,WD,E);
    syhf(K1,IMP,WM,B,KSS,KG,SM,EX,P,EM,WU[i-1]+WL[i-1]+WD[i-1],R,S[i-1],SDCount,RS,RSS,RG,S1);
    S :=S1;
  end; //end if

  SetLength(IU,N+1);   SetLength(QRS,N+1);    SetLength(QRG,N+1);  
  for i :=1 to N do IU :=RS+RSS;//地面径流和壤中流之和
  //-------流域汇流计算
  m :=30;
  if SDCount=8 then Hour :=3 else if SDCount=1 then Hour:=24;
  SetLength(UH,m+1);
  UnitLine(NN,NK,Hour,UH); //计算单位线
  Confluence(IU,UH,QRS);//地面径流汇流计算
  QJL :=25;
  //地下汇流
  QRG[0] :=QRG0;
  for i :=1 to N do QRG :=QRG[i-1]*KKG+RG*(1-KKG)*U0;
  for i :=1 to N do Q :=QRG+QRS*U0+QJL;
  QRG0 :=Q[High(Q)]; //最后时刻的流量
end;

//-------------------单位线为数组形式
procedure Tfrm_XinAnJiangModel.XinAnJiang2(rq: array of TDateTime; K, K4,
  K5, K6, K7, K8, WM, B, IMP, WUM, WLM, C, SM, EX, KSS, KG, KKG: Double;
  UH: array of double; N: Integer; P, EM: array of double;
  SDCount: Integer; U0,QJL: Double; var QRG0, FR, S, WU, WL, WD: Double;
  var Q: array of double);
var
  Year,Month,Day :Word;
  i:Integer;
  R,RS,RSS,RG :array of double;
  WU0,WL0,WD0,E :array of double;
  IU,QRS,QRG :array of double; //IU=RS+RSS--净雨
  Hour,S0,S1,K1 :Double; //QJL-基流
begin
  //----计算各时段的RS、RSS、RG
  SetLength(R,N+1);     SetLength(RS,N+1);    SetLength(RSS,N+1);   SetLength(RG,N+1);
  SetLength(WU0,N+2);   SetLength(WL0,N+2);   SetLength(WD0,N+2);   SetLength(E,N+1);
  SetLength(IU,N+1);    SetLength(QRS,N+1);   SetLength(QRG,N+1);
  S0:=S; WU0[1]:=WU; WL0[1]:=WL; WD0[1]:=WD;
  for i :=1 to N do begin
    DecodeDate(rq,Year,Month,Day);
    case Month of
      4 :K1:=K4;
      5 :K1 :=K5;
      6 :K1:=K6;
      7 :K1:=K7;
      8 :K1:=K8;
      else K1:=K;
    end;
    xmcl(K1,WM,B,IMP,P,EM,WU0+WL0+WD0,R);
    zsf(K1,WM,WUM,WLM,C,P,EM,WU0,WL0,WD0,R,WU0[i+1],WL0[i+1],WD0[i+1],E);
    syhf(K1,IMP,WM,B,KSS,KG,SM,EX,P,EM,WU0+WL0+WD0,R,S0,SDCount,RS,RSS,RG,S1);
    S0 :=S1;
  end; //end if
  S :=S1; WU :=WU0[N+1]; WL :=WL0[N+1]; WD :=WD0[N+1]; FR :=FR;
  for i :=1 to N do IU :=RS+RSS;//地面径流和壤中流之和
  //-------流域汇流计算
  Confluence(IU,UH,QRS);//地面径流汇流计算

  //地下汇流
  QRG[0] :=QRG0;
  for i :=1 to N do QRG :=QRG[i-1]*KKG+RG*(1-KKG)*U0;
  for i :=1 to N do Q :=QRG+QRS*U0+QJL;
  QRG0 :=Q[High(Q)]; //最后时刻的流量
end;

procedure Tfrm_XinAnJiangModel.Confluence(I, U: array of Double;
  var Q: array of double);
var
  j,k :Integer;
begin
  if High(Q)>=High(U) then begin //当净雨时段数大于单位线数值的个数
    for j :=1 to High(Q) do begin
      for k :=1 to j do begin
        if k>High(U)+1 then break;
        Q[j] :=I[j-k+1]*U[k]+Q[j];
        //Form1.Memo5.Lines.Add(FormatFloat('0.00',I[j-k+1]));
        //Form1.Memo6.Lines.Add(FormatFloat('0.00',U[k]));
      end;
    end;//end for j
  end
  else begin  //当净雨时段数小于单位线数值的个数
    for j :=1 to High(Q) do begin
      for k :=1 to j do begin
        Q[j] :=I[k]*U[j-k+1]+Q[j];
        if k>High(I) then break;
      end;
    end;
  end; //end if
end;

procedure Tfrm_XinAnJiangModel.UnitLine(n, K, Hour: Double;
  var U: array of double);
var
  i :Integer;
  Sum1 :Double;
begin
  Sum1 :=0;
  U[1] :=0;  //单位线第一个设为0
  for i :=1 to High(U)-1 do begin
    if i<=High(U)-2 then begin
      U[i+1] :=GammaP(n,i*Hour/K)-GammaP(n,(i-1)*Hour/K);
      Sum1 :=Sum1+U[i+1];
    end
    else U[i+1] :=1-Sum1;
  end;//end for
  for i :=1 to High(U)-1 do begin
    U :=U[i+1];
    //Form1.Memo1.Lines.Add('U['+inttostr(i)+']='+FormatFloat('0.00',U));
  end;
  //U[1]:=0.1; U[2]:=0.8; U[3]:=0.1;
  //U[1]:=0.17; U[2]:=0.4; U[3]:=0.2; U[4]:=0.12; U[5]:=0.06; U[6]:=0.03;   U[7]:=0.02;
end;

function Tfrm_XinAnJiangModel.Gamma(XX: Double): Double;
begin
  Result :=Exp(GammaLn(XX));
end;

procedure Tfrm_XinAnJiangModel.GammaConFrac(var GammaCF: Double; A,
  X: Double; var GamLN: Double);
const
  IMax =10; EPS=3.0E-7;
var
  N :Integer;
  Gold,G,Fac,B0,B1,A0,A1,ANF,ANA,AN :Double;
begin
  GamLN :=GammaLn(A);
  if X<0 then begin
    MessageDlg('X<0，退出程序，请核实！',mtError,[mbOK],0);
    GammaCF :=Null;
    exit;
  end
  else if X=0 then GammaCF :=0
  else begin
    Gold :=0; Fac :=1;
    A0 :=1;  A1 :=X;
    B0 :=0;  B1 :=1;
    for N :=1 to IMax do begin
      AN :=N;
      ANA :=AN-A;
      A0 :=(A1+A0*ANA)*Fac;
      B0 :=(B1+B0*ANA)*Fac;
      ANF :=AN*Fac;
      A1 :=X*A0+ANF*A1;
      B1 :=X*B0+ANF*B1;
      if A1<>0 then begin
        Fac :=1.0/A1;
        G :=B1*Fac;
        if ABS((G-Gold)/G)<EPS then begin
          GammaCF :=G*Exp(-X+A*Ln(X)-GamLN);
          break;
        end;
        Gold :=G;
      end;//end if
    end;//end for
    if N>IMax then begin
      GammaCF :=Null;
      Showmessage('A too large,IMax too small');
    end;
  end;
end;

function Tfrm_XinAnJiangModel.GammaLn(XX: Double): Double;
const
  Half=0.5; One=1.0; FPF=5.5;
var
  j :Integer;
  Stp,X,Tmp,Ser :Double;//X为XX-1
  Cof :array[1..6] of Double;//Lanczos逼近公式的常系数Ci
begin
  if XX<=0 then begin
    MessageDlg('珈玛函数的α值不能小于等于0，请核实！',mtError,[mbOK],0);
    GammaLn :=Null;
    exit;
  end;
  if XX=1 then Result :=0
  else begin
    Stp :=Power(2*Pi,0.5);
    Cof[1]:=76.18009173;  Cof[2]:=-86.50532033;
    Cof[3]:=24.01409822;  Cof[4]:=-1.231739516;
    Cof[5]:=1.20858003E-3;Cof[6]:=-5.6382E-6;
    X :=XX-One;   //Z=X
    Tmp :=X+FPF;  //Z+5.5;
    Tmp :=(X+Half)*Ln(Tmp)-Tmp;
    Ser :=One;
    For j :=1 To 6 do
    begin
      X :=X+One;
      Ser :=Ser+Cof[j]/X;
    end;
    Result :=Tmp+Ln(Stp*Ser);
  end;//end if
end;

function Tfrm_XinAnJiangModel.GammaP(A, X: Double): Double;
var
  GammaCF,GamLn :double;
begin
  if (X<0) or (A<=0) then
    MessageDlg('A或X<0，退出程序，请核实！',mtError,[mbOK],0);
  if X<A+1 then begin
    GammaSeries(GammaCF,A,X,GamLn);
    GammaP :=GammaCF;
  end else begin
    GammaConFrac(GammaCF,A,X,GamLN);
    GammaP :=1-GammaCF;
  end;
end;

function Tfrm_XinAnJiangModel.GammaQ(A, X: Double): Double;
var
  GammaCF,GamLn :Double;
begin
  if (X<0) or (A<=0) then
    MessageDlg('A或X<0，退出程序，请核实！',mtError,[mbOK],0);
  if X<A+1 then begin
    GammaSeries(GammaCF,A,X,GamLn);
    GammaQ :=1-GammaCF;
  end else begin
    GammaConFrac(GammaCF,A,X,GamLN);
    GammaQ :=GammaCF;
  end;
end;

procedure Tfrm_XinAnJiangModel.GammaSeries(var GammaSer: Double; A,
  X: Double; var GamLN: Double);
const
  IMax=100; EPS=0.3E-6;
var
  N :Integer;
  Sum,Del,AP :Double;
begin
  GamLN :=GammaLn(A);
  if X<0 then begin
    MessageDlg('X<0，退出程序，请核实！',mtError,[mbOK],0);
    GammaSer :=Null;
    exit;
  end
  else if X=0 then GammaSer :=0
  else begin
    AP :=A;
    Sum :=1.0/A;
    Del :=Sum;
    for N :=1 To IMax do begin
      AP :=AP+1;
      Del :=Del*X/AP;
      Sum :=Sum+Del;
      if abs(Del)<abs(Sum)*EPS then begin
        GammaSer :=Sum*Exp(-X+A*Ln(X)-GamLn);
        break;
      end;
    end;//end for
  end;//end if
end;