/************************************************************************* * 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. * * * * C++ Release By J-P Moreau, Paris (07/01/2005). * * (www.jpmoreau.fr) * *************************************************************************/ // Functions using only module COMPLEX.cpp #include #include #include #include "complex.h" REAL Sign(REAL a, REAL b) { if (b <0.0) return (-ABS(a)); else return ABS(a); } void ZUCHK(REAL YR, REAL YI, int *NZ, REAL ASCLE, REAL TOL) { /***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 */ REAL SS, ST, WR, WI; *NZ = 0; WR = ABS(YR); WI = ABS(YI); ST = DMIN(WR,WI); if (ST > ASCLE) return; SS = DMAX(WR,WI); ST = ST/TOL; if (SS < ST) *NZ = 1; } //ZUCHK() REAL DGAMLN(REAL Z, int *IERR) { /***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 REAL PRECISION ROUTINE **** ! DGAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR ! Z > 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 <= 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 > 0 ! ! OUTPUT DGAMLN IS REAL PRECISION ! DGAMLN - NATURAL LOG OF THE GAMMA FUNCTION AT Z <> 0 ! 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 module COMPLEX). !***END PROLOGUE DGAMLN */ //Labels: e10,e20,e40,e50,e70 REAL CON, FLN, FZ, RLN, S, TLG, TRM, TST, T1, WDTOL, ZDMY, ZINC, ZM, ZMIN, ZP, ZSQ; int I, I1M, K, MZ, NZ; REAL *CF; REAL *GLN; void *vmblock = NULL; // LNGAMMA(N), N=1,100 vmblock = vminit(); GLN = (REAL *) vmalloc(vmblock, VEKTOR, 101, 0); //index 0 not used CF = (REAL *) vmalloc(vmblock, VEKTOR, 23, 0); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return 0; } 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) goto e70; if (Z > 101.0) goto e10; NZ = (int) floor(Z); FZ = Z - 1.0*NZ; if (FZ > 0.0) goto e10; if (NZ > 100) goto e10; vmfree(vmblock); return (GLN[NZ]); e10: 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 = (int) (floor(ZM) + 1); ZMIN = 1.0*MZ; ZDMY = Z; ZINC = 0.0; if (Z >= ZMIN) goto e20; ZINC = ZMIN - 1.0*NZ; ZDMY = Z + ZINC; e20: ZP = 1.0/ZDMY; T1 = CF[1]*ZP; S = T1; if (ZP < WDTOL) goto e40; ZSQ = ZP*ZP; TST = T1*WDTOL; for (K=2; K<=22; K++) { ZP = ZP*ZSQ; TRM = CF[K]*ZP; if (ABS(TRM) < TST) goto e40; S += TRM; } e40: if (ZINC != 0.0) goto e50; TLG = log(Z); vmfree(vmblock); return (Z*(TLG-1.0) + 0.5*(CON-TLG) + S); e50: ZP = 1.0; NZ = (int) floor(ZINC); for (I=1; I<=NZ; I++) ZP *= (Z+1.0*(I-1)); TLG = log(ZDMY); vmfree(vmblock); return (ZDMY*(TLG-1.0) - log(ZP) + 0.5*(CON-TLG) + S); e70: *IERR=1; vmfree(vmblock); return 0; } //DGAMLN() void ZUNIK(REAL ZRR, REAL ZRI, REAL FNU, int IKFLG, int IPMTR, REAL TOL, int INIT, REAL *PHIR, REAL *PHII, REAL *ZETA1R, REAL *ZETA1I, REAL *ZETA2R, REAL *ZETA2I, REAL *SUMR, REAL *SUMI, REAL *CWRKR, REAL *CWRKI) { /***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 */ //Labels: e30, e40, e60 REAL *C; void *vmblock = NULL; REAL AC, CONEI, CONER, CRFNI, CRFNR, RFN, SI, SR, SRI, SRR, STI, STR, TEST, TI, TR, T2I, T2R, ZEROI, ZEROR, ZNI, ZNR; int I, IDUM, J, K, L; REAL CON[3]; //index 0 not used //*** First executable statement ZUNIK ZEROR=0.0; ZEROI= 0.0; CONER=1.0; CONEI=0.0; CON[0]=0.0; CON[1]=3.98942280401432678e-01; CON[2]=1.25331413731550025; //Initialize table C vmblock = vminit(); C = (REAL *) vmalloc(vmblock, VEKTOR, 121, 0); //index 0 not used if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return; } 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) goto e40; /*---------------------------------------------------------------------- ! 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) { vmfree(vmblock); 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; K<=15; K++) { SR = ZEROR; SI = ZEROI; for (J=1; J<=K; J++) { L = L + 1; STR = SR*T2R - SI*T2I + C[L]; SI = SR*T2I + SI*T2R; SR = STR; } 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 && TEST < TOL) goto e30; } K = 15; e30: INIT = K; e40: if (IKFLG == 2) goto e60; /*---------------------------------------------------------------------- ! COMPUTE SUM FOR THE I FUNCTION !---------------------------------------------------------------------*/ SR = ZEROR; SI = ZEROI; for (I=1; I<=INIT; I++) { SR = SR + CWRKR[I]; SI = SI + CWRKI[I]; } *SUMR = SR; *SUMI = SI; *PHIR = CWRKR[16]*CON[1]; *PHII = CWRKI[16]*CON[1]; vmfree(vmblock); return; /*---------------------------------------------------------------------- ! COMPUTE SUM FOR THE K FUNCTION !---------------------------------------------------------------------*/ e60: SR = ZEROR; SI = ZEROI; TR = CONER; for (I=1; I<=INIT; I++) { SR = SR + TR*CWRKR[I]; SI = SI + TR*CWRKI[I]; TR = -TR; } *SUMR = SR; *SUMI = SI; *PHIR = CWRKR[16]*CON[2]; *PHII = CWRKI[16]*CON[2]; vmfree(vmblock); } //ZUNIK() void ZUNHJ(REAL ZR, REAL ZI, REAL FNU, int IPMTR, REAL TOL, REAL *PHIR, REAL *PHII, REAL *ARGR, REAL *ARGI, REAL *ZETA1R, REAL *ZETA1I, REAL *ZETA2R, REAL *ZETA2I, REAL *ASUMR, REAL *ASUMI, REAL *BSUMR, REAL *BSUMI) { /***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 */ //Labels: e20, e50,e60, e80,e90,e110,e120,e130,e140,e180,e200,e220 REAL ANG, ATOL, AW2, AZTH, BTOL, CONEI, CONER, EX1, EX2, FN13, FN23, GXPI, 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; int IAS, IBS, IS, J, JR, JU, K, KMAX, KP1, KS, L, LR, LRP1, L1, L2, M, IDUM; REAL *AR, *BR, *UPR, *UXPI, *CRR, *CRI, *DRR, *DRI; //size 0..16 REAL *C; //size 0..105 REAL *ALFA; //size 0..180 REAL *BETA; //size 0..210 REAL *GAMA, *AP, *PR, *XPI; //size 0..30 void *vmblock = NULL; //Initialize dynamic tables vmblock = vminit(); C = (REAL *) vmalloc(vmblock, VEKTOR, 106, 0); //index 0 not used ALFA = (REAL *) vmalloc(vmblock, VEKTOR, 181, 0); BETA = (REAL *) vmalloc(vmblock, VEKTOR, 211, 0); GAMA = (REAL *) vmalloc(vmblock, VEKTOR, 31, 0); AP = (REAL *) vmalloc(vmblock, VEKTOR, 31, 0); PR = (REAL *) vmalloc(vmblock, VEKTOR, 31, 0); XPI = (REAL *) vmalloc(vmblock, VEKTOR, 31, 0); AR = (REAL *) vmalloc(vmblock, VEKTOR, 17, 0); BR = (REAL *) vmalloc(vmblock, VEKTOR, 17, 0); UPR = (REAL *) vmalloc(vmblock, VEKTOR, 17, 0); UXPI = (REAL *) vmalloc(vmblock, VEKTOR, 17, 0); CRR = (REAL *) vmalloc(vmblock, VEKTOR, 17, 0); CRI = (REAL *) vmalloc(vmblock, VEKTOR, 17, 0); DRR = (REAL *) vmalloc(vmblock, VEKTOR, 17, 0); DRI = (REAL *) vmalloc(vmblock, VEKTOR, 17, 0); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return; } AR[0]= 0.0; 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[0]= 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[0] = 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[0]= 0.0; 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[0]= 0.0; 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[0]= 0.0; 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; GXPI=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 = pow(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) goto e130; /*---------------------------------------------------------------------- ! POWER SERIES FOR CABS(W2) <= 0.25 !---------------------------------------------------------------------*/ K = 1; PR[1] = CONER; XPI[1] = CONEI; SUMAR = GAMA[1]; SUMAI = ZEROI; AP[1] = 1.0; if (AW2 < TOL) goto e20; for (K=2; K<=30; K++) { PR[K] = PR[K-1]*W2R - XPI[K-1]*W2I; XPI[K] = PR[K-1]*W2I + XPI[K-1]*W2R; SUMAR = SUMAR + PR[K]*GAMA[K]; SUMAI = SUMAI + XPI[K]*GAMA[K]; AP[K] = AP[K-1]*AW2; if (AP[K] < TOL) goto e20; } K = 30; e20: 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) goto e120; /*---------------------------------------------------------------------- ! SUM SERIES FOR ASUM AND BSUM !---------------------------------------------------------------------*/ SUMBR = ZEROR; SUMBI = ZEROI; for (K=1; K<=KMAX; K++) { SUMBR = SUMBR + PR[K]*BETA[K]; SUMBI = SUMBI + XPI[K]*BETA[K]; } *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) goto e110; for (IS=2; IS<8; IS++) { ATOL = ATOL/RFNU2; PP = PP*RFNU2; if (IAS == 1) goto e60; SUMAR = ZEROR; SUMAI = ZEROI; for (K=1; K<=KMAX; K++) { M = L1 + K; SUMAR = SUMAR + PR[K]*ALFA[M]; SUMAI = SUMAI + XPI[K]*ALFA[M]; if (AP[K] < ATOL) goto e50; } e50: *ASUMR = *ASUMR + SUMAR*PP; *ASUMI = *ASUMI + SUMAI*PP; if (PP < TOL) IAS = 1; e60: if (IBS == 1) goto e90; SUMBR = ZEROR; SUMBI = ZEROI; for (K=1; K<=KMAX; K++) { M = L2 + K; SUMBR = SUMBR + PR[K]*BETA[M]; SUMBI = SUMBI + XPI[K]*BETA[M]; if (AP[K] < ATOL) goto e80; } e80: *BSUMR = *BSUMR + SUMBR*PP; *BSUMI = *BSUMI + SUMBI*PP; if (PP < BTOL) IBS = 1; e90: if (IAS == 1 && IBS == 1) goto e110; L1 += 30; L2 += 30; } //IS loop e110: *ASUMR = *ASUMR + CONER; PP = RFNU*RFN13; *BSUMR = *BSUMR*PP; *BSUMI = *BSUMI*PP; e120: vmfree(vmblock); return; /*---------------------------------------------------------------------- ! CABS(W2) > 0.25 !---------------------------------------------------------------------*/ e130: ZSQRT(W2R, W2I, &WR, &WI); if (WR < 0.0) WR = 0.0; if (WI < 0.0) 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) ZCI = 0.0; if (ZCI > HPI) ZCI = HPI; if (ZCR < 0.0) 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 && ZTHI < 0.0) goto e140; ANG = HPI; if (ZTHR == 0.0) goto e140; ANG = atan(ZTHI/ZTHR); if (ZTHR < 0.0) ANG += GXPI; e140: PP = pow(AZTH,EX2); ANG = ANG*EX2; ZETAR = PP*cos(ANG); ZETAI = PP*sin(ANG); if (ZETAI < 0.0) 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) goto e120; 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; UXPI[2] = STR*TFNI + STI*TFNR; *BSUMR = UPR[2] + ZCR; *BSUMI = UXPI[2] + ZCI; *ASUMR = ZEROR; *ASUMI = ZEROI; if (RFNU < TOL) goto e220; PRZTHR = RZTHR; PRZTHI = RZTHI; PTFNR = TFNR; PTFNI = TFNI; UPR[1] = CONER; UXPI[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) { LRP1 = LR + 1; /*---------------------------------------------------------------------- ! COMPUTE TWO ADDITIONAL CR, DR, AND UP FOR TWO MORE TERMS IN ! NEXT SUMA AND SUMB !---------------------------------------------------------------------*/ for (K=LR; K<=LRP1; K++) { KS++; KP1++; L++; ZAR = C[L]; ZAI = ZEROI; for (J=2; J<=KP1; J++) { L++; STR = ZAR*T2R - T2I*ZAI + C[L]; ZAI = ZAR*T2I + ZAI*T2R; ZAR = STR; } STR = PTFNR*TFNR - PTFNI*TFNI; PTFNI = PTFNR*TFNI + PTFNI*TFNR; PTFNR = STR; UPR[KP1] = PTFNR*ZAR - PTFNI*ZAI; UXPI[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]; } //K loop PP = PP*RFNU2; if (IAS == 1) goto e180; SUMAR = UPR[LRP1]; SUMAI = UXPI[LRP1]; JU = LRP1; for (JR=1; JR<=LR; JR++) { JU = JU - 1; SUMAR = SUMAR + CRR[JR]*UPR[JU] - CRI[JR]*UXPI[JU]; SUMAI = SUMAI + CRR[JR]*UXPI[JU] + CRI[JR]*UPR[JU]; } *ASUMR = *ASUMR + SUMAR; *ASUMI = *ASUMI + SUMAI; TEST = ABS(SUMAR) + ABS(SUMAI); if (PP < TOL && TEST < TOL) IAS = 1; e180: if (IBS == 1) goto e200; SUMBR = UPR[LR+2] + UPR[LRP1]*ZCR - UXPI[LRP1]*ZCI; SUMBI = UXPI[LR+2] + UPR[LRP1]*ZCI + UXPI[LRP1]*ZCR; JU = LRP1; for (JR=1; JR<=LR; JR++) { JU = JU - 1; SUMBR = SUMBR + DRR[JR]*UPR[JU] - DRI[JR]*UXPI[JU]; SUMBI = SUMBI + DRR[JR]*UXPI[JU] + DRI[JR]*UPR[JU]; } *BSUMR = *BSUMR + SUMBR; *BSUMI = *BSUMI + SUMBI; TEST = ABS(SUMBR) + ABS(SUMBI); if (PP < BTOL && TEST < BTOL) IBS = 1; e200: if (IAS == 1 && IBS == 1) goto e220; LR += 2; } //while LR e220: *ASUMR = *ASUMR + CONER; STR = -(*BSUMR)*RFN13; STI = -(*BSUMI)*RFN13; ZDIV(STR, STI, RTZTR, RTZTI, BSUMR, BSUMI); goto e120; } //ZUNHJ() void ZUOIK(REAL ZR, REAL ZI, REAL FNU, int KODE, int IKFLG, int N, REAL *YR, REAL *YI, int *NUF, REAL TOL, REAL ELIM, REAL ALIM) { /***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 < 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 > 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 = N MEANS ALL Y VALUES WERE SET TO ZERO ! IKFLG=2 AND 0 < NUF < 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 */ //Labels: e10,e20,e30,e40,e50,e60,e70,e80,e90,e110,e120,e130,e140,e150,e160,e170, // e180,e190,e200,e210 REAL 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; int I, IDUM, IFORM, INIT, NN, NW; REAL *CWRKR, *CWRKI; void *vmblock = NULL; ZEROR=0.0; ZEROI=0.0; // Initialize CWRKR, CWRKI vmblock = vminit(); CWRKR = (REAL *) vmalloc(vmblock, VEKTOR, 17, 0); //index 0 not used CWRKI = (REAL *) vmalloc(vmblock, VEKTOR, 17, 0); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return; } AIC=1.265512123484645396; *NUF = 0; NN = N; ZRR = ZR; ZRI = ZI; if (ZR >= 0.0) goto e10; ZRR = -ZR; ZRI = -ZI; e10: ZBR = ZRR; ZBI = ZRI; AX = ABS(ZR)*1.7321; AY = ABS(ZI); IFORM = 1; if (AY > AX) IFORM = 2; GNU = DMAX(FNU,1.0); if (IKFLG == 1) goto e20; 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. !---------------------------------------------------------------------*/ e20: if (IFORM == 2) goto e30; 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 e50; e30: ZNR = ZRI; ZNI = -ZRR; if (ZI > 0.0) goto e40; ZNR = -ZNR; e40: 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); e50: if (KODE == 1) goto e60; CZR -= ZBR; CZI -= ZBI; e60: if (IKFLG == 1) goto e70; CZR = -CZR; CZI = -CZI; e70: APHI = ZABS(PHIR,PHII); RCZ = CZR; /*---------------------------------------------------------------------- ! OVERFLOW TEST !---------------------------------------------------------------------*/ if (RCZ > ELIM) goto e210; if (RCZ < ALIM) goto e80; RCZ += log(APHI); if (IFORM == 2) RCZ -= (0.25*log(AARG) - AIC); if (RCZ > ELIM) goto e210; goto e130; /*---------------------------------------------------------------------- ! UNDERFLOW TEST !---------------------------------------------------------------------*/ e80: if (RCZ < -ELIM) goto e90; if (RCZ > -ALIM) goto e130; RCZ += log(APHI); if (IFORM == 2) RCZ -= (0.25*log(AARG) - AIC); if (RCZ > -ELIM) goto e110; e90: for (I=1; I<=NN; I++) { YR[I] = ZEROR; YI[I] = ZEROI; } *NUF = NN; vmfree(vmblock); return; e110: ASCLE = 1000*D1MACH(1)/TOL; ZLOG(PHIR, PHII, &STR, &STI, &IDUM); CZR = CZR + STR; CZI = CZI + STI; if (IFORM == 1) goto e120; ZLOG(ARGR, ARGI, &STR, &STI, &IDUM); CZR = CZR - 0.25*STR - AIC; CZI = CZI - 0.25*STI; e120: AX = EXP(RCZ)/TOL; AY = CZI; CZR = AX*COS(AY); CZI = AX*SIN(AY); ZUCHK(CZR, CZI, &NW, ASCLE, TOL); if (NW != 0) goto e90; e130: if (IKFLG == 2) { vmfree(vmblock); return; } if (N == 1) { vmfree(vmblock); return; } /*---------------------------------------------------------------------- ! SET UNDERFLOWS ON I SEQUENCE !---------------------------------------------------------------------*/ e140: GNU = FNU + 1.0*(NN-1); if (IFORM == 2) goto e150; 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 e160; e150: 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); e160: if (KODE == 1) goto e170; CZR = CZR - ZBR; CZI = CZI - ZBI; e170: APHI = ZABS(PHIR,PHII); RCZ = CZR; if (RCZ < -ELIM) goto e180; if (RCZ > -ALIM) { vmfree(vmblock); return; } RCZ += log(APHI); if (IFORM == 2) RCZ -= (0.25*log(AARG) - AIC); if (RCZ > -ELIM) goto e190; e180: YR[NN] = ZEROR; YI[NN] = ZEROI; NN = NN - 1; *NUF = *NUF + 1; if (NN == 0) { vmfree(vmblock); return; } goto e140; e190: ASCLE = 1000*D1MACH(1)/TOL; ZLOG(PHIR, PHII, &STR, &STI, &IDUM); CZR = CZR + STR; CZI = CZI + STI; if (IFORM == 1) goto e200; ZLOG(ARGR, ARGI, &STR, &STI, &IDUM); CZR = CZR - 0.25*STR - AIC; CZI = CZI - 0.25*STI; e200: AX = EXP(RCZ)/TOL; AY = CZI; CZR = AX*COS(AY); CZI = AX*SIN(AY); ZUCHK(CZR, CZI, &NW, ASCLE, TOL); if (NW != 0) goto e180; vmfree(vmblock); return; e210: *NUF = -1; vmfree(vmblock); } //ZUOIK() void ZS1S2(REAL *ZRR, REAL *ZRI, REAL *S1R, REAL *S1I, REAL *S2R, REAL *S2I, int *NZ, REAL ASCLE, REAL ALIM, int *IUF) { /***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 e10 REAL AA, ALN, AS1, AS2, C1I, C1R, S1DI, S1DR, ZEROI, ZEROR; int IDUM; ZEROR=0.0; ZEROI=0.0; *NZ = 0; AS1 = ZABS(*S1R,*S1I); AS2 = ZABS(*S2R,*S2I); if (*S1R == 0.0 && *S1I == 0.0) goto e10; if (AS1 == 0.0) goto e10; ALN = -(*ZRR) -(*ZRR) + log(AS1); S1DR = *S1R; S1DI = *S1I; *S1R = ZEROR; *S1I = ZEROI; AS1 = ZEROR; if (ALN < -ALIM) goto e10; 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; e10: AA = DMAX(AS1,AS2); if (AA > ASCLE) return; *S1R = ZEROR; *S1I = ZEROI; *S2R = ZEROR; *S2I = ZEROI; *NZ = 1; *IUF = 0; } //ZS1S2() //for debug only void test(REAL ZR, REAL ZI, REAL FNU, int KODE, int N, REAL *YR, REAL *YI, int *NZ, REAL TOL, REAL ELIM, REAL ALIM) { printf(" %f %f %f %d %d %f %f %d %f %f %f\n", ZR, ZI, FNU, KODE, N, YR[1], YI[1], *NZ, TOL, ELIM, ALIM); } void ZSERI(REAL ZR, REAL ZI, REAL FNU, int KODE, int N, REAL *YR, REAL *YI, int *NZ, REAL TOL, REAL ELIM, REAL ALIM) { /***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 */ //Labels: e10,e20,e30,e40,e50,e60,e70,e80,e90,e100,e120,e140,e150,e160, // e170,e190 REAL 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; int I, IB, idebug,IDUM, IFLAG, IL, K, L, M, NN, NW; REAL *WR, *WI; void *vmblock = NULL; FILE *fp; *NZ=0; // Initialize WR, WI vmblock = vminit(); WR = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0); //index 0 not used WI = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return; } ZEROR=0.0; ZEROI=0.0; CONER=1.0; CONEI=0.0; idebug=1; // Note: In C++ version, results are not correct with idebug=0 (any idea?). if (idebug) fp=fopen("cbess0.txt","w"); AZ = ZABS(ZR,ZI); if (AZ == 0.0) goto e160; ARM = 1000*D1MACH(1); RTR1 = SQRT(ARM); CRSCR = 1.0; IFLAG = 0; if (AZ < ARM) goto e150; HZR = 0.5*ZR; HZI = 0.5*ZI; CZR = ZEROR; CZI = ZEROI; if (AZ <= RTR1) goto e10; ZMLT(HZR, HZI, HZR, HZI, &CZR, &CZI); e10: ACZ = ZABS(CZR,CZI); NN = N; ZLOG(HZR, HZI, &CKR, &CKI, &IDUM); e20: DFNU = FNU + 1.0*(NN-1); FNUP = DFNU + 1.0; /*---------------------------------------------------------------------- ! UNDERFLOW TEST !---------------------------------------------------------------------*/ AK1R = CKR*DFNU; AK1I = CKI*DFNU; AK = DGAMLN(FNUP,&IDUM); AK1R -= AK; if (KODE == 2) AK1R -= ZR; if (AK1R > -ELIM) goto e40; e30: *NZ = *NZ + 1; YR[NN] = ZEROR; YI[NN] = ZEROI; if (ACZ > DFNU) goto e190; NN--; if (NN == 0) { vmfree(vmblock); return; } goto e20; e40: if (AK1R > -ALIM) goto e50; IFLAG = 1; SS = 1.0/TOL; CRSCR = TOL; ASCLE = ARM*SS; e50: AA = EXP(AK1R); if (IFLAG == 1) AA *= SS; COEFR = AA*COS(AK1I); COEFI = AA*SIN(AK1I); ATOL = TOL*ACZ/FNUP; IL = IMIN(2,NN); for (I=1; I<=IL; I++) { DFNU = FNU + 1.0*(NN-I); FNUP = DFNU + 1.0; S1R = CONER; S1I = CONEI; if (ACZ < TOL*FNUP) goto e70; AK1R = CONER; AK1I = CONEI; AK = FNUP + 2.0; S = FNUP; AA = 2.0; e60: RS = 1.0/S; if (idebug) fprintf(fp," ** RS=%e\n", RS); 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) goto e60; e70: S2R = S1R*COEFR - S1I*COEFI; S2I = S1R*COEFI + S1I*COEFR; WR[I] = S2R; WI[I] = S2I; if (IFLAG == 0) goto e80; ZUCHK(S2R, S2I, &NW, ASCLE, TOL); if (NW != 0) goto e30; e80: M = NN - I + 1; YR[M] = S2R*CRSCR; YI[M] = S2I*CRSCR; if (I == IL) goto e90; ZDIV(COEFR, COEFI, HZR, HZI, &STR, &STI); COEFR = STR*DFNU; COEFI = STI*DFNU; e90: ;} // I loop if (idebug) fclose(fp); if (NN <= 2) { vmfree(vmblock); 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) goto e120; IB = 3; e100: for (I=IB; I<=NN; I++) { 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 -= 1.0; K--; } vmfree(vmblock); return; /*---------------------------------------------------------------------- ! RECUR BACKWARD WITH SCALED VALUES !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE ! UNDERFLOW LIMIT = ASCLE = D1MACH(1)*SS*1000. !---------------------------------------------------------------------*/ e120: S1R = WR[1]; S1I = WI[1]; S2R = WR[2]; S2I = WI[2]; for (L=3; L<=NN; L++) { 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) goto e140; } vmfree(vmblock); return; e140: IB = L + 1; if (IB > NN) { vmfree(vmblock); return; } goto e100; e150: *NZ = N; if (FNU == 0.0) *NZ = *NZ - 1; e160: YR[1] = ZEROR; YI[1] = ZEROI; if (FNU != 0.0) goto e170; YR[1] = CONER; YI[1] = CONEI; e170: if (N == 1) { vmfree(vmblock); return; } for (I=2; I<=N; I++) { YR[I] = ZEROR; YI[I] = ZEROI; } vmfree(vmblock); return; /*---------------------------------------------------------------------- ! RETURN WITH NZ < 0 IF CABS(Z*Z/4) > FNU+N-NZ-1 COMPLETE ! THE CALCULATION IN CBINU WITH N=N-ABS(NZ) !---------------------------------------------------------------------*/ e190: *NZ = -(*NZ); vmfree(vmblock); } //ZSERI() void ZASYI(REAL ZR, REAL ZI, REAL FNU, int KODE, int N, REAL *YR, REAL *YI, int *NZ, REAL RL, REAL TOL, REAL ELIM, REAL ALIM) { /***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 */ //Labels: e10,e20,e30,e50,e60, e100,e110 REAL 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; int I, IB, IL, INU, J, JL, K, KODED, M, NN; 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) goto e10; CZR = ZEROR; CZI = ZI; e10: if (ABS(CZR) > ELIM) goto e100; DNU2 = DFNU + DFNU; KODED = 1; if (ABS(CZR) > ALIM && N > 2) goto e20; KODED = 0; ZEXP(CZR, CZI, &STR, &STI); ZMLT(AK1R, AK1I, STR, STI, &AK1R, &AK1I); e20: FDN = 0.0; if (DNU2 > RTR1) 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 = (int) (floor(RL+RL) + 2); P1R = ZEROR; P1I = ZEROI; if (ZI == 0.0) goto e30; /*---------------------------------------------------------------------- ! CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF ! SIGNIFICANCE WHEN FNU OR N IS LARGE !---------------------------------------------------------------------*/ INU = (int) floor(FNU); ARG = (FNU-1.0*INU)*PI; INU = INU + N - IL; AK = -SIN(ARG); BK = COS(ARG); if (ZI < 0.0) BK = -BK; P1R = AK; P1I = BK; if ((INU % 2)==0) goto e30; P1R = -P1R; P1I = -P1I; e30: for (K=1; K<=IL; K++) { 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; J<=JL; J++) { 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) goto e50; } goto e110; e50: S2R = CS1R; S2I = CS1I; if (ZR+ZR >= ELIM) goto e60; 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; e60: 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; } // K loop if (N <= 2) 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; I<=NN; I++) { 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 -= 1.0; K--; if (KODED == 0) return; } ZEXP(CZR, CZI, &CKR, &CKI); for (I=1; I<=NN; I++) { STR = YR[I]*CKR - YI[I]*CKI; YI[I] = YR[I]*CKI + YI[I]*CKR; YR[I] = STR; } return; e100: *NZ = -1; return; e110: *NZ=-2; } //ZASYI() void ZMLRI(REAL ZR, REAL ZI, REAL FNU, int KODE, int N, REAL *YR, REAL *YI, int * NZ, REAL TOL) { /***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 */ //Labels: e20,e30,e40,e50,e70,e90, e110 REAL 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; int I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, M; ZEROR=0.0; ZEROI=0.0; CONER=1.0; CONEI=0.0; SCLE = D1MACH(1)/TOL; *NZ=0; AZ = ZABS(ZR,ZI); IAZ = (int) floor(AZ); IFNU = (int) floor(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; I<81; I++) { 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) goto e20; AK += 1.0; } goto e110; e20: I++; K = 0; if (INU < IAZ) goto e40; /*---------------------------------------------------------------------- ! 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; K<81; K++) { 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) goto e30; if (ITIME == 2) goto e40; 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; e30: ;} goto e110; /*---------------------------------------------------------------------- ! BACKWARD RECURRENCE AND SUM NORMALIZING RELATION !---------------------------------------------------------------------*/ e40: K++; 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; I<=KM; I++) { 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 -= 1.0; } YR[N] = P2R; YI[N] = P2I; if (N == 1) goto e70; for (I=2; I<=N; I++) { 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 -= 1.0; M = N - I + 1; YR[M] = P2R; YI[M] = P2I; } e70: if (IFNU <= 0) goto e90; for (I=1; I<=IFNU; I++) { 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 -= 1.0; } e90: PTR = ZR; PTI = ZI; if (KODE == 2) 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; I<=N; I++) { STR = YR[I]*CNORMR - YI[I]*CNORMI; YI[I] = YR[I]*CNORMI + YI[I]*CNORMR; YR[I] = STR; } return; e110: *NZ=-2; } //ZMLRI() //COSH is defined in basis.h REAL SINH(REAL r) { return (1.0+COSH(r)*COSH(r)); } void ZSHCH(REAL ZR, REAL ZI, REAL *CSHR, REAL *CSHI, REAL *CCHR, REAL *CCHI) { /***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 */ REAL CH, CN, SH, SN; SH = SINH(ZR); CH = COSH(ZR); SN = SIN(ZI); CN = COS(ZI); *CSHR = SH*CN; *CSHI = CH*SN; *CCHR = CH*CN; *CCHI = SH*SN; } void ZRATI(REAL ZR, REAL ZI, REAL FNU, int N, REAL *CYR, REAL *CYI, REAL TOL) { /***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 */ //Labels: e10, e20, e40, e50 REAL 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; int I,ID, IDNU, INU, ITIME, K, KK, MAGZ; CZEROR=0.0; CZEROI=0.0; CONER=1.0; CONEI=0.0; RT2=1.41421356237309505; AZ = ZABS(ZR,ZI); INU = (int) floor(FNU); IDNU = INU + N - 1; MAGZ = (int) floor(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) 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; e10: K++; 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) goto e10; if (ITIME == 2) goto e20; 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 e10; e20: 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; I<=KK; I++) { 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 -= CONER; } if (P1R != CZEROR || P1I != CZEROI) goto e40; P1R = TOL; P1I = TOL; e40: ZDIV(P2R, P2I, P1R, P1I, &CYR[N], &CYI[N]); if (N == 1) return; K = N - 1; AK = 1.0*K; T1R = AK; T1I = CZEROI; CDFNUR = FNU*RZR; CDFNUI = FNU*RZI; for (I=2; I<=N; I++) { 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) goto e50; PTR = TOL; PTI = TOL; AK = TOL*RT2; e50: RAK = CONER/AK; CYR[K] = RAK*PTR*RAK; CYI[K] = -RAK*PTI*RAK; T1R = T1R - CONER; K--; } } //ZRATI() //end of file Cbess0.cpp