C ALGORITHM 831, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 30, NO. 2, June, 2004, P. 159--164. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Doc/ # Fortran77/ # Fortran77/Drivers/ # Fortran77/Drivers/Res # Fortran77/Drivers/wktest.f # Fortran77/Src/ # Fortran77/Src/airy.f # Fortran77/Src/d1mach.f # Fortran77/Src/src.f # This archive created: Wed Aug 25 13:06:12 2004 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' cd .. if test ! -d 'Fortran77' then mkdir 'Fortran77' fi cd 'Fortran77' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'Res' then echo shar: will not over-write existing file "'Res'" else cat << "SHAR_EOF" > 'Res' ERR MEANS RELATIVE ERROR ******************************************** TEST OF THE VALUES OF THE UNSCALED FUNCTIONS ******************************************** X A ERR(KIA) ERR(KIAP) ERR(LIA) ERR(LIAP) 1 150 0.155431223D-14 0.111022302D-13 0.721644966D-14 0.333066907D-14 2 140 0.111022302D-15 0.544009282D-14 0.222044605D-15 0.111022302D-14 4 130 0.111022302D-14 0.488498131D-14 0.133226763D-14 0.321964677D-14 6 120 0.250910404D-13 0.888178420D-15 0.577315973D-14 0.244249065D-14 8 110 0.247579734D-13 0.164313008D-13 0.355271368D-14 0.843769499D-14 10 100 0.233146835D-14 0.208721929D-13 0.111022302D-14 0.455191440D-14 20 90 0.921485110D-14 0.810462808D-14 0.555111512D-15 0.167643677D-13 30 80 0.233146835D-14 0.440758541D-13 0.111022302D-14 0.754951657D-14 40 70 0.821565038D-14 0.999200722D-15 0.432986980D-14 0.266453526D-14 50 60 0.754951657D-14 0.122124533D-14 0.255351296D-14 0.477395901D-14 60 50 0.139888101D-13 0.421884749D-14 0.643929354D-14 0.555111512D-15 70 40 0.410782519D-14 0.677236045D-14 0.666133815D-14 0.284217094D-13 80 30 0.000000000D+00 0.111022302D-14 0.133226763D-14 0.421884749D-14 90 20 0.168753900D-13 0.378586051D-13 0.821565038D-14 0.532907052D-14 100 10 0.145439216D-13 0.810462808D-14 0.135447209D-13 0.117683641D-13 110 8 0.199840144D-14 0.190958360D-13 0.666133815D-15 0.666133815D-14 120 6 0.532907052D-14 0.610622664D-14 0.666133815D-14 0.621724894D-14 130 4 0.255351296D-14 0.943689571D-14 0.865973959D-14 0.438538095D-13 140 2 0.321964677D-14 0.281996648D-13 0.166533454D-13 0.310862447D-14 150 1 0.244249065D-14 0.488498131D-14 0.102140518D-13 0.444089210D-15 ****************************************** TEST OF THE VALUES OF THE SCALED FUNCTIONS ****************************************** X A ERR(KIA) ERR(KIAP) ERR(LIA) ERR(LIAP) 10 1000 0.222044605D-14 0.432986980D-14 0.162092562D-13 0.133226763D-14 30 600 0.444089210D-14 0.176525461D-13 0.153210777D-13 0.444089210D-15 50 550 0.399680289D-14 0.488498131D-14 0.244249065D-14 0.283106871D-13 70 500 0.777156117D-14 0.421884749D-14 0.262012634D-13 0.107691633D-13 90 450 0.122124533D-14 0.643929354D-14 0.111022302D-15 0.111022302D-14 110 400 0.155431223D-14 0.599520433D-14 0.128785871D-13 0.842548253D-12 130 350 0.268673972D-13 0.999200722D-15 0.355271368D-13 0.111022302D-15 150 300 0.599520433D-14 0.259792188D-13 0.621724894D-14 0.122124533D-13 170 230 0.333066907D-13 0.599520433D-14 0.140998324D-13 0.932587341D-14 190 210 0.186517468D-13 0.250910404D-13 0.643929354D-14 0.365263375D-13 210 190 0.308642001D-13 0.466293670D-14 0.865973959D-14 0.799360578D-14 230 170 0.444089210D-15 0.621724894D-14 0.117683641D-13 0.521804822D-14 300 150 0.155431223D-14 0.599520433D-14 0.145439216D-13 0.732747196D-14 350 130 0.643929354D-14 0.510702591D-14 0.152100554D-13 0.175415238D-13 400 110 0.177635684D-14 0.666133815D-15 0.199840144D-13 0.135447209D-13 450 90 0.444089210D-14 0.888178420D-15 0.102140518D-13 0.147659662D-13 500 70 0.222044605D-15 0.377475828D-14 0.193178806D-13 0.310862447D-14 550 50 0.832667268D-14 0.588418203D-14 0.290878432D-13 0.377475828D-14 600 30 0.910382880D-14 0.421884749D-14 0.177635684D-14 0.133226763D-13 1000 10 0.124344979D-13 0.444089210D-15 0.368594044D-13 0.710542736D-14 ********************************************* TEST OF THE WRONSKIAN RELATION FOR DKIA, DLIA ********************************************* MAXIMUM VALUE OF THE WRONSKIAN CHECK = 6.814548925149211E-13 SHAR_EOF fi # end of overwriting check if test -f 'wktest.f' then echo shar: will not over-write existing file "'wktest.f'" else cat << "SHAR_EOF" > 'wktest.f' PROGRAM WKTEST CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC THIS IS A TEST PROGRAM FOR THE ROUTINES DKIA AND DLIA. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC FIRST, WE COMPARE THE VALUES OBTAINED FROM CCC THE ROUTINES WITH 20(NON-SCALED)+20(SCALED) CCC PRE-COMPUTED VALUES OF THE K_IA(X), K_IA'(X), CCC L_IA(X), L_IA'(X) FUNCTIONS. CCC RELATIVE ERRORS ARE SHOWN. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION A,A1(20),A2(20),ERR1,ERR2,ERR3,ERR4, + KI(20),KIA,KIAP,KIP(20),KIPS(20),KIS(20),LI(20),LIA, + LIAP,LIP(20),LIPS(20),LIS(20),WKILI,WMAX,X,X1(20), + X2(20) INTEGER I,IERR1,IERR2,IFAC,J DATA X1/1.D0,2.D0,4.D0,6.D0,8.D0,10.D0,20.D0,30.D0, + 40.D0,50.D0,60.D0,70.D0,80.D0,90.D0,100.D0,110.D0, + 120.D0,130.D0,140.D0,150.D0/ DATA A1/150.D0,140.D0,130.D0,120.D0,110.D0,100.D0,90.D0, + 80.D0,70.D0,60.D0,50.D0,40.D0,30.D0,20.D0,10.D0, + 8.D0,6.D0,4.D0,2.D0,1.D0/ DATA X2/10.D0,30.D0,50.D0,70.D0,90.D0,110.D0,130.D0,150.D0, + 170.D0,190.D0,210.D0,230.D0,300.D0,350.D0,400.D0, + 450.D0,500.D0,550.D0,600.D0,1000.D0/ DATA A2/1000.D0,600.D0,550.D0,500.D0,450.D0,400.D0, + 350.D0,300.D0,230.D0,210.D0,190.D0,170.D0,150.D0, + 130.D0,110.D0,90.D0,70.D0,50.D0,30.D0,10.D0/ DATA KI/4.6459642814487D-104,-1.9411952078515D-97, + -4.2454441778618D-90,1.8459505138282D-84, + -1.8586595362018D-76,-6.6567482660523D-70, + 3.2023919540906D-63,6.5889532942282D-56, + -3.6683946835532D-49,-4.9527096312644D-42, + 3.4304632381987D-37,5.0365831118869D-37, + 8.8239913220317D-39,1.1752645784103D-41, + 2.8302336906411D-45,1.5088757838133D-49, + 7.5473164457079D-54,3.5955363786557D-58, + 1.6489206017421D-62,7.3120381225095D-67/ DATA KIP/1.2621961347941D-101,-4.4169741448774D-95, + -5.2927180566489D-89,6.2650669669425D-82, + 1.5597363053925D-75,-1.3578086419107D-68, + 4.4928562943487D-62,-1.0038564433216D-55, + -6.5451759692658D-49,-8.8936889262989D-43, + -1.9827493875810D-37,-4.1859420782458D-37, + -8.2438494849568D-39,-1.1527224542051D-41, + -2.8303037451092D-45,-1.5117590940217D-49, + -7.5693363356202D-54,-3.6076494841867D-58, + -1.6546320846520D-62,-7.3362098087182D-67/ DATA LI/-6.0718133124582D+100,1.0341432633441D+94, + 6.0637749124694D+86,-2.6510722091492D+80, + -2.1834147196919D+73,5.9492950832761D+66, + -1.0138118922006D+60,9.1476313015052D+52, + 2.3193283928606D+46,1.3580216912150D+39, + 4.4179418706912D+34,1.7284668020496D+34, + 7.6408995646371D+35,4.8483988671162D+38, + 1.7755623619757D+42,3.0205086366265D+46, + 5.5276887787140D+50,1.0702170837946D+55, + 2.1661539486827D+59,4.5588191419541D+63/ DATA LIP/5.0284086667481D+102,-2.2265454909523D+95, + -5.1327055521212D+88,3.1161042355098D+80, + -4.8930148931703D+74,-2.8872891775895D+67, + 1.3898637806421D+60,3.6652900137178D+53, + -2.6768078097588D+46,-3.7943306091860D+39, + 2.3049351003502D+34,1.3998483121471D+34, + 7.0273838443872D+35,4.6987316485906D+38, + 1.7576708289761D+42,2.9986747348300D+46, + 5.4976358957640D+50,1.0655822048506D+55, + 2.1581811248591D+59,4.5434955896932D+63/ DATA KIS/7.8472200752424D-02,-2.9337081559928D-02, + 8.8250602497608D-02,-3.8547847571986D-02, + 8.9831637526943D-02,-6.4291284802431D-03, + 0.11613787400357D0,-0.12881102538026D0, + 0.13421837026894D0,0.14445491788195D0, + 0.13106736280859D0,1.0045472166660D-01, + 7.7698096235925D-02,6.9491690653934D-02, + 6.3886304416143D-02,5.9669771294809D-02, + 5.6313293442941D-02,5.3540106335192D-02, + 5.1187650274291D-02,3.9629311299016D-02/ DATA KIPS/-1.1207165096700D0,-1.9596302401387D0, + 0.66483564964473D0,0.74865993938252D0, + -0.38512843465564D0,-0.44631514237633D0, + -0.19107259534381D0,0.15077089591635D0, + 0.13730072565776D0,-1.0296790701385D-01, + -5.7446432538801D-02,-6.8138584839609D-02, + -6.7460748750760D-02,-6.4635364385151D-02, + -6.1509430896961D-02,-5.8533212149322D-02, + -5.5816099440864D-02,-5.3367462361202D-02, + -5.1166370921143D-02,-3.9647141469294D-02/ DATA LIS/1.7837655381280D-03,1.5613761576696D-02, + -9.6591519368687D-03,-1.6847454253603D-02, + 1.2512482829322D-02,2.0317388853924D-02, + 1.2170039362900D-02,-1.3867160855031D-02, + -2.3894808167256D-02,3.5321796899890D-02, + 4.2742676164314D-02,3.2131359207386D-02, + 2.4769121304795D-02,2.2141459308109D-02, + 2.0350669922650D-02,1.9004999051914D-02, + 1.7934433775710D-02,1.7050233762627D-02, + 1.6300362441923D-02,1.2617556248941D-02/ DATA LIPS/1.2488614257325D0,-9.3265377462138D-02, + 0.15386015576951D0,-4.3392830231001D-02, + 7.0044344716655D-02,-3.5698138442038D-03, + 4.6211855795097D-02,-3.5524151657274D-02, + 1.9383177095253D-02,1.1257120351664D-02, + 1.7597748590334D-02,2.1486702723200D-02, + 2.1395503168649D-02,2.0520749363341D-02, + 1.9538489291148D-02,1.8599008454981D-02, + 1.7739503408225D-02,1.6963995250943D-02, + 1.6266350796578D-02,1.2610614366489D-02/ OPEN(UNIT=1,FILE='TESTF.DAT',STATUS='UNKNOWN') CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC WE CHECK THE VALUES OF THE UNSCALED FUNCTIONS C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC WRITE(1,*)'ERR MEANS RELATIVE ERROR' WRITE(1,*)'********************************************' WRITE(1,*)'TEST OF THE VALUES OF THE UNSCALED FUNCTIONS' WRITE(1,*)'********************************************' WRITE(1,30)'X','A','ERR(KIA)','ERR(KIAP)','ERR(LIA)','ERR(LIAP)' DO 10 I=1,20 X=X1(I) A=A1(I) IERR1=0 IERR2=0 IFAC=1 CALL DKIA(IFAC,X,A,KIA,KIAP,IERR1) CALL DLIA(IFAC,X,A,LIA,LIAP,IERR2) ERR1=ABS(1.D0-KIA/KI(I)) ERR2=ABS(1.D0-KIAP/KIP(I)) ERR3=ABS(1.D0-LIA/LI(I)) ERR4=ABS(1.D0-LIAP/LIP(I)) WRITE(1,31)INT(X),INT(A),ERR1,ERR2,ERR3,ERR4 10 CONTINUE WRITE(1,*)'******************************************' WRITE(1,*)'TEST OF THE VALUES OF THE SCALED FUNCTIONS' WRITE(1,*)'******************************************' WRITE(1,30)'X','A','ERR(KIA)','ERR(KIAP)','ERR(LIA)','ERR(LIAP)' DO 20 I=1,20 IFAC=2 X=X2(I) A=A2(I) CALL DKIA(IFAC,X,A,KIA,KIAP,IERR1) CALL DLIA(IFAC,X,A,LIA,LIAP,IERR2) ERR1=ABS(1.D0-KIA/KIS(I)) ERR2=ABS(1.D0-KIAP/KIPS(I)) ERR3=ABS(1.D0-LIA/LIS(I)) ERR4=ABS(1.D0-LIAP/LIPS(I)) WRITE(1,31)INT(X),INT(A),ERR1,ERR2,ERR3,ERR4 20 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC TEST OF THE WRONSKIAN RELATION: W[KIA(X),LIA(X)]=1/X CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC WE TEST THE RELATION FOR 0 'airy.f' SUBROUTINE AIZ(IFUN,IFAC,X0,Y0,GAIR,GAII,IERRO) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COMPUTATION OF THE AIRY FUNCTION AI(Z) OR ITS DERIVATIVE AI'(Z) C THE CODE USES: C 1. MACLAURIN SERIES FOR |Y|<3 AND -2.515. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INPUTS: C IFUN: C * IFUN=1, THE CODE COMPUTES AI(Z) C * IFUN=2, THE CODE COMPUTES AI'(Z) C IFAC: C * IFAC=1, THE CODE COMPUTES AI(Z) OR AI'(Z) C * IFAC=2, THE CODE COMPUTES NORMALIZED AI(Z) OR AI'(Z) C X0: REAL PART OF THE ARGUMENT Z C Y0: IMAGINARY PART OF THE ARGUMENT Z C C OUTPUTS: C GAIR: REAL PART OF AI(Z) OR AI'(Z) C GAII: IMAGINARY PART OF AI(Z) OR AI'(Z) C C IERRO: ERROR FLAG C * IERRO=0, SUCCESSFUL COMPUTATION C * IERRO=1, COMPUTATION OUT OF RANGE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C MACHINE DEPENDENT CONSTANTS: FUNCTION D1MACH CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ACCURACY: C C 1) SCALED AIRY FUNCTIONS: C RELATIVE ACCURACY BETTER THAN 10**(-13) EXCEPT CLOSE TO C THE ZEROS, WHERE 10**(-13) IS THE ABSOLUTE PRECISION. C GRADUAL LOSS OF PRECISION TAKES PLACE FOR |Z|>1000 C (REACHING 10**(-8) ABSOLUTE ACCURACY FOR |Z| CLOSE C TO 10**(6)) IN THE CASE OF PHASE(Z) CLOSE TO PI. C 2) UNSCALED AIRY FUNCTIONS: C THE FUNCTION OVERFLOWS/UNDERFLOWS FOR C 3/2*|Z|**(3/2)>LOG(OVER). C FOR |Z|<30: C A) RELATIVE ACCURACY FOR THE MODULUS (EXCEPT AT THE C ZEROS) BETTER THAN 10**(-13). C B) ABSOLUTE ACCURACY FOR MIN(R(Z),1/R(Z)) BETTER C THAN 10**(-13), WHERE R(Z)=REAL(AI)/IMAG(AI) C OR R(Z)=REAL(AI')/IMAG(AI'). C FOR |Z|>30, GRADUAL LOSS OF PRECISION TAKES PLACE C AS |Z| INCREASES. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C AUTHORS: C AMPARO GIL (U. AUTONOMA DE MADRID, MADRID, SPAIN). C E-MAIL: AMPARO.GIL@UAM.ES C JAVIER SEGURA (U. CARLOS III DE MADRID, MADRID, SPAIN). C E-MAIL: JSEGURA@MATH.UC3M.ES C NICO M. TEMME (CWI, AMSTERDAM, THE NETHERLANDS). C E-MAIL: NICO.TEMME@CWI.NL CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C REFERENCES: C COMPUTING AIRY FUNCTIONS BY NUMERICAL QUADRATURE. C NUMERICAL ALGORITHMS (2001). C A. GIL, J. SEGURA, N.M. TEMME CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION X0,Y0,GAIR,GAII,X,W,XD,WD DOUBLE PRECISION OVER,UNDER,DL1,DL2,COVER,D1MACH DOUBLE PRECISION PI,PIHAL,PIH3,PISR,A,ALF,THET,R,TH15,S1,C1, * R32,FACTO,TH025,S3,C3,F23,PI23,SQRT3,XA,YA,F23R,DF1,DF2, * S11,C11,DEX,DRE,DIMA,GAR,GAI,C,S,U,V,V0,AR,AI,AR1,AI1, * RO,COE1,COE2,REX,DFR,DFI,AR11,AI11,PHASE INTEGER IFUN,IFAC,IERRO,IEXPF,IEXPF2,N DIMENSION X(25),W(25) DIMENSION XD(25),WD(25) COMMON/PARAM1/PI,PIHAL COMMON/PARAM2/PIH3,PISR,A,ALF COMMON/PARAM3/THET,R,TH15,S1,C1,R32 COMMON/PARAM4/FACTO,TH025,S3,C3 SAVE X,W SAVE XD,WD DATA X,W/.283891417994567679D-1,.170985378860034935D0, *.435871678341770460D0,.823518257913030858D0,1.33452543254227372D0, *1.96968293206435071D0,2.72998134002859938D0,3.61662161916100897D0, *4.63102611052654146D0,5.77485171830547694D0,7.05000568630218682D0, *8.45866437513237792D0,10.0032955242749393D0,11.6866845947722423D0, *13.5119659344693551D0,15.4826596959377140D0,17.6027156808069112D0, *19.8765656022785451D0,22.3091856773962780D0,24.9061720212974207D0, *27.6738320739497190D0,30.6192963295084111D0,33.7506560850239946D0, *37.0771349708391198D0,40.6093049694341322D0,.143720408803313866D0, *.230407559241880881D0,.242253045521327626D0,.203636639103440807D0, *.143760630622921410D0,.869128834706078120D-1,.4541750018329 * 15883D-1,.206118031206069497D-1,.814278821268606972D-2,.280266 *075663377634D-2,.840337441621719716D-3,.219303732907765020D-3, *.497401659009257760D-4,.978508095920717661D-5,.166542824603725 *563D-5,.244502736801316287D-6,.308537034236207072D-7,.3332960 *72940112245D-8,.306781892316295828D-9,.239331309885375719D-10, *.157294707710054952D-11,.864936011664392267D-13,.394819815 *638647111D-14,.148271173082850884D-15,.453390377327054458D-17/ DATA XD,WD/.435079659953445D-1,.205779160144678D0, *.489916161318751D0,.896390483211727D0,1.42582496737580D0, *2.07903190767599D0,2.85702335104978D0,3.76102058198275D0, *4.79246521225895D0,5.95303247470003D0,7.24464710774066D0, *8.66950223642504D0,10.2300817341775D0,11.9291866622602D0, *13.7699665302828D0,15.7559563095946D0,17.8911203751898D0, *20.1799048700978D0,22.6273004064466D0,25.2389175786164D0, *28.0210785229929D0,30.9809287996116D0,34.1265753192057D0, *37.4672580871163D0,41.0135664833476D0,.576354557898966D-1, *.139560003272262D0,.187792315011311D0,.187446935256946D0, *.150716717316301D0,.101069904453380D0,.575274105486025D-1, *.280625783448681D-1,.117972164134041D-1,.428701743297432D-2, *.134857915232883D-2,.367337337105948D-3,.865882267841931D-4, *.176391622890609D-4,.309929190938078D-5,.468479653648208D-6, *.607273267228907D-7,.672514812555074D-8,.633469931761606D-9, *.504938861248542D-10,.338602527895834D-11,.189738532450555D-12, *.881618802142698D-14,.336676636121976D-15,.104594827170761D-16/ CC CONSTANTS CCCCCCCCCCCCCCCCCCCCCCC PI=3.1415926535897932385D0 PIHAL=1.5707963267948966192D0 PIH3=4.71238898038469D0 F23=.6666666666666666D0 PI23=2.09439510239320D0 PISR=1.77245385090552D0 SQRT3=1.7320508075688772935D0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC YA=Y0 XA=X0 IERRO=0 IEXPF=0 IEXPF2=0 IF (YA.LT.0.D0) YA=-YA R=SQRT(XA*XA+YA*YA) R32=R*SQRT(R) THET=PHASE(XA,YA) COVER=2.D0/3.D0*R32*ABS(COS(1.5D0*THET)) CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER) OVER=D1MACH(2)*1.D-3 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER) UNDER=D1MACH(1)*1.D+3 DL1=LOG(OVER) DL2=-LOG(UNDER) IF (DL1.GT.DL2) OVER=1/UNDER IF (IFAC.EQ.1) THEN IF (COVER.GE.LOG(OVER)) THEN CCC OVERFLOW/UNDERFLOW PROBLEMS. CCC CALCULATION ABORTED IERRO=1 GAIR=0 GAII=0 ENDIF IF (COVER.GE.(LOG(OVER)*0.2)) IEXPF2=1 ELSE IF (COVER.GE.(LOG(OVER)*0.2)) IEXPF=1 ENDIF IF (IERRO.EQ.0) THEN IF (IFUN.EQ.1) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC CALCULATION OF AI(Z) CCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCC SERIES, INTEGRALS OR EXPANSIONS CCCCCCCCCCCCCCCCCCCCCCCC IF ((YA.LT.3.D0).AND.(XA.LT.1.3D0).AND.(XA.GT.-2.5D0)) THEN CCC SERIES CCC CALL SERAI(XA,YA,GAR,GAI) IF (IFAC.EQ.2) THEN THET=PHASE(XA,YA) TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) F23R=F23*R32 DF1=F23R*C1 DF2=F23R*S1 S11=SIN(DF2) C11=COS(DF2) DEX=EXP(DF1) DRE=DEX*C11 DIMA=DEX*S11 GAIR=DRE*GAR-DIMA*GAI GAII=DRE*GAI+DIMA*GAR ELSE GAIR=GAR GAII=GAI IF (Y0.EQ.0.) GAII=0.D0 ENDIF ELSE IF (R.GT.15.D0) THEN CCC ASYMPTOTIC EXPANSIONS CCC THET=PHASE(XA,YA) FACTO=0.5D0/PISR*R**(-0.25D0) IF (THET.GT.PI23) THEN CCCCCCCCCCC CONNECTION FORMULAE CCCCCCCCCCCCCCCCCCCCCCCCCC CCC N= 1: TRANSFORM Z TO W= U+IV=Z EXP( 2 PI I/3) N=1 C=-0.5D0 S=N*0.5*SQRT3 U=XA*C-YA*S V=XA*S+YA*C V0=V IF (V.LT.0.D0) V=-V THET=PHASE(U,V) TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL EXPAI(AR1,AI1) IF (V0.LT.0.D0) AI1=-AI1 AR=-(C*AR1-S*AI1) AI=-(S*AR1+C*AI1) IF (IEXPF.EQ.0) THEN IF (IEXPF2.EQ.0) THEN CCC N=-1: TRANSFORM Z TO W= U+IV=Z EXP(-2 PI I/3) N=-1 C=-0.5D0 S=N*0.5*SQRT3 U=XA*C-YA*S V=XA*S+YA*C V0=V IF (V.LT.0.D0) V=-V THET=PHASE(U,V) TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL EXPAI(AR1,AI1) IF (V0.LT.0.D0) AI1=-AI1 THET=PHASE(XA,YA) TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) RO=1.333333333333333D0*R32 COE1=RO*C1 COE2=RO*S1 REX=EXP(COE1) DFR=REX*COS(COE2) DFI=REX*SIN(COE2) AR11=DFR*AR1-DFI*AI1 AI11=DFR*AI1+DFI*AR1 GAIR=AR-(C*AR11-S*AI11) GAII=AI-(S*AR11+C*AI11) ELSE THET=PHASE(XA,YA) TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) GAIR=AR GAII=AI ENDIF ELSE GAIR=AR GAII=AI ENDIF ELSE CCCCCCC ASYMPTOTIC EXPANSION CCCCCCCCCCCCCCC THET=PHASE(XA,YA) TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL EXPAI(GAIR,GAII) ENDIF ELSE CCC INTEGRALS A=0.1666666666666666D0 ALF=-A FACTO=0.280514117723058D0*R**(-0.25D0) THET=PHASE(XA,YA) IF (THET.LE.PIHAL) THEN TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL AIRY1(X,W,GAIR,GAII) ENDIF IF ((THET.GT.PIHAL).AND.(THET.LE.PI23)) THEN TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL AIRY2(X,W,GAIR,GAII) ENDIF IF (THET.GT.PI23) THEN N=1 C=-0.5D0 S=N*0.5*SQRT3 U=XA*C-YA*S V=XA*S+YA*C V0=V IF (V.LT.0.D0) V=-V THET=PHASE(U,V) IF (THET.LE.PIHAL) THEN TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL AIRY1(X,W,AR1,AI1) ENDIF IF ((THET.GT.PIHAL).AND.(THET.LE.PI23)) THEN TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL AIRY2(X,W,AR1,AI1) ENDIF IF (V0.LT.0.D0) AI1=-AI1 AR=-(C*AR1-S*AI1) AI=-(S*AR1+C*AI1) N=-1 C=-0.5D0 S=N*0.5*SQRT3 U=XA*C-YA*S V=XA*S+YA*C V0=V IF (V.LT.0.D0) V=-V THET=PHASE(U,V) IF (THET.LE.PIHAL) THEN TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL AIRY1(X,W,AR1,AI1) ENDIF IF ((THET.GT.PIHAL).AND.(THET.LE.PI23)) THEN TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL AIRY2(X,W,AR1,AI1) ENDIF IF (V0.LT.0.D0) AI1=-AI1 THET=PHASE(XA,YA) TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) RO=1.333333333333333D0*R32 COE1=RO*C1 COE2=RO*S1 REX=EXP(COE1) DFR=REX*COS(COE2) DFI=REX*SIN(COE2) AR11=DFR*AR1-DFI*AI1 AI11=DFR*AI1+DFI*AR1 GAIR=AR-(C*AR11-S*AI11) GAII=AI-(S*AR11+C*AI11) ENDIF ENDIF IF (IFAC.EQ.1) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC CALCULATION OF THE UNSCALED AI(Z) CCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC F23R=F23*R32 DF1=F23R*C1 DF2=F23R*S1 S11=SIN(DF2) C11=COS(DF2) DEX=EXP(-DF1) DRE=DEX*C11 DIMA=-DEX*S11 GAR=DRE*GAIR-DIMA*GAII GAI=DRE*GAII+DIMA*GAIR GAIR=GAR GAII=GAI IF (Y0.EQ.0.) GAII=0.D0 ENDIF ENDIF ELSE CCCC CALCULATION OF AIŽ(Z) CCCCCCCCCCC ALF=0.1666666666666666D0 FACTO=-0.270898621247918D0*R**0.25D0 CCCCCCCCCCCCCCC SERIES OR INTEGRALS CCCCCCCCCCCCCCCCCCCCCCCCCC IF ((YA.LT.3.D0).AND.(XA.LT.1.3D0).AND.(XA.GT.-2.5D0)) THEN CCC SERIES CALL SERAID(XA,YA,GAR,GAI) IF (IFAC.EQ.2) THEN THET=PHASE(XA,YA) TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) F23R=F23*R32 DF1=F23R*C1 DF2=F23R*S1 S11=SIN(DF2) C11=COS(DF2) DEX=EXP(DF1) DRE=DEX*C11 DIMA=DEX*S11 GAIR=DRE*GAR-DIMA*GAI GAII=DRE*GAI+DIMA*GAR ELSE GAIR=GAR GAII=GAI IF (Y0.EQ.0.) GAII=0.D0 ENDIF ELSE IF (R.GT.15.D0) THEN CCC ASYMPTOTIC EXPANSIONS CCCCCCCCCCCCC THET=PHASE(XA,YA) FACTO=0.5D0/PISR*R**0.25D0 IF (THET.GT.PI23) THEN CCCCCCC CONNECTION FORMULAE CCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC N= 1: TRANSFORM Z TO W= U+IV=Z EXP( 2 PI I/3) N=1 C=-0.5D0 S=N*0.5*SQRT3 U=XA*C-YA*S V=XA*S+YA*C V0=V IF (V.LT.0.D0) V=-V THET=PHASE(U,V) TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL EXPAID(AR1,AI1) IF (V0.LT.0.D0) AI1=-AI1 AR=-(C*AR1+S*AI1) AI=-(-S*AR1+C*AI1) IF (IEXPF.EQ.0) THEN IF (IEXPF2.EQ.0) THEN CCC N=-1: TRANSFORM Z TO W= U+IV=Z EXP(-2 PI I/3) N=-1 C=-0.5D0 S=N*0.5*SQRT3 U=XA*C-YA*S V=XA*S+YA*C V0=V IF (V.LT.0.D0) V=-V THET=PHASE(U,V) TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL EXPAID(AR1,AI1) IF (V0.LT.0.D0) AI1=-AI1 THET=PHASE(XA,YA) TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) RO=1.333333333333333D0*R32 COE1=RO*C1 COE2=RO*S1 REX=EXP(COE1) DFR=REX*COS(COE2) DFI=REX*SIN(COE2) AR11=DFR*AR1-DFI*AI1 AI11=DFR*AI1+DFI*AR1 GAIR=AR-(C*AR11+S*AI11) GAII=AI-(-S*AR11+C*AI11) ELSE THET=PHASE(XA,YA) TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) GAIR=AR GAII=AI ENDIF ELSE GAIR=AR GAII=AI ENDIF ELSE TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL EXPAID(GAIR,GAII) ENDIF ELSE CCC INTEGRALS CCCCCCCCCCCCCCCC THET=PHASE(XA,YA) IF (THET.LE.PIHAL) THEN TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL AIRY1D(XD,WD,GAIR,GAII) ENDIF IF ((THET.GT.PIHAL).AND.(THET.LE.PI23)) THEN TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL AIRY2D(XD,WD,GAIR,GAII) ENDIF IF (THET.GT.PI23) THEN N=1 C=-0.5D0 S=N*0.5*SQRT3 U=XA*C-YA*S V=XA*S+YA*C V0=V IF (V.LT.0.D0) V=-V THET=PHASE(U,V) IF (THET.LE.PIHAL) THEN TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL AIRY1D(XD,WD,AR1,AI1) ENDIF IF ((THET.GT.PIHAL).AND.(THET.LE.PI23)) THEN TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL AIRY2D(XD,WD,AR1,AI1) ENDIF IF (V0.LT.0.D0) AI1=-AI1 AR=-(C*AR1+S*AI1) AI=-(-S*AR1+C*AI1) N=-1 C=-0.5D0 S=N*0.5*SQRT3 U=XA*C-YA*S V=XA*S+YA*C V0=V IF (V.LT.0.D0) V=-V THET=PHASE(U,V) IF (THET.LE.PIHAL) THEN TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL AIRY1D(XD,WD,AR1,AI1) ENDIF IF ((THET.GT.PIHAL).AND.(THET.LE.PI23)) THEN TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) TH025=THET*0.25D0 S3=SIN(TH025) C3=COS(TH025) CALL AIRY2D(XD,WD,AR1,AI1) ENDIF IF (V0.LT.0.D0) AI1=-AI1 THET=PHASE(XA,YA) TH15=1.5D0*THET S1=SIN(TH15) C1=COS(TH15) RO=1.333333333333333D0*R32 COE1=RO*C1 COE2=RO*S1 REX=EXP(COE1) DFR=REX*COS(COE2) DFI=REX*SIN(COE2) AR11=DFR*AR1-DFI*AI1 AI11=DFR*AI1+DFI*AR1 GAIR=AR-(C*AR11+S*AI11) GAII=AI-(-S*AR11+C*AI11) ENDIF ENDIF IF (IFAC.EQ.1) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC CALCULATION OF THE UNSCALED AI'(z) CCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC F23R=F23*R32 DF1=F23R*C1 DF2=F23R*S1 S11=SIN(DF2) C11=COS(DF2) DEX=EXP(-DF1) DRE=DEX*C11 DIMA=-DEX*S11 GAR=DRE*GAIR-DIMA*GAII GAI=DRE*GAII+DIMA*GAIR GAIR=GAR GAII=GAI IF (Y0.EQ.0) GAII=0.D0 ENDIF ENDIF ENDIF ENDIF IF (Y0.LT.0.D0) GAII=-GAII RETURN END SUBROUTINE AIRY1(X,W,GAIR,GAII) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC COMPUTES AI(Z) BY GAUSS-LAGUERRE QUADRATURE IN THE SECTOR C CCC 0 <= PHASE(Z) <= PI/2 C CCC C CCC INPUTS: C CCC X,W, NODES AND WEIGHTS FOR THE GAUSSIAN C CCC QUADRATURE C CCC OUTPUTS: C CCC GAIR, GAII, REAL AND IMAGINARY PARTS OF AI(Z) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION X,W,GAIR,GAII DOUBLE PRECISION PIH3,PISR,A,ALF,THET,R,TH15,S1,C1, * R32,FACTO,TH025,S3,C3,SUMAR,SUMAI,DF1,DF1C1,PHI,PHI6, * S2,C2,DMODU,DMODU2,FUNR,FUNI,FAC1,FAC2,PHASE INTEGER I DIMENSION X(25),W(25) COMMON/PARAM2/PIH3,PISR,A,ALF COMMON/PARAM3/THET,R,TH15,S1,C1,R32 COMMON/PARAM4/FACTO,TH025,S3,C3 SUMAR=0.D0 SUMAI=0.D0 DO 1 I=1,25 DF1=1.5D0*X(I)/R32 DF1C1=DF1*C1 PHI=PHASE(2.D0+DF1C1,DF1*S1) PHI6=PHI/6.D0 S2=SIN(PHI6) C2=COS(PHI6) DMODU=SQRT(4.D0+DF1*DF1+4.D0*DF1C1) DMODU2=DMODU**ALF FUNR=DMODU2*C2 FUNI=DMODU2*S2 SUMAR=SUMAR+W(I)*FUNR SUMAI=SUMAI+W(I)*FUNI 1 CONTINUE FAC1=FACTO*C3 FAC2=FACTO*S3 GAIR=FAC1*SUMAR+FAC2*SUMAI GAII=FAC1*SUMAI-FAC2*SUMAR RETURN END SUBROUTINE AIRY2(X,W,GAIR,GAII) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC COMPUTES AI(Z) BY GAUSS-LAGUERRE QUADRATURE IN THE SECTOR C CCC PI/2 < PHASE(Z) <= 2PI/3 C CCC C CCC INPUTS: C CCC X,W, NODES AND WEIGHTS FOR THE GAUSSIAN C CCC QUADRATURE C CCC OUTPUTS: C CCC GAIR, GAII, REAL AND IMAGINARY PARTS OF AI(Z) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION X,W,GAIR,GAII DOUBLE PRECISION PIH3,PISR,A,ALF,THET,R,TH15,S1,C1, * R32,FACTO,TH025,S3,C3,SUMAR,SUMAI,DF1,DF1C1,PHI,PHI6, * S2,C2,DMODU,DMODU2,FUNR,FUNI,FAC1,FAC2,PHASE DOUBLE PRECISION SQR2,SQR2I,TAU,TGTAU,B,ANG,CTAU,CFAC,CT,ST, * SUMR,SUMI,TTAU,BETA INTEGER I DIMENSION X(25),W(25) COMMON/PARAM2/PIH3,PISR,A,ALF COMMON/PARAM3/THET,R,TH15,S1,C1,R32 COMMON/PARAM4/FACTO,TH025,S3,C3 SQR2=1.41421356237310D0 SQR2I=0.707106781186548D0 TAU=TH15-PIH3*0.5D0 TGTAU=DTAN(TAU) B=5.D0*A ANG=TAU*B CTAU=COS(TAU) CFAC=CTAU**(-B) CT=COS(ANG) ST=SIN(ANG) SUMR=0.D0 SUMI=0.D0 DO 2 I=1,25 DF1=3.D0*X(I)/(CTAU*R32) DF1C1=DF1*SQR2I*0.5D0 PHI=PHASE(2.D0-DF1C1,DF1C1) PHI6=PHI/6.D0 TTAU=X(I)*TGTAU BETA=PHI6-TTAU S2=SIN(BETA) C2=COS(BETA) DMODU=SQRT(4.D0+DF1*DF1*0.25D0-SQR2*DF1) DMODU2=DMODU**ALF FUNR=DMODU2*C2 FUNI=DMODU2*S2 SUMR=SUMR+W(I)*FUNR SUMI=SUMI+W(I)*FUNI 2 CONTINUE SUMAR=CFAC*(CT*SUMR-ST*SUMI) SUMAI=CFAC*(CT*SUMI+ST*SUMR) FAC1=FACTO*C3 FAC2=FACTO*S3 GAIR=FAC1*SUMAR+FAC2*SUMAI GAII=FAC1*SUMAI-FAC2*SUMAR RETURN END SUBROUTINE AIRY1D(X,W,GAIR,GAII) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC COMPUTES AI'(Z) BY GAUSS-LAGUERRE QUADRATURE IN THE SECTOR C CCC 0 <= PHASE(Z) <= PI/2 C CCC C CCC INPUTS: C CCC X,W, NODES AND WEIGHTS FOR THE GAUSSIAN C CCC QUADRATURE C CCC OUTPUTS: C CCC GAIR,GAII, REAL AND IMAGINARY PARTS OF AI'(Z) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION X,W,GAIR,GAII DOUBLE PRECISION PIH3,PISR,A,ALF,THET,R,TH15,S1,C1, * R32,FACTO,TH025,S3,C3,SUMAR,SUMAI,DF1,DF1C1,PHI,PHI6, * S2,C2,DMODU,DMODU2,FUNR,FUNI,FAC1,FAC2,PHASE INTEGER I DIMENSION X(25),W(25) COMMON/PARAM2/PIH3,PISR,A,ALF COMMON/PARAM3/THET,R,TH15,S1,C1,R32 COMMON/PARAM4/FACTO,TH025,S3,C3 SUMAR=0.D0 SUMAI=0.D0 DO 3 I=1,25 DF1=1.5D0*X(I)/R32 DF1C1=DF1*C1 PHI=PHASE(2.D0+DF1C1,DF1*S1) PHI6=-PHI*ALF S2=SIN(PHI6) C2=COS(PHI6) DMODU=SQRT(4.D0+DF1*DF1+4.D0*DF1C1) DMODU2=DMODU**ALF FUNR=DMODU2*C2 FUNI=DMODU2*S2 SUMAR=SUMAR+W(I)*FUNR SUMAI=SUMAI+W(I)*FUNI 3 CONTINUE FAC1=FACTO*C3 FAC2=FACTO*S3 GAIR=FAC1*SUMAR-FAC2*SUMAI GAII=FAC1*SUMAI+FAC2*SUMAR RETURN END SUBROUTINE AIRY2D(X,W,GAIR,GAII) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC COMPUTES AI'(Z) BY GAUSS-LAGUERRE QUADRATURE IN THE SECTOR C CCC PI/2 < PHASE(Z) <= 3PI/2 C CCC C CCC INPUTS: C CCC X,W, NODES AND WEIGHTS FOR THE GAUSSIAN C CCC QUADRATURE C CCC OUTPUTS: C CCC GAIR,GAII, REAL AND IMAGINARY PARTS OF AI'(Z) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION X,W,GAIR,GAII DOUBLE PRECISION PIH3,PISR,A,ALF,THET,R,TH15,S1,C1, * R32,FACTO,TH025,S3,C3,SUMAR,SUMAI,DF1,DF1C1,PHI,PHI6, * S2,C2,DMODU,DMODU2,FUNR,FUNI,FAC1,FAC2,PHASE DOUBLE PRECISION SQR2,SQR2I,TAU,TGTAU,B,ANG,CTAU,CFAC,CT,ST, * SUMR,SUMI,TTAU,BETA INTEGER I DIMENSION X(25),W(25) COMMON/PARAM2/PIH3,PISR,A,ALF COMMON/PARAM3/THET,R,TH15,S1,C1,R32 COMMON/PARAM4/FACTO,TH025,S3,C3 SQR2=1.41421356237310D0 SQR2I=0.707106781186548D0 TAU=TH15-PIH3*0.5D0 TGTAU=DTAN(TAU) B=7.D0*ALF ANG=TAU*B CTAU=COS(TAU) CFAC=CTAU**(-B) CT=COS(ANG) ST=SIN(ANG) SUMR=0.D0 SUMI=0.D0 DO 4 I=1,25 DF1=3.D0*X(I)/(CTAU*R32) DF1C1=DF1*SQR2I*0.5D0 PHI=PHASE(2.D0-DF1C1,DF1C1) PHI6=-PHI/6.D0 TTAU=X(I)*TGTAU BETA=PHI6-TTAU S2=SIN(BETA) C2=COS(BETA) DMODU=SQRT(4.D0+DF1*DF1*0.25D0-SQR2*DF1) DMODU2=DMODU**ALF FUNR=DMODU2*C2 FUNI=DMODU2*S2 SUMR=SUMR+W(I)*FUNR SUMI=SUMI+W(I)*FUNI 4 CONTINUE SUMAR=CFAC*(CT*SUMR-ST*SUMI) SUMAI=CFAC*(CT*SUMI+ST*SUMR) FAC1=FACTO*C3 FAC2=FACTO*S3 GAIR=FAC1*SUMAR-FAC2*SUMAI GAII=FAC1*SUMAI+FAC2*SUMAR RETURN END DOUBLE PRECISION FUNCTION PHASE(X,Y) DOUBLE PRECISION PI,PIHAL,X,Y,AY,P CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC COMPUTES THE PHASE OF Z = X + IY, IN (-PI,PI] CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON/PARAM1/PI,PIHAL IF ((X.EQ.0).AND.(Y.EQ.0)) THEN P=0.D0 ELSE AY=ABS(Y) IF (X.GE.AY) THEN P=ATAN(AY/X) ELSEIF ((X+AY).GE.0.D0) THEN P=PIHAL-ATAN(X/AY) ELSE P=PI+ATAN(AY/X) ENDIF IF (Y.LT.0.D0) P=-P ENDIF PHASE=P END SUBROUTINE FGP(X,Y,EPS,FR,FI,GR,GI) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COMPUTES THE FUNCTIONS F AND G FOR THE SERIES C C OF AI'(Z). C C THIS ROUTINE IS CALLED BY SERAID. C C C C INPUTS: C C X,Y, REAL AND IMAGINARY PARTS OF Z C C EPS, PRECISION FOR THE COMPUTATION OF C C THE SERIES C C OUTPUTS: C C FR,FI, REAL AND IMAGINARY PARTS OF F C C GR,GI, REAL AND IMAGINARY PARTS OF G C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION X,Y,EPS,FR,FI,GR,GI INTEGER A,B,K3 DOUBLE PRECISION X2,Y2,U,V,P,Q,CR,CI,DR,DI X2=X*X Y2=Y*Y K3=0 U=X*(X2-3*Y2) V=Y*(3*X2-Y2) CR=0.5D0 CI=0.D0 DR=1.D0 DI=0.D0 FR=0.5D0 FI=0.D0 GR=1.D0 GI=0.D0 70 A=(K3+5)*(K3+3) B=(K3+1)*(K3+3) P=(U*CR-V*CI)/A Q=(V*CR+U*CI)/A CR=P CI=Q P=(U*DR-V*DI)/B Q=(V*DR+U*DI)/B DR=P DI=Q FR=FR+CR FI=FI+CI GR=GR+DR GI=GI+DI K3=K3+3 IF ((ABS(CR)+ABS(DR)+ABS(CI)+ABS(DI)).GE.EPS) GOTO 70 U=X2-Y2 V=2.D0*X*Y P=U*FR-V*FI Q=U*FI+V*FR FR=P FI=Q RETURN END SUBROUTINE FG(X,Y,EPS,FR,FI,GR,GI) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COMPUTES THE FUNCTIONS F AND G IN EXPRESSION C C 10.4.2 OF ABRAMOWITZ & STEGUN FOR THE SERIES C C OF AI(Z). C C THIS ROUTINE IS CALLED BY SERAI. C C C C INPUTS: C C X,Y, REAL AND IMAGINARY PARTS OF Z C C EPS, PRECISION FOR THE COMPUTATION C C OF THE SERIES. C C OUTPUTS: C C FR,FI, REAL AND IMAGINARY PARTS OF F C C GR,GI, REAL AND IMAGINARY PARTS OF G C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER A,B,K3 DOUBLE PRECISION X2,Y2,U,V,P,Q,CR,CI,DR,DI DOUBLE PRECISION X,Y,EPS,FR,FI,GR,GI X2=X*X Y2=Y*Y K3=0 U=X*(X2-3.D0*Y2) V=Y*(3.D0*X2-Y2) CR=1.D0 CI=0.D0 DR=1.D0 DI=0.D0 FR=1.D0 FI=0.D0 GR=1.D0 GI=0.D0 71 A=(K3+2)*(K3+3) B=(K3+4)*(K3+3) P=(U*CR-V*CI)/A Q=(V*CR+U*CI)/A CR=P CI=Q P=(U*DR-V*DI)/B Q=(V*DR+U*DI)/B DR=P DI=Q FR=FR+CR FI=FI+CI GR=GR+DR GI=GI+DI K3=K3+3 IF ((ABS(CR)+ABS(DR)+ABS(CI)+ABS(DI)).GE.EPS) GOTO 71 P=X*GR-Y*GI Q=X*GI+Y*GR GR=P GI=Q RETURN END SUBROUTINE SERAI(X,Y,AIR,AII) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC AIRY AI(Z), TAYLOR, COMPLEX Z CCC CCC CCC CCC INPUTS: CCC CCC X,Y, REAL AND IMAGINARY CCC CCC PARTS OF Z CCC CCC OUTPUTS: CCC CCC AIR,AII, REAL AND IMAGINARY CCC CCC PARTS OF AI(Z) CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION X,Y,EPS,AIR,AII DOUBLE PRECISION FZR,FZI,GZR,GZI,CONS1,CONS2 DOUBLE PRECISION D1MACH EPS=D1MACH(3) IF (EPS.LT.1.D-15) EPS=1.D-15 CONS1=0.355028053887817239260D0 CONS2=0.258819403792806798405D0 CALL FG(X,Y,EPS,FZR,FZI,GZR,GZI) AIR=CONS1*FZR-CONS2*GZR AII=CONS1*FZI-CONS2*GZI RETURN END SUBROUTINE SERAID(X,Y,AIR,AII) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC AIRY AI'(Z), TAYLOR, COMPLEX Z CCC CCC CCC CCC INPUTS: CCC CCC X,Y, REAL AND IMAGINARY CCC CCC PARTS OF Z CCC CCC OUTPUTS: CCC CCC AIR,AII, REAL AND IMAGINARY CCC CCC PARTS OF AI'(Z) CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION X,Y,EPS,AIR,AII DOUBLE PRECISION FZR,FZI,GZR,GZI,CONS1,CONS2 DOUBLE PRECISION D1MACH EPS=D1MACH(3) IF (EPS.LT.1.D-15) EPS=1.D-15 CONS1=0.355028053887817239260D0 CONS2=0.258819403792806798405D0 CALL FGP(X,Y,EPS,FZR,FZI,GZR,GZI) AIR=CONS1*FZR-CONS2*GZR AII=CONS1*FZI-CONS2*GZI RETURN END SUBROUTINE EXPAI(GAIR,GAII) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC AIRY AI(Z), ASYMPTOTIC EXPANSION, COMPLEX Z CCC CCC CCC CCC OUTPUTS: CCC CCC GAIR, GAII, REAL AND IMAGINARY CCC CCC PARTS OF AI(Z) CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION EPS,GAIR,GAII DOUBLE PRECISION THET,R,TH15,S1, * C1,R32,FACTO,TH025,S3,C3 DOUBLE PRECISION DF1,PSIIR,PSIII,CK,DFRR,DFII,SUMAR,SUMAI, * DFR,DFI,DELTAR,DELTAI,FAC1,FAC2 DOUBLE PRECISION CO,DF DOUBLE PRECISION D1MACH INTEGER K DIMENSION CO(20) COMMON/PARAM3/THET,R,TH15,S1,C1,R32 COMMON/PARAM4/FACTO,TH025,S3,C3 SAVE CO DATA CO/-6.944444444444445D-2,3.713348765432099D-2, * -3.799305912780064D-2,5.764919041266972D-2,-0.116099064025515D0, * 0.291591399230751D0,-0.877666969510017D0,3.07945303017317D0, * -12.3415733323452D0,55.6227853659171D0,-278.465080777603D0, * 1533.16943201280D0,-9207.20659972641D0,59892.5135658791D0, * -419524.875116551D0,3148257.41786683D0,-25198919.8716024D0, * 214288036.963680D0,-1929375549.18249D0,18335766937.8906D0/ EPS=D1MACH(3) IF (EPS.LT.1.D-15) EPS=1.D-15 DF1=1.5D0/R32 PSIIR=C1 PSIII=-S1 K=0 CK=1.D0 DF=1.D0 DFRR=1.D0 DFII=0.D0 SUMAR=1.D0 SUMAI=0.D0 80 DF=DF*DF1 CK=CO(K+1)*DF DFR=DFRR DFI=DFII DFRR=DFR*PSIIR-DFI*PSIII DFII=DFR*PSIII+DFI*PSIIR DELTAR=DFRR*CK DELTAI=DFII*CK SUMAR=SUMAR+DELTAR SUMAI=SUMAI+DELTAI K=K+1 IF (SUMAR.NE.0) THEN IF (ABS(DELTAR/SUMAR).GT.EPS) GOTO 80 ENDIF IF (SUMAI.NE.0) THEN IF (ABS(DELTAI/SUMAI).GT.EPS) GOTO 80 ENDIF FAC1=FACTO*C3 FAC2=FACTO*S3 GAIR=FAC1*SUMAR+FAC2*SUMAI GAII=FAC1*SUMAI-FAC2*SUMAR RETURN END SUBROUTINE EXPAID(GAIR,GAII) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC AIRY AI'(Z), ASYMPTOTIC EXPANSION, COMPLEX Z CCC CCC CCC CCC OUTPUTS: CCC CCC GAIR, GAII, REAL AND IMAGINARY CCC CCC PARTS OF AI'(Z) CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION EPS,GAIR,GAII DOUBLE PRECISION THET,R,TH15,S1, * C1,R32,FACTO,TH025,S3,C3 DOUBLE PRECISION DF1,PSIIR,PSIII,VK,DFRR,DFII,SUMAR,SUMAI, * DFR,DFI,DELTAR,DELTAI,FAC1,FAC2 DOUBLE PRECISION CO,DF DOUBLE PRECISION D1MACH INTEGER K DIMENSION CO(20) COMMON/PARAM3/THET,R,TH15,S1,C1,R32 COMMON/PARAM4/FACTO,TH025,S3,C3 SAVE CO DATA CO/9.722222222222222D-2,-4.388503086419753D-2, * 4.246283078989484D-2,-6.266216349203230D-2, * 0.124105896027275D0,-0.308253764901079D0, * 0.920479992412945D0,-3.21049358464862D0, * 12.8072930807356D0,-57.5083035139143D0, * 287.033237109221D0,-1576.35730333710D0, * 9446.35482309593D0,-61335.7066638521D0, * 428952.400400069D0,-3214536.52140086D0, * 25697908.3839113D0,-218293420.832160D0, * 1963523788.99103D0,-18643931088.1072D0/ EPS=D1MACH(3) IF (EPS.LT.1.D-15) EPS=1.D-15 DF1=1.5D0/R32 PSIIR=C1 PSIII=-S1 K=0 DF=1.D0 DFRR=1.D0 DFII=0.D0 SUMAR=1.D0 SUMAI=0.D0 81 DF=DF*DF1 VK=CO(K+1)*DF DFR=DFRR DFI=DFII DFRR=DFR*PSIIR-DFI*PSIII DFII=DFR*PSIII+DFI*PSIIR DELTAR=DFRR*VK DELTAI=DFII*VK SUMAR=SUMAR+DELTAR SUMAI=SUMAI+DELTAI K=K+1 IF (SUMAR.NE.0) THEN IF (ABS(DELTAR/SUMAR).GT.EPS) GOTO 81 ENDIF IF (SUMAI.NE.0) THEN IF (ABS(DELTAI/SUMAI).GT.EPS) GOTO 81 ENDIF FAC1=FACTO*C3 FAC2=FACTO*S3 GAIR=-(FAC1*SUMAR-FAC2*SUMAI) GAII=-(FAC1*SUMAI+FAC2*SUMAR) RETURN END SUBROUTINE BIZ(IFUN,IFAC,X0,Y0,GBIR,GBII,IERRO) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COMPUTATION OF THE AIRY FUNCTION BI(Z) OR ITS DERIVATIVE BI'(Z) C THE CODE USES THE CONNECTION OF BI(Z) WITH AI(Z). C BIZ CALLS THE ROUTINE AIZ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INPUTS: C IFUN: C * IFUN=1, THE CODE COMPUTES BI(Z) C * IFUN=2, THE CODE COMPUTES BI'(Z) C IFAC: C * IFAC=1, THE CODE COMPUTES BI(Z) OR BI'(Z) C * IFAC=2, THE CODE COMPUTES NORMALIZED BI(Z) OR BI'(Z) C X0: REAL PART OF THE ARGUMENT Z C Y0: IMAGINARY PART OF THE ARGUMENT Z C C OUTPUTS: C GBIR: REAL PART OF BI(Z) OR BI'(Z) C GBII: IMAGINARY PART OF BI(Z) OR BI'(Z) C C IERRO: ERROR FLAG C * IERRO=0, SUCCESSFUL COMPUTATION C * IERRO=1, COMPUTATION OUT OF RANGE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C MACHINE DEPENDENT CONSTANTS: FUNCTION D1MACH CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ACCURACY: C C 1) SCALED AIRY FUNCTIONS: C RELATIVE ACCURACY BETTER THAN 10**(-13) EXCEPT CLOSE TO C THE ZEROS, WHERE 10**(-13) IS THE ABSOLUTE PRECISION. C GRADUAL LOSS OF PRECISION TAKES PLACE FOR |Z|>1000 C IN THE CASE OF PHASE(Z) CLOSE TO +3*PI/2 OR -3*PI/2. C 2) UNSCALED AIRY FUNCTIONS: C THE FUNCTION OVERFLOWS/UNDERFLOWS FOR C 3/2*|Z|**(3/2)>LOG(OVER). C FOR |Z|<30: C A) RELATIVE ACCURACY FOR THE MODULUS (EXCEPT AT THE C ZEROS) BETTER THAN 10**(-13). C B) ABSOLUTE ACCURACY FOR MIN(R(Z),1/R(Z)) BETTER C THAN 10**(-13), WHERE R(Z)=REAL(BI)/IMAG(BI) C OR R(Z)=REAL(BI')/IMAG(BI'). C FOR |Z|>30, GRADUAL LOSS OF PRECISION TAKES PLACE C AS |Z| INCREASES. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C AUTHORS: C AMPARO GIL (U. AUTONOMA DE MADRID, MADRID, SPAIN). C E-MAIL: AMPARO.GIL@UAM.ES C JAVIER SEGURA (U. CARLOS III DE MADRID, MADRID, SPAIN). C E-MAIL: JSEGURA@MATH.UC3M.ES C NICO M. TEMME (CWI, AMSTERDAM, THE NETHERLANDS). C E-MAIL: NICO.TEMME@CWI.NL CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C REFERENCES: C COMPUTING AIRY FUNCTIONS BY NUMERICAL QUADRATURE. C NUMERICAL ALGORITHMS (2001). C A. GIL, J. SEGURA, N.M. TEMME CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION X0,Y0,GBIR,GBII DOUBLE PRECISION OVER,UNDER,DL1,DL2,COVER,D1MACH DOUBLE PRECISION PI,PI3,PI23,SQRT3,C,S,C1,S1,U,V,X,Y,AR,AI, * APR,API,BR,BI,BPR,BPI,BBR,BBI,BBPR,BBPI,PHASE DOUBLE PRECISION THET,R,R32,THET32,A1,B1,DF1,EXPO,EXPOI, * ETAR,ETAI,ETAGR,ETAGI,PIHAL INTEGER IFUN,IFAC,IEXPF,IERR,IERRO COMMON/PARAM1/PI,PIHAL SQRT3=1.7320508075688772935D0 PI=3.1415926535897932385D0 PIHAL=1.5707963267948966192D0 PI3=PI/3.D0 PI23=2.D0*PI3 X=X0 C=0.5D0*SQRT3 S=0.5D0 IERRO=0 IEXPF=0 IF (Y0.LT.0.D0) THEN Y=-Y0 ELSE Y=Y0 ENDIF R=SQRT(X*X+Y*Y) R32=R*SQRT(R) THET=PHASE(X,Y) COVER=2.D0/3.D0*R32*ABS(COS(1.5D0*THET)) CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER) OVER=D1MACH(2)*1.D-3 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER) UNDER=D1MACH(1)*1.D+3 DL1=LOG(OVER) DL2=-LOG(UNDER) IF (DL1.GT.DL2) OVER=1/UNDER IF (IFAC.EQ.1) THEN IF (COVER.GE.LOG(OVER)) THEN CCC OVERFLOW/UNDERFLOW PROBLEMS. CCC CALCULATION ABORTED IERRO=1 GBIR=0 GBII=0 ENDIF ELSE IF (COVER.GE.(LOG(OVER)*0.2)) IEXPF=1 ENDIF IF (IERRO.EQ.0) THEN IF (IFAC.EQ.1) THEN IF (Y.EQ.0.D0) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC TAKE TWICE THE REAL PART OF EXP(-PI I/6) AI_(1)(Z) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C1=-0.5D0 S1=-0.5D0*SQRT3 U=X*C1-Y*S1 V=X*S1+Y*C1 CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR) IF (IFUN.EQ.1) THEN BR=SQRT3*AR+AI BI=0.D0 ELSE U=AR*C1-AI*S1 V=AR*S1+AI*C1 APR=U API=V BPR=SQRT3*APR+API BPI=0.D0 ENDIF ELSE IF ((X.LT.0.D0).AND.(Y.LT.-X*SQRT3)) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC 2 PI/3 < PHASE(Z) < PI CCC BI(Z)=EXP(I PI/6) AI_(-1)(Z) + EXP(-I PI/6) AI_(1)(Z) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C1=-0.5D0 S1=0.5D0*SQRT3 U=X*C1-Y*S1 V=X*S1+Y*C1 CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR) IF (IFUN.EQ.1) THEN BR=C*AR-S*AI BI=C*AI+S*AR ELSE U=AR*C1-AI*S1 V=AR*S1+AI*C1 APR=U API=V BPR=C*APR-S*API BPI=C*API+S*APR ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC WE NEED ALSO AI_(1)(Z) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C1=-0.5D0 S1=-0.5D0*SQRT3 U=X*C1-Y*S1 V=X*S1+Y*C1 CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR) IF (IFUN.EQ.1) THEN S=-S BR=BR+C*AR-S*AI BI=BI+C*AI+S*AR ELSE U=AR*C1-AI*S1 V=AR*S1+AI*C1 APR=U API=V S=-S BPR=BPR+C*APR-S*API BPI=BPI+C*API+S*APR ENDIF ELSE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC BI(Z) = I AI(Z) + 2 EXP(-I PI/6) AI_(1)(Z) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C1=-0.5D0 S1=-0.5D0*SQRT3 U=X*C1-Y*S1 V=X*S1+Y*C1 CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR) IF (IFUN.EQ.1) THEN BR=SQRT3*AR+AI BI=-AR+SQRT3*AI ELSE U=AR*C1-AI*S1 V=AR*S1+AI*C1 APR=U API=V BPR=SQRT3*APR+API BPI=-APR+SQRT3*API ENDIF CALL AIZ(IFUN,IFAC,X,Y,AR,AI,IERR) IF (IFUN.EQ.1) THEN BR=BR-AI BI=BI+AR ELSE BPR=BPR-AI BPI=BPI+AR ENDIF ENDIF ENDIF ELSE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC SCALED BI AIRY FUNCTIONS C CCC WE USE THE FOLLOWING NORMALIZATION: C CCC LET ARGZ=ARG(Z), THEN: C CCC A) IF 0 <= ARGZ <= PI/3 C CCC BI=EXP(-2/3Z^3/2)BI C CCC B) IF PI/3 <= ARGZ <= PI C CCC BI=EXP(2/3Z^3/2)BI C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC THET=PHASE(X,Y) IF (THET.LE.PI3) THEN C1=-0.5D0 S1=-0.5D0*SQRT3 U=X*C1-Y*S1 V=X*S1+Y*C1 CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR) IF (IFUN.EQ.1) THEN BR=SQRT3*AR+AI BI=-AR+SQRT3*AI ELSE U=AR*C1-AI*S1 V=AR*S1+AI*C1 APR=U API=V BPR=SQRT3*APR+API BPI=-APR+SQRT3*API ENDIF IF (IEXPF.EQ.0) THEN R=SQRT(X*X+Y*Y) R32=R*SQRT(R) THET32=THET*1.5D0 A1=COS(THET32) B1=SIN(THET32) DF1=4.D0/3.D0*R32 EXPO=EXP(DF1*A1) EXPOI=1.D0/EXPO ETAR=EXPO*COS(DF1*B1) ETAI=EXPO*SIN(DF1*B1) ETAGR=EXPOI*COS(-DF1*B1) ETAGI=EXPOI*SIN(-DF1*B1) CALL AIZ(IFUN,IFAC,X,Y,AR,AI,IERR) IF (IFUN.EQ.1) THEN BR=BR-AR*ETAGI-ETAGR*AI BI=BI+AR*ETAGR-ETAGI*AI ELSE BPR=BPR-AR*ETAGI-ETAGR*AI BPI=BPI+AR*ETAGR-ETAGI*AI ENDIF ENDIF ENDIF IF ((THET.GT.PI3).AND.(THET.LE.PI23)) THEN IF (IEXPF.EQ.0) THEN R=SQRT(X*X+Y*Y) R32=R*SQRT(R) THET32=THET*1.5D0 A1=COS(THET32) B1=SIN(THET32) DF1=4.D0/3.D0*R32 EXPO=EXP(DF1*A1) EXPOI=1.D0/EXPO ETAR=EXPO*COS(DF1*B1) ETAI=EXPO*SIN(DF1*B1) ETAGR=EXPOI*COS(-DF1*B1) ETAGI=EXPOI*SIN(-DF1*B1) C1=-0.5D0 S1=-0.5D0*SQRT3 U=X*C1-Y*S1 V=X*S1+Y*C1 CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR) IF (IFUN.EQ.1) THEN BBR=SQRT3*AR+AI BBI=-AR+SQRT3*AI BR=BBR*ETAR-BBI*ETAI BI=BBI*ETAR+BBR*ETAI ELSE U=AR*C1-AI*S1 V=AR*S1+AI*C1 APR=U API=V BBPR=SQRT3*APR+API BBPI=-APR+SQRT3*API BPR=BBPR*ETAR-BBPI*ETAI BPI=BBPI*ETAR+BBPR*ETAI ENDIF ELSE IF (IFUN.EQ.1) THEN BR=0.D0 BI=0.D0 ELSE BPR=0.D0 BPI=0.D0 ENDIF ENDIF CALL AIZ(IFUN,IFAC,X,Y,AR,AI,IERR) IF (IFUN.EQ.1) THEN BR=BR-AI BI=BI+AR ELSE BPR=BPR-AI BPI=BPI+AR ENDIF ENDIF IF (THET.GT.PI23) THEN C1=-0.5D0 S1=0.5D0*SQRT3 U=X*C1-Y*S1 V=X*S1+Y*C1 CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR) IF (IFUN.EQ.1) THEN BR=C*AR-S*AI BI=C*AI+S*AR ELSE U=AR*C1-AI*S1 V=AR*S1+AI*C1 APR=U API=V BPR=C*APR-S*API BPI=C*API+S*APR ENDIF IF (IEXPF.EQ.0) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC WE NEED ALSO AI_(1)(Z) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC R=SQRT(X*X+Y*Y) R32=R*SQRT(R) THET32=THET*1.5D0 A1=COS(THET32) B1=SIN(THET32) DF1=4.D0/3.D0*R32 EXPO=EXP(DF1*A1) EXPOI=1.D0/EXPO ETAR=EXPO*COS(DF1*B1) ETAI=EXPO*SIN(DF1*B1) ETAGR=EXPOI*COS(-DF1*B1) ETAGI=EXPOI*SIN(-DF1*B1) C1=-0.5D0 S1=-0.5D0*SQRT3 U=X*C1-Y*S1 V=X*S1+Y*C1 CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR) IF (IFUN.EQ.1) THEN S=-S BBR=C*AR-S*AI BBI=C*AI+S*AR BR=BR+ETAR*BBR-ETAI*BBI BI=BI+BBI*ETAR+ETAI*BBR ELSE U=AR*C1-AI*S1 V=AR*S1+AI*C1 APR=U API=V S=-S BBPR=C*APR-S*API BBPI=C*API+S*APR BPR=BPR+ETAR*BBPR-ETAI*BBPI BPI=BPI+BBPI*ETAR+ETAI*BBPR ENDIF ENDIF ENDIF ENDIF IF (Y0.LT.0) THEN BI=-BI BPI=-BPI ENDIF IF (IFUN.EQ.1) THEN GBIR=BR GBII=BI ELSE GBIR=BPR GBII=BPI ENDIF ENDIF RETURN END SHAR_EOF fi # end of overwriting check if test -f 'd1mach.f' then echo shar: will not over-write existing file "'d1mach.f'" else cat << "SHAR_EOF" > 'd1mach.f' DOUBLE PRECISION FUNCTION D1MACH(I) INTEGER I C C DOUBLE-PRECISION MACHINE CONSTANTS C D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. C D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. C D1MACH( 5) = LOG10(B) C INTEGER SMALL(2) INTEGER LARGE(2) INTEGER RIGHT(2) INTEGER DIVER(2) INTEGER LOG10(2) INTEGER SC, CRAY1(38), J COMMON /D9MACH/ CRAY1 SAVE SMALL, LARGE, RIGHT, DIVER, LOG10, SC DOUBLE PRECISION DMACH(5) EQUIVALENCE (DMACH(1),SMALL(1)) EQUIVALENCE (DMACH(2),LARGE(1)) EQUIVALENCE (DMACH(3),RIGHT(1)) EQUIVALENCE (DMACH(4),DIVER(1)) EQUIVALENCE (DMACH(5),LOG10(1)) C THIS VERSION ADAPTS AUTOMATICALLY TO MOST CURRENT MACHINES. C R1MACH CAN HANDLE AUTO-DOUBLE COMPILING, BUT THIS VERSION OF C D1MACH DOES NOT, BECAUSE WE DO NOT HAVE QUAD CONSTANTS FOR C MANY MACHINES YET. C TO COMPILE ON OLDER MACHINES, ADD A C IN COLUMN 1 C ON THE NEXT LINE DATA SC/0/ C AND REMOVE THE C FROM COLUMN 1 IN ONE OF THE SECTIONS BELOW. C CONSTANTS FOR EVEN OLDER MACHINES CAN BE OBTAINED BY C mail netlib@research.bell-labs.com C send old1mach from blas C PLEASE SEND CORRECTIONS TO dmg OR ehg@bell-labs.com. C C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. C DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /, SC/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 32-BIT INTEGERS. C DATA SMALL(1),SMALL(2) / 8388608, 0 / C DATA LARGE(1),LARGE(2) / 2147483647, -1 / C DATA RIGHT(1),RIGHT(2) / 612368384, 0 / C DATA DIVER(1),DIVER(2) / 620756992, 0 / C DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /, SC/987/ C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /, SC/987/ C C ON FIRST CALL, IF NO DATA UNCOMMENTED, TEST MACHINE TYPES. IF (SC .NE. 987) THEN DMACH(1) = 1.D13 IF ( SMALL(1) .EQ. 1117925532 * .AND. SMALL(2) .EQ. -448790528) THEN * *** IEEE BIG ENDIAN *** SMALL(1) = 1048576 SMALL(2) = 0 LARGE(1) = 2146435071 LARGE(2) = -1 RIGHT(1) = 1017118720 RIGHT(2) = 0 DIVER(1) = 1018167296 DIVER(2) = 0 LOG10(1) = 1070810131 LOG10(2) = 1352628735 ELSE IF ( SMALL(2) .EQ. 1117925532 * .AND. SMALL(1) .EQ. -448790528) THEN * *** IEEE LITTLE ENDIAN *** SMALL(2) = 1048576 SMALL(1) = 0 LARGE(2) = 2146435071 LARGE(1) = -1 RIGHT(2) = 1017118720 RIGHT(1) = 0 DIVER(2) = 1018167296 DIVER(1) = 0 LOG10(2) = 1070810131 LOG10(1) = 1352628735 ELSE IF ( SMALL(1) .EQ. -2065213935 * .AND. SMALL(2) .EQ. 10752) THEN * *** VAX WITH D_FLOATING *** SMALL(1) = 128 SMALL(2) = 0 LARGE(1) = -32769 LARGE(2) = -1 RIGHT(1) = 9344 RIGHT(2) = 0 DIVER(1) = 9472 DIVER(2) = 0 LOG10(1) = 546979738 LOG10(2) = -805796613 ELSE IF ( SMALL(1) .EQ. 1267827943 * .AND. SMALL(2) .EQ. 704643072) THEN * *** IBM MAINFRAME *** SMALL(1) = 1048576 SMALL(2) = 0 LARGE(1) = 2147483647 LARGE(2) = -1 RIGHT(1) = 856686592 RIGHT(2) = 0 DIVER(1) = 873463808 DIVER(2) = 0 LOG10(1) = 1091781651 LOG10(2) = 1352628735 ELSE IF ( SMALL(1) .EQ. 1120022684 * .AND. SMALL(2) .EQ. -448790528) THEN * *** CONVEX C-1 *** SMALL(1) = 1048576 SMALL(2) = 0 LARGE(1) = 2147483647 LARGE(2) = -1 RIGHT(1) = 1019215872 RIGHT(2) = 0 DIVER(1) = 1020264448 DIVER(2) = 0 LOG10(1) = 1072907283 LOG10(2) = 1352628735 ELSE IF ( SMALL(1) .EQ. 815547074 * .AND. SMALL(2) .EQ. 58688) THEN * *** VAX G-FLOATING *** SMALL(1) = 16 SMALL(2) = 0 LARGE(1) = -32769 LARGE(2) = -1 RIGHT(1) = 15552 RIGHT(2) = 0 DIVER(1) = 15568 DIVER(2) = 0 LOG10(1) = 1142112243 LOG10(2) = 2046775455 ELSE DMACH(2) = 1.D27 + 1 DMACH(3) = 1.D27 LARGE(2) = LARGE(2) - RIGHT(2) IF (LARGE(2) .EQ. 64 .AND. SMALL(2) .EQ. 0) THEN CRAY1(1) = 67291416 DO 10 J = 1, 20 CRAY1(J+1) = CRAY1(J) + CRAY1(J) 10 CONTINUE CRAY1(22) = CRAY1(21) + 321322 DO 20 J = 22, 37 CRAY1(J+1) = CRAY1(J) + CRAY1(J) 20 CONTINUE IF (CRAY1(38) .EQ. SMALL(1)) THEN * *** CRAY *** CALL I1MCRY(SMALL(1), J, 8285, 8388608, 0) SMALL(2) = 0 CALL I1MCRY(LARGE(1), J, 24574, 16777215, 16777215) CALL I1MCRY(LARGE(2), J, 0, 16777215, 16777214) CALL I1MCRY(RIGHT(1), J, 16291, 8388608, 0) RIGHT(2) = 0 CALL I1MCRY(DIVER(1), J, 16292, 8388608, 0) DIVER(2) = 0 CALL I1MCRY(LOG10(1), J, 16383, 10100890, 8715215) CALL I1MCRY(LOG10(2), J, 0, 16226447, 9001388) ELSE WRITE(*,9000) STOP 779 END IF ELSE WRITE(*,9000) STOP 779 END IF END IF SC = 987 END IF * SANITY CHECK IF (DMACH(4) .GE. 1.0D0) STOP 778 IF (I .LT. 1 .OR. I .GT. 5) THEN WRITE(*,*) 'D1MACH(I): I =',I,' is out of bounds.' STOP END IF D1MACH = DMACH(I) RETURN 9000 FORMAT(/' Adjust D1MACH by uncommenting data statements'/ *' appropriate for your machine.') * /* Standard C source for D1MACH -- remove the * in column 1 */ *#include *#include *#include *double d1mach_(long *i) *{ * switch(*i){ * case 1: return DBL_MIN; * case 2: return DBL_MAX; * case 3: return DBL_EPSILON/FLT_RADIX; * case 4: return DBL_EPSILON; * case 5: return log10((double)FLT_RADIX); * } * fprintf(stderr, "invalid argument: d1mach(%ld)\n", *i); * exit(1); return 0; /* some compilers demand return values */ *} END SUBROUTINE I1MCRY(A, A1, B, C, D) **** SPECIAL COMPUTATION FOR OLD CRAY MACHINES **** INTEGER A, A1, B, C, D A1 = 16777216*B + C A = 16777216*A1 + D END SHAR_EOF fi # end of overwriting check if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << "SHAR_EOF" > 'src.f' SUBROUTINE DKIA(IFAC,X,A,DKI,DKID,IERRO) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCC THE PURPOUSE OF THE ROUTINE IS THE CALCULATION OF THE CCCC MODIFIED BESSEL FUNCTIONS K_IA(X) AND K'_IA(X), CCCC WHERE I IS THE IMAGINARY UNIT AND X IS THE ARGUMENT OF THE CCCC FUNCTIONS. WE WILL REFER TO A AS THE ORDER OF THE FUNCTIONS. CCCC CCCC THE ROUTINE HAS THE OPTION OF COMPUTING SCALED FUNCTIONS. CCCC THIS SCALING CAN BE USED TO ENLARGE THE RANGE OF CCCC COMPUTATION. CCCC THE SCALED FUNCTIONS ARE DEFINED AS FOLLOWS CCCC (S STANDS FOR SCALED AND CCCC L=SQRT{X**2-A**2} + A*ARCSIN(A/X)): CCCC CCCC EXP(L)*K_IA(X), IF X >=ABS(A) CCCC SK_IA(X) = CCCC EXP(ABS(A)*PI/2)*K_IA(X), IF X < ABS(A) CCCC CCCC CCCC EXP(L)*K'_IA(X), IF X >=ABS(A) CCCC SK'_IA(X) = CCCC EXP(ABS(A)*PI/2)*K'_IA(X), IF X < ABS(A) CCCC CCCC THE RANGE OF THE PARAMETERS (X,A) FOR THE COMPUTATION OF CCCC SCALED FUNCTIONS IS: CCCC 0 < X <= 1500, -1500 <= A <= 1500. CCCC FOR FUNCTIONS WITHOUT SCALING, THE RANGE IS USUALLY LARGER CCCC THAN CCCC 0 < X <= 500, -400 <= A <= 400, CCCC DEPENDING ON THE MACHINE OVERFLOW/UNDERFLOW PARAMETERS, WHICH CCCC ARE SET UP BY THE ROUTINE D1MACH. CCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCC METHODS OF COMPUTATION: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCC 1) SERIES, IF ABS(A) > 0.044*ABS(X-3.1)**1.9+(X-3.1) CCCC 2) CONTINUED FRACTION, CCCC IF X>3 AND ABS(A) < 380*(ABS((X-3)/2300))**0.572 CCCC 3) AIRY-TYPE ASYMPTOTIC EXPANSION, CCCC IF ABS(A) > 0.4*X+7.5 AND ABS(A) < 3.7*X-10 CCCC 4) QUADRATURES, CCCC IN THE REST OF CASES CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCC DESCRIPTION OF INPUT/OUTPUT VARIABLES: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCC INPUTS: CCCC X : ARGUMENT OF THE FUNCTIONS CCC A : ORDER OF THE FUNCTIONS CCCC IFAC : INTEGER FLAG FOR THE SCALING CCCC * IF IFAC=1, THE CODE COMPUTES KIA(X), KIA'(X) CCCC * OTHERWISE, THE CODE COMPUTES SCALED KIA(X), KIA'(X) CCCC OUTPUTS CCCC DKI : KIA(X) FUNCTION CCCC DKID : DERIVATIVE OF THE KIA(X) FUNCTION CCCC IERRO: ERROR FLAG CCCC * IF IERRO=0, COMPUTATION SUCCESSFUL. CCCC * IF IERRO=1, COMPUTATION OUT OF RANGE. CCCC * IF IERRO=2, VARIABLES X AND/OR A, OUT OF RANGE. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCC ACCURACY: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCC THE RELATIVE ACCURACY IS: CCCC BETTER THAN 10**(-13) FOR (X,A) IN (0,200]X[-200,200]; CCCC BETTER THAN 5.10**(-13) FOR (X,A) IN (0,500]X[-500,500]; CCCC CLOSE TO 10**(-12) FOR (X,A) IN (0,1500]X[-1500,1500]. CCCC NEAR THE ZEROS OF THE FUNCTIONS (THERE ARE INFINITELY CCCC MANY OF THEM IN THE OSCILLATORY REGION) RELATIVE PRECISION CCCC LOOSES MEANING AND ONLY ABSOLUTE PRECISION MAKES SENSE. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C AUTHORS: C AMPARO GIL (U. CANTABRIA, SANTANDER, SPAIN). C E-MAIL: AMPARO.GIL@UNICAN.ES C JAVIER SEGURA (U. CANTABRIA, SANTANDER, SPAIN). C E-MAIL: SEGURAJJ@UNICAN.ES C NICO M. TEMME (CWI, AMSTERDAM, THE NETHERLANDS). C E-MAIL: NICO.TEMME@CWI.NL CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C REFERENCES: C THIS IS THE COMPANION SOFTWARE OF THE ARTICLES C 1)'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL DIFFERENTIAL C EQUATION FOR IMAGINARY ORDERS AND POSITIVE ARGUMENTS', C A. GIL, J. SEGURA, N.M. TEMME C ACM TRANS. MATH. SOFT. (2004) C 2)'MODIFIED BESSEL FUNCTIONS OF IMAGINARY ORDER AND C POSITIVE ARGUMENT', C A. GIL, J. SEGURA, N.M. TEMME C ACM TRANS. MATH. SOFT. (2004) C ADDITIONAL REFERENCES: C - 'COMPUTATION OF THE MODIFIED BESSEL FUNCTION OF THE C THIRD KIND FOR IMAGINARY ORDERS' C A. GIL, J. SEGURA, N.M. TEMME C J. COMPUT. PHYS. 175 (2002) 398-411 C - 'COMPUTATION OF THE MODIFIED BESSEL FUNCTIONS OF THE C THIRD KIND OF IMAGINARY ORDERS: C UNIFORM AIRY-TYPE ASYMPTOTIC EXPANSION' C A. GIL, J. SEGURA, N.M. TEMME, PROCEEDINGS OPSFA 2001 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC LOCAL VARIABLES: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C D: X-3.1 C DF1: 0.044*ABS(D)**1.9+(X-3.1) C DF2: 380*(ABS((X-3)/2300))**0.572 C DF3: 0.4*X+7.5 C DF4: 3.7*X-10 C PNU: ORDER OF THE FUNCTION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION A,D,DF1,DF2,DF3,DF4,DKI, + DKID,PNU,X INTEGER IERRO,IFAC IERRO=0 PNU=A IF ((X.GT.1500.D0).OR.(X.LE.0.D0)) THEN IERRO=2 DKI=0.D0 DKID=0.D0 ENDIF IF (ABS(PNU).GT.1500.D0) THEN IERRO=2 DKI=0.D0 DKID=0.D0 ENDIF IF (IERRO.EQ.0) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC THESE FUNCTIONS ARE EVEN FUNCTIONS IN THE C CCC PARAMETER A (=PNU) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF (PNU.LT.0.D0) PNU=-PNU D=X-3.1D0 DF1=0.044D0*ABS(D)**1.9D0+D DF2=380.D0*(ABS((X-3.D0)/2300.D0))**0.572D0 DF3=0.4D0*X+7.5D0 DF4=3.7D0*X-10.D0 IF (PNU.GT.DF1) THEN CALL SERKIA(IFAC,X,PNU,DKI,DKID,IERRO) ELSEIF ((X.GT.3.D0).AND.(PNU.LT.DF2)) THEN CALL FRAKIA(IFAC,X,PNU,DKI,DKID,IERRO) ELSEIF ((PNU.GT.DF3).AND.(PNU.LT.DF4)) THEN CALL AIEXKI(IFAC,X,PNU,DKI,DKID,IERRO) ELSE CALL DKINT(IFAC,X,PNU,DKI,DKID,IERRO) ENDIF ENDIF RETURN END SUBROUTINE DKINT(IFAC,XX,PNUA,DKINF,DKIND,IERRO) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC CALCULATION OF THE FUNCTIONS K,K' USING NON-OSCILLATING C CCC INTEGRAL REPRESENTATIONS C CCC C CCC INPUTS: C CCC XX: ARGUMENT OF THE FUNCTIONS C CCC PNUA: ORDER OF THE FUNCTIONS C CCC IFAC: C CCC =1, NON SCALED FUNCTIONS C CCC OTYHERWISE, SCALED FUNCTIONS C CCC OUTPUTS: C CCC DKINF: K FUNCTION C CCC DKIND: DERIVATIVE OF THE K FUNCTION C CCC IERRO: ERROR FLAG C CCC * IF IERRO=0, COMPUTATION SUCCESSFUL. C CCC * IF IERRO=1, COMPUTATION OUT OF RANGE. C CCC * IF IERRO=2, ARGUMENT AND/OR ORDER, OUT OF RANGE. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC METHOD OF COMPUTATION: CCC * IF XX>=PNUA, THE NON-OSCILLATING INTEGRAL REPRESENTATIONS CCC FOR THE MONOTONIC REGION ARE USED CCC * IF XX=1, 2/3*(-PSI)**3/2=SQRT(Z**2-1)-ARCCOS(1/Z) C PSI12 : SQRT(PSI) C PSI2 : PSI*PSI C PSI3 : PSI**3 C SAS : ACCUMULATES THE CONTRIBUTION OF A_S(PSI) C SBS : ACCUMULATES THE CONTRIBUTION OF B_S(PSI) C SCS : ACCUMULATES THE CONTRIBUTION OF C_S(PSI) C SDS : ACCUMULATES THE CONTRIBUTION OF D_S(PSI) C SIG : (-)**K C SINTH : SIN(THET) C THET : ASIN(A/X) C UNDER : UNDERFLOW NUMBER C Y : Z-1 C Z : X/A C Z2 : Z**2 C Z21M : 1-Z**2 C ZDS : SQRT(1-Z**2) C ZMASF : EXPANSION IN TERMS OF (Z-1) FOR THE C CALCULATION OF PSI(Z) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION A,A0EX,A13,A2,A23,A2K,AC(0:5,0:20), + AII,AIR,APIHAL,ARG,AS,ASP,B0EX,B0PEX,BC(0:5,0:20), + BS,BSO,BSP,BSPO,C0EX,CHI,CHIEX,CHIN(0:26),COSTH, + COZMAS(0:12),D0EX,D1MACH,DAII,DAIR,DF,DKAI,DKAID DOUBLE PRECISION DZZ,ETA,ETAJ,ETAK,ETAL,EXPAM, + EXPAPI,F2,F4,FAC,FACD,FDOMIN,PHI(0:35),PHIEX,PHIEX2, + PHIS,PI,PIHALF,PSI,PSI12,PSI2,PSI3,SAS,SBS,SCS,SDS, + SIG,SINTH,THET,UNDER,X,Y,Z,Z2,Z21M,ZDS,ZMASF INTEGER IERRO,IERROK,IFAC,IFACA,IFUN,INDA(0:5),INDB(0:5), + J,K,L SAVE PHI,CHIN,AC,BC CCCC VALUES OF THE COEFFICIENTS DATA PHI/1.D0,0.2D0,.25714285714285714286D-1, + -.56507936507936507937D-2,-.39367965367965367965D-2, + -.5209362066504924D-3,.3708541264731741D-3, + .2123827840293627D-3,.2150629788753145D-4, + -.2636904062981799D-4,-.1405469826493129D-4, + -.1149328899029441D-5,.1972641193938624D-5, + .1014324305593329D-5,.7034331100192200D-7, + -.1525044777392676D-6,-.7677866256900572D-7, + -.4669842638693018D-8,.1206673645965329D-7, + .5990877668092344D-8,.3269102150077715D-9, + -.97138350244428335095D-9,-.47745934295232233834D-9, + -.23750318642839155779D-10,.79244598109106655567D-10, + .38653584230817865528D-10,.17734105846426807873D-11, + -.65332110030864751956D-11,-.31673512094772575686D-11, + -.13524195177121030660D-12,.54324103217903338951D-12, + .26204918647967626464D-12,.10488134973664584898D-13, + -.45490420121539001435D-13,-.21851238232690420878D-13, + -.82456080260379042800D-15/ DATA CHIN/.2D0,.1142857142857142857142857D-1, + -.2438095238095238095238095D-1,-.1003471449185734900020614D-1, + .8811404468547325690182833D-3,.2318339655809043564145605D-2, + .7794494413564441575646057D-3,-.1373952504077228744949558D-3, + -.2162301322540308393022684D-3,-.6585004634375583559042795D-4, + .1502851367097217578058824D-4,.1991904617871647675455262D-4, + .5731972706719818525615427D-5,-.1496427612747891044606885D-5, + -.1821231830428939670133992D-5,-.5054254875930882534538217D-6, + .1433283859497625931203415D-6,.1657681319078678321113361D-6, + .4485592642702540575627044D-7,-.1346138242826094098161839D-7, + -.1504601862773535117708677D-7,-.3995660198654805921651406D-8, + .1250124931952495738882300D-8,.1363187174221864073749532D-8, + .3567608459777506132029204D-9,-.1152749290139859167119863D-9, + -.1233547289799408517691696D-9/ DATA COZMAS/.9428090415820634D0,-.4242640687119285D0, + .2904188565587606D0,-.2234261009999160D0, + .1821754153944745D0,-.1539625154198624D0, + .1333737583085222D0,-.1176617834148007D0, + .1052687194772381D0,-.9524025714638822D-1, + .8695738197073783D-1,-.8000034897653656D-1, + .7407421797273086D-1/ DATA AC(0,0)/1.D0/ DATA AC(1,0)/-.44444444444444444445D-2/ DATA AC(1,1)/-.18441558441558441558D-2/ DATA AC(1,2)/.11213675213675213675D-2/ DATA AC(1,3)/.13457752124418791086D-2/ DATA AC(1,4)/.0003880626562979504D0/ DATA AC(1,5)/-.0001830686723781799D0/ DATA AC(1,6)/-.0001995460887806733D0/ DATA AC(1,7)/-.00005256191234041587D0/ DATA AC(1,8)/.00002460619652459158D0/ DATA AC(1,9)/.00002519246780924541D0/ DATA AC(1,10)/.6333157376533242D-5/ DATA AC(1,11)/-.2957485733830202D-5/ DATA AC(1,12)/-.2925255920564838D-5/ DATA AC(1,13)/-.7159702610502009D-6/ DATA AC(1,14)/.3331510720390949D-6/ DATA AC(1,15)/.3227670475692310D-6/ DATA AC(1,16)/.7767729381664199D-7/ DATA AC(1,17)/-.3600954237921120D-7/ DATA AC(1,18)/-.3441724449034226D-7/ DATA AC(1,19)/-.8188194356398772D-8/ DATA AC(1,20)/.3783148485152038D-8/ DATA AC(2,0)/.69373554135458897365D-3/ DATA AC(2,1)/.46448349036584330703D-3/ DATA AC(2,2)/-.42838130171535112460D-3/ DATA AC(2,3)/-.0007026702868771135D0/ DATA AC(2,4)/-.0002632580046778811D0/ DATA AC(2,5)/.0001663853666288703D0/ DATA AC(2,6)/.0002212087687818584D0/ DATA AC(2,7)/.00007020345615329662D0/ DATA AC(2,8)/-.00004000421782540614D0/ DATA AC(2,9)/-.00004786324966453962D0/ DATA AC(2,10)/-.00001394600741473631D0/ DATA AC(2,11)/.7536186591273727D-5/ DATA AC(2,12)/.8478502161067667D-5/ DATA AC(2,13)/.2345355228453912D-5/ DATA AC(2,14)/-.1225943294710883D-5/ DATA AC(2,15)/-.1325082343401027D-5/ DATA AC(2,16)/-.3539954776569997D-6/ DATA AC(2,17)/.1808291719376674D-6/ DATA AC(2,18)/.1900383515233655D-6/ DATA AC(3,0)/-.35421197145774384076D-3/ DATA AC(3,1)/-.31232252789031883276D-3/ DATA AC(3,2)/.3716442237502298D-3/ DATA AC(3,3)/.0007539269155977733D0/ DATA AC(3,4)/.0003408300059444739D0/ DATA AC(3,5)/-.0002634968172069594D0/ DATA AC(3,6)/-.0004089275726648432D0/ DATA AC(3,7)/-.0001501108759563460D0/ DATA AC(3,8)/.00009964015205538056D0/ DATA AC(3,9)/.0001352492955751283D0/ DATA AC(3,10)/.00004443117087272903D0/ DATA AC(3,11)/-.00002713205071914117D0/ DATA AC(3,12)/-.00003396796969771860D0/ DATA AC(3,13)/-.00001040708865273043D0/ DATA AC(3,14)/.6024639065414414D-5/ DATA AC(3,15)/.7143919607846883D-5/ DATA AC(4,0)/.378194199201773D-3/ DATA AC(4,1)/.000404943905523634D0/ DATA AC(4,2)/-.000579130526946453D0/ DATA AC(4,3)/-.00138017901171011D0/ DATA AC(4,4)/-.000722520056780091D0/ DATA AC(4,5)/.000651265924036825D0/ DATA AC(4,6)/.00114674563328389D0/ DATA AC(4,7)/.000474423189340405D0/ DATA AC(4,8)/-.000356495172735468D0/ DATA AC(4,9)/-.000538157791035111D0/ DATA AC(4,10)/-.000195687390661225D0/ DATA AC(4,11)/.000132563525210293D0/ DATA AC(4,12)/.000181949256267291D0/ DATA AC(5,0)/-.69114139728829416760D-3/ DATA AC(5,1)/-.00085995326611774383285D0/ DATA AC(5,2)/.0014202335568143511489D0/ DATA AC(5,3)/.0038535426995603052443D0/ DATA AC(5,4)/.0022752811642901374595D0/ DATA AC(5,5)/-.0023219572034556988366D0/ DATA AC(5,6)/-.0045478643394434635622D0/ DATA AC(5,7)/-.0020824431758272449919D0/ DATA AC(5,8)/.0017370443573195808719D0/ DATA BC(0,0)/.14285714285714285714D-1/ DATA BC(0,1)/.88888888888888888889D-2/ DATA BC(0,2)/.20482374768089053803D-2/ DATA BC(0,3)/-.57826617826617826618D-3/ DATA BC(0,4)/-.60412089799844901886D-3/ DATA BC(0,5)/-.0001472685745626922D0/ DATA BC(0,6)/.00005324102148009784D0/ DATA BC(0,7)/.00005206561006583416D0/ DATA BC(0,8)/.00001233115050894939D0/ DATA BC(0,9)/-.4905932728531366D-5/ DATA BC(0,10)/-.4632230987136350D-5/ DATA BC(0,11)/-.1077174523455235D-5/ DATA BC(0,12)/.4475963978932822D-6/ DATA BC(0,13)/.4152586188464624D-6/ DATA BC(0,14)/.9555819293589234D-7/ DATA BC(0,15)/-.4060599208403059D-7/ DATA BC(0,16)/-.3731367187988482D-7/ DATA BC(0,17)/-.8532670645553778D-8/ DATA BC(0,18)/.3673017245573624D-8/ DATA BC(0,19)/.3355960460784536D-8/ DATA BC(0,20)/.7643107095110475D-9/ DATA BC(1,0)/-.11848595848595848596D-2/ DATA BC(1,1)/-.13940630797773654917D-2/ DATA BC(1,2)/-.48141005586383737645D-3/ DATA BC(1,3)/.26841705366016142958D-3/ DATA BC(1,4)/.0003419706982709903D0/ DATA BC(1,5)/.0001034548234902078D0/ DATA BC(1,6)/-.00005418191982095504D0/ DATA BC(1,7)/-.00006202184830690167D0/ DATA BC(1,8)/-.00001724885886056087D0/ DATA BC(1,9)/.8744675992887053D-5/ DATA BC(1,10)/.9420684216180929D-5/ DATA BC(1,11)/.2494922112085850D-5/ DATA BC(1,12)/-.1238458608836357D-5/ DATA BC(1,13)/-.1285461713809769D-5/ DATA BC(1,14)/-.3299710862537507D-6/ DATA BC(1,15)/.1613441105788315D-6/ DATA BC(1,16)/.1633623194402374D-6/ DATA BC(1,17)/.4104252949605779D-7/ DATA BC(1,18)/-.1984317042326989D-7/ DATA BC(1,19)/-.1973948142769706D-7/ DATA BC(1,20)/-.4882194808588752D-8/ DATA BC(2,0)/.43829180944898810994D-3/ DATA BC(2,1)/.71104865116708668943D-3/ DATA BC(2,2)/.31858383945387580576D-3/ DATA BC(2,3)/-.0002404809426804458D0/ DATA BC(2,4)/-.0003722966038621536D0/ DATA BC(2,5)/-.0001352752059595618D0/ DATA BC(2,6)/.00008691694372704142D0/ DATA BC(2,7)/.0001158750753591377D0/ DATA BC(2,8)/.00003724965927846447D0/ DATA BC(2,9)/-.00002198334949606935D0/ DATA BC(2,10)/-.00002686449633870452D0/ DATA BC(2,11)/-.8023061612032524D-5/ DATA BC(2,12)/.4494756592180126D-5/ DATA BC(2,13)/.5193504763856015D-5/ DATA BC(2,14)/.1477156191529617D-5/ DATA BC(2,15)/-.7988793826096784D-6/ DATA BC(3,0)/-.37670439477105454219D-3/ DATA BC(3,1)/-.75856271658798642365D-3/ DATA BC(3,2)/-.0004103253968775048D0/ DATA BC(3,3)/.0003791263310429010D0/ DATA BC(3,4)/.0006850981673903450D0/ DATA BC(3,5)/.0002878310571932216D0/ DATA BC(3,6)/-.0002157010636115705D0/ DATA BC(3,7)/-.0003260863991373500D0/ DATA BC(3,8)/-.0001181317008748678D0/ DATA BC(3,9)/.00007887526841582582D0/ DATA BC(3,10)/.0001072081833420685D0/ DATA BC(3,11)/.00003544595251288735D0/ DATA BC(3,12)/-.00002201447920733824D0/ DATA BC(3,13)/-.00002789336359620813D0/ DATA BC(4,0)/.00058453330122076187255D0/ DATA BC(4,1)/.0013854690422372401251D0/ DATA BC(4,2)/.00086830374184946900245D0/ DATA BC(4,3)/-.00093502904801345951693D0/ DATA BC(4,4)/-.0019175486005525492059D0/ DATA BC(4,5)/-.00090795047113308137941D0/ DATA BC(4,6)/.00077050429806392235104D0/ DATA BC(4,7)/.0012953100128255983488D0/ DATA BC(4,8)/.00051933869471899550762D0/ DATA BC(4,9)/-.00038482631948759834653D0/ DATA BC(4,10)/-.00057335393099012476502D0/ DATA BC(5,0)/-.0014301070053470410656D0/ DATA BC(5,1)/-.0038637811942002539408D0/ DATA BC(5,2)/-.0027322816261168328245D0/ DATA BC(5,3)/.0033294980346743452748D0/ DATA BC(5,4)/.0075972237660887795911D0/ DATA BC(5,5)/.0039816655673062060620D0/ DATA BC(5,6)/-.0037510180460986006595D0/ DATA INDA/0,20,18,15,12,8/ DATA INDB/20,20,15,13,10,6/ CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER) UNDER=D1MACH(1)*1.D+8 CCCC CONSTANTS CCCCCCCCCCCCCC PI=ACOS(-1.D0) PIHALF=PI*0.5D0 IERROK=0 IF (IFAC.EQ.1) THEN IF (X.GE.A) THEN SINTH=A/X THET=ASIN(SINTH) COSTH=COS(THET) FDOMIN=X*(COSTH+THET*SINTH) IF (-FDOMIN.LE.LOG(UNDER)) IERROK=1 ELSE IF ((-PI*A*0.5D0).LE.LOG(UNDER)) IERROK=1 ENDIF ENDIF IF (IERROK.EQ.0) THEN DF=2.D0**(1.D0/3.D0) CCCC VARIABLES CCCCCCCCCCCCCCC Z=X/A A2=A*A A13=A**(1.D0/3.D0) A23=A13*A13 IF (IFAC.EQ.1) THEN APIHAL=A*PIHALF EXPAPI=EXP(-APIHAL) FAC=PI*EXPAPI/A13 FACD=2.D0*PI*EXPAPI/Z/A23 ELSE IF (X.LT.A) THEN FAC=PI/A13 FACD=2.D0*PI/Z/A23 ELSE SINTH=A/X THET=ASIN(SINTH) COSTH=COS(THET) FDOMIN=X*(COSTH+THET*SINTH) EXPAM=EXP(-A*PIHALF+FDOMIN) FAC=PI*EXPAM/A13 FACD=2.D0*PI*EXPAM/Z/A23 ENDIF ENDIF F4=A13**4 IF (Z.LE.0.9D0) THEN ZDS=SQRT((1.D0-Z)*(1.D0+Z)) PSI=(1.5D0*(LOG((1.D0+ZDS)/Z)-ZDS))**(2.D0/3.D0) ELSEIF (Z.GT.1.1D0) THEN ZDS=SQRT((Z-1.D0)*(1.D0+Z)) PSI=-(1.5D0*(ZDS-ACOS(1.D0/Z)))**(2.D0/3.D0) ELSE Y=Z-1.D0 ZMASF=COZMAS(12) DO 8 K=0,11 J=11-K ZMASF=COZMAS(J)+Y*ZMASF 8 CONTINUE ZMASF=ABS(Y)**(1.5D0)*ZMASF IF (Z.LT.1.D0) THEN PSI=(1.5D0*ZMASF)**(2.D0/3.D0) ELSE PSI=-(1.5D0*ZMASF)**(2.D0/3.D0) ENDIF ENDIF ETA=PSI/DF ARG=-A23*PSI PSI2=PSI*PSI PSI3=PSI2*PSI IF ((Z.GT.0.8D0).AND.(Z.LT.1.2D0)) THEN PHIS=0.D0 CHI=0.D0 SAS=0.D0 SBS=0.D0 SDS=0.D0 SCS=1.D0 BS=0.D0 BSPO=0.D0 DO 41 L=0,20 IF (L.EQ.0) THEN ETAL=1.D0 ELSE ETAL=ETAL*ETA ENDIF CHI=CHI+CHIN(L)*ETAL 41 CONTINUE CHI=CHI/DF DO 10 K=0,20 BSO=BS BSPO=BSP AS=0.D0 BS=0.D0 ASP=0.D0 BSP=0.D0 IF (K.EQ.0) THEN ETAK=1.D0 A2K=1.D0 SIG=1.D0 ELSE ETAK=ETAK*ETA A2K=A2K*A2 SIG=-1.D0*SIG ENDIF PHIS=PHIS+PHI(K)*ETAK F2=SIG/A2K IF (K.LE.5) THEN DO 20 J=0,20 IF (J.EQ.0) THEN ETAJ=1.D0 ELSE ETAJ=ETAJ*ETA ENDIF IF (J.LE.INDA(K)) THEN AS=AS+AC(K,J)*ETAJ ENDIF IF (J.LE.INDB(K)) THEN BS=BS+BC(K,J)*ETAJ ENDIF IF ((J+1).LE.INDA(K)) THEN ASP=ASP+(J+1)*AC(K,J+1)*ETAJ ENDIF IF ((J+1).LE.INDB(K)) THEN BSP=BSP+(J+1)*BC(K,J+1)*ETAJ ENDIF 20 CONTINUE ASP=1.D0/DF*ASP BS=BS*DF SAS=SAS+AS*F2 SBS=SBS+BS*F2/DF SDS=SDS-(CHI*AS+ASP+PSI*BS)*F2 ENDIF IF ((K.GT.0).AND.(K.LE.6)) THEN SCS=SCS+(AS+CHI*BSO+BSPO)*F2 ENDIF 10 CONTINUE PHIS=DF*PHIS SBS=DF*SBS ELSE CCCC EXACT VALUES CCCCCCCCCCCCCCCCCCCCC Z2=Z*Z Z21M=1.D0-Z2 PHIEX=(4.D0*PSI/Z21M)**0.25D0 PHIEX2=PHIEX*PHIEX A0EX=1.D0 B0EX=-5.D0/48.D0/(PSI*PSI)+PHIEX2/48.D0/PSI* + (5.D0/Z21M-3.D0) CHIEX=0.25D0/PSI*(1.D0-Z2*PHIEX2**3*0.25D0) CHI=CHIEX D0EX=-(7.D0/48.D0/PSI+PHIEX2/48.D0* + (9.D0-7.D0/Z21M)) IF (PSI.GT.0.D0) THEN PSI12=SQRT(PSI) DZZ=SQRT(Z21M) ELSE PSI12=SQRT(-PSI) DZZ=SQRT(ABS(Z21M)) ENDIF B0PEX=5.D0/24.D0/PSI3+PHIEX2/48.D0*((2.D0*CHI*PSI-1.D0)/ + PSI2*(5.D0/Z21M-3.D0)-10.D0*Z2*PSI12/ + Z21M**2/DZZ/PSI) C0EX=1.D0 SAS=A0EX SBS=B0EX SCS=C0EX SDS=D0EX BS=0.D0 BSP=0.D0 A2K=1.D0 SIG=1.D0 DO 30 K=1,6 BSO=BS BSPO=BSP AS=0.D0 BS=0.D0 ASP=0.D0 BSP=0.D0 A2K=A2K*A2 SIG=-1.D0*SIG F2=SIG/A2K IF (K.LE.5) THEN DO 40 J=0,20 IF (J.EQ.0) THEN ETAJ=1.D0 ELSE ETAJ=ETAJ*ETA ENDIF IF (J.LE.INDA(K)) THEN AS=AS+AC(K,J)*ETAJ ENDIF IF (J.LE.INDB(K)) THEN BS=BS+BC(K,J)*ETAJ ENDIF IF ((J+1).LE.INDA(K)) THEN ASP=ASP+(J+1)*AC(K,J+1)*ETAJ ENDIF IF ((J+1).LE.INDB(K)) THEN BSP=BSP+(J+1)*BC(K,J+1)*ETAJ ENDIF 40 CONTINUE BS=DF*BS ASP=ASP/DF SAS=SAS+AS*F2 SBS=SBS+BS*F2 SDS=SDS-(CHI*AS+ASP+PSI*BS)*F2 ENDIF IF (K.GE.2) THEN SCS=SCS+(AS+CHI*BSO+BSPO)*F2 ELSE SCS=SCS+(AS+CHI*B0EX+B0PEX)*F2 ENDIF 30 CONTINUE PHIS=PHIEX ENDIF CCCCC CALL THE AIRY ROUTINE CCCCCCCCC IFUN=1 IFACA=1 CALL AIZ(IFUN,IFACA,ARG,0.D0,AIR,AII,IERRO) IFUN=2 IFACA=1 CALL AIZ(IFUN,IFACA,ARG,0.D0,DAIR,DAII,IERRO) DKAI=FAC*PHIS*(AIR*SAS+DAIR*SBS/F4) DKAID=FACD/PHIS*(DAIR*SCS+AIR*SDS/A23) ELSE DKAI=0.D0 DKAID=0.D0 ENDIF RETURN END SUBROUTINE DLIA(IFAC,X,A,DLI,DLID,IERRO) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCC THE PURPOUSE OF THE ROUTINE IS THE CALCULATION OF THE CCCC MODIFIED BESSEL FUNCTIONS L_IA(X) AND L'_IA(X), CCCC WHERE I IS THE IMAGINARY UNIT AND X IS THE ARGUMENT OF THE CCCC FUNCTIONS. WE WILL REFER TO A AS THE ORDER OF THE FUNCTIONS. CCCC CCCC THE ROUTINE HAS THE OPTION OF COMPUTING SCALED FUNCTIONS. CCCC THIS SCALING CAN BE USED TO ENLARGE THE RANGE OF CCCC COMPUTATION. CCCC THE SCALED FUNCTIONS ARE DEFINED AS FOLLOWS: CCCC (S STANDS FOR SCALED AND CCCC L=SQRT{X**2-A**2} + A*ARCSIN(A/X)): CCCC CCCC EXP(-L)*L_IA(X), IF X >=ABS(A) CCCC SL_IA(X) = CCCC EXP(-ABS(A)*PI/2)*L_IA(X), IF X < ABS(A) CCCC CCCC CCCC EXP(-L)*L'_IA(X), IF X >=ABS(A) CCCC SL'_IA(X) = CCCC EXP(-ABS(A)*PI/2)*L'_IA(X), IF X 0.044*ABS(X-3.1)**1.9+(X-3.1) OR CCCC ABS(A) <= 25 AND X <= 60 CCCC 2) ASYMPTOTIC EXPANSION, IF ABS(A) < 0.7*X-10 CCCC 3) AIRY-TYPE ASYMPTOTIC EXPANSION, CCCC IF ABS(A) < 3.7*X-10 CCCC 4) QUADRATURES, CCCC IN THE REST OF THE PLANE (X,A) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCC DESCRIPTION OF INPUT/OUTPUT VARIABLES: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCC INPUTS: CCCC X : ARGUMENT OF THE FUNCTIONS CCC A : ORDER OF THE FUNCTIONS CCCC IFAC: INTEGER FLAG FOR THE SCALING CCCC * IFAC=1, THE CODE COMPUTES L_IA(X), L'_IA(X) CCCC * OTHERWISE, THE CODE COMPUTES SCALED L_IA(X), L'_IA(X) CCCC OUTPUTS: CCCC DLI : LIA(X) FUNCTION CCCC DLID : DERIVATIVE OF THE LIA(X) FUNCTION CCCC IERRO: ERROR FLAG CCCC * IF IERRO=0, COMPUTATION SUCCESSFUL. CCCC * IF IERRO=1, COMPUTATION OUT OF RANGE. CCCC * IF IERRO=2, ARGUMENT X AND/OR ORDER A OUT OF RANGE. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCC ACCURACY: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCC THE RELATIVE ACCURACY IS: CCCC BETTER THAN 10**(-13) FOR (X,A) IN (0,200]X[-200,200]; CCCC BETTER THAN 5.10**(-13) FOR (X,A) IN (0,500]X[-500,500]; CCCC CLOSE TO 10**(-12) FOR (X,A) IN (0,1500]X[-1500,1500]. CCCC NEAR THE ZEROS OF THE FUNCTIONS (THERE ARE INFINITELY CCCC MANY OF THEM IN THE OSCILLATORY REGION) RELATIVE PRECISION CCCC LOOSES MEANING AND ONLY ABSOLUTE PRECISION MAKES SENSE. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C AUTHORS: C AMPARO GIL (U. CANTABRIA, SANTANDER, SPAIN). C E-MAIL: AMPARO.GIL@UNICAN.ES C JAVIER SEGURA (U. CANTABRIA, SANTANDER, SPAIN). C E-MAIL: SEGURAJJ@UNICAN.ES C NICO M. TEMME (CWI, AMSTERDAM, THE NETHERLANDS). C E-MAIL: NICO.TEMME@CWI.NL CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C REFERENCES: C THIS IS THE COMPANION SOFTWARE OF THE ARTICLES C 1)'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL DIFFERENTIAL C EQUATION FOR IMAGINARY ORDERS AND POSITIVE ARGUMENTS', C A. GIL, J. SEGURA, N.M. TEMME C ACM TRANS. MATH. SOFT. (2004) C 2)'MODIFIED BESSEL FUNCTIONS OF IMAGINARY ORDER AND C POSITIVE ARGUMENT', C A. GIL, J. SEGURA, N.M. TEMME C ACM TRANS. MATH. SOFT. (2004) C ADDITIONAL REFERENCES: C - 'COMPUTATION OF THE MODIFIED BESSEL FUNCTION OF THE C THIRD KIND FOR IMAGINARY ORDERS' C A. GIL, J. SEGURA, N.M. TEMME C J. COMPUT. PHYS. 175 (2002) 398-411 C - 'COMPUTATION OF THE MODIFIED BESSEL FUNCTIONS OF THE C THIRD KIND OF IMAGINARY ORDERS: C UNIFORM AIRY-TYPE ASYMPTOTIC EXPANSION' C A. GIL, J. SEGURA, N.M. TEMME, PROCEEDINGS OPSFA 2001 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC LOCAL VARIABLES: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C D: X-3.1 C DF1: 0.044*ABS(D)**1.9+(X-3.1) C DF2: 0.7*X-10 C DF3: 3.7*X-10 C PNU: ORDER OF THE FUNCTION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION A,D,DF1,DF2,DF3,DLI,DLID,PNU,X INTEGER IERRO,IFAC IERRO=0 PNU=A IF ((X.GT.1500.D0).OR.(X.LE.0.D0)) THEN IERRO=2 DLI=0.D0 DLID=0.D0 ENDIF IF (ABS(PNU).GT.1500.D0) THEN IERRO=2 DLI=0.D0 DLID=0.D0 ENDIF IF (IERRO.EQ.0) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC THESE FUNCTIONS ARE EVEN FUNCTIONS IN THE C CCC PARAMETER A (=PNU) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF (PNU.LT.0.D0) PNU=-PNU D=X-3.1D0 DF1=0.044D0*ABS(D)**1.9D0+D DF2=0.7D0*X-10.D0 DF3=3.7D0*X-10.D0 IF ((PNU.GT.DF1).OR.(((PNU.LE.25.D0).AND.(X.LE.60.D0)))) + THEN CALL SERLIA(IFAC,X,PNU,DLI,DLID,IERRO) ELSEIF (PNU.LT.DF2) THEN CALL EXPLIA(IFAC,X,PNU,DLI,DLID,IERRO) ELSEIF (PNU.LT.DF3) THEN CALL AIEXLI(IFAC,X,PNU,DLI,DLID,IERRO) ELSE CALL DLINT(IFAC,X,PNU,DLI,DLID,IERRO) ENDIF ENDIF RETURN END SUBROUTINE DLINT(IFAC,XX,PNUA,DLINF,DLIND,IERRO) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC CALCULATION OF L, L' USING NON-OSCILLATING INTEGRAL C CCC REPRESENTATIONS C CCC C CCC INPUT: C CCC XX: ARGUMENT OF THE FUNCTIONS C CCC PNUA: ORDER OF THE FUNCTIONS C CCC IFAC: C CCC =1, NON SCALED FUNCTIONS C CCC OTHERWISE, SCALED FUNCTIONS C CCC OUTPUT: C CCC DLINF: L FUNCTION C CCC DLIND: DERIVATIVE OF THE L FUNCTION C CCC IERRO: ERROR FLAG C CCC * IF IERRO=0, COMPUTATION SUCCESSFUL. C CCC * IF IERRO=1, COMPUTATION OUT OF RANGE. C CCC * IF IERRO=2, ARGUMENT AND/OR ORDER, OUT OF RANGE. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC METHOD OF COMPUTATION: CCC * XX=1, 2/3*(-PSI)**3/2=SQRT(Z**2-1)-ARCCOS(1/Z) C PSI12 : SQRT(PSI) C PSI2 : PSI*PSI C PSI3 : PSI**3 C SAS : ACCUMULATES THE CONTRIBUTION OF A_S(PSI) C SBS : ACCUMULATES THE CONTRIBUTION OF B_S(PSI) C SCS : ACCUMULATES THE CONTRIBUTION OF C_S(PSI) C SDS : ACCUMULATES THE CONTRIBUTION OF D_S(PSI) C SIG : (-)**K C SINTH : SIN(THET) C THET : ASIN(A/X) C Y : Z-1 C Z : X/A C Z2 : Z**2 C Z21M : 1-Z**2 C ZDS : SQRT(1-Z**2) C ZMASF : EXPANSION IN TERMS OF (Z-1) FOR THE C CALCULATION OF PSI(Z) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION A,A0EX,A13,A2,A23,A2K,AC(0:5,0:20), + APIHAL,ARG,AS,ASP,B0EX,B0PEX,BC(0:5,0:20),BII,BIR, + BS,BSO,BSP,BSPO,C0EX,CHI,CHIEX,CHIN(0:26),COSTH, + COZMAS(0:12),D0EX,D1MACH,DBII,DBIR,DF,DLAI,DLAID DOUBLE PRECISION DZZ,ETA,ETAJ,ETAK,ETAL,EXPAM,EXPAPI, + F2,F4,FAC,FACD,FDOMIN,OVER,PHI(0:35),PHIEX,PHIEX2, + PHIS,PI,PIHALF,PSI,PSI12,PSI2,PSI3,SAS,SBS,SCS,SDS, + SIG,SINTH,THET,X,Y,Z,Z2,Z21M,ZDS,ZMASF INTEGER IERRO,IERROL,IFAC,IFACA,IFUN,INDA(0:5),INDB(0:5), + J,K,L SAVE PHI,CHIN,AC,BC CCCC VALUES OF THE COEFFICIENTS DATA PHI/1.D0,0.2D0,.25714285714285714286D-1, + -.56507936507936507937D-2,-.39367965367965367965D-2, + -.5209362066504924D-3,.3708541264731741D-3, + .2123827840293627D-3,.2150629788753145D-4, + -.2636904062981799D-4,-.1405469826493129D-4, + -.1149328899029441D-5,.1972641193938624D-5, + .1014324305593329D-5,.7034331100192200D-7, + -.1525044777392676D-6,-.7677866256900572D-7, + -.4669842638693018D-8,.1206673645965329D-7, + .5990877668092344D-8,.3269102150077715D-9, + -.97138350244428335095D-9,-.47745934295232233834D-9, + -.23750318642839155779D-10,.79244598109106655567D-10, + .38653584230817865528D-10,.17734105846426807873D-11, + -.65332110030864751956D-11,-.31673512094772575686D-11, + -.13524195177121030660D-12,.54324103217903338951D-12, + .26204918647967626464D-12,.10488134973664584898D-13, + -.45490420121539001435D-13,-.21851238232690420878D-13, + -.82456080260379042800D-15/ DATA CHIN/.2D0,.1142857142857142857142857D-1, + -.2438095238095238095238095D-1,-.1003471449185734900020614D-1, + .8811404468547325690182833D-3,.2318339655809043564145605D-2, + .7794494413564441575646057D-3,-.1373952504077228744949558D-3, + -.2162301322540308393022684D-3,-.6585004634375583559042795D-4, + .1502851367097217578058824D-4,.1991904617871647675455262D-4, + .5731972706719818525615427D-5,-.1496427612747891044606885D-5, + -.1821231830428939670133992D-5,-.5054254875930882534538217D-6, + .1433283859497625931203415D-6,.1657681319078678321113361D-6, + .4485592642702540575627044D-7,-.1346138242826094098161839D-7, + -.1504601862773535117708677D-7,-.3995660198654805921651406D-8, + .1250124931952495738882300D-8,.1363187174221864073749532D-8, + .3567608459777506132029204D-9,-.1152749290139859167119863D-9, + -.1233547289799408517691696D-9/ DATA COZMAS/.9428090415820634D0,-.4242640687119285D0, + .2904188565587606D0,-.2234261009999160D0, + .1821754153944745D0,-.1539625154198624D0, + .1333737583085222D0,-.1176617834148007D0, + .1052687194772381D0,-.9524025714638822D-1, + .8695738197073783D-1,-.8000034897653656D-1, + .7407421797273086D-1/ DATA AC(0,0)/1.D0/ DATA AC(1,0)/-.44444444444444444445D-2/ DATA AC(1,1)/-.18441558441558441558D-2/ DATA AC(1,2)/.11213675213675213675D-2/ DATA AC(1,3)/.13457752124418791086D-2/ DATA AC(1,4)/.0003880626562979504D0/ DATA AC(1,5)/-.0001830686723781799D0/ DATA AC(1,6)/-.0001995460887806733D0/ DATA AC(1,7)/-.00005256191234041587D0/ DATA AC(1,8)/.00002460619652459158D0/ DATA AC(1,9)/.00002519246780924541D0/ DATA AC(1,10)/.6333157376533242D-5/ DATA AC(1,11)/-.2957485733830202D-5/ DATA AC(1,12)/-.2925255920564838D-5/ DATA AC(1,13)/-.7159702610502009D-6/ DATA AC(1,14)/.3331510720390949D-6/ DATA AC(1,15)/.3227670475692310D-6/ DATA AC(1,16)/.7767729381664199D-7/ DATA AC(1,17)/-.3600954237921120D-7/ DATA AC(1,18)/-.3441724449034226D-7/ DATA AC(1,19)/-.8188194356398772D-8/ DATA AC(1,20)/.3783148485152038D-8/ DATA AC(2,0)/.69373554135458897365D-3/ DATA AC(2,1)/.46448349036584330703D-3/ DATA AC(2,2)/-.42838130171535112460D-3/ DATA AC(2,3)/-.0007026702868771135D0/ DATA AC(2,4)/-.0002632580046778811D0/ DATA AC(2,5)/.0001663853666288703D0/ DATA AC(2,6)/.0002212087687818584D0/ DATA AC(2,7)/.00007020345615329662D0/ DATA AC(2,8)/-.00004000421782540614D0/ DATA AC(2,9)/-.00004786324966453962D0/ DATA AC(2,10)/-.00001394600741473631D0/ DATA AC(2,11)/.7536186591273727D-5/ DATA AC(2,12)/.8478502161067667D-5/ DATA AC(2,13)/.2345355228453912D-5/ DATA AC(2,14)/-.1225943294710883D-5/ DATA AC(2,15)/-.1325082343401027D-5/ DATA AC(2,16)/-.3539954776569997D-6/ DATA AC(2,17)/.1808291719376674D-6/ DATA AC(2,18)/.1900383515233655D-6/ DATA AC(3,0)/-.35421197145774384076D-3/ DATA AC(3,1)/-.31232252789031883276D-3/ DATA AC(3,2)/.3716442237502298D-3/ DATA AC(3,3)/.0007539269155977733D0/ DATA AC(3,4)/.0003408300059444739D0/ DATA AC(3,5)/-.0002634968172069594D0/ DATA AC(3,6)/-.0004089275726648432D0/ DATA AC(3,7)/-.0001501108759563460D0/ DATA AC(3,8)/.00009964015205538056D0/ DATA AC(3,9)/.0001352492955751283D0/ DATA AC(3,10)/.00004443117087272903D0/ DATA AC(3,11)/-.00002713205071914117D0/ DATA AC(3,12)/-.00003396796969771860D0/ DATA AC(3,13)/-.00001040708865273043D0/ DATA AC(3,14)/.6024639065414414D-5/ DATA AC(3,15)/.7143919607846883D-5/ DATA AC(4,0)/.378194199201773D-3/ DATA AC(4,1)/.000404943905523634D0/ DATA AC(4,2)/-.000579130526946453D0/ DATA AC(4,3)/-.00138017901171011D0/ DATA AC(4,4)/-.000722520056780091D0/ DATA AC(4,5)/.000651265924036825D0/ DATA AC(4,6)/.00114674563328389D0/ DATA AC(4,7)/.000474423189340405D0/ DATA AC(4,8)/-.000356495172735468D0/ DATA AC(4,9)/-.000538157791035111D0/ DATA AC(4,10)/-.000195687390661225D0/ DATA AC(4,11)/.000132563525210293D0/ DATA AC(4,12)/.000181949256267291D0/ DATA AC(5,0)/-.69114139728829416760D-3/ DATA AC(5,1)/-.00085995326611774383285D0/ DATA AC(5,2)/.0014202335568143511489D0/ DATA AC(5,3)/.0038535426995603052443D0/ DATA AC(5,4)/.0022752811642901374595D0/ DATA AC(5,5)/-.0023219572034556988366D0/ DATA AC(5,6)/-.0045478643394434635622D0/ DATA AC(5,7)/-.0020824431758272449919D0/ DATA AC(5,8)/.0017370443573195808719D0/ DATA BC(0,0)/.14285714285714285714D-1/ DATA BC(0,1)/.88888888888888888889D-2/ DATA BC(0,2)/.20482374768089053803D-2/ DATA BC(0,3)/-.57826617826617826618D-3/ DATA BC(0,4)/-.60412089799844901886D-3/ DATA BC(0,5)/-.0001472685745626922D0/ DATA BC(0,6)/.00005324102148009784D0/ DATA BC(0,7)/.00005206561006583416D0/ DATA BC(0,8)/.00001233115050894939D0/ DATA BC(0,9)/-.4905932728531366D-5/ DATA BC(0,10)/-.4632230987136350D-5/ DATA BC(0,11)/-.1077174523455235D-5/ DATA BC(0,12)/.4475963978932822D-6/ DATA BC(0,13)/.4152586188464624D-6/ DATA BC(0,14)/.9555819293589234D-7/ DATA BC(0,15)/-.4060599208403059D-7/ DATA BC(0,16)/-.3731367187988482D-7/ DATA BC(0,17)/-.8532670645553778D-8/ DATA BC(0,18)/.3673017245573624D-8/ DATA BC(0,19)/.3355960460784536D-8/ DATA BC(0,20)/.7643107095110475D-9/ DATA BC(1,0)/-.11848595848595848596D-2/ DATA BC(1,1)/-.13940630797773654917D-2/ DATA BC(1,2)/-.48141005586383737645D-3/ DATA BC(1,3)/.26841705366016142958D-3/ DATA BC(1,4)/.0003419706982709903D0/ DATA BC(1,5)/.0001034548234902078D0/ DATA BC(1,6)/-.00005418191982095504D0/ DATA BC(1,7)/-.00006202184830690167D0/ DATA BC(1,8)/-.00001724885886056087D0/ DATA BC(1,9)/.8744675992887053D-5/ DATA BC(1,10)/.9420684216180929D-5/ DATA BC(1,11)/.2494922112085850D-5/ DATA BC(1,12)/-.1238458608836357D-5/ DATA BC(1,13)/-.1285461713809769D-5/ DATA BC(1,14)/-.3299710862537507D-6/ DATA BC(1,15)/.1613441105788315D-6/ DATA BC(1,16)/.1633623194402374D-6/ DATA BC(1,17)/.4104252949605779D-7/ DATA BC(1,18)/-.1984317042326989D-7/ DATA BC(1,19)/-.1973948142769706D-7/ DATA BC(1,20)/-.4882194808588752D-8/ DATA BC(2,0)/.43829180944898810994D-3/ DATA BC(2,1)/.71104865116708668943D-3/ DATA BC(2,2)/.31858383945387580576D-3/ DATA BC(2,3)/-.0002404809426804458D0/ DATA BC(2,4)/-.0003722966038621536D0/ DATA BC(2,5)/-.0001352752059595618D0/ DATA BC(2,6)/.00008691694372704142D0/ DATA BC(2,7)/.0001158750753591377D0/ DATA BC(2,8)/.00003724965927846447D0/ DATA BC(2,9)/-.00002198334949606935D0/ DATA BC(2,10)/-.00002686449633870452D0/ DATA BC(2,11)/-.8023061612032524D-5/ DATA BC(2,12)/.4494756592180126D-5/ DATA BC(2,13)/.5193504763856015D-5/ DATA BC(2,14)/.1477156191529617D-5/ DATA BC(2,15)/-.7988793826096784D-6/ DATA BC(3,0)/-.37670439477105454219D-3/ DATA BC(3,1)/-.75856271658798642365D-3/ DATA BC(3,2)/-.0004103253968775048D0/ DATA BC(3,3)/.0003791263310429010D0/ DATA BC(3,4)/.0006850981673903450D0/ DATA BC(3,5)/.0002878310571932216D0/ DATA BC(3,6)/-.0002157010636115705D0/ DATA BC(3,7)/-.0003260863991373500D0/ DATA BC(3,8)/-.0001181317008748678D0/ DATA BC(3,9)/.00007887526841582582D0/ DATA BC(3,10)/.0001072081833420685D0/ DATA BC(3,11)/.00003544595251288735D0/ DATA BC(3,12)/-.00002201447920733824D0/ DATA BC(3,13)/-.00002789336359620813D0/ DATA BC(4,0)/.00058453330122076187255D0/ DATA BC(4,1)/.0013854690422372401251D0/ DATA BC(4,2)/.00086830374184946900245D0/ DATA BC(4,3)/-.00093502904801345951693D0/ DATA BC(4,4)/-.0019175486005525492059D0/ DATA BC(4,5)/-.00090795047113308137941D0/ DATA BC(4,6)/.00077050429806392235104D0/ DATA BC(4,7)/.0012953100128255983488D0/ DATA BC(4,8)/.00051933869471899550762D0/ DATA BC(4,9)/-.00038482631948759834653D0/ DATA BC(4,10)/-.00057335393099012476502D0/ DATA BC(5,0)/-.0014301070053470410656D0/ DATA BC(5,1)/-.0038637811942002539408D0/ DATA BC(5,2)/-.0027322816261168328245D0/ DATA BC(5,3)/.0033294980346743452748D0/ DATA BC(5,4)/.0075972237660887795911D0/ DATA BC(5,5)/.0039816655673062060620D0/ DATA BC(5,6)/-.0037510180460986006595D0/ DATA INDA/0,20,18,15,12,8/ DATA INDB/20,20,15,13,10,6/ CCCC CONSTANTS CCCCCCCCCCCCCC PI=ACOS(-1.D0) CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER) OVER=D1MACH(2)*1.D-8 IERROL=0 IF (IFAC.EQ.1) THEN IF (X.GE.A) THEN SINTH=A/X THET=ASIN(SINTH) COSTH=COS(THET) FDOMIN=X*(COSTH+THET*SINTH) IF (FDOMIN.GE.LOG(OVER)) IERROL=1 ELSE IF ((PI*A*0.5D0).GE.LOG(OVER)) IERROL=1 ENDIF ENDIF IF (IERROL.EQ.0) THEN PIHALF=PI*0.5D0 DF=2.D0**(1.D0/3.D0) CCCC VARIABLES CCCCCCCCCCCCCCC Z=X/A A2=A*A A13=A**(1.D0/3.D0) A23=A13*A13 IF (IFAC.EQ.1) THEN APIHAL=A*PIHALF EXPAPI=EXP(APIHAL) FAC=EXPAPI*0.5D0/A13 FACD=EXPAPI/Z/A23 ELSE IF (X.LT.A) THEN FAC=0.5D0/A13 FACD=1.D0/Z/A23 ELSE SINTH=A/X THET=ASIN(SINTH) COSTH=COS(THET) FDOMIN=X*(COSTH+THET*SINTH) EXPAM=EXP(A*PIHALF-FDOMIN) FAC=EXPAM*0.5D0/A13 FACD=EXPAM/Z/A23 ENDIF ENDIF F4=A13**4 IF (Z.LE.0.9D0) THEN ZDS=SQRT((1.D0-Z)*(1.D0+Z)) PSI=(1.5D0*(LOG((1.D0+ZDS)/Z)-ZDS))**(2.D0/3.D0) ELSEIF (Z.GT.1.1D0) THEN ZDS=SQRT((Z-1.D0)*(1.D0+Z)) PSI=-(1.5D0*(ZDS-ACOS(1.D0/Z)))**(2.D0/3.D0) ELSE Y=Z-1.D0 ZMASF=COZMAS(12) DO 8 K=0,11 J=11-K ZMASF=COZMAS(J)+Y*ZMASF 8 CONTINUE ZMASF=ABS(Y)**(1.5D0)*ZMASF IF (Z.LT.1.D0) THEN PSI=(1.5D0*ZMASF)**(2.D0/3.D0) ELSE PSI=-(1.5D0*ZMASF)**(2.D0/3.D0) ENDIF ENDIF ETA=PSI/DF ARG=-A23*PSI PSI2=PSI*PSI PSI3=PSI2*PSI IF ((Z.GT.0.8D0).AND.(Z.LT.1.2D0)) THEN PHIS=0.D0 CHI=0.D0 SAS=0.D0 SBS=0.D0 SDS=0.D0 SCS=1.D0 BS=0.D0 BSPO=0.D0 DO 41 L=0,20 IF (L.EQ.0) THEN ETAL=1.D0 ELSE ETAL=ETAL*ETA ENDIF CHI=CHI+CHIN(L)*ETAL 41 CONTINUE CHI=CHI/DF DO 10 K=0,20 BSO=BS BSPO=BSP AS=0.D0 BS=0.D0 ASP=0.D0 BSP=0.D0 IF (K.EQ.0) THEN ETAK=1.D0 A2K=1.D0 SIG=1.D0 ELSE ETAK=ETAK*ETA A2K=A2K*A2 SIG=-1.D0*SIG ENDIF PHIS=PHIS+PHI(K)*ETAK F2=SIG/A2K IF (K.LE.5) THEN DO 20 J=0,20 IF (J.EQ.0) THEN ETAJ=1.D0 ELSE ETAJ=ETAJ*ETA ENDIF IF (J.LE.INDA(K)) THEN AS=AS+AC(K,J)*ETAJ ENDIF IF (J.LE.INDB(K)) THEN BS=BS+BC(K,J)*ETAJ ENDIF IF ((J+1).LE.INDA(K)) THEN ASP=ASP+(J+1)*AC(K,J+1)*ETAJ ENDIF IF ((J+1).LE.INDB(K)) THEN BSP=BSP+(J+1)*BC(K,J+1)*ETAJ ENDIF 20 CONTINUE ASP=1.D0/DF*ASP BS=BS*DF SAS=SAS+AS*F2 SBS=SBS+BS*F2/DF SDS=SDS-(CHI*AS+ASP+PSI*BS)*F2 ENDIF IF ((K.GT.0).AND.(K.LE.6)) THEN SCS=SCS+(AS+CHI*BSO+BSPO)*F2 ENDIF 10 CONTINUE PHIS=DF*PHIS SBS=DF*SBS ELSE CCCC EXACT VALUES CCCCCCCCCCCCCCCCCCCCC Z2=Z*Z Z21M=1.D0-Z2 PHIEX=(4.D0*PSI/Z21M)**0.25D0 PHIEX2=PHIEX*PHIEX A0EX=1.D0 B0EX=-5.D0/48.D0/(PSI*PSI)+PHIEX2/48.D0/PSI* + (5.D0/Z21M-3.D0) CHIEX=0.25D0/PSI*(1.D0-Z2*PHIEX2**3*0.25D0) CHI=CHIEX D0EX=-(7.D0/48.D0/PSI+PHIEX2/48.D0* + (9.D0-7.D0/Z21M)) IF (PSI.GT.0.D0) THEN PSI12=SQRT(PSI) DZZ=SQRT(Z21M) ELSE PSI12=SQRT(-PSI) DZZ=SQRT(ABS(Z21M)) ENDIF B0PEX=5.D0/24.D0/PSI3+PHIEX2/48.D0*((2.D0*CHI*PSI-1.D0)/ + PSI2*(5.D0/Z21M-3.D0)-10.D0*Z2*PSI12/ + Z21M**2/DZZ/PSI) C0EX=1.D0 SAS=A0EX SBS=B0EX SCS=C0EX SDS=D0EX BS=0.D0 BSP=0.D0 A2K=1.D0 SIG=1.D0 DO 30 K=1,6 BSO=BS BSPO=BSP AS=0.D0 BS=0.D0 ASP=0.D0 BSP=0.D0 A2K=A2K*A2 SIG=-1.D0*SIG F2=SIG/A2K IF (K.LE.5) THEN DO 40 J=0,20 IF (J.EQ.0) THEN ETAJ=1.D0 ELSE ETAJ=ETAJ*ETA ENDIF IF (J.LE.INDA(K)) THEN AS=AS+AC(K,J)*ETAJ ENDIF IF (J.LE.INDB(K)) THEN BS=BS+BC(K,J)*ETAJ ENDIF IF ((J+1).LE.INDA(K)) THEN ASP=ASP+(J+1)*AC(K,J+1)*ETAJ ENDIF IF ((J+1).LE.INDB(K)) THEN BSP=BSP+(J+1)*BC(K,J+1)*ETAJ ENDIF 40 CONTINUE BS=DF*BS ASP=ASP/DF SAS=SAS+AS*F2 SBS=SBS+BS*F2 SDS=SDS-(CHI*AS+ASP+PSI*BS)*F2 ENDIF IF (K.GE.2) THEN SCS=SCS+(AS+CHI*BSO+BSPO)*F2 ELSE SCS=SCS+(AS+CHI*B0EX+B0PEX)*F2 ENDIF 30 CONTINUE PHIS=PHIEX ENDIF CCCCC CALL THE AIRY ROUTINE CCCCCCCCC IFUN=1 IFACA=1 CALL BIZ(IFUN,IFACA,ARG,0.D0,BIR,BII,IERRO) IFUN=2 IFACA=1 CALL BIZ(IFUN,IFACA,ARG,0.D0,DBIR,DBII,IERRO) DLAI=FAC*PHIS*(BIR*SAS+DBIR*SBS/F4) DLAID=FACD/PHIS*(DBIR*SCS+BIR*SDS/A23) ELSE DLAI=0.D0 DLAID=0.D0 ENDIF RETURN END DOUBLE PRECISION FUNCTION V(U) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC AUXILIARY FUNCTION: CCC THIS FUNCTION COMPUTES THE LAST EXPRESSION CCC IN EQ.(21) OF THE REFERENCE CCC 'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL CCC DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS CCC AND POSITIVE ARGUMENTS', CCC A. GIL, J. SEGURA, N.M. TEMME CCC ACM TRANS. MATH. SOFT. (2004) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC INPUT VARIABLE: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C U : ARGUMENT OF THE FUNCTION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC LOCAL VARIABLES: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COSTH : COS(THET) C FDOMIN : X*(COS(THET)+THET*SIN(THET)) C SINTH : SIN(THET) C THET : ASIN(PNU/X) C UNDER : UNDERFLOW NUMBER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION COSTH,D1MACH,FDOMIN,SINTH, + THET,U,UNDER COMMON/PARMON/THET,SINTH,COSTH,FDOMIN CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER) UNDER=D1MACH(1)*1.D+6 IF (ABS(U).LT.UNDER) THEN V=ASIN(SINTH) ELSE V=ASIN(U/SINH(U)*SINTH) ENDIF RETURN END DOUBLE PRECISION FUNCTION PHIR(U) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC AUXILIARY FUNCTION: CCC THIS FUNCTION COMPUTES THE EXPONENT IN THE CCC INTEGRAND IN EQ.(20) OF THE REFERENCE CCC 'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL CCC DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS CCC AND POSITIVE ARGUMENTS', CCC A. GIL, J. SEGURA, N.M. TEMME CCC ACM TRANS. MATH. SOFT. (2004) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC INPUT VARIABLE: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C U : ARGUMENT OF THE FUNCTION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC LOCAL VARIABLES: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COSTH : COS(THET) C COSVU : COS(V(U)) C FDOMIN : X*(COS(THET)+THET*SIN(THET)) C FU1 : 2*SINH(U/2)**2 C FU2 : -2*SIN(FU3/2)*SIN(0.5*(V(U)+THET)) C FU3 : ASIN(-SIN(THET)/(COS(THET)*U/SINH(U) C +COS(V(U))*(SINH(U)+U)*FUAC/(SINH(U)*SINH(U))) C FUAC : U**3/6+U**5/120+U**7/5040 C OVER : OVERFLOW NUMBER C PNU : ORDER OF THE FUNCTION C SINHU : SINH(U) C SINHUH : SINH(U/2) C SINTH : SIN(THET) C THET : ASIN(PNU/X) C U2 : U**2 C U3 : U**3 C U5 : U**5 C U7 : U**7 C UH : U/2 C UNDER : UNDERFLOW NUMBER C V(U) : ASIN(U/SINH(U)*SIN(THET)) C VU : V(U) C X : ARGUMENT OF THE FUNCTIONS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION COSTH,COSVU,D1MACH,FDOMIN,FU1,FU2, + FU3,FUAC,OVER,PNU,SINHU,SINHUH,SINTH,THET,U,U2, + U3,U5,U7,UH,UNDER,V,VU,X COMMON/ARGU/X,PNU COMMON/PARMON/THET,SINTH,COSTH,FDOMIN CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER) UNDER=D1MACH(1)*1.D+6 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER) OVER=D1MACH(2)*1.D-6 IF (U.LT.200.D0) THEN IF (ABS(U).LT.0.1D0) THEN IF (ABS(U).LT.UNDER) THEN PHIR=0.D0 ELSE UH=U*0.5D0 SINHUH=SINH(UH) FU1=2.D0*SINHUH*SINHUH U2=U*U U3=U2*U U5=U3*U2 U7=U5*U2 FUAC=U3/6.D0+U5/120.D0+U7/5040.D0 SINHU=SINH(U) VU=V(U) COSVU=COS(VU) FU3=ASIN(-SINTH/(COSTH*U/SINHU+COSVU)* + (SINHU+U)*FUAC/(SINHU*SINHU)) FU2=-2.D0*SIN(FU3*0.5D0)*SIN(0.5D0*(VU+THET)) PHIR=X*(FU1*COSVU+FU2+SINTH*FU3) ENDIF ELSE VU=V(U) COSVU=COS(VU) PHIR=X*COSH(U)*COSVU+PNU*VU-FDOMIN ENDIF ELSE PHIR=OVER ENDIF RETURN END DOUBLE PRECISION FUNCTION SIGMA(U) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC AUXILIARY FUNCTION: CCC THIS FUNCTION COMPUTES THE FUNCTION CCC SIGMA FROM EQ.(31) OF THE REFERENCE CCC 'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL CCC DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS CCC AND POSITIVE ARGUMENTS', CCC A. GIL, J. SEGURA, N.M. TEMME CCC ACM TRANS. MATH. SOFT. (2004) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC INPUT VARIABLE: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C U : ARGUMENT OF THE FUNCTION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC LOCAL VARIABLES: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ARGU : (COSH(MU)*U-DMUFAC)/SINH(U) C COSCHI : COS(CHI) C COSHM : COSH(MU) C D1 : COSH(MU)*U-DMUFAC C DMU : SOLUTION MU OF COSH(MU)=PNU/X C DMUFAC : MU*COSH(MU)-SINH(MU) C SINCHI : SIN(CHI) C SINHM : SINH(MU) C SINHU : SINH(U) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION ARGU,COSCHI,COSHM,D1,DMU,DMUFAC, + PI,SINCHI,SINHM,SINHU,U COMMON/PAROS1/DMU,COSHM,SINHM,DMUFAC,COSCHI,SINCHI PI=ACOS(-1.D0) D1=COSHM*U-DMUFAC SINHU=SINH(U) ARGU=D1/SINH(U) IF (ABS(ARGU).GT.1.D0) THEN ARGU=1.D0 ENDIF SIGMA=ASIN(ARGU) IF (U.LT.DMU) SIGMA=PI-SIGMA RETURN END DOUBLE PRECISION FUNCTION PHIB(U) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC THIS FUNCTION COMPUTES THE FIRST EXPRESSION CCC OF EQ.(30) OF THE REFERENCE CCC 'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL CCC DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS CCC AND POSITIVE ARGUMENTS', CCC A. GIL, J. SEGURA, N.M. TEMME CCC ACM TRANS. MATH. SOFT. (2004) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC INPUT VARIABLE: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C U : ARGUMENT OF THE FUNCTION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC LOCAL VARIABLES: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C OVER : OVERFLOW NUMBER C PNU : ORDER OF THE FUNCTIONS C SIGMA(U) : ASIN((COSH(MU)*U-DMUFAC)/SINH(U)) C SIGMAU : SIGMA(U) C X : ARGUMENT OF THE FUNCTIONS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION U,SIGMA,SIGMAU,OVER,D1MACH, + X,PNU COMMON/ARGU/X,PNU CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER) OVER=D1MACH(2)*1.D-8 IF (U.LT.200.D0) THEN SIGMAU=SIGMA(U) PHIB=X*COSH(U)*COS(SIGMAU)+PNU*SIGMAU ELSE PHIB=OVER ENDIF RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. # End of shell archive exit 0