24小时热门版块排行榜

返回列表

艾雨0802

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 16
帖子: 4
在线: 9.2小时
虫号: 1062419
注册: 2010-07-22
性别: MM
专业: 遥感机理与方法

[求助] 用1stopt迭代求参数，求高手帮我看看代码哪里有问题，谢谢

我是菜鸟一个，昨天才接触1Stopt，更加不会pascal语言。想迭代求解参数，使方程值最小。但是，点了“执行”后确没有啥反应。我用的1.5pro版本。跪求高手指点！谢谢

algorithm=5par_retrieval ;
constant theta(1:9)=[30,32.5,35,37.5,40,42.5,45,47.5,50] ;
constant TBh_obs(1:9)=[241.17748,240.65136,239.83198,238.66749,237.10444,235.08808,232.56285,229.47270,225.76159] ;
constant p=0.01744,mus=0.1;
parameter eb[3,30],tau[0,1.0],omega[0,0.5],rms[0.0025,0.03],cl[0.05,0.5];
minimum;
startprogram [pascal];
procedure mainmodel;
var i:integer;
tt:array[1..9] of double;
Ah:array[1..9] of double;
Bh:array[1..9] of double;
Ch:array[1..9] of double;
rh:array[1..9] of double;
eh:array[1..9] of double;
E_eff:array[1..9] of double;
cf,temp:double;
Begin
cf:=0.0
tt[1]:=exp(-tau/COS(theta[1]*p)) ;
rh[1]:=ABS((mus*COS(theta[1]*p)-SQRT(mus*eb-SIN(theta[1]*p)*SIN(theta[1]*p)))/(mus*COS(theta[1]*p)+SQRT(mus*eb-SIN(theta[1]*p)*SIN(theta[1]*p))))^2;
Ah[1]:=0.068502*(theta[1]*p)^2 -0.058486*(theta[1]*p) +0.976321;
Bh[1]:=-0.051377*(theta[1]*p)^2 +0.014978*(theta[1]*p) +0.045456 ;
Ch[1]:=0.601618*(theta[1]*p)^2 -0.151848*(theta[1]*p) -0.607679;
eh[1]:=1-rh[1]*Ah[1]*exp(Bh[1]*(rms^2/cl)^2+Ch[1]*(rms^2/cl)) ;
E_eff[1]:=eh[1]*tt[1]+(1-tt[1])*(1-omega)+(1-eh[1])*(1-omega)*(1-tt[1])*tt[1];
for i:=2 to 9 do begin
      tt:=exp(-tau/COS(theta*p)) ;
      rh:=ABS((mus*COS(theta*p)-SQRT(mus*eb-SIN(theta*p)*SIN(theta*p)))/(mus*COS(theta*p)+SQRT(mus*eb-SIN(theta*p)*SIN(theta*p))))^2;
      Ah:=0.068502*(theta*p)^2 -0.058486*(theta*p) +0.976321;
      Bh:=-0.051377*(theta*p)^2 +0.014978*(theta*p) +0.045456 ;
      Ch:=0.601618*(theta*p)^2 -0.151848*(theta*p) -0.607679;
      eh:=1-rh*Ah*exp(Bh*(rms^2/cl)^2+Ch*(rms^2/cl)) ;
      E_eff:=eh*tt+(1-tt)*(1-omega)+(1-eh)*(1-omega)*(1-tt)*tt;
      cf:=cf+(E_eff/E_eff[1]-TBh_obs/TBh_obs[1])^2 ;
      temp:=cf;
end ;
Functionresult:=temp  ;
End;
Endprogram;