code 33050; procedure EXPONENTIALLY FITTED TAYLOR(T,TE,M0,M,U,SIGMA,PHI,DIAMETER, DERIVATIVE,I,K,ALFA,NORM,AETA,RETA,ETA,RHO,HMIN,HSTART,OUTPUT); integer M0,M,I,K,NORM; real T,TE,SIGMA,PHI,DIAMETER,ALFA,AETA,RETA,ETA,RHO,HMIN,HSTART; array U; procedure DERIVATIVE,OUTPUT; begin integer KL; real Q,EC0,EC1,EC2,H,HI,H0,H1,H2,BETAN,T2,SIGMAL,PHIL; real array C,RO[M0:M],BETA,BETHA[1:3]; boolean LAST,START; procedure COEFFICIENT; begin real B,B1,B2,BB,E,BETA2,BETA3; B:=H*SIGMAL; B1:=B*COS(PHIL); BB:=B*B; if ABS(B)<"-3 then begin BETA2:=.5-BB/24; BETA3:=1/6+B1/12; BETHA[3]:=.5+B1/3 end else if B1<-40 then begin BETA2:=(-2*B1-4*B1*B1/BB+1)/BB; BETA3:=(1+2*B1/BB)/BB; BETHA[3]:=1/BB end else begin E:=EXP(B1)/BB; B2:=B*SIN(PHIL); BETA2:=(-2*B1-4*B1*B1/BB+1)/BB; BETA3:=(1+2*B1/BB)/BB; if ABS(B2/B)<"-5 then begin BETA2:=BETA2-E*(B1-3); BETA3:=BETA3+E*(B1-2)/B1; BETHA[3]:=1/BB+E*(B1-1) end else begin BETA2:=BETA2-E*SIN(B2-3*PHIL)/B2*B; BETA3:=BETA3+E*SIN(B2-2*PHIL)/B2; BETHA[3]:=1/BB+E*SIN(B2-PHIL)/B2*B; end end; BETA[1]:=BETHA[1]:=1; BETA[2]:=BETA2; BETA[3]:=BETA3; BETHA[2]:=1-BB*BETA3; B:=ABS(B); Q:=if B<1.5 then 4-2*B/3 else if B<6 then (30-2*B)/9 else 2; real procedure NORMFUNCTION(NORM,W); integer NORM; array W; begin integer J; real S,X; S:=0; if NORM=1 then begin for J:=M0 step 1 until M do begin X:=ABS(W[J]); if X>S then S:=X end end else S:=SQRT(VECVEC(M0,M,0,W,W)); NORMFUNCTION:=S; end; procedure LOCAL ERROR BOUND; ETA:=AETA+RETA * NORMFUNCTION(NORM,U); procedure LOCAL ERROR CONSTRUCTION(I); integer I; begin if I=1 then INIVEC(M0,M,RO,0); if I<4 then ELMVEC(M0,M,0,RO,C,BETHA[I]*HI); if I=4 then begin ELMVEC(M0,M,0,RO,C,-H); RHO:=NORMFUNCTION(NORM,RO); EC0:=EC1;EC1:=EC2;EC2:=RHO/H**Q; end end; procedure STEPSIZE; begin real HACC,HSTAB,HCR,HMAX,A,B,C; if not START then LOCAL ERROR BOUND; if START then begin H1:=H2:=HACC:=HSTART; EC2:=EC1:=1; KL:=1; START:=false end else if KL<3 then begin HACC:=(ETA/RHO)**(1/Q)*H2; if HACC>10*H2 then HACC:=10*H2 else KL:=KL+1 end else begin A:=(H0*(EC2-EC1)-H1*(EC1-EC0))/(H2*H0-H1*H1); H:=H2*(if ETA<RHO then (ETA/RHO)**(1/Q) else ALFA); if A>0 then begin B:=(EC2-EC1-A*(H2-H1))/H1; C:=EC2-A*H2-B*T2; HACC:=0; HMAX:=H; if ^ZEROIN(HACC,H,HACC**Q*(A*HACC+B*T+C)-ETA, "-3*H2) then HACC:=HMAX end else HACC:=H; if HACC<.5*H2 then HACC:=.5*H2; end; if HACC<HMIN then HACC:=HMIN; H:=HACC; if H*SIGMAL>1 then begin A:=ABS(DIAMETER/SIGMAL+"-14)/2; B:=2*ABS(SIN(PHIL)); BETAN:=(if A>B then 1/A else 1/B)/A; HSTAB:=ABS(BETAN/SIGMAL); if HSTAB<"-14*T then goto ENDOFEFT; if H>HSTAB then H:=HSTAB end; HCR:=H2*H2/H1; if KL>2 and ABS(H-HCR)<"-6*HCR then H:=if H<HCR then HCR*(1-"-7) else HCR*(1+"-7); if T+H>TE then begin LAST:=true; HSTART:=H; H:=TE-T end; H0:=H1;H1:=H2;H2:=H; end; procedure DIFFERENCE SCHEME; begin HI:=1; SIGMAL:=SIGMA; PHIL:=PHI; STEPSIZE; COEFFICIENT; for I:=1,2,3 do begin HI:=HI*H; if I>1 then DERIVATIVE(I,C); LOCALERRORCONSTRUCTION(I); ELMVEC(M0,M,0,U,C,BETA[I]*HI) end; T2:=T; K:=K+1; if LAST then begin LAST:=false; T:=TE; START:=true end else T:=T+H; DUPVEC(M0,M,0,C,U); DERIVATIVE(1,C); LOCALERRORCONSTRUCTION(4); OUTPUT; end; START:=true; LAST:=false; DUPVEC(M0,M,0,C,U); DERIVATIVE(1,C); if K=0 then begin LOCAL ERROR BOUND; HSTART:=ETA/NORMFUNCTION(NORM,C) end; NEXT LEVEL: DIFFERENCE SCHEME; if T^=TE then goto NEXT LEVEL; ENDOFEFT: end EXPONENTIAL FITTED TAYLOR; eop