{************************************************************************* * Procedures and Functions used By programs TZBESJ, TZBESK, TZBESY * * (Evalute Bessel Functions with complex argument, 1st to 3rd kind) * * ---------------------------------------------------------------------- * * Reference: AMOS, DONALD E., SANDIA NATIONAL LABORATORIES, 1983. * * * * Pascal Release By J-P Moreau, Paris (07/01/2005). * * (www.jpmoreau.fr) * *************************************************************************} UNIT CBESS0; {Procedures or functions using only unit COMPLEX} INTERFACE Uses COMPLEX; Function Power(y,x:Double): Double; Function Sign(a,b : Double) : Double; Procedure ZUCHK(YR, YI:double; var NZ:integer; ASCLE, TOL:double); Function DGAMLN(Z:Double; IERR:Integer): Double; Procedure ZUNHJ(ZR, ZI, FNU:double; IPMTR:integer; Var TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI:double); Procedure ZUNIK(ZRR, ZRI, FNU:double; IKFLG, IPMTR:Integer; TOL:double; INIT:integer; PHIR, PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI: double; Var CWRKR, CWRKI:VEC16); Procedure ZUOIK(ZR, ZI, FNU:Double; KODE, IKFLG, N: Integer; Var YR, YI:VEC; Var NUF:Integer; Var TOL, ELIM, ALIM:Double); Procedure ZS1S2(Var ZRR, ZRI, S1R, S1I, S2R, S2I:double; NZ:Integer; Var ASCLE, ALIM: double; Var IUF:Integer); Procedure ZSERI(ZR, ZI, FNU:Double; KODE, N: Integer; Var YR, YI: VEC; Var NZ: Integer; Var TOL, ELIM, ALIM:Double); Procedure ZASYI(ZR, ZI, FNU:Double; KODE, N: Integer; Var YR, YI: VEC; Var NZ: Integer; Var RL, TOL, ELIM, ALIM:Double); Procedure ZMLRI(ZR, ZI, FNU:Double; KODE, N:Integer; Var YR, YI:VEC; Var NZ:Integer; Var TOL:Double); Procedure ZSHCH(ZR, ZI: Double; Var CSHR, CSHI, CCHR, CCHI:Double); Procedure ZRATI(ZR, ZI, FNU:Double; Var N:Integer; Var CYR, CYI:VEC; TOL:Double); IMPLEMENTATION Function Power(y,x:Double): Double; Begin IF x<0 THEN EXIT; Power:=Exp(x*Ln(y)) End; Function Sign(a,b : Double) : Double; Begin if (b <0.0) then Sign := - Abs(a) else Sign := Abs(a) End; Procedure ZUCHK(YR, YI:double; var NZ:integer; ASCLE, TOL:double); {***BEGIN PROLOGUE ZUCHK !***REFER TO ZSERI,ZUOIK,ZUNK1,ZUNK2,ZUNI1,ZUNI2,ZKSCL ! ! Y ENTERS AS A SCALED QUANTITY WHOSE MAGNITUDE IS GREATER THAN ! EXP(-ALIM)=ASCLE=1000*D1MACH(1)/TOL. THE TEST IS MADE TO SEE ! IF THE MAGNITUDE OF THE REAL OR IMAGINARY PART WOULD UNDERFLOW ! WHEN Y IS SCALED (BY TOL) TO ITS PROPER VALUE. Y IS ACCEPTED ! IF THE UNDERFLOW IS AT LEAST ONE PRECISION BELOW THE MAGNITUDE ! OF THE LARGEST COMPONENT; OTHERWISE THE PHASE ANGLE DOES NOT HAVE ! ABSOLUTE ACCURACY AND AN UNDERFLOW IS ASSUMED. ! !***ROUTINES CALLED (NONE) !***END PROLOGUE ZUCHK ! ! COMPLEX Y } Label Return; Var SS, ST, WR, WI:double; Begin NZ := 0; WR := ABS(YR); WI := ABS(YI); ST := DMIN(WR,WI); IF ST > ASCLE THEN GOTO RETURN; SS := DMAX(WR,WI); ST := ST/TOL; IF SS < ST THEN NZ := 1; Return:End; {ZUCHK} FUNCTION DGAMLN(Z:Double; IERR:Integer): Double; {***BEGIN PROLOGUE DGAMLN !***DATE WRITTEN 830501 (YYMMDD) !***REVISION DATE 830501 (YYMMDD) !***CATEGORY NO. B5F !***KEYWORDS GAMMA FUNCTION,LOGARITHM OF GAMMA FUNCTION !***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES !***PURPOSE TO COMPUTE THE LOGARITHM OF THE GAMMA FUNCTION !***DESCRIPTION ! ! **** A DOUBLE PRECISION ROUTINE **** ! DGAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR ! Z.GT.0. THE ASYMPTOTIC EXPANSION IS USED TO GENERATE VALUES ! GREATER THAN ZMIN WHICH ARE ADJUSTED BY THE RECURSION ! G(Z+1)=Z*G(Z) FOR Z.LE.ZMIN. THE FUNCTION WAS MADE AS ! PORTABLE AS POSSIBLE BY COMPUTIMG ZMIN FROM THE NUMBER OF BASE ! 10 DIGITS IN A WORD, RLN=MAX(-ALOG10(R1MACH(4)),0.5E-18) ! LIMITED TO 18 DIGITS OF (RELATIVE) ACCURACY. ! ! SINCE INTEGER ARGUMENTS ARE COMMON, A TABLE LOOK UP ON 100 ! VALUES IS USED FOR SPEED OF EXECUTION. ! ! DESCRIPTION OF ARGUMENTS ! ! INPUT Z IS D0UBLE PRECISION ! Z - ARGUMENT, Z.GT.0.0D0 ! ! OUTPUT DGAMLN IS DOUBLE PRECISION ! DGAMLN - NATURAL LOG OF THE GAMMA FUNCTION AT Z.NE.0.0D0 ! IERR - ERROR FLAG ! IERR=0, NORMAL RETURN, COMPUTATION COMPLETED ! IERR=1, Z <= 0.0, NO COMPUTATION ! ! !***REFERENCES COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT ! BY D. E. AMOS, SAND83-0083, MAY, 1983. !***ROUTINES CALLED I1MACH,D1MACH (of unit COMPLEX). !***END PROLOGUE DGAMLN } Label 10,20,40,50,70, Return; Type pTGLN = ^TGLN; TGLN = Array[1..100] of double; Var CON, FLN, FZ, RLN, S, TLG, TRM, TST, T1, WDTOL, ZDMY, ZINC, ZM, ZMIN, ZP, ZSQ: Double; I, I1M, K, MZ, NZ: Integer; CF:array[1..22] of double; GLN: pTGLN; Begin { LNGAMMA(N), N:=1,100 } New(GLN); GLN^[1]:= 0.00000000000000000; GLN^[2]:= 0.00000000000000000; GLN^[3]:= 6.93147180559945309E-01; GLN^[4]:= 1.79175946922805500; GLN^[5]:= 3.17805383034794562; GLN^[6]:= 4.78749174278204599; GLN^[7]:= 6.57925121201010100; GLN^[8]:= 8.52516136106541430; GLN^[9]:= 1.06046029027452502E+01; GLN^[10]:=1.28018274800814696E+01; GLN^[11]:=1.51044125730755153E+01; GLN^[12]:=1.75023078458738858E+01; GLN^[13]:=1.99872144956618861E+01; GLN^[14]:=2.25521638531234229E+01; GLN^[15]:=2.51912211827386815E+01; GLN^[16]:=2.78992713838408916E+01; GLN^[17]:=3.06718601060806728E+01; GLN^[18]:=3.35050734501368889E+01; GLN^[19]:=3.63954452080330536E+01; GLN^[20]:=3.93398841871994940E+01; GLN^[21]:=4.23356164607534850E+01; GLN^[22]:=4.53801388984769080E+01; GLN^[23]:=4.84711813518352239E+01; GLN^[24]:=5.16066755677643736E+01; GLN^[25]:=5.47847293981123192E+01; GLN^[26]:=5.80036052229805199E+01; GLN^[27]:=6.12617017610020020E+01; GLN^[28]:=6.45575386270063311E+01; GLN^[29]:=6.78897431371815350E+01; GLN^[30]:=7.12570389671680090E+01; GLN^[31]:=7.46582363488301644E+01; GLN^[32]:=7.80922235533153106E+01; GLN^[33]:=8.15579594561150372E+01; GLN^[34]:=8.50544670175815174E+01; GLN^[35]:=8.85808275421976788E+01; GLN^[36]:=9.21361756036870925E+01; GLN^[37]:=9.57196945421432025E+01; GLN^[38]:=9.93306124547874269E+01; GLN^[39]:=1.02968198614513813E+02; GLN^[40]:=1.06631760260643459E+02; GLN^[41]:=1.10320639714757395E+02; GLN^[42]:=1.14034211781461703E+02; GLN^[43]:=1.17771881399745072E+02; GLN^[44]:=1.21533081515438634E+02; GLN^[45]:=1.25317271149356895E+02; GLN^[46]:=1.29123933639127215E+02; GLN^[47]:=1.32952575035616310E+02; GLN^[48]:=1.36802722637326368E+02; GLN^[49]:=1.40673923648234259E+02; GLN^[50]:=1.44565743946344886E+02; GLN^[51]:=1.48477766951773032E+02; GLN^[52]:=1.52409592584497358E+02; GLN^[53]:=1.56360836303078785E+02; GLN^[54]:=1.60331128216630907E+02; GLN^[55]:=1.64320112263195181E+02; GLN^[56]:=1.68327445448427652E+02; GLN^[57]:=1.72352797139162802E+02; GLN^[58]:=1.76395848406997352E+02; GLN^[59]:=1.80456291417543771E+02; GLN^[60]:=1.84533828861449491E+02; GLN^[61]:=1.88628173423671591E+02; GLN^[62]:=1.92739047287844902E+02; GLN^[63]:=1.96866181672889994E+02; GLN^[64]:=2.01009316399281527E+02; GLN^[65]:=2.05168199482641199E+02; GLN^[66]:=2.09342586752536836E+02; GLN^[67]:=2.13532241494563261E+02; GLN^[68]:=2.17736934113954227E+02; GLN^[69]:=2.21956441819130334E+02; GLN^[70]:=2.26190548323727593E+02; GLN^[71]:=2.30439043565776952E+02; GLN^[72]:=2.34701723442818268E+02; GLN^[73]:=2.38978389561834323E+02; GLN^[74]:=2.43268849002982714E+02; GLN^[75]:=2.47572914096186884E+02; GLN^[76]:=2.51890402209723194E+02; GLN^[77]:=2.56221135550009525E+02; GLN^[78]:=2.60564940971863209E+02; GLN^[79]:=2.64921649798552801E+02; GLN^[80]:=2.69291097651019823E+02; GLN^[81]:=2.73673124285693704E+02; GLN^[82]:=2.78067573440366143E+02; GLN^[83]:=2.82474292687630396E+02; GLN^[84]:=2.86893133295426994E+02; GLN^[85]:=2.91323950094270308E+02; GLN^[86]:=2.95766601350760624E+02; GLN^[87]:=3.00220948647014132E+02; GLN^[88]:=3.04686856765668715E+02; GLN^[89]:=3.09164193580146922E+02; GLN^[90]:=3.13652829949879062E+02; GLN^[91]:=3.18152639620209327E+02; GLN^[92]:=3.22663499126726177E+02; GLN^[93]:=3.27185287703775217E+02; GLN^[94]:=3.31717887196928473E+02; GLN^[95]:=3.36261181979198477E+02; GLN^[96]:=3.40815058870799018E+02; GLN^[97]:=3.45379407062266854E+02; GLN^[98]:=3.49954118040770237E+02; GLN^[99]:=3.54539085519440809E+02; GLN^[100]:=3.59134205369575399E+02; { COEFFICIENTS OF ASYMPTOTIC EXPANSION } CF[1]:= 8.33333333333333333E-02; CF[2]:= -2.77777777777777778E-03; CF[3]:= 7.93650793650793651E-04; CF[4]:= -5.95238095238095238E-04; CF[5]:= 8.41750841750841751E-04; CF[6]:= -1.91752691752691753E-03; CF[7]:= 6.41025641025641026E-03; CF[8]:= -2.95506535947712418E-02; CF[9]:= 1.79644372368830573E-01; CF[10]:=-1.39243221690590112E+00; CF[11]:=1.34028640441683920E+01; CF[12]:=-1.56848284626002017E+02; CF[13]:=2.19310333333333333E+03; CF[14]:=-3.61087712537249894E+04; CF[15]:=6.91472268851313067E+05; CF[16]:=-1.52382215394074162E+07; CF[17]:=3.82900751391414141E+08; CF[18]:=-1.08822660357843911E+10; CF[19]:=3.47320283765002252E+11; CF[20]:=-1.23696021422692745E+13; CF[21]:=4.88788064793079335E+14; CF[22]:=-2.13203339609193739E+16; { LN(2*PI) } CON:=1.83787706640934548; IERR:=0; IF Z <= 0.0 THEN GOTO 70; IF Z > 101.0 THEN GOTO 10; NZ := Round(Z); FZ := Z - 1.0*NZ; IF FZ > 0.0 THEN GOTO 10; IF NZ > 100 THEN GOTO 10; DGAMLN := GLN^[NZ]; GOTO RETURN; 10: WDTOL := D1MACH(4); WDTOL := DMAX(WDTOL,0.5E-18); I1M := I1MACH(14); RLN := D1MACH(5)*I1M; FLN := DMIN(RLN,20.0); FLN := DMAX(FLN,3.0); FLN := FLN - 3.0; ZM := 1.8000 + 0.3875*FLN; MZ := Round(ZM) + 1; ZMIN := 1.0*MZ; ZDMY := Z; ZINC := 0.0; IF Z >= ZMIN THEN GOTO 20; ZINC := ZMIN - 1.0*NZ; ZDMY := Z + ZINC; 20: ZP := 1.0/ZDMY; T1 := CF[1]*ZP; S := T1; IF ZP < WDTOL THEN GOTO 40; ZSQ := ZP*ZP; TST := T1*WDTOL; For K:=2 to 22 do begin ZP := ZP*ZSQ; TRM := CF[K]*ZP; IF ABS(TRM) < TST THEN GOTO 40; S := S + TRM end; 40: IF ZINC <> 0.0 THEN GOTO 50; TLG := Ln(Z); DGAMLN := Z*(TLG-1.0) + 0.5*(CON-TLG) + S; GOTO RETURN; 50: ZP := 1.0; NZ := Round(ZINC); For I:=1 to NZ do ZP := ZP*(Z+1.0*(I-1)); TLG := Ln(ZDMY); DGAMLN := ZDMY*(TLG-1.0) - Ln(ZP) + 0.5*(CON-TLG) + S; GOTO RETURN; 70: IERR:=1; RETURN:Dispose(GLN) End; {DGAMLN} Procedure ZUNIK(ZRR, ZRI, FNU:double; IKFLG, IPMTR:Integer; TOL:double; INIT:integer; PHIR, PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI: double; Var CWRKR, CWRKI:VEC16); {***BEGIN PROLOGUE ZUNIK !***REFER TO ZBESI,ZBESK ! ! ZUNIK COMPUTES PARAMETERS FOR THE UNIFORM ASYMPTOTIC ! EXPANSIONS OF THE I AND K FUNCTIONS ON IKFLG= 1 OR 2 ! RESPECTIVELY BY: ! ! W(FNU,ZR) = PHI*EXP(ZETA)*SUM ! ! WHERE ZETA=-ZETA1 + ZETA2 OR ! ZETA1 - ZETA2 ! ! THE FIRST CALL MUST HAVE INIT=0. SUBSEQUENT CALLS WITH THE ! SAME ZR AND FNU WILL RETURN THE I OR K FUNCTION ON IKFLG= ! 1 OR 2 WITH NO CHANGE IN INIT. CWRK IS A COMPLEX WORK ! ARRAY. IPMTR=0 COMPUTES ALL PARAMETERS. IPMTR=1 COMPUTES PHI, ! ZETA1,ZETA2. ! !***ROUTINES CALLED ZDIV,ZLOG,ZSQRT !***END PROLOGUE ZUNIK ! COMPLEX CFN,CON,CONE,CRFN,CWRK,CZERO,P,HI,S,SR,SUM,T,T2,ZETA1, ! ZETA2,ZN,ZR } Label 30, 40, 60, Return; Type pTC = ^TC; TC = Array[1..120] of Double; Var AC, CONEI, CONER, CRFNI, CRFNR, RFN, SI, SR, SRI, SRR, STI, STR, TEST, TI, TR, T2I, T2R, ZEROI, ZEROR, ZNI, ZNR: Double; I, IDUM, J, K, L: Integer; C: pTC; CON: array[1..2] of Double; Begin ZEROR:=0.0; ZEROI:= 0.0; CONER:=1.0; CONEI:=0.0; CON[1]:=3.98942280401432678E-01; CON[2]:=1.25331413731550025; {Initialize table C} New(C); C^[1]:= 1.00000000000000000E+00; C^[2]:= -2.08333333333333333E-01; C^[3]:= 1.25000000000000000E-01; C^[4]:= 3.34201388888888889E-01; C^[5]:= -4.01041666666666667E-01; C^[6]:= 7.03125000000000000E-02; C^[7]:= -1.02581259645061728E+00; C^[8]:= 1.84646267361111111E+00; C^[9]:= -8.91210937500000000E-01; C^[10]:= 7.32421875000000000E-02; C^[11]:= 4.66958442342624743E+00; C^[12]:=-1.12070026162229938E+01; C^[13]:= 8.78912353515625000E+00; C^[14]:=-2.36408691406250000E+00; C^[15]:= 1.12152099609375000E-01; C^[16]:=-2.82120725582002449E+01; C^[17]:= 8.46362176746007346E+01; C^[18]:=-9.18182415432400174E+01; C^[19]:= 4.25349987453884549E+01; C^[20]:=-7.36879435947963170E+00; C^[21]:= 2.27108001708984375E-01; C^[22]:= 2.12570130039217123E+02; C^[23]:=-7.65252468141181642E+02; C^[24]:= 1.05999045252799988E+03; C^[25]:=-6.99579627376132541E+02; C^[26]:= 2.18190511744211590E+02; C^[27]:=-2.64914304869515555E+01; C^[28]:= 5.72501420974731445E-01; C^[29]:=-1.91945766231840700E+03; C^[30]:= 8.06172218173730938E+03; C^[31]:=-1.35865500064341374E+04; C^[32]:= 1.16553933368645332E+04; C^[33]:=-5.30564697861340311E+03; C^[34]:= 1.20090291321635246E+03; C^[35]:=-1.08090919788394656E+02; C^[36]:= 1.72772750258445740E+00; C^[37]:= 2.02042913309661486E+04; C^[38]:=-9.69805983886375135E+04; C^[39]:= 1.92547001232531532E+05; C^[40]:=-2.03400177280415534E+05; C^[41]:= 1.22200464983017460E+05; C^[42]:=-4.11926549688975513E+04; C^[43]:= 7.10951430248936372E+03; C^[44]:=-4.93915304773088012E+02; C^[45]:= 6.07404200127348304E+00; C^[46]:=-2.42919187900551333E+05; C^[47]:= 1.31176361466297720E+06; C^[48]:=-2.99801591853810675E+06; C^[49]:= 3.76327129765640400E+06; C^[50]:=-2.81356322658653411E+06; C^[51]:= 1.26836527332162478E+06; C^[52]:=-3.31645172484563578E+05; C^[53]:= 4.52187689813627263E+04; C^[54]:=-2.49983048181120962E+03; C^[55]:= 2.43805296995560639E+01; C^[56]:= 3.28446985307203782E+06; C^[57]:=-1.97068191184322269E+07; C^[58]:= 5.09526024926646422E+07; C^[59]:=-7.41051482115326577E+07; C^[60]:= 6.63445122747290267E+07; C^[61]:=-3.75671766607633513E+07; C^[62]:= 1.32887671664218183E+07; C^[63]:=-2.78561812808645469E+06; C^[64]:= 3.08186404612662398E+05; C^[65]:=-1.38860897537170405E+04; C^[66]:= 1.10017140269246738E+02; C^[67]:=-4.93292536645099620E+07; C^[68]:= 3.25573074185765749E+08; C^[69]:=-9.39462359681578403E+08; C^[70]:= 1.55359689957058006E+09; C^[71]:=-1.62108055210833708E+09; C^[72]:= 1.10684281682301447E+09; C^[73]:=-4.95889784275030309E+08; C^[74]:= 1.42062907797533095E+08; C^[75]:=-2.44740627257387285E+07; C^[76]:= 2.24376817792244943E+06; C^[77]:=-8.40054336030240853E+04; C^[78]:= 5.51335896122020586E+02; C^[79]:= 8.14789096118312115E+08; C^[80]:=-5.86648149205184723E+09; C^[81]:= 1.86882075092958249E+10; C^[82]:=-3.46320433881587779E+10; C^[83]:= 4.12801855797539740E+10; C^[84]:=-3.30265997498007231E+10; C^[85]:= 1.79542137311556001E+10; C^[86]:=-6.56329379261928433E+09; C^[87]:= 1.55927986487925751E+09; C^[88]:=-2.25105661889415278E+08; C^[89]:= 1.73951075539781645E+07; C^[90]:=-5.49842327572288687E+05; C^[91]:= 3.03809051092238427E+03; C^[92]:=-1.46792612476956167E+10; C^[93]:= 1.14498237732025810E+11; C^[94]:=-3.99096175224466498E+11; C^[95]:= 8.19218669548577329E+11; C^[96]:=-1.09837515608122331E+12; C^[97]:= 1.00815810686538209E+12; C^[98]:= -6.45364869245376503E+11; C^[99]:= 2.87900649906150589E+11; C^[100]:=-8.78670721780232657E+10; C^[101]:= 1.76347306068349694E+10; C^[102]:=-2.16716498322379509E+09; C^[103]:= 1.43157876718888981E+08; C^[104]:=-3.87183344257261262E+06; C^[105]:= 1.82577554742931747E+04; C^[106]:= 2.86464035717679043E+11; C^[107]:=-2.40629790002850396E+12; C^[108]:= 9.10934118523989896E+12; C^[109]:=-2.05168994109344374E+13; C^[110]:= 3.05651255199353206E+13; C^[111]:=-3.16670885847851584E+13; C^[112]:= 2.33483640445818409E+13; C^[113]:=-1.23204913055982872E+13; C^[114]:= 4.61272578084913197E+12; C^[115]:=-1.19655288019618160E+12; C^[116]:= 2.05914503232410016E+11; C^[117]:=-2.18229277575292237E+10; C^[118]:= 1.24700929351271032E+09; C^[119]:=-2.91883881222208134E+07; C^[120]:= 1.18838426256783253E+05; IF INIT <> 0 THEN GOTO 40; {----------------------------------------------------------------------- ! INITIALIZE ALL VARIABLES !----------------------------------------------------------------------} RFN := 1.0/FNU; TR := ZRR*RFN; TI := ZRI*RFN; SR := CONER + (TR*TR-TI*TI); SI := CONEI + (TR*TI+TI*TR); ZSQRT(SR, SI, SRR, SRI); STR := CONER + SRR; STI := CONEI + SRI; ZDIV(STR, STI, TR, TI, ZNR, ZNI); ZLOG(ZNR, ZNI, STR, STI, IDUM); ZETA1R := FNU*STR; ZETA1I := FNU*STI; ZETA2R := FNU*SRR; ZETA2I := FNU*SRI; ZDIV(CONER, CONEI, SRR, SRI, TR, TI); SRR := TR*RFN; SRI := TI*RFN; ZSQRT(SRR, SRI, CWRKR[16], CWRKI[16]); PHIR := CWRKR[16]*CON[IKFLG]; PHII := CWRKI[16]*CON[IKFLG]; IF IPMTR <> 0 THEN GOTO RETURN; ZDIV(CONER, CONEI, SR, SI, T2R, T2I); CWRKR[1] := CONER; CWRKI[1] := CONEI; CRFNR := CONER; CRFNI := CONEI; AC := 1.0; L := 1; For K:=2 to 15 do begin SR := ZEROR; SI := ZEROI; For J:=1 to K do begin L := L + 1; STR := SR*T2R - SI*T2I + C^[L]; SI := SR*T2I + SI*T2R; SR := STR end; STR := CRFNR*SRR - CRFNI*SRI; CRFNI := CRFNR*SRI + CRFNI*SRR; CRFNR := STR; CWRKR[K] := CRFNR*SR - CRFNI*SI; CWRKI[K] := CRFNR*SI + CRFNI*SR; AC := AC*RFN; TEST := ABS(CWRKR[K]) + ABS(CWRKI[K]); IF (AC < TOL) AND (TEST < TOL) THEN GOTO 30; end; K := 15; 30: INIT := K; 40: IF IKFLG = 2 THEN GOTO 60; {----------------------------------------------------------------------- ! COMPUTE SUM FOR THE I FUNCTION !----------------------------------------------------------------------} SR := ZEROR; SI := ZEROI; For I:=1 to INIT do begin SR := SR + CWRKR[I]; SI := SI + CWRKI[I] end; SUMR := SR; SUMI := SI; PHIR := CWRKR[16]*CON[1]; PHII := CWRKI[16]*CON[1]; GOTO RETURN; {----------------------------------------------------------------------- ! COMPUTE SUM FOR THE K FUNCTION !----------------------------------------------------------------------} 60: SR := ZEROR; SI := ZEROI; TR := CONER; For I:=1 to INIT do begin SR := SR + TR*CWRKR[I]; SI := SI + TR*CWRKI[I]; TR := -TR end; SUMR := SR; SUMI := SI; PHIR := CWRKR[16]*CON[2]; PHII := CWRKI[16]*CON[2]; Return: Dispose(C); End; {ZUNIK} Procedure ZUNHJ(ZR, ZI, FNU:double; IPMTR:integer; Var TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI:double); {***BEGIN PROLOGUE ZUNHJ !***REFER TO ZBESI,ZBESK ! ! REFERENCES ! HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ AND I.A. ! STEGUN, AMS55, NATIONAL BUREAU OF STANDARDS, 1965, CHAPTER 9. ! ! ASYMPTOTICS AND SPECIAL FUNCTIONS BY F.W.J. OLVER, ACADEMIC ! PRESS, N.Y., 1974, PAGE 420 ! ! ABSTRACT ! ZUNHJ COMPUTES PARAMETERS FOR BESSEL FUNCTIONS C(FNU,Z) = ! J(FNU,Z), Y(FNU,Z) OR H(I,FNU,Z) I=1,2 FOR LARGE ORDERS FNU ! BY MEANS OF THE UNIFORM ASYMPTOTIC EXPANSION ! ! C(FNU,Z)=C1*PHI*( ASUM*AIRY(ARG) + C2*BSUM*DAIRY(ARG) ) ! ! FOR PROPER CHOICES OF C1, C2, AIRY AND DAIRY WHERE AIRY IS ! AN AIRY FUNCTION AND DAIRY IS ITS DERIVATIVE. ! ! (2/3)*FNU*ZETA^1.5 = ZETA1-ZETA2, ! ! ZETA1=0.5*FNU*CLOG((1+W)/(1-W)), ZETA2=FNU*W FOR SCALING ! PURPOSES IN AIRY FUNCTIONS FROM CAIRY OR CBIRY. ! ! MCONJ=SIGN OF AIMAG(Z), BUT IS AMBIGUOUS WHEN Z IS REAL AND ! MUST BE SPECIFIED. IPMTR=0 RETURNS ALL PARAMETERS. IPMTR=1 ! COMPUTES ALL EXCEPT ASUM AND BSUM. ! !***ROUTINES CALLED ZABS,ZDIV,ZLOG,ZSQRT !***END PROLOGUE ZUNHJ ! COMPLEX ARG,ASUM,BSUM,CFNU,CONE,CR,CZERO,DR,P,PHI,PRZTH,PTFN, ! RFN13,RTZTA,RZTH,SUMA,SUMB,TFN,T2,UP,W,W2,Z,ZA,ZB,ZC,ZETA,ZETA1, ! ZETA2,ZTH } Label 20, 50,60, 80,90,110,120,130,140,180,200,220, Return; Type pVEC30 = ^VEC30; VEC30 = Array[1..30] of Double; pVEC105 = ^VEC105; VEC105 = Array[1..105] of Double; pVEC210 = ^VEC210; VEC210 = Array[1..210] of Double; Var ANG, ATOL, AW2, AZTH, BTOL, CONEI, CONER, EX1, EX2, FN13, FN23, GPI, HPI, PP, PRZTHI, PRZTHR, PTFNI, PTFNR, RAW, RAW2, RAZTH, RFNU, RFNU2, RFN13, RTZTI, RTZTR, RZTHI, RZTHR, STI, STR, SUMAI, SUMAR, SUMBI, SUMBR, TEST, TFNI, TFNR, THPI, TZAI, TZAR, T2I, T2R, WI, WR, W2I, W2R, ZAI, ZAR, ZBI, ZBR, ZCI, ZCR, ZEROI, ZEROR, ZETAI, ZETAR, ZTHI, ZTHR: Double; IAS, IBS, IS, J, JR, JU, K, KMAX, KP1, KS, L, LR, LRP1, L1, L2, M, IDUM: Integer; AR, BR, UPR, UPI, CRR, CRI, DRR, DRI: VEC16; C: pVEC105; ALFA, BETA: pVEC210; GAMA, AP, PR, PI: pVEC30; Begin New(C); New(ALFA); New(BETA); New(GAMA); New(AP); New(PR); New(PI); AR[1]:= 1.00000000000000000E+00; AR[2]:= 1.04166666666666667E-01; AR[3]:= 8.35503472222222222E-02; AR[4]:= 1.28226574556327160E-01; AR[5]:= 2.91849026464140464E-01; AR[6]:= 8.81627267443757652E-01; AR[7]:= 3.32140828186276754E+00; AR[8]:= 1.49957629868625547E+01; AR[9]:= 7.89230130115865181E+01; AR[10]:=4.74451538868264323E+02; AR[11]:=3.20749009089066193E+03; AR[12]:=2.40865496408740049E+04; AR[13]:=1.98923119169509794E+05; AR[14]:=1.79190200777534383E+06; AR[15]:=0.0; AR[16]:=0.0; BR[1]:= 1.00000000000000000E+00; BR[2]:= -1.45833333333333333E-01; BR[3]:= -9.87413194444444444E-02; BR[4]:= -1.43312053915895062E-01; BR[5]:= -3.17227202678413548E-01; BR[6]:= -9.42429147957120249E-01; BR[7]:= -3.51120304082635426E+00; BR[8]:= -1.57272636203680451E+01; BR[9]:= -8.22814390971859444E+01; BR[10]:=-4.92355370523670524E+02; BR[11]:=-3.31621856854797251E+03; BR[12]:=-2.48276742452085896E+04; BR[13]:=-2.04526587315129788E+05; BR[14]:=-1.83844491706820990E+06; BR[15]:=0.0; BR[16]:=0.0; C^[1]:= 1.00000000000000000E+00; C^[2]:= -2.08333333333333333E-01; C^[3]:= 1.25000000000000000E-01; C^[4]:= 3.34201388888888889E-01; C^[5]:= -4.01041666666666667E-01; C^[6]:= 7.03125000000000000E-02; C^[7]:= -1.02581259645061728E+00; C^[8]:= 1.84646267361111111E+00; C^[9]:= -8.91210937500000000E-01; C^[10]:= 7.32421875000000000E-02; C^[11]:= 4.66958442342624743E+00; C^[12]:=-1.12070026162229938E+01; C^[13]:= 8.78912353515625000E+00; C^[14]:=-2.36408691406250000E+00; C^[15]:= 1.12152099609375000E-01; C^[16]:=-2.82120725582002449E+01; C^[17]:= 8.46362176746007346E+01; C^[18]:=-9.18182415432400174E+01; C^[19]:= 4.25349987453884549E+01; C^[20]:=-7.36879435947963170E+00; C^[21]:= 2.27108001708984375E-01; C^[22]:= 2.12570130039217123E+02; C^[23]:=-7.65252468141181642E+02; C^[24]:= 1.05999045252799988E+03; C^[25]:=-6.99579627376132541E+02; C^[26]:= 2.18190511744211590E+02; C^[27]:=-2.64914304869515555E+01; C^[28]:= 5.72501420974731445E-01; C^[29]:=-1.91945766231840700E+03; C^[30]:= 8.06172218173730938E+03; C^[31]:=-1.35865500064341374E+04; C^[32]:= 1.16553933368645332E+04; C^[33]:=-5.30564697861340311E+03; C^[34]:= 1.20090291321635246E+03; C^[35]:=-1.08090919788394656E+02; C^[36]:= 1.72772750258445740E+00; C^[37]:= 2.02042913309661486E+04; C^[38]:=-9.69805983886375135E+04; C^[39]:= 1.92547001232531532E+05; C^[40]:=-2.03400177280415534E+05; C^[41]:= 1.22200464983017460E+05; C^[42]:=-4.11926549688975513E+04; C^[43]:= 7.10951430248936372E+03; C^[44]:=-4.93915304773088012E+02; C^[45]:= 6.07404200127348304E+00; C^[46]:=-2.42919187900551333E+05; C^[47]:= 1.31176361466297720E+06; C^[48]:=-2.99801591853810675E+06; C^[49]:= 3.76327129765640E+0006; C^[50]:=-2.81356322658653E+0006; C^[51]:= 1.26836527332162E+0006; C^[52]:=-3.31645172484564E+0005; C^[53]:= 4.52187689813627E+0004; C^[54]:=-2.49983048181121E+0003; C^[55]:= 2.43805296995561E+0001; C^[56]:= 3.28446985307204E+0006; C^[57]:=-1.97068191184322E+0007; C^[58]:= 5.09526024926646E+0007; C^[59]:=-7.41051482115327E+0007; C^[60]:= 6.63445122747290E+0007; C^[61]:=-3.75671766607634E+0007; C^[62]:= 1.32887671664218E+0007; C^[63]:=-2.78561812808645E+0006; C^[64]:= 3.08186404612662E+0005; C^[65]:=-1.38860897537170E+0004; C^[66]:= 1.10017140269247E+0002; C^[67]:=-4.93292536645100E+0007; C^[68]:= 3.25573074185766E+0008; C^[69]:=-9.39462359681578E+0008; C^[70]:= 1.55359689957058E+0009; C^[71]:=-1.62108055210834E+0009; C^[72]:= 1.10684281682301E+0009; C^[73]:=-4.95889784275030E+0008; C^[74]:= 1.42062907797533E+0008; C^[75]:=-2.44740627257387E+0007; C^[76]:= 2.24376817792245E+0006; C^[77]:=-8.40054336030241E+0004; C^[78]:= 5.51335896122021E+0002; C^[79]:= 8.14789096118312E+0008; C^[80]:=-5.86648149205185E+0009; C^[81]:= 1.86882075092958E+0010; C^[82]:=-3.46320433881588E+0010; C^[83]:= 4.12801855797540E+0010; C^[84]:=-3.30265997498007E+0010; C^[85]:= 1.79542137311556E+0010; C^[86]:=-6.56329379261928E+0009; C^[87]:= 1.55927986487926E+0009; C^[88]:=-2.25105661889415E+0008; C^[89]:= 1.73951075539782E+0007; C^[90]:=-5.49842327572289E+0005; C^[91]:= 3.03809051092238E+0003; C^[92]:=-1.46792612476956E+0010; C^[93]:= 1.14498237732026E+0011; C^[94]:=-3.99096175224466E+0011; C^[95]:= 8.19218669548577E+0011; C^[96]:=-1.09837515608122E+0012; C^[97]:= 1.00815810686538E+0012; C^[98]:=-6.45364869245376E+0011; C^[99]:= 2.87900649906151E+0011; C^[100]:=-8.78670721780233E+0010; C^[101]:= 1.76347306068350E+0010; C^[102]:=-2.16716498322379E+0009; C^[103]:= 1.43157876718889E+0008; C^[104]:=-3.87183344257261E+0006; C^[105]:= 1.82577554742932E+0004; ALFA^[1]:=-4.44444444444444E-0003; ALFA^[2]:=-9.22077922077922E-0004; ALFA^[3]:=-8.84892884892885E-0005; ALFA^[4]:= 1.65927687832450E-0004; ALFA^[5]:= 2.46691372741793E-0004; ALFA^[6]:= 2.65995589346255E-0004; ALFA^[7]:= 2.61824297061501E-0004; ALFA^[8]:= 2.48730437344656E-0004; ALFA^[9]:= 2.32721040083232E-0004; ALFA^[10]:= 2.16362485712365E-0004; ALFA^[11]:= 2.00738858762752E-0004; ALFA^[12]:= 1.86267636637545E-0004; ALFA^[13]:= 1.73060775917876E-0004; ALFA^[14]:= 1.61091705929016E-0004; ALFA^[15]:= 1.50274774160908E-0004; ALFA^[16]:= 1.40503497391270E-0004; ALFA^[17]:= 1.31668816545923E-0004; ALFA^[18]:= 1.23667445598253E-0004; ALFA^[19]:= 1.16405271474738E-0004; ALFA^[20]:= 1.09798298372713E-0004; ALFA^[21]:= 1.03772410422993E-0004; ALFA^[22]:= 9.82626078369363E-0005; ALFA^[23]:= 9.32120517249503E-0005; ALFA^[24]:= 8.85710852478712E-0005; ALFA^[25]:= 8.42963105715700E-0005; ALFA^[26]:= 8.03497548407791E-0005; ALFA^[27]:= 7.66981345359207E-0005; ALFA^[28]:= 7.33122157481778E-0005; ALFA^[29]:= 7.01662625163141E-0005; ALFA^[30]:= 6.72375633790160E-0005; ALFA^[31]:= 6.93735541354589E-0004; ALFA^[32]:= 2.32241745182922E-0004; ALFA^[33]:=-1.41986273556691E-0005; ALFA^[34]:=-1.16444931672049E-0004; ALFA^[35]:=-1.50803558053049E-0004; ALFA^[36]:=-1.55121924918096E-0004; ALFA^[37]:=-1.46809756646466E-0004; ALFA^[38]:=-1.33815503867491E-0004; ALFA^[39]:=-1.19744975684254E-0004; ALFA^[40]:=-1.06184319207974E-0004; ALFA^[41]:=-9.37699549891194E-0005; ALFA^[42]:=-8.26923045588193E-0005; ALFA^[43]:=-7.29374348155221E-0005; ALFA^[44]:=-6.44042357721016E-0005; ALFA^[45]:=-5.69611566009369E-0005; ALFA^[46]:=-5.04731044303562E-0005; ALFA^[47]:=-4.48134868008883E-0005; ALFA^[48]:=-3.98688727717599E-0005; ALFA^[49]:=-3.55400532972043E-0005; ALFA^[50]:=-3.17414256609023E-0005; ALFA^[51]:=-2.83996793904175E-0005; ALFA^[52]:=-2.54522720634871E-0005; ALFA^[53]:=-2.28459297164725E-0005; ALFA^[54]:=-2.05352753106481E-0005; ALFA^[55]:=-1.84816217627666E-0005; ALFA^[56]:=-1.66519330021394E-0005; ALFA^[57]:=-1.50179412980119E-0005; ALFA^[58]:=-1.35554031379041E-0005; ALFA^[59]:=-1.22434746473858E-0005; ALFA^[60]:=-1.10641884811308E-0005; ALFA^[61]:=-3.54211971457744E-0004; ALFA^[62]:=-1.56161263945159E-0004; ALFA^[63]:= 3.04465503594936E-0005; ALFA^[64]:= 1.30198655773243E-0004; ALFA^[65]:= 1.67471106699712E-0004; ALFA^[66]:= 1.70222587683593E-0004; ALFA^[67]:= 1.56501427608595E-0004; ALFA^[68]:= 1.36339170977445E-0004; ALFA^[69]:= 1.14886692029825E-0004; ALFA^[70]:= 9.45869093034688E-0005; ALFA^[71]:= 7.64498419250898E-0005; ALFA^[72]:= 6.07570334965197E-0005; ALFA^[73]:= 4.74394299290509E-0005; ALFA^[74]:= 3.62757512005344E-0005; ALFA^[75]:= 2.69939714979225E-0005; ALFA^[76]:= 1.93210938247939E-0005; ALFA^[77]:= 1.30056674793963E-0005; ALFA^[78]:= 7.82620866744497E-0006; ALFA^[79]:= 3.59257485819352E-0006; ALFA^[80]:= 1.44040049814252E-0007; ALFA^[81]:=-2.65396769697939E-0006; ALFA^[82]:=-4.91346867098486E-0006; ALFA^[83]:=-6.72739296091248E-0006; ALFA^[84]:=-8.17269379678658E-0006; ALFA^[85]:=-9.31304715093561E-0006; ALFA^[86]:=-1.02011418798016E-0005; ALFA^[87]:=-1.08805962510593E-0005; ALFA^[88]:=-1.13875481509604E-0005; ALFA^[89]:=-1.17519675674556E-0005; ALFA^[90]:=-1.19987364870944E-0005; ALFA^[91]:= 3.78194199201773E-0004; ALFA^[92]:= 2.02471952761816E-0004; ALFA^[93]:=-6.37938506318862E-0005; ALFA^[94]:=-2.38598230603006E-0004; ALFA^[95]:=-3.10916256027362E-0004; ALFA^[96]:=-3.13680115247576E-0004; ALFA^[97]:=-2.78950273791323E-0004; ALFA^[98]:=-2.28564082619141E-0004; ALFA^[99]:=-1.75245280340847E-0004; ALFA^[100]:=-1.25544063060690E-0004; ALFA^[101]:=-8.22982872820208E-0005; ALFA^[102]:=-4.62860730588116E-0005; ALFA^[103]:=-1.72334302366962E-0005; ALFA^[104]:= 5.60690482304602E-0006; ALFA^[105]:= 2.31395443148287E-0005; ALFA^[106]:= 3.62642745856794E-0005; ALFA^[107]:= 4.58006124490189E-0005; ALFA^[108]:= 5.24595294959114E-0005; ALFA^[109]:= 5.68396208545815E-0005; ALFA^[110]:= 5.94349820393104E-0005; ALFA^[111]:= 6.06478527578422E-0005; ALFA^[112]:= 6.08023907788436E-0005; ALFA^[113]:= 6.01577894539460E-0005; ALFA^[114]:= 5.89199657344698E-0005; ALFA^[115]:= 5.72515823777593E-0005; ALFA^[116]:= 5.52804375585853E-0005; ALFA^[117]:= 5.31063773802880E-0005; ALFA^[118]:= 5.08069302012326E-0005; ALFA^[119]:= 4.84418647620095E-0005; ALFA^[120]:= 4.60568581607475E-0005; ALFA^[121]:=-6.91141397288294E-0004; ALFA^[122]:=-4.29976633058872E-0004; ALFA^[123]:= 1.83067735980039E-0004; ALFA^[124]:= 6.60088147542014E-0004; ALFA^[125]:= 8.75964969951186E-0004; ALFA^[126]:= 8.77335235958236E-0004; ALFA^[127]:= 7.49369585378991E-0004; ALFA^[128]:= 5.63832329756981E-0004; ALFA^[129]:= 3.68059319971443E-0004; ALFA^[130]:= 1.88464535514456E-0004; ALFA^[131]:= 3.70663057664904E-0005; ALFA^[132]:=-8.28520220232137E-0005; ALFA^[133]:=-1.72751952869173E-0004; ALFA^[134]:=-2.36314873605873E-0004; ALFA^[135]:=-2.77966150694907E-0004; ALFA^[136]:=-3.02079514155457E-0004; ALFA^[137]:=-3.12594712643820E-0004; ALFA^[138]:=-3.12872558758067E-0004; ALFA^[139]:=-3.05678038466324E-0004; ALFA^[140]:=-2.93226470614557E-0004; ALFA^[141]:=-2.77255655582935E-0004; ALFA^[142]:=-2.59103928467032E-0004; ALFA^[143]:=-2.39784014396480E-0004; ALFA^[144]:=-2.20048260045423E-0004; ALFA^[145]:=-2.00443911094971E-0004; ALFA^[146]:=-1.81358692210971E-0004; ALFA^[147]:=-1.63057674478657E-0004; ALFA^[148]:=-1.45712672175206E-0004; ALFA^[149]:=-1.29425421983925E-0004; ALFA^[150]:=-1.14245691942446E-0004; ALFA^[151]:= 1.92821964248776E-0003; ALFA^[152]:= 1.35592576302022E-0003; ALFA^[153]:=-7.17858090421303E-0004; ALFA^[154]:=-2.58084802575270E-0003; ALFA^[155]:=-3.49271130826168E-0003; ALFA^[156]:=-3.46986299340961E-0003; ALFA^[157]:=-2.82285233351310E-0003; ALFA^[158]:=-1.88103076404891E-0003; ALFA^[159]:=-8.89531718383948E-0004; ALFA^[160]:= 3.87912102631035E-0006; ALFA^[161]:= 7.28688540119691E-0004; ALFA^[162]:= 1.26566373053458E-0003; ALFA^[163]:= 1.62518158372674E-0003; ALFA^[164]:= 1.83203153216373E-0003; ALFA^[165]:= 1.91588388990528E-0003; ALFA^[166]:= 1.90588846755546E-0003; ALFA^[167]:= 1.82798982421826E-0003; ALFA^[168]:= 1.70389506421122E-0003; ALFA^[169]:= 1.55097127171098E-0003; ALFA^[170]:= 1.38261421852276E-0003; ALFA^[171]:= 1.20881424230065E-0003; ALFA^[172]:= 1.03676532638345E-0003; ALFA^[173]:= 8.71437918068619E-0004; ALFA^[174]:= 7.16080155297701E-0004; ALFA^[175]:= 5.72637002558129E-0004; ALFA^[176]:= 4.42089819465802E-0004; ALFA^[177]:= 3.24724948503091E-0004; ALFA^[178]:= 2.20342042730247E-0004; ALFA^[179]:= 1.28412898401354E-0004; ALFA^[180]:= 4.82005924552095E-0005; BETA^[1]:= 1.79988721413553E-0002; BETA^[2]:= 5.59964911064388E-0003; BETA^[3]:= 2.88501402231133E-0003; BETA^[4]:= 1.80096606761054E-0003; BETA^[5]:= 1.24753110589199E-0003; BETA^[6]:= 9.22878876572938E-0004; BETA^[7]:= 7.14430421727287E-0004; BETA^[8]:= 5.71787281789705E-0004; BETA^[9]:= 4.69431007606482E-0004; BETA^[10]:= 3.93232835462917E-0004; BETA^[11]:= 3.34818889318298E-0004; BETA^[12]:= 2.88952148495752E-0004; BETA^[13]:= 2.52211615549573E-0004; BETA^[14]:= 2.22280580798883E-0004; BETA^[15]:= 1.97541838033063E-0004; BETA^[16]:= 1.76836855019718E-0004; BETA^[17]:= 1.59316899661821E-0004; BETA^[18]:= 1.44347930197334E-0004; BETA^[19]:= 1.31448068119965E-0004; BETA^[20]:= 1.20245444949303E-0004; BETA^[21]:= 1.10449144504599E-0004; BETA^[22]:= 1.01828770740567E-0004; BETA^[23]:= 9.41998224204238E-0005; BETA^[24]:= 8.74130545753834E-0005; BETA^[25]:= 8.13466262162801E-0005; BETA^[26]:= 7.59002269646219E-0005; BETA^[27]:= 7.09906300634154E-0005; BETA^[28]:= 6.65482874842468E-0005; BETA^[29]:= 6.25146958969275E-0005; BETA^[30]:= 5.88403394426252E-0005; BETA^[31]:=-1.49282953213429E-0003; BETA^[32]:=-8.78204709546389E-0004; BETA^[33]:=-5.02916549572035E-0004; BETA^[34]:=-2.94822138512746E-0004; BETA^[35]:=-1.75463996970783E-0004; BETA^[36]:=-1.04008550460816E-0004; BETA^[37]:=-5.96141953046458E-0005; BETA^[38]:=-3.12038929076098E-0005; BETA^[39]:=-1.26089735980230E-0005; BETA^[40]:=-2.42892608575730E-0007; BETA^[41]:= 8.05996165414274E-0006; BETA^[42]:= 1.36507009262147E-0005; BETA^[43]:= 1.73964125472926E-0005; BETA^[44]:= 1.98672978842134E-0005; BETA^[45]:= 2.14463263790823E-0005; BETA^[46]:= 2.23954659232457E-0005; BETA^[47]:= 2.28967783814713E-0005; BETA^[48]:= 2.30785389811178E-0005; BETA^[49]:= 2.30321976080909E-0005; BETA^[50]:= 2.28236073720349E-0005; BETA^[51]:= 2.25005881105292E-0005; BETA^[52]:= 2.20981015361991E-0005; BETA^[53]:= 2.16418427448104E-0005; BETA^[54]:= 2.11507649256221E-0005; BETA^[55]:= 2.06388749782171E-0005; BETA^[56]:= 2.01165241997082E-0005; BETA^[57]:= 1.95913450141179E-0005; BETA^[58]:= 1.90689367910437E-0005; BETA^[59]:= 1.85533719641637E-0005; BETA^[60]:= 1.80475722259674E-0005; BETA^[61]:= 5.52213076721293E-0004; BETA^[62]:= 4.47932581552385E-0004; BETA^[63]:= 2.79520653992021E-0004; BETA^[64]:= 1.52468156198447E-0004; BETA^[65]:= 6.93271105657044E-0005; BETA^[66]:= 1.76258683069991E-0005; BETA^[67]:=-1.35744996343269E-0005; BETA^[68]:=-3.17972413350427E-0005; BETA^[69]:=-4.18861861696693E-0005; BETA^[70]:=-4.69004889379141E-0005; BETA^[71]:=-4.87665447413787E-0005; BETA^[72]:=-4.87010031186735E-0005; BETA^[73]:=-4.74755620890087E-0005; BETA^[74]:=-4.55813058138628E-0005; BETA^[75]:=-4.33309644511266E-0005; BETA^[76]:=-4.09230193157750E-0005; BETA^[77]:=-3.84822638603221E-0005; BETA^[78]:=-3.60857167535411E-0005; BETA^[79]:=-3.37793306123367E-0005; BETA^[80]:=-3.15888560772110E-0005; BETA^[81]:=-2.95269561750807E-0005; BETA^[82]:=-2.75978914828336E-0005; BETA^[83]:=-2.58006174666884E-0005; BETA^[84]:=-2.41308356761280E-0005; BETA^[85]:=-2.25823509518346E-0005; BETA^[86]:=-2.11479656768913E-0005; BETA^[87]:=-1.98200638885295E-0005; BETA^[88]:=-1.85909870801065E-0005; BETA^[89]:=-1.74532699844210E-0005; BETA^[90]:=-1.63997823854498E-0005; BETA^[91]:=-4.74617796559960E-0004; BETA^[92]:=-4.77864567147322E-0004; BETA^[93]:=-3.20390228067038E-0004; BETA^[94]:=-1.61105016119962E-0004; BETA^[95]:=-4.25778101285435E-0005; BETA^[96]:= 3.44571294294967E-0005; BETA^[97]:= 7.97092684075675E-0005; BETA^[98]:= 1.03138236708272E-0004; BETA^[99]:= 1.12466775262204E-0004; BETA^[100]:= 1.13103642108481E-0004; BETA^[101]:= 1.08651634848774E-0004; BETA^[102]:= 1.01437951597662E-0004; BETA^[103]:= 9.29298396593364E-0005; BETA^[104]:= 8.40293133016090E-0005; BETA^[105]:= 7.52727991349134E-0005; BETA^[106]:= 6.69632521975731E-0005; BETA^[107]:= 5.92564547323195E-0005; BETA^[108]:= 5.22169308826976E-0005; BETA^[109]:= 4.58539485165361E-0005; BETA^[110]:= 4.01445513891487E-0005; BETA^[111]:= 3.50481730031328E-0005; BETA^[112]:= 3.05157995034347E-0005; BETA^[113]:= 2.64956119950516E-0005; BETA^[114]:= 2.29363633690998E-0005; BETA^[115]:= 1.97893056664022E-0005; BETA^[116]:= 1.70091984636413E-0005; BETA^[117]:= 1.45547428261524E-0005; BETA^[118]:= 1.23886640995878E-0005; BETA^[119]:= 1.04775876076583E-0005; BETA^[120]:= 8.79179954978479E-0006; BETA^[121]:= 7.36465810572578E-0004; BETA^[122]:= 8.72790805146194E-0004; BETA^[123]:= 6.22614862573135E-0004; BETA^[124]:= 2.85998154194304E-0004; BETA^[125]:= 3.84737672879366E-0006; BETA^[126]:=-1.87906003636972E-0004; BETA^[127]:=-2.97603646594555E-0004; BETA^[128]:=-3.45998126832656E-0004; BETA^[129]:=-3.53382470916038E-0004; BETA^[130]:=-3.35715635775049E-0004; BETA^[131]:=-3.04321124789040E-0004; BETA^[132]:=-2.66722723047613E-0004; BETA^[133]:=-2.27654214122820E-0004; BETA^[134]:=-1.89922611854562E-0004; BETA^[135]:=-1.55058918599094E-0004; BETA^[136]:=-1.23778240761874E-0004; BETA^[137]:=-9.62926147717644E-0005; BETA^[138]:=-7.25178327714425E-0005; BETA^[139]:=-5.22070028895634E-0005; BETA^[140]:=-3.50347750511901E-0005; BETA^[141]:=-2.06489761035552E-0005; BETA^[142]:=-8.70106096849767E-0006; BETA^[143]:= 1.13698686675100E-0006; BETA^[144]:= 9.16426474122779E-0006; BETA^[145]:= 1.56477785428873E-0005; BETA^[146]:= 2.08223629482467E-0005; BETA^[147]:= 2.48923381004595E-0005; BETA^[148]:= 2.80340509574146E-0005; BETA^[149]:= 3.03987774629862E-0005; BETA^[150]:= 3.21156731406701E-0005; BETA^[151]:=-1.80182191963886E-0003; BETA^[152]:=-2.43402962938043E-0003; BETA^[153]:=-1.83422663549857E-0003; BETA^[154]:=-7.62204596354010E-0004; BETA^[155]:= 2.39079475256927E-0004; BETA^[156]:= 9.49266117176881E-0004; BETA^[157]:= 1.34467449701540E-0003; BETA^[158]:= 1.48457495259449E-0003; BETA^[159]:= 1.44732339830618E-0003; BETA^[160]:= 1.30268261285657E-0003; BETA^[161]:= 1.10351597375643E-0003; BETA^[162]:= 8.86047440419792E-0004; BETA^[163]:= 6.73073208165665E-0004; BETA^[164]:= 4.77603872856582E-0004; BETA^[165]:= 3.05991926358789E-0004; BETA^[166]:= 1.60315694594722E-0004; BETA^[167]:= 4.00749555270613E-0005; BETA^[168]:=-5.66607461635252E-0005; BETA^[169]:=-1.32506186772983E-0004; BETA^[170]:=-1.90296187989614E-0004; BETA^[171]:=-2.32811450376937E-0004; BETA^[172]:=-2.62628811464669E-0004; BETA^[173]:=-2.82050469867599E-0004; BETA^[174]:=-2.93081563192861E-0004; BETA^[175]:=-2.97435962176317E-0004; BETA^[176]:=-2.96557334239348E-0004; BETA^[177]:=-2.91647363312091E-0004; BETA^[178]:=-2.83696203837734E-0004; BETA^[179]:=-2.73512317095673E-0004; BETA^[180]:=-2.61750155806769E-0004; BETA^[181]:= 6.38585891212051E-0003; BETA^[182]:= 9.62374215806378E-0003; BETA^[183]:= 7.61878061207001E-0003; BETA^[184]:= 2.83219055545628E-0003; BETA^[185]:=-2.09841352012720E-0003; BETA^[186]:=-5.73826764216626E-0003; BETA^[187]:=-7.70804244495415E-0003; BETA^[188]:=-8.21011692264844E-0003; BETA^[189]:=-7.65824520346905E-0003; BETA^[190]:=-6.47209729391045E-0003; BETA^[191]:=-4.99132412004966E-0003; BETA^[192]:=-3.45612289713133E-0003; BETA^[193]:=-2.01785580014171E-0003; BETA^[194]:=-7.59430686781961E-0004; BETA^[195]:= 2.84173631523859E-0004; BETA^[196]:= 1.10891667586337E-0003; BETA^[197]:= 1.72901493872729E-0003; BETA^[198]:= 2.16812590802685E-0003; BETA^[199]:= 2.45357710494540E-0003; BETA^[200]:= 2.61281821058335E-0003; BETA^[201]:= 2.67141039656277E-0003; BETA^[202]:= 2.65203073395980E-0003; BETA^[203]:= 2.57411652877287E-0003; BETA^[204]:= 2.45389126236094E-0003; BETA^[205]:= 2.30460058071796E-0003; BETA^[206]:= 2.13684837686713E-0003; BETA^[207]:= 1.95896528478871E-0003; BETA^[208]:= 1.77737008679454E-0003; BETA^[209]:= 1.59690280765839E-0003; BETA^[210]:= 1.42111975664439E-0003; GAMA^[1]:= 6.29960524947437E-0001; GAMA^[2]:= 2.51984209978975E-0001; GAMA^[3]:= 1.54790300415656E-0001; GAMA^[4]:= 1.10713062416159E-0001; GAMA^[5]:= 8.57309395527395E-0002; GAMA^[6]:= 6.97161316958684E-0002; GAMA^[7]:= 5.86085671893714E-0002; GAMA^[8]:= 5.04698873536311E-0002; GAMA^[9]:= 4.42600580689155E-0002; GAMA^[10]:= 3.93720661543510E-0002; GAMA^[11]:= 3.54283195924455E-0002; GAMA^[12]:= 3.21818857502098E-0002; GAMA^[13]:= 2.94646240791158E-0002; GAMA^[14]:= 2.71581677112934E-0002; GAMA^[15]:= 2.51768272973862E-0002; GAMA^[16]:= 2.34570755306079E-0002; GAMA^[17]:= 2.19508390134907E-0002; GAMA^[18]:= 2.06210828235646E-0002; GAMA^[19]:= 1.94388240897881E-0002; GAMA^[20]:= 1.83810633800683E-0002; GAMA^[21]:= 1.74293213231963E-0002; GAMA^[22]:= 1.65685837786612E-0002; GAMA^[23]:= 1.57865285987918E-0002; GAMA^[24]:= 1.50729501494096E-0002; GAMA^[25]:= 1.44193250839955E-0002; GAMA^[26]:= 1.38184805735342E-0002; GAMA^[27]:= 1.32643378994277E-0002; GAMA^[28]:= 1.27517121970499E-0002; GAMA^[29]:= 1.22761545318763E-0002; GAMA^[30]:= 1.18338262398482E-0002; EX1:=3.33333333333333333E-01; EX2:=6.66666666666666667E-01; HPI:=1.57079632679489662E+00; GPI:=3.14159265358979324E+00; THPI:=4.71238898038468986; ZEROR:=0.0; ZEROI:=0.0; CONER:=1.0; CONEI:=0.0; RFNU := 1.0/FNU; ZBR := ZR*RFNU; ZBI := ZI*RFNU; RFNU2 := RFNU*RFNU; {----------------------------------------------------------------------- ! COMPUTE IN THE FOURTH QUADRANT !----------------------------------------------------------------------} FN13 := Power(FNU,EX1); FN23 := FN13*FN13; RFN13 := 1.0/FN13; W2R := CONER - ZBR*ZBR + ZBI*ZBI; W2I := CONEI - ZBR*ZBI - ZBR*ZBI; AW2 := ZABS(W2R,W2I); IF AW2 > 0.25 THEN GOTO 130; {----------------------------------------------------------------------- ! POWER SERIES FOR CABS(W2).LE.0.25D0 !----------------------------------------------------------------------} K := 1; PR^[1] := CONER; PI^[1] := CONEI; SUMAR := GAMA^[1]; SUMAI := ZEROI; AP^[1] := 1.0; IF AW2 < TOL THEN GOTO 20; For K:=2 to 30 do begin PR^[K] := PR^[K-1]*W2R - PI^[K-1]*W2I; PI^[K] := PR^[K-1]*W2I + PI^[K-1]*W2R; SUMAR := SUMAR + PR^[K]*GAMA^[K]; SUMAI := SUMAI + PI^[K]*GAMA^[K]; AP^[K] := AP^[K-1]*AW2; IF AP^[K] < TOL THEN GOTO 20 end; K := 30; 20: KMAX := K; ZETAR := W2R*SUMAR - W2I*SUMAI; ZETAI := W2R*SUMAI + W2I*SUMAR; ARGR := ZETAR*FN23; ARGI := ZETAI*FN23; ZSQRT(SUMAR, SUMAI, ZAR, ZAI); ZSQRT(W2R, W2I, STR, STI); ZETA2R := STR*FNU; ZETA2I := STI*FNU; STR := CONER + EX2*(ZETAR*ZAR-ZETAI*ZAI); STI := CONEI + EX2*(ZETAR*ZAI+ZETAI*ZAR); ZETA1R := STR*ZETA2R - STI*ZETA2I; ZETA1I := STR*ZETA2I + STI*ZETA2R; ZAR := ZAR + ZAR; ZAI := ZAI + ZAI; ZSQRT(ZAR, ZAI, STR, STI); PHIR := STR*RFN13; PHII := STI*RFN13; IF IPMTR = 1 THEN GOTO 120; {----------------------------------------------------------------------- ! SUM SERIES FOR ASUM AND BSUM !----------------------------------------------------------------------} SUMBR := ZEROR; SUMBI := ZEROI; For K:=1 to KMAX do begin SUMBR := SUMBR + PR^[K]*BETA^[K]; SUMBI := SUMBI + PI^[K]*BETA^[K]; end; ASUMR := ZEROR; ASUMI := ZEROI; BSUMR := SUMBR; BSUMI := SUMBI; L1 := 0; L2 := 30; BTOL := TOL*(ABS(BSUMR)+ABS(BSUMI)); ATOL := TOL; PP := 1.0; IAS := 0; IBS := 0; IF RFNU2 < TOL THEN GOTO 110; For IS:=2 to 7 do begin ATOL := ATOL/RFNU2; PP := PP*RFNU2; IF IAS = 1 THEN GOTO 60; SUMAR := ZEROR; SUMAI := ZEROI; For K:=1 to KMAX do begin M := L1 + K; SUMAR := SUMAR + PR^[K]*ALFA^[M]; SUMAI := SUMAI + PI^[K]*ALFA^[M]; IF AP^[K] < ATOL THEN GOTO 50; end; 50: ASUMR := ASUMR + SUMAR*PP; ASUMI := ASUMI + SUMAI*PP; IF PP < TOL THEN IAS := 1; 60: IF IBS = 1 THEN GOTO 90; SUMBR := ZEROR; SUMBI := ZEROI; For K:=1 to KMAX do begin M := L2 + K; SUMBR := SUMBR + PR^[K]*BETA^[M]; SUMBI := SUMBI + PI^[K]*BETA^[M]; IF AP^[K] < ATOL THEN GOTO 80; end; 80: BSUMR := BSUMR + SUMBR*PP; BSUMI := BSUMI + SUMBI*PP; IF PP < BTOL THEN IBS := 1; 90: IF (IAS = 1) AND (IBS = 1) THEN GOTO 110; L1 := L1 + 30; L2 := L2 + 30 end; {IS loop} 110: ASUMR := ASUMR + CONER; PP := RFNU*RFN13; BSUMR := BSUMR*PP; BSUMI := BSUMI*PP; 120: GOTO RETURN; {----------------------------------------------------------------------- ! CABS(W2).GT.0.25D0 !----------------------------------------------------------------------} 130: ZSQRT(W2R, W2I, WR, WI); IF WR < 0.0 THEN WR := 0.0; IF WI < 0.0 THEN WI := 0.0; STR := CONER + WR; STI := WI; ZDIV(STR, STI, ZBR, ZBI, ZAR, ZAI); ZLOG(ZAR, ZAI, ZCR, ZCI, IDUM); IF ZCI < 0.0 THEN ZCI := 0.0; IF ZCI > HPI THEN ZCI := HPI; IF ZCR < 0.0 THEN ZCR := 0.0; ZTHR := (ZCR-WR)*1.5; ZTHI := (ZCI-WI)*1.5; ZETA1R := ZCR*FNU; ZETA1I := ZCI*FNU; ZETA2R := WR*FNU; ZETA2I := WI*FNU; AZTH := ZABS(ZTHR,ZTHI); ANG := THPI; IF (ZTHR >= 0.0) AND (ZTHI < 0.0) THEN GOTO 140; ANG := HPI; IF ZTHR = 0.0 THEN GOTO 140; ANG := ArcTan(ZTHI/ZTHR); IF ZTHR < 0.0 THEN ANG := ANG + GPI; 140: PP := Power(AZTH,EX2); ANG := ANG*EX2; ZETAR := PP*COS(ANG); ZETAI := PP*SIN(ANG); IF ZETAI < 0.0 THEN ZETAI := 0.0; ARGR := ZETAR*FN23; ARGI := ZETAI*FN23; ZDIV(ZTHR, ZTHI, ZETAR, ZETAI, RTZTR, RTZTI); ZDIV(RTZTR, RTZTI, WR, WI, ZAR, ZAI); TZAR := ZAR + ZAR; TZAI := ZAI + ZAI; ZSQRT(TZAR, TZAI, STR, STI); PHIR := STR*RFN13; PHII := STI*RFN13; IF IPMTR = 1 THEN GOTO 120; RAW := 1.0/SQRT(AW2); STR := WR*RAW; STI := -WI*RAW; TFNR := STR*RFNU*RAW; TFNI := STI*RFNU*RAW; RAZTH := 1.0/AZTH; STR := ZTHR*RAZTH; STI := -ZTHI*RAZTH; RZTHR := STR*RAZTH*RFNU; RZTHI := STI*RAZTH*RFNU; ZCR := RZTHR*AR[2]; ZCI := RZTHI*AR[2]; RAW2 := 1.0/AW2; STR := W2R*RAW2; STI := -W2I*RAW2; T2R := STR*RAW2; T2I := STI*RAW2; STR := T2R*C^[2] + C^[3]; STI := T2I*C^[2]; UPR[2] := STR*TFNR - STI*TFNI; UPI[2] := STR*TFNI + STI*TFNR; BSUMR := UPR[2] + ZCR; BSUMI := UPI[2] + ZCI; ASUMR := ZEROR; ASUMI := ZEROI; IF RFNU < TOL THEN GOTO 220; PRZTHR := RZTHR; PRZTHI := RZTHI; PTFNR := TFNR; PTFNI := TFNI; UPR[1] := CONER; UPI[1] := CONEI; PP := 1.0; BTOL := TOL*(ABS(BSUMR)+ABS(BSUMI)); KS := 0; KP1 := 2; L := 3; IAS := 0; IBS := 0; LR:=2; While LR<=12 do begin LRP1 := LR + 1; {----------------------------------------------------------------------- ! COMPUTE TWO ADDITIONAL CR, DR, AND UP FOR TWO MORE TERMS IN ! NEXT SUMA AND SUMB !----------------------------------------------------------------------} For K:=LR to LRP1 do begin KS := KS + 1; KP1 := KP1 + 1; L := L + 1; ZAR := C^[L]; ZAI := ZEROI; For J:=2 to KP1 do begin L := L + 1; STR := ZAR*T2R - T2I*ZAI + C^[L]; ZAI := ZAR*T2I + ZAI*T2R; ZAR := STR end; STR := PTFNR*TFNR - PTFNI*TFNI; PTFNI := PTFNR*TFNI + PTFNI*TFNR; PTFNR := STR; UPR[KP1] := PTFNR*ZAR - PTFNI*ZAI; UPI[KP1] := PTFNI*ZAR + PTFNR*ZAI; CRR[KS] := PRZTHR*BR[KS+1]; CRI[KS] := PRZTHI*BR[KS+1]; STR := PRZTHR*RZTHR - PRZTHI*RZTHI; PRZTHI := PRZTHR*RZTHI + PRZTHI*RZTHR; PRZTHR := STR; DRR[KS] := PRZTHR*AR[KS+2]; DRI[KS] := PRZTHI*AR[KS+2] end; {K loop} PP := PP*RFNU2; IF IAS = 1 THEN GOTO 180; SUMAR := UPR[LRP1]; SUMAI := UPI[LRP1]; JU := LRP1; For JR:=1 to LR do begin JU := JU - 1; SUMAR := SUMAR + CRR[JR]*UPR[JU] - CRI[JR]*UPI[JU]; SUMAI := SUMAI + CRR[JR]*UPI[JU] + CRI[JR]*UPR[JU] end; ASUMR := ASUMR + SUMAR; ASUMI := ASUMI + SUMAI; TEST := ABS(SUMAR) + ABS(SUMAI); IF (PP < TOL) AND (TEST < TOL) THEN IAS := 1; 180: IF IBS = 1 THEN GOTO 200; SUMBR := UPR[LR+2] + UPR[LRP1]*ZCR - UPI[LRP1]*ZCI; SUMBI := UPI[LR+2] + UPR[LRP1]*ZCI + UPI[LRP1]*ZCR; JU := LRP1; For JR:=1 to LR do begin JU := JU - 1; SUMBR := SUMBR + DRR[JR]*UPR[JU] - DRI[JR]*UPI[JU]; SUMBI := SUMBI + DRR[JR]*UPI[JU] + DRI[JR]*UPR[JU] end; BSUMR := BSUMR + SUMBR; BSUMI := BSUMI + SUMBI; TEST := ABS(SUMBR) + ABS(SUMBI); IF (PP < BTOL) AND (TEST < BTOL) THEN IBS := 1; 200: IF (IAS = 1) AND (IBS = 1) THEN GOTO 220; Inc(LR,2) end; {while LR} 220: ASUMR := ASUMR + CONER; STR := -BSUMR*RFN13; STI := -BSUMI*RFN13; ZDIV(STR, STI, RTZTR, RTZTI, BSUMR, BSUMI); GOTO 120; Return: End; {ZUNHJ} Procedure ZUOIK(ZR, ZI, FNU:Double; KODE, IKFLG, N: Integer; Var YR, YI:VEC; Var NUF:Integer; Var TOL, ELIM, ALIM:Double); {***BEGIN PROLOGUE ZUOIK !***REFER TO ZBESI,ZBESK,ZBESH ! ! ZUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC ! EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM ! (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW ! WHERE ALIM.LT.ELIM. IF THE MAGNITUDE, BASED ON THE LEADING ! EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN ! THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER ! MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE ! EXP(-ELIM):=SMALLEST MACHINE NUMBER*1E+03 AND EXP(-ALIM)= ! EXP(-ELIM)/TOL. ! ! IKFLG=1 MEANS THE I SEQUENCE IS TESTED ! =2 MEANS THE K SEQUENCE IS TESTED ! NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE ! =-1 MEANS AN OVERFLOW WOULD OCCUR ! IKFLG=1 AND NUF.GT.0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO ! THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE ! IKFLG=2 AND NUF.EQ.N MEANS ALL Y VALUES WERE SET TO ZERO ! IKFLG=2 AND 0.LT.NUF.LT.N NOT CONSIDERED. Y MUST BE SET BY ! ANOTHER ROUTINE. ! !***ROUTINES CALLED ZUCHK,ZUNHJ,ZUNIK,D1MACH,ZABS,ZLOG !***END PROLOGUE ZUOIK ! COMPLEX ARG,ASUM,BSUM,CWRK,CZ,CZERO,PHI,SUM,Y,Z,ZB,ZETA1,ZETA2,ZN, ! *ZR } Label 10,20,30,40,50,60,70,80,90,110,120,130,140,150,160,170,180,190,200,210,Return; Var AARG, AIC, APHI, ARGI, ARGR, ASUMI, ASUMR, ASCLE, AX, AY, BSUMI, BSUMR, CZI, CZR, FNN, GNN, GNU, PHII, PHIR, RCZ, STR, STI, SUMI, SUMR, ZBI, ZBR, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZNI, ZNR, ZRI, ZRR: Double; I, IDUM, IFORM, INIT, NN, NW: Integer; CWRKR,CWRKI: VEC16; Begin ZEROR:=0.0; ZEROI:=0.0; AIC:=1.265512123484645396; NUF := 0; NN := N; ZRR := ZR; ZRI := ZI; IF ZR >= 0.0 THEN GOTO 10; ZRR := -ZR; ZRI := -ZI; 10: ZBR := ZRR; ZBI := ZRI; AX := ABS(ZR)*1.7321; AY := ABS(ZI); IFORM := 1; IF AY > AX THEN IFORM := 2; GNU := DMAX(FNU,1.0); IF IKFLG = 1 THEN GOTO 20; FNN := 1.0*NN; GNN := FNU + FNN - 1.0; GNU := DMAX(GNN,FNN); {----------------------------------------------------------------------- ! ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE ! REAL PARTS OF ZETA1, ZETA2 AND ZB. NO ATTEMPT IS MADE TO GET ! THE SIGN OF THE IMAGINARY PART CORRECT. !----------------------------------------------------------------------} 20: IF IFORM = 2 THEN GOTO 30; INIT := 0; ZUNIK(ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, PHIR, PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI); CZR := -ZETA1R + ZETA2R; CZI := -ZETA1I + ZETA2I; GOTO 50; 30: ZNR := ZRI; ZNI := -ZRR; IF ZI > 0.0 THEN GOTO 40; ZNR := -ZNR; 40: ZUNHJ(ZNR, ZNI, GNU, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI); CZR := -ZETA1R + ZETA2R; CZI := -ZETA1I + ZETA2I; AARG := ZABS(ARGR,ARGI); 50: IF KODE = 1 THEN GOTO 60; CZR := CZR - ZBR; CZI := CZI - ZBI; 60: IF IKFLG = 1 THEN GOTO 70; CZR := -CZR; CZI := -CZI; 70: APHI := ZABS(PHIR,PHII); RCZ := CZR; {----------------------------------------------------------------------- ! OVERFLOW TEST !----------------------------------------------------------------------} IF RCZ > ELIM THEN GOTO 210; IF RCZ < ALIM THEN GOTO 80; RCZ := RCZ + Ln(APHI); IF IFORM = 2 THEN RCZ := RCZ - 0.25*Ln(AARG) - AIC; IF RCZ > ELIM THEN GOTO 210; GOTO 130; {----------------------------------------------------------------------- ! UNDERFLOW TEST !----------------------------------------------------------------------} 80: IF RCZ < -ELIM THEN GOTO 90; IF RCZ > -ALIM THEN GOTO 130; RCZ := RCZ + Ln(APHI); IF IFORM = 2 THEN RCZ := RCZ - 0.25*Ln(AARG) - AIC; IF RCZ > -ELIM THEN GOTO 110; 90: For I:=1 to NN do begin YR[I] := ZEROR; YI[I] := ZEROI end; NUF := NN; GOTO RETURN; 110: ASCLE := 1000*D1MACH(1)/TOL; ZLOG(PHIR, PHII, STR, STI, IDUM); CZR := CZR + STR; CZI := CZI + STI; IF IFORM = 1 THEN GOTO 120; ZLOG(ARGR, ARGI, STR, STI, IDUM); CZR := CZR - 0.25*STR - AIC; CZI := CZI - 0.25*STI; 120: AX := EXP(RCZ)/TOL; AY := CZI; CZR := AX*COS(AY); CZI := AX*SIN(AY); ZUCHK(CZR, CZI, NW, ASCLE, TOL); IF NW <> 0 THEN GOTO 90; 130: IF IKFLG = 2 THEN GOTO RETURN; IF N = 1 THEN GOTO RETURN; {----------------------------------------------------------------------- ! SET UNDERFLOWS ON I SEQUENCE !----------------------------------------------------------------------} 140: GNU := FNU + 1.0*(NN-1); IF IFORM = 2 THEN GOTO 150; INIT := 0; ZUNIK(ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, PHIR, PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI); CZR := -ZETA1R + ZETA2R; CZI := -ZETA1I + ZETA2I; GOTO 160; 150: ZUNHJ(ZNR, ZNI, GNU, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI); CZR := -ZETA1R + ZETA2R; CZI := -ZETA1I + ZETA2I; AARG := ZABS(ARGR,ARGI); 160: IF KODE = 1 THEN GOTO 170; CZR := CZR - ZBR; CZI := CZI - ZBI; 170: APHI := ZABS(PHIR,PHII); RCZ := CZR; IF RCZ < -ELIM THEN GOTO 180; IF RCZ > -ALIM THEN GOTO RETURN; RCZ := RCZ + Ln(APHI); IF IFORM = 2 THEN RCZ := RCZ - 0.25*Ln(AARG) - AIC; IF RCZ > -ELIM THEN GOTO 190; 180: YR[NN] := ZEROR; YI[NN] := ZEROI; NN := NN - 1; NUF := NUF + 1; IF NN = 0 THEN GOTO RETURN; GOTO 140; 190: ASCLE := 1000*D1MACH(1)/TOL; ZLOG(PHIR, PHII, STR, STI, IDUM); CZR := CZR + STR; CZI := CZI + STI; IF IFORM = 1 THEN GOTO 200; ZLOG(ARGR, ARGI, STR, STI, IDUM); CZR := CZR - 0.25*STR - AIC; CZI := CZI - 0.25*STI; 200: AX := EXP(RCZ)/TOL; AY := CZI; CZR := AX*COS(AY); CZI := AX*SIN(AY); ZUCHK(CZR, CZI, NW, ASCLE, TOL); IF NW <> 0 THEN GOTO 180; GOTO RETURN; 210: NUF := -1; Return: End; {ZUOIK} Procedure ZS1S2(Var ZRR, ZRI, S1R, S1I, S2R, S2I:double; NZ:Integer; Var ASCLE, ALIM: double; Var IUF:Integer); {***BEGIN PROLOGUE ZS1S2 !***REFER TO ZBESK,ZAIRY ! ! ZS1S2 TESTS FOR A POSSIBLE UNDERFLOW RESULTING FROM THE ! ADDITION OF THE I AND K FUNCTIONS IN THE ANALYTIC CON- ! TINUATION FORMULA WHERE S1=K FUNCTION AND S2=I FUNCTION. ! ON KODE=1 THE I AND K FUNCTIONS ARE DIFFERENT ORDERS OF ! MAGNITUDE, BUT FOR KODE=2 THEY CAN BE OF THE SAME ORDER ! OF MAGNITUDE AND THE MAXIMUM MUST BE AT LEAST ONE ! PRECISION ABOVE THE UNDERFLOW LIMIT. ! !***ROUTINES CALLED ZABS,ZEXP,ZLOG !***END PROLOGUE ZS1S2 ! COMPLEX CZERO,C1,S1,S1D,S2,ZR } Label 10, Return; Var AA, ALN, AS1, AS2, C1I, C1R, S1DI, S1DR, ZEROI, ZEROR: Double; IDUM: Integer; Begin ZEROR:=0.0; ZEROI:=0.0; NZ := 0; AS1 := ZABS(S1R,S1I); AS2 := ZABS(S2R,S2I); IF (S1R = 0.0) AND (S1I = 0.0) THEN GOTO 10; IF AS1 = 0.0 THEN GOTO 10; ALN := -ZRR - ZRR + Ln(AS1); S1DR := S1R; S1DI := S1I; S1R := ZEROR; S1I := ZEROI; AS1 := ZEROR; IF ALN < -ALIM THEN GOTO 10; ZLOG(S1DR, S1DI, C1R, C1I, IDUM); C1R := C1R - ZRR - ZRR; C1I := C1I - ZRI - ZRI; ZEXP(C1R, C1I, S1R, S1I); AS1 := ZABS(S1R,S1I); IUF := IUF + 1; 10: AA := DMAX(AS1,AS2); IF AA > ASCLE THEN GOTO RETURN; S1R := ZEROR; S1I := ZEROI; S2R := ZEROR; S2I := ZEROI; NZ := 1; IUF := 0; Return:End; {ZS1S2} Procedure ZSERI(ZR, ZI, FNU:Double; KODE, N: Integer; Var YR, YI: VEC; Var NZ: Integer; Var TOL, ELIM, ALIM:Double); {***BEGIN PROLOGUE ZSERI !***REFER TO ZBESI,ZBESK ! ! ZSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z) >= 0 BY ! MEANS OF THE POWER SERIES FOR LARGE CABS(Z) IN THE ! REGION CABS(Z) <= 2*SQRT(FNU+1). NZ=0 IS A NORMAL RETURN. ! NZ > 0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO ! DUE TO UNDERFLOW. NZ < 0 MEANS UNDERFLOW OCCURRED, BUT THE ! CONDITION CABS(Z) <= 2*SQRT(FNU+1) WAS VIOLATED AND THE ! COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ). ! !***ROUTINES CALLED DGAMLN,D1MACH,ZUCHK,ZABS,ZDIV,ZLOG,ZMLT !***END PROLOGUE ZSERI ! COMPLEX AK1,CK,COEF,CONE,CRSC,CSCL,CZ,CZERO,HZ,RZ,S1,S2,Y,Z } Label 10,20,30,40,50,60,70,80,90,100,120,140,150,160,170,190, Return; Var AA, ACZ, AK, AK1I, AK1R, ARM, ASCLE, ATOL,AZ, CKI, CKR, COEFI, COEFR, CONEI, CONER, CRSCR, CZI, CZR, DFNU, FNUP, HZI, HZR, RAZ, RS, RTR1, RZI, RZR, S, SS, STI, STR, S1I, S1R, S2I, S2R, ZEROI, ZEROR: Double; I, IB, IDUM, IFLAG, IL, K, L, M, NN, NW: Integer; WR, WI: VEC; Begin ZEROR:=0.0; ZEROI:=0.0; CONER:=1.0; CONEI:=0.0; NZ := 0; AZ := ZABS(ZR,ZI); IF AZ = 0.0 THEN GOTO 160; ARM := 1000*D1MACH(1); RTR1 := SQRT(ARM); CRSCR := 1.0; IFLAG := 0; IF AZ < ARM THEN GOTO 150; HZR := 0.5*ZR; HZI := 0.5*ZI; CZR := ZEROR; CZI := ZEROI; IF AZ <= RTR1 THEN GOTO 10; ZMLT(HZR, HZI, HZR, HZI, CZR, CZI); 10: ACZ := ZABS(CZR,CZI); NN := N; ZLOG(HZR, HZI, CKR, CKI, IDUM); 20: DFNU := FNU + 1.0*(NN-1); FNUP := DFNU + 1.0; {----------------------------------------------------------------------- ! UNDERFLOW TEST !----------------------------------------------------------------------} AK1R := CKR*DFNU; AK1I := CKI*DFNU; AK := DGAMLN(FNUP,IDUM); AK1R := AK1R - AK; IF KODE = 2 THEN AK1R := AK1R - ZR; IF AK1R > -ELIM THEN GOTO 40; 30: NZ := NZ + 1; YR[NN] := ZEROR; YI[NN] := ZEROI; IF ACZ > DFNU THEN GOTO 190; NN := NN - 1; IF NN = 0 THEN GOTO RETURN; GOTO 20; 40: IF AK1R > -ALIM THEN GOTO 50; IFLAG := 1; SS := 1.0/TOL; CRSCR := TOL; ASCLE := ARM*SS; 50: AA := EXP(AK1R); IF IFLAG = 1 THEN AA := AA*SS; COEFR := AA*COS(AK1I); COEFI := AA*SIN(AK1I); ATOL := TOL*ACZ/FNUP; IL := IMIN(2,NN); For I:=1 to IL do begin DFNU := FNU + 1.0*(NN-I); FNUP := DFNU + 1.0; S1R := CONER; S1I := CONEI; IF ACZ < TOL*FNUP THEN GOTO 70; AK1R := CONER; AK1I := CONEI; AK := FNUP + 2.0; S := FNUP; AA := 2.0; 60: RS := 1.0/S; STR := AK1R*CZR - AK1I*CZI; STI := AK1R*CZI + AK1I*CZR; AK1R := STR*RS; AK1I := STI*RS; S1R := S1R + AK1R; S1I := S1I + AK1I; S := S + AK; AK := AK + 2.0; AA := AA*ACZ*RS; IF AA > ATOL THEN GOTO 60; 70: S2R := S1R*COEFR - S1I*COEFI; S2I := S1R*COEFI + S1I*COEFR; WR[I] := S2R; WI[I] := S2I; IF IFLAG = 0 THEN GOTO 80; ZUCHK(S2R, S2I, NW, ASCLE, TOL); IF NW <> 0 THEN GOTO 30; 80: M := NN - I + 1; YR[M] := S2R*CRSCR; YI[M] := S2I*CRSCR; IF I = IL THEN GOTO 90; ZDIV(COEFR, COEFI, HZR, HZI, STR, STI); COEFR := STR*DFNU; COEFI := STI*DFNU; 90: end; IF NN <= 2 THEN GOTO RETURN; K := NN - 2; AK := 1.0*K; RAZ := 1.0/AZ; STR := ZR*RAZ; STI := -ZI*RAZ; RZR := (STR+STR)*RAZ; RZI := (STI+STI)*RAZ; IF IFLAG = 1 THEN GOTO 120; IB := 3; 100: For I:=IB to NN do begin YR[K] := (AK+FNU)*(RZR*YR[K+1]-RZI*YI[K+1]) + YR[K+2]; YI[K] := (AK+FNU)*(RZR*YI[K+1]+RZI*YR[K+1]) + YI[K+2]; AK := AK - 1.0; K := K - 1 end; GOTO RETURN; {----------------------------------------------------------------------- ! RECUR BACKWARD WITH SCALED VALUES !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! EXP(-ALIM):=EXP(-ELIM)/TOL:=APPROX. ONE PRECISION ABOVE THE ! UNDERFLOW LIMIT = ASCLE := D1MACH(1)*SS*1000. !----------------------------------------------------------------------} 120: S1R := WR[1]; S1I := WI[1]; S2R := WR[2]; S2I := WI[2]; For L:=3 to NN do begin CKR := S2R; CKI := S2I; S2R := S1R + (AK+FNU)*(RZR*CKR-RZI*CKI); S2I := S1I + (AK+FNU)*(RZR*CKI+RZI*CKR); S1R := CKR; S1I := CKI; CKR := S2R*CRSCR; CKI := S2I*CRSCR; YR[K] := CKR; YI[K] := CKI; AK := AK - 1.0; K := K - 1; IF ZABS(CKR,CKI) > ASCLE THEN GOTO 140; end; GOTO RETURN; 140: IB := L + 1; IF IB > NN THEN GOTO RETURN; GOTO 100; 150: NZ := N; IF FNU = 0.0 THEN NZ := NZ - 1; 160: YR[1] := ZEROR; YI[1] := ZEROI; IF FNU <> 0.0 THEN GOTO 170; YR[1] := CONER; YI[1] := CONEI; 170: IF N = 1 THEN GOTO RETURN; For I:=2 to N do begin YR[I] := ZEROR; YI[I] := ZEROI end; GOTO RETURN; {----------------------------------------------------------------------- ! RETURN WITH NZ.LT.0 IF CABS(Z*Z/4).GT.FNU+N-NZ-1 COMPLETE ! THE CALCULATION IN CBINU WITH N:=N-IABS(NZ) !----------------------------------------------------------------------} 190: NZ := -NZ; Return:End; {ZSERI} Procedure ZASYI(ZR, ZI, FNU:Double; KODE, N: Integer; Var YR, YI: VEC; Var NZ: Integer; Var RL, TOL, ELIM, ALIM:Double); {***BEGIN PROLOGUE ZASYI !***REFER TO ZBESI,ZBESK ! ZASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z) >= 0 BY ! MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE ! REGION CABS(Z) > MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN. ! NZ < 0 INDICATES AN OVERFLOW ON KODE=1. ! !***ROUTINES CALLED D1MACH,ZABS,ZDIV,ZEXP,ZMLT,ZSQRT !***END PROLOGUE ZASYI ! COMPLEX AK1,CK,CONE,CS1,CS2,CZ,CZERO,DK,EZ,P1,RZ,S2,Y,Z } Label 10,20,30,50,60, 100,110, Return; Var AA, AEZ, AK, AK1I, AK1R, ARG, ARM, ATOL, AZ, BB, BK, CKI, CKR, CONEI, CONER, CS1I, CS1R, CS2I, CS2R, CZI, CZR, DFNU, DKI, DKR, DNU2, EZI, EZR, FDN, PI, P1I, P1R, RAZ, RTPI, RTR1, RZI, RZR, S, SGN, SQK, STI, STR, S2I, S2R, TZI, TZR, ZEROI, ZEROR: Double; I, IB, IL, INU, J, JL, K, KODED, M, NN: Integer; Begin RTPI:=0.159154943091895336; ZEROR:=0.0; ZEROI:=0.0; CONER:=1.0; CONEI:=0.0; NZ := 0; AZ := ZABS(ZR,ZI); ARM := 1000*D1MACH(1); RTR1 := SQRT(ARM); IL := IMIN(2,N); DFNU := FNU + 1.0*(N-IL); {----------------------------------------------------------------------- ! OVERFLOW TEST !----------------------------------------------------------------------} RAZ := 1.0/AZ; STR := ZR*RAZ; STI := -ZI*RAZ; AK1R := RTPI*STR*RAZ; AK1I := RTPI*STI*RAZ; ZSQRT(AK1R, AK1I, AK1R, AK1I); CZR := ZR; CZI := ZI; IF KODE <> 2 THEN GOTO 10; CZR := ZEROR; CZI := ZI; 10: IF ABS(CZR) > ELIM THEN GOTO 100; DNU2 := DFNU + DFNU; KODED := 1; IF (ABS(CZR) > ALIM) AND (N > 2) THEN GOTO 20; KODED := 0; ZEXP(CZR, CZI, STR, STI); ZMLT(AK1R, AK1I, STR, STI, AK1R, AK1I); 20: FDN := 0.0; IF DNU2 > RTR1 THEN FDN := DNU2*DNU2; EZR := ZR*8.0; EZI := ZI*8.0; {----------------------------------------------------------------------- ! WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE ! FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE ! EXPANSION FOR THE IMAGINARY PART. !----------------------------------------------------------------------} AEZ := 8.0*AZ; S := TOL/AEZ; JL := Round(RL+RL) + 2; P1R := ZEROR; P1I := ZEROI; IF ZI = 0.0 THEN GOTO 30; {----------------------------------------------------------------------- ! CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF ! SIGNIFICANCE WHEN FNU OR N IS LARGE !----------------------------------------------------------------------} INU := Round(FNU); ARG := (FNU-1.0*INU)*PI; INU := INU + N - IL; AK := -SIN(ARG); BK := COS(ARG); IF ZI < 0.0 THEN BK := -BK; P1R := AK; P1I := BK; IF (INU Mod 2)=0 THEN GOTO 30; P1R := -P1R; P1I := -P1I; 30: For K:=1 to IL do begin SQK := FDN - 1.0; ATOL := S*ABS(SQK); SGN := 1.0; CS1R := CONER; CS1I := CONEI; CS2R := CONER; CS2I := CONEI; CKR := CONER; CKI := CONEI; AK := 0.0; AA := 1.0; BB := AEZ; DKR := EZR; DKI := EZI; For J:=1 to JL do begin ZDIV(CKR, CKI, DKR, DKI, STR, STI); CKR := STR*SQK; CKI := STI*SQK; CS2R := CS2R + CKR; CS2I := CS2I + CKI; SGN := -SGN; CS1R := CS1R + CKR*SGN; CS1I := CS1I + CKI*SGN; DKR := DKR + EZR; DKI := DKI + EZI; AA := AA*ABS(SQK)/BB; BB := BB + AEZ; AK := AK + 8.0; SQK := SQK - AK; IF AA <= ATOL THEN GOTO 50 end; GOTO 110; 50: S2R := CS1R; S2I := CS1I; IF ZR+ZR >= ELIM THEN GOTO 60; TZR := ZR + ZR; TZI := ZI + ZI; ZEXP(-TZR, -TZI, STR, STI); ZMLT(STR, STI, P1R, P1I, STR, STI); ZMLT(STR, STI, CS2R, CS2I, STR, STI); S2R := S2R + STR; S2I := S2I + STI; 60: FDN := FDN + 8.0*DFNU + 4.0; P1R := -P1R; P1I := -P1I; M := N - IL + K; YR[M] := S2R*AK1R - S2I*AK1I; YI[M] := S2R*AK1I + S2I*AK1R end; {K loop} IF N <= 2 THEN GOTO RETURN; NN := N; K := NN - 2; AK := 1.0*K; STR := ZR*RAZ; STI := -ZI*RAZ; RZR := (STR+STR)*RAZ; RZI := (STI+STI)*RAZ; IB := 3; For I:=IB to NN do begin YR[K] := (AK+FNU)*(RZR*YR[K+1]-RZI*YI[K+1]) + YR[K+2]; YI[K] := (AK+FNU)*(RZR*YI[K+1]+RZI*YR[K+1]) + YI[K+2]; AK := AK - 1.0; K := K - 1; IF KODED = 0 THEN GOTO RETURN end; ZEXP(CZR, CZI, CKR, CKI); For I:=1 to NN do begin STR := YR[I]*CKR - YI[I]*CKI; YI[I] := YR[I]*CKI + YI[I]*CKR; YR[I] := STR end; GOTO RETURN; 100: NZ := -1; GOTO RETURN; 110: NZ:=-2; Return:End; {ZASYI} Procedure ZMLRI(ZR, ZI, FNU:Double; KODE, N:Integer; Var YR, YI:VEC; Var NZ:Integer; Var TOL:Double); {***BEGIN PROLOGUE ZMLRI !***REFER TO ZBESI,ZBESK ! ! ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z) >= 0 BY THE ! MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES. ! !***ROUTINES CALLED DGAMLN,D1MACH,ZABS,ZEXP,ZLOG,ZMLT !***END PROLOGUE ZMLRI ! COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z } Label 20,30,40,50,70,90, 110, Return; Var ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI, CNORMR, CONEI, CONER, FKAP, FKK, FLAM, FNF, PTI, PTR, P1I, P1R, P2I, P2R, RAZ, RHO, RHO2, RZI, RZR, SCLE, STI, STR, SUMI, SUMR, TFNF, TST, ZEROI, ZEROR: Double; I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, M: Integer; Begin ZEROR:=0.0; ZEROI:=0.0; CONER:=1.0; CONEI:=0.0; SCLE := D1MACH(1)/TOL; NZ:=0; AZ := ZABS(ZR,ZI); IAZ := Round(AZ); IFNU := Round(FNU); INU := IFNU + N - 1; AT := 1.0*IAZ + 1.0; RAZ := 1.0/AZ; STR := ZR*RAZ; STI := -ZI*RAZ; CKR := STR*AT*RAZ; CKI := STI*AT*RAZ; RZR := (STR+STR)*RAZ; RZI := (STI+STI)*RAZ; P1R := ZEROR; P1I := ZEROI; P2R := CONER; P2I := CONEI; ACK := (AT+1.0)*RAZ; RHO := ACK + SQRT(ACK*ACK-1.0); RHO2 := RHO*RHO; TST := (RHO2+RHO2)/((RHO2-1.0)*(RHO-1.0)); TST := TST/TOL; {----------------------------------------------------------------------- ! COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES !----------------------------------------------------------------------} AK := AT; For I:=1 to 80 do begin PTR := P2R; PTI := P2I; P2R := P1R - (CKR*PTR-CKI*PTI); P2I := P1I - (CKI*PTR+CKR*PTI); P1R := PTR; P1I := PTI; CKR := CKR + RZR; CKI := CKI + RZI; AP := ZABS(P2R,P2I); IF AP > TST*AK*AK THEN GOTO 20; AK := AK + 1.0 end; GOTO 110; 20: I := I + 1; K := 0; IF INU < IAZ THEN GOTO 40; {----------------------------------------------------------------------- ! COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS !----------------------------------------------------------------------} P1R := ZEROR; P1I := ZEROI; P2R := CONER; P2I := CONEI; AT := 1.0*INU + 1.0; STR := ZR*RAZ; STI := -ZI*RAZ; CKR := STR*AT*RAZ; CKI := STI*AT*RAZ; ACK := AT*RAZ; TST := SQRT(ACK/TOL); ITIME := 1; For K:=1 to 80 do begin PTR := P2R; PTI := P2I; P2R := P1R - (CKR*PTR-CKI*PTI); P2I := P1I - (CKR*PTI+CKI*PTR); P1R := PTR; P1I := PTI; CKR := CKR + RZR; CKI := CKI + RZI; AP := ZABS(P2R,P2I); IF AP < TST THEN GOTO 30; IF ITIME = 2 THEN GOTO 40; ACK := ZABS(CKR,CKI); FLAM := ACK + SQRT(ACK*ACK-1.0); FKAP := AP/ZABS(P1R,P1I); RHO := DMIN(FLAM,FKAP); TST := TST*SQRT(RHO/(RHO*RHO-1.0)); ITIME := 2; 30: end; GOTO 110; {----------------------------------------------------------------------- ! BACKWARD RECURRENCE AND SUM NORMALIZING RELATION !----------------------------------------------------------------------} 40: K := K + 1; KK := IMAX(I+IAZ,K+INU); FKK := 1.0*KK; P1R := ZEROR; P1I := ZEROI; {----------------------------------------------------------------------- ! SCALE P2 AND SUM BY SCLE !----------------------------------------------------------------------} P2R := SCLE; P2I := ZEROI; FNF := FNU - 1.0*IFNU; TFNF := FNF + FNF; BK := DGAMLN(FKK+TFNF+1.0,IDUM)-DGAMLN(FKK+1.0,IDUM)-DGAMLN(TFNF+1.0,IDUM); BK := EXP(BK); SUMR := ZEROR; SUMI := ZEROI; KM := KK - INU; For I:=1 to KM do begin PTR := P2R; PTI := P2I; P2R := P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI); P2I := P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI); P1R := PTR; P1I := PTI; AK := 1.0 - TFNF/(FKK+TFNF); ACK := BK*AK; SUMR := SUMR + (ACK+BK)*P1R; SUMI := SUMI + (ACK+BK)*P1I; BK := ACK; FKK := FKK - 1.0 end; YR[N] := P2R; YI[N] := P2I; IF N = 1 THEN GOTO 70; For I:=2 to N do begin PTR := P2R; PTI := P2I; P2R := P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI); P2I := P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI); P1R := PTR; P1I := PTI; AK := 1.0 - TFNF/(FKK+TFNF); ACK := BK*AK; SUMR := SUMR + (ACK+BK)*P1R; SUMI := SUMI + (ACK+BK)*P1I; BK := ACK; FKK := FKK - 1.0; M := N - I + 1; YR[M] := P2R; YI[M] := P2I end; 70: IF IFNU <= 0 THEN GOTO 90; For I:=1 to IFNU do begin PTR := P2R; PTI := P2I; P2R := P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI); P2I := P1I + (FKK+FNF)*(RZR*PTI+RZI*PTR); P1R := PTR; P1I := PTI; AK := 1.0 - TFNF/(FKK+TFNF); ACK := BK*AK; SUMR := SUMR + (ACK+BK)*P1R; SUMI := SUMI + (ACK+BK)*P1I; BK := ACK; FKK := FKK - 1.0 end; 90: PTR := ZR; PTI := ZI; IF KODE = 2 THEN PTR := ZEROR; ZLOG(RZR, RZI, STR, STI, IDUM); P1R := -FNF*STR + PTR; P1I := -FNF*STI + PTI; AP := DGAMLN(1.0+FNF,IDUM); PTR := P1R - AP; PTI := P1I; {----------------------------------------------------------------------- ! THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW ! IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES !----------------------------------------------------------------------} P2R := P2R + SUMR; P2I := P2I + SUMI; AP := ZABS(P2R,P2I); P1R := 1.0/AP; ZEXP(PTR, PTI, STR, STI); CKR := STR*P1R; CKI := STI*P1R; PTR := P2R*P1R; PTI := -P2I*P1R; ZMLT(CKR, CKI, PTR, PTI, CNORMR, CNORMI); For I:=1 to N do begin STR := YR[I]*CNORMR - YI[I]*CNORMI; YI[I] := YR[I]*CNORMI + YI[I]*CNORMR; YR[I] := STR end; GOTO RETURN; 110: NZ:=-2; Return:End; {ZMLRI} Procedure ZSHCH(ZR, ZI: Double; Var CSHR, CSHI, CCHR, CCHI:Double); {***BEGIN PROLOGUE ZSHCH !***REFER TO ZBESK,ZBESH ! ! ZSHCH COMPUTES THE COMPLEX HYPERBOLIC FUNCTIONS CSH=SINH(X+I*Y) ! AND CCH=COSH(X+I*Y), WHERE I^2:=-1. ! !***ROUTINES CALLED (NONE) !***END PROLOGUE ZSHCH } Var CH, CN, SH, SN: Double; Function cosh(r:double):double; begin cosh:=(exp(r)+exp(-r))/2.0 end; Function sinh(r:double):double; begin sinh:=(exp(r)-exp(-r))/2.0 end; Begin SH := SINH(ZR); CH := COSH(ZR); SN := SIN(ZI); CN := COS(ZI); CSHR := SH*CN; CSHI := CH*SN; CCHR := CH*CN; CCHI := SH*SN End; Procedure ZRATI(ZR, ZI, FNU:Double; Var N:Integer; Var CYR, CYI:VEC; TOL:Double); {***BEGIN PROLOGUE ZRATI !***REFER TO ZBESI,ZBESK,ZBESH ! ! ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD ! RECURRENCE. THE STARTING INDEX IS DETERMINED BY FORWARD ! RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B, ! MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973, ! BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER, ! BY D. J. SOOKNE. ! !***ROUTINES CALLED ZABS,ZDIV !***END PROLOGUE ZRATI ! COMPLEX Z,CY(1),CONE,CZERO,P1,P2,T1,RZ,PT,CDFNU } Label 10,20,40,50, Return; Var AK, AMAGZ, AP1, AP2, ARG, AZ, CDFNUI, CDFNUR, CONEI, CONER, CZEROI, CZEROR, DFNU, FDNU, FLAM, FNUP, PTI, PTR, P1I, P1R, P2I, P2R, RAK, RAP1, RHO, RT2, RZI, RZR, TEST, TEST1, TTI, TTR, T1I, T1R: Double; I,ID, IDNU, INU, ITIME, K, KK, MAGZ: Integer; Begin CZEROR:=0.0; CZEROI:=0.0; CONER:=1.0; CONEI:=0.0; RT2:=1.41421356237309505; AZ := ZABS(ZR,ZI); INU := Round(FNU); IDNU := INU + N - 1; MAGZ := Round(AZ); AMAGZ := 1.0*(MAGZ+1); FDNU := 1.0*IDNU; FNUP := DMAX(AMAGZ,FDNU); ID := IDNU - MAGZ - 1; ITIME := 1; K := 1; PTR := 1.0/AZ; RZR := PTR*(ZR+ZR)*PTR; RZI := -PTR*(ZI+ZI)*PTR; T1R := RZR*FNUP; T1I := RZI*FNUP; P2R := -T1R; P2I := -T1I; P1R := CONER; P1I := CONEI; T1R := T1R + RZR; T1I := T1I + RZI; IF ID > 0 THEN ID := 0; AP2 := ZABS(P2R,P2I); AP1 := ZABS(P1R,P1I); {----------------------------------------------------------------------- ! THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU ! GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT ! P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR ! PREMATURELY. !----------------------------------------------------------------------} ARG := (AP2+AP2)/(AP1*TOL); TEST1 := SQRT(ARG); TEST := TEST1; RAP1 := 1.0/AP1; P1R := P1R*RAP1; P1I := P1I*RAP1; P2R := P2R*RAP1; P2I := P2I*RAP1; AP2 := AP2*RAP1; 10: K := K + 1; AP1 := AP2; PTR := P2R; PTI := P2I; P2R := P1R - (T1R*PTR-T1I*PTI); P2I := P1I - (T1R*PTI+T1I*PTR); P1R := PTR; P1I := PTI; T1R := T1R + RZR; T1I := T1I + RZI; AP2 := ZABS(P2R,P2I); IF AP1 <= TEST THEN GOTO 10; IF ITIME = 2 THEN GOTO 20; AK := ZABS(T1R,T1I)*0.5; FLAM := AK + SQRT(AK*AK-1.0); RHO := DMIN(AP2/AP1,FLAM); TEST := TEST1*SQRT(RHO/(RHO*RHO-1.0)); ITIME := 2; GOTO 10; 20: KK := K + 1 - ID; AK := 1.0*KK; T1R := AK; T1I := CZEROI; DFNU := FNU + 1.0*(N-1); P1R := 1.0/AP2; P1I := CZEROI; P2R := CZEROR; P2I := CZEROI; For I:=1 to KK do begin PTR := P1R; PTI := P1I; RAP1 := DFNU + T1R; TTR := RZR*RAP1; TTI := RZI*RAP1; P1R := (PTR*TTR-PTI*TTI) + P2R; P1I := (PTR*TTI+PTI*TTR) + P2I; P2R := PTR; P2I := PTI; T1R := T1R - CONER end; IF (P1R <> CZEROR) OR (P1I <> CZEROI) THEN GOTO 40; P1R := TOL; P1I := TOL; 40: ZDIV(P2R, P2I, P1R, P1I, CYR[N], CYI[N]); IF N = 1 THEN GOTO RETURN; K := N - 1; AK := 1.0*K; T1R := AK; T1I := CZEROI; CDFNUR := FNU*RZR; CDFNUI := FNU*RZI; For I:=2 to N do begin PTR := CDFNUR + (T1R*RZR-T1I*RZI) + CYR[K+1]; PTI := CDFNUI + (T1R*RZI+T1I*RZR) + CYI[K+1]; AK := ZABS(PTR,PTI); IF AK <> CZEROR THEN GOTO 50; PTR := TOL; PTI := TOL; AK := TOL*RT2; 50: RAK := CONER/AK; CYR[K] := RAK*PTR*RAK; CYI[K] := -RAK*PTI*RAK; T1R := T1R - CONER; K := K - 1 end; Return: End; {ZRATI} END. {end of file Cbess0.pas}