#! /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/ # Doc/airy_README # Doc/makefile # Src/ # Src/Fortran90/ # Src/Fortran90/Drivers/ # Src/Fortran90/Drivers/test_install.f90 # Src/Fortran90/Src/ # Src/Fortran90/Src/airy_complex # Src/Fortran90/Src/airy_functions.f90 # Src/Fortran90/Src/airy_head # Src/Fortran90/Src/airy_parameters # Src/Fortran90/Src/airy_real # Src/Fortran90/Src/airyerr.log # Src/Fortran90/Src/res # Src/Fortran90/Src/test_install.f90 # This archive created: Mon Mar 7 16:41:34 2005 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'airy_README' then echo shar: will not over-write existing file "'airy_README'" else cat << "SHAR_EOF" > 'airy_README' User's Guide for TOMS Algorithm XXX: Airy Functions Bruce Fabijonas Southern Methodist University bfabi@smu.edu Introduction This guide is a companion to B.R. Fabijonas. Algorithm XXX: Airy Functions. ACM Transactions on Mathematical Software Vol(No) xx-xx, 2004, based on the methods explained in B.R. Fabijonas, D.W. Lozier and F.W.J. Olver. Computation of Complex Airy functions and zeros using asymptotics and the differential equation. ACM Transactions on Mathematical Software Vol(No) xx-xx, 2004. Installation and usage instructions for the software described in this paper are given here. The functions computed are the Airy functions Ai(z), Ai(x), Bi(x) where z is complex and x is real, the derivatives of these func- tions, and the auxiliary modulus and phase functions when x is negative. Also, the zeros of Ai(x), Ai'(x), Bi(x), Bi'(x), Bi(z) and Bi'(z), are computed. The programs, written in Fortran-90, are designed to accommodate precisions up to 64 decimals. Single and double precision are obtained automatically by compiling the package as distributed. Quadruple precision is available for compilers that support it by "uncommenting" certain well-marked sections of code. This package has been tested with the following compilers: Compaq Fortran Compiler V5.5-1877-48BBF GNU Fortran 95(GCC 3.5-tree-ssa 20021226 (G95)) NAG Fortran 95 Release 5.0(298) Files The distributed files are airy_README this document airy_functions.f90 TOMS Algorithm XXX airy_head INCLUDE file airy_parameters INCLUDE file airy_real INCLUDE file airy_complex INCLUDE file test_install.f90 program to test installation parameter_info.f90 program to return parameters specific to the compiler and computer being used sample_airy.f90 program to demonstrate versatility When the TOMS algorithm is compiled, for example by the command f90 -c airy_functions.f90 the resulting files are airy_functions.o airy_functions_real_single.mod airy_functions_real_double.mod airy_functions_complex_single.mod airy_functions_complex_double.mod These may be gathered together into a library, for example by the command ar crv airy_lib.a *.o *.mod Main Program Example A program foo.f90 using any or all of these modules must include the corresponding USE statements, e.g. program foo use airy_functions_real_single and may be compiled with the appropriate object files or li- brary, for example by f90 -o foo foo.f90 airy_functions.o or f90 -o foo foo.f90 airy_lib.a Fortran-90's INTERFACE capability has been used in the coding of airy_functions.f90. This allows the same subroutine call to be used with different Fortran-90 KINDs and TYPEs. For example, program foo use airy_functions_real_single use airy_functions_complex_double real(kind(0.0e0)) xr, air, dair complex(kind(0.0d0)) xc, aic, daic integer ierr, ierc call airy_ai( xr, air, dair, ierr ) call airy_ai( xc, aic, daic, ierc ) end program foo returns the real-valued single precision Airy function Ai(xr) in the variable air, its derivative Ai'(xr) in the variable dair, the complex-valued double precision Airy function Ai(xc) in the variable aic, and its derivative Ai'(xc) in the variable daic. Usage of the Real Subroutines The single and double precision real subroutines are contained in the modules airy_functions_real_single.mod and airy_functions_re- al_double.mod, respectively. The calling Fortran 90 program must contain an appropriate USE statement; see Main Program Example above. The subroutines are airy_ai( x, ai, dai, ierr, modify_switch= ) airy_bi( x, bi, dbi, ierr, modify_switch= ) airy_aux( x, m_mod, theta, ierr, n_mod =, phi= ) airy_ai_zero( n, ai_zero, ierr, ai_assoc=, dai_zero=, dai_assoc= ) airy_bi_zero( n, bi_zero, ierr, bi_assoc=, dbi_zero=, dbi_assoc= ) airy_info( radius, max_step_size=, n_terms_taylor=, n_terms_asymp=, n_partitions= ) airy_aux_info( radius, max_step_size=, n_terms_taylor=, n_terms_asympc, n_terms_asymp_phase=, n_terms_asymp_mod=, n_partitions= ) where x argument TYPE(REAL), INTENT(IN), KIND an allowable real kind ai, dai, bi, dbi s(x)*Ai(x), s(x)*Ai'(x), s(x)*Bi(x), s(x)*Bi'(x) where s(x) is determined by MODIFY_SWITCH (described below) TYPE(REAL), INTENT(OUT), KIND same as x ierr error flag (described below) TYPE(INTEGER), INTENT(OUT) modify_switch default value: .FALSE. if .FALSE., s(x) = 1 if .TRUE. and x>0, s(x) = exp((2/3)*x^(3/2)) for airy_ai and s(x) = exp((-2/3)*x^(3/2)) for airy_bi if .TRUE. and x<=0, ierr is set to -1 TYPE(LOGICAL), INTENT(IN), OPTIONAL m_mod if x<=0, modulus function M(x) = sqrt(Ai^2(x)+Bi^2(x)) if x>0, ierr is set to -2 TYPE(REAL), INTENT(OUT) theta if x<=0, phase function theta(x) = arctan(Ai(x)/Bi(x)) if x>0, ierr is set to -2 TYPE(REAL), INTENT(OUT) n_mod if x<=0, modulus function N(x) = sqrt((Ai'(x))^2+(Bi'(x))^2) if x>0, ierr is set to -2 TYPE(REAL), INTENT(OUT), OPTIONAL phi if x<=0, phase function phi(x) = arctan(Ai'(x)/Bi'(x)) if x>0, ierr is set to -2 TYPE(REAL), INTENT(OUT), OPTIONAL n if n>=1, the index of the desired zero if n<=0, ierr is set to -3 TYPE(INTEGER), INTENT(IN) ai_zero if a scalar, the nth zero of Ai(x) TYPE(REAL), INTENT(OUT) if a vector, the nth through the (n+k-1)st zeros of Ai(x) TYPE(REAL), DIMENSION(k), INTENT(OUT) dai_zero if a scalar, the nth zero of Ai'(x) TYPE(REAL), INTENT(OUT) if a vector, the nth through the (n+k-1)st zeros of Ai'(x) TYPE(REAL), DIMENSION(k), INTENT(OUT) ai_assoc if a scalar, Ai'(a_n) where a_n is the nth zero of Ai(x) TYPE(REAL), INTENT(OUT) if a vector, values of Ai'(a) at zeros n to n+k-1 of Ai(x) TYPE(REAL), DIMENSION(k), INTENT(OUT) dai_assoc if a scalar, Ai(a_n) where a_n is the nth zero of Ai'(x) TYPE(REAL), INTENT(OUT) if a vector, values of Ai(a) at zeros n to n+k-1 of Ai'(x) TYPE(REAL), DIMENSION(k), INTENT(OUT) bi_zero, dbi_zero, bi_assoc, dbi_assoc analogous to ai_zero, dai_zero, ai_assoc, dai_assoc for Bi(x) radius radius of the circle separating the integration and asymptotic expansion regions of the mathematical algorithm TYPE(REAL), INTENT(OUT) max_step_size maximum distance between partition points on an integration ray TYPE(REAL), INTENT(OUT), OPTIONAL n_terms_taylor number of terms in the Taylor series integration TYPE(INTEGER), INTENT(OUT), OPTIONAL n_terms_asymp number of asymptotic expansion terms for the Airy functions TYPE(INTEGER), INTENT(OUT), OPTIONAL n_terms_asymp_phase number of asymptotic expansion terms for the phase functions TYPE(INTEGER), INTENT(OUT), OPTIONAL n_terms_asymp_mod number of asymptotic expansion terms for the modulus functions TYPE(INTEGER), INTENT(OUT), OPTIONAL n_partitions number of partitions used in the integration routine TYPE(INTEGER), INTENT(OUT), OPTIONAL The arguments of airy_info are dynamically determined upon the first call to any of the subroutines and are here for diag- nostic purposes only. Usage of the Complex Subroutines The single and double precision complex subroutines are contained in the modules airy_functions_complex_single.mod and airy_func- tions_complex_double.mod, respectively. The calling Fortran 90 program must contain an appropriate USE statement; see Main Pro- gram Example above. The subroutines are airy_ai( z, ai, dai, ierr, argument_switch=, modify_switch= ) airy_ai( r, alpha, ai, dai, ierr, argument_switch=, modify_switch= ) airy_bi_zero( n, bi_zero, ierr, bi_assoc=, dbi_zero=, dbi_assoc= ) where z argument (see ARGUMENT_SWITCH below) TYPE(COMPLEX), INTENT(IN), KIND an allowable complex kind ai, dai if a scalar, s(z)*Ai(z), s(z)*Ai'(z) where s(z) is determined by MODIFY_SWITCH (described below) TYPE(COMPLEX), INTENT(OUT), KIND same as z if a vector, s(z)*Ai(z), s(z)*Ai'(z) where s(z) is determined by MODIFY_SWITCH (described below) note: ARGUMENT_SWITCH must be set to .TRUE. TYPE(COMPLEX), DIMENSION(k), INTENT(OUT), KIND and SIZE the same as r (see below) ierr error flag (described below) TYPE(INTEGER), INTENT(OUT) argument_switch default value: .FALSE. if .FALSE., z = (x,y) = x+i*y if .TRUE. and x>=0 and -pi=1, the index of the desired complex zero if n<=0, ierr is set to -3 TYPE(INTEGER), INTENT(IN) bi_zero if a scalar, nth complex zero of Bi(x) in the upper half-plane TYPE(COMPLEX), INTENT(OUT) if a vector, the nth through the (n+k-1)st zeros of Bi(x) TYPE(COMPLEX), DIMENSION(k), INTENT(OUT) dbi_zero if a scalar, the nth complex zero of Bi'(x) TYPE(COMPLEX), INTENT(OUT) if a vector, the nth through the (n+k-1)st zeros of Bi'(x) TYPE(COMPLEX), DIMENSION(k), INTENT(OUT) bi_assoc if a scalar, Bi'(b_n) where b_n is the nth complex zero of Bi(x) TYPE(COMPLEX), INTENT(OUT) if a vector, the values of Bi'(a) at zeros n to n+k-1 of Bi(x) TYPE(COMPLEX), DIMENSION(k), INTENT(OUT) dbi_assoc if a scalar, Bi(b_n) where b_n is the nth complex zero of Bi'(x) TYPE(COMPLEX), INTENT(OUT) if a vector, the values of Bi(a) at zeros n to n+k-1 of Bi'(x) TYPE(COMPLEX), DIMENSION(k), INTENT(OUT) Since Ai(z) and Ai'(z) have no zeros in the complex plane other than those that lie on the negative real axis, calling the sub- routine airy_ai_zero with complex arguments is an error that results in ierr being set to -3. Error and Diagnostic Flags Error flags are negative integers. They indicate that a computa- tion was requested that falls outside the valid range of the software or that invalid input was provided. Returned values are likely to be greatly in error. Error flags are listed below and noted also in the Subroutine Usage sections above. Diagnostic flags are positive integers. They indicate that a special algorithmic step was taken, often to avoid vagaries of floating-point arithmetic such as underflow or overflow. Never- theless, the computed values retain full accuracy. The flags are ierr =-50 the largest zero requested exceeds the largest computable zero; the noncomputable ones were returned as zero ierr =-10 z = r*exp(i*alpha) was provided in polar form with r < 0 and/or abs(alpha) > pi ierr = -7 the real part of (2/3)*z^(3/2) was too large in the region where the function is exponentially small; function values were set to zero to avoid underflow ierr = -6 the real part of (2/3)*z^(3/2) was too large in the region where the function is exponentially large; function values were set to zero to avoid overflow ierr = -4 machine's largest integer was exceeded ierr = -3 zeros are given for n => 1 only; or Ai(z) has no complex zeros ierr = -2 modulus and phase functions are generated for x real and x <= 0 only ierr = -1 in the real subroutines, scaling is for x > 0 only. ierr = 0 normal computation, no errors or diagnostics The vectorized subroutine airy_ai allows the user to fix a phase angle alpha along which Ai(z), Ai'(z) are wanted at a sequence of points z = r*exp(i*alpha). When applicable, the vectorized version is much faster than the linear version. Note that ARGUMENT_SWITCH must be set to TRUE, and that r, ai, dai must be vectors of the same size SHAR_EOF fi # end of overwriting check if test -f 'makefile' then echo shar: will not over-write existing file "'makefile'" else cat << "SHAR_EOF" > 'makefile' ASOURCE = . AFUNS = $(ASOURCE)/airy_functions.f90 FC = f95 #FC = g95 #FC = nag_f95 all: test_install parameter_info sample_airy airy_functions.o: $(ASOURCE)/airy_functions.f90 $(ASOURCE)/airy_real $(ASOURCE)/airy_complex $(ASOURCE)/airy_head $(ASOURCE)/airy_parameters $(FC) -c -o airy_functions.o $(ASOURCE)/airy_functions.f90 test_install: test_install.f90 airy_functions.o $(FC) -o test_install test_install.f90 airy_functions.o clean: rm -f test_install airy_functions.o rm -f *.mod airyerr.log %.o: %.f90 $(FC) -c $< -o $@ SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test ! -d 'Fortran90' then mkdir 'Fortran90' fi cd 'Fortran90' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'test_install.f90' then echo shar: will not over-write existing file "'test_install.f90'" else cat << "SHAR_EOF" > 'test_install.f90' ! updated 04/26/2004 ! This is a program which tests the installation ! The Airy function pacakge by B.R.Fabijonas. ! reference: B.R. Fabijonas, Algorithm XXX: Airy functions, ! ACM Trans. Math. Softw., Vol. xxx, pp. xxx-xxx. ! Please report all bug reports to bfabi@mail.smu.edu ! ! This programs tests all featuers of the package: ! Airy functions of real argument computed by integration of ODE ! Airy functions or real argument computed by summation of ! asymptotic expansions ! Airy functions of complex argument computed by integration of ODE ! Airy functions or complex argument computed by summation of ! asymptotic expansions ! Real and complex zeros of the Airy functions computed by table lookup ! Real and complex zeros of the Airy functions computed by ! summation of asymptotic expansions ! Vectorized version of zeros of Ai(x) and its derivative with ! associated values ! Vectorized version of zeros of Bi(x) and its derivative with ! associated values ! Vectorized version of zeros of Bi(z) and its derivative with ! associated values ! Auxiliary functions for negative real arguments computed by ! integration of ODE ! Auxiliary functions for negative real arguments computed by ! summation of asymptotic expansions ! ! A module to interface the error computation and switch ! between relative precision and absolute error. ! the endings: ! s = real single ! d = real double ! c = complex single ! z = complex double ! module error implicit none private public errg integer, parameter :: prs = kind(0.0e0) integer, parameter :: prd = kind(0.0d0) real(prd), parameter :: zerod = 0.0_prd, & pid = 3.1415926535897932384626433832795028841971_prd real(prs), parameter :: zeros = 0.0_prs, pis = pid interface errg module procedure errds module procedure errs module procedure errd module procedure errzc module procedure errz module procedure errc end interface interface psi module procedure psis module procedure psid end interface contains ! ! real double and single real(prd) function errds(arg1,arg2) real(prd) arg1 real(prs) arg2 errds = abs(psi(abs(arg1)) - psi(abs(arg2))) end function errds ! ! real single real(prs) function errs(arg1,arg2) real(prs) arg1, arg2 errs = abs(psi(arg1) - psi(arg2)) end function errs ! ! real double real(prd) function errd(arg1,arg2) real(prd) arg1, arg2 errd = abs(psi(arg1) - psi(arg2)) end function errd ! ! complex double and single real(prd) function errzc(arg1,arg2) complex(prd) arg1, arg3 complex(prs) arg2 real(prd) x1,x2 arg3 = arg2 if ( abs(arg1) < 1.0_prd .or. abs(arg3) < 1.0_prd ) then errzc = abs(arg1-arg2) else x1 = atan2(aimag(arg1),real(arg1,prd)) x2 = atan2(aimag(arg3),real(arg3,prd)) x2 = mod(abs(x1-x2),2.0_prd*pid) errzc = abs(cmplx(log(abs(arg1/arg3)),min(x2,2.0_prd*pid-x2),prd)) end if end function errzc ! ! complex single real(prs) function errc(arg1,arg2) complex(prs) arg1, arg2 real(prs) x1,x2 if ( abs(arg1) < 1.0_prs .or. abs(arg2) < 1.0_prs ) then errc = abs(arg1-arg2) else x1 = atan2(aimag(arg1),real(arg1)) x2 = atan2(aimag(arg2),real(arg2)) x2 = mod(abs(x1-x2),2.0_prs*pis) errc = abs(cmplx(log(abs(arg1/arg2)),min(x2,2.0_prs*pis-x2))) end if end function errc ! ! complex double real(prd) function errz(arg1,arg2) complex(prd) arg1, arg2 real(prd) x1,x2 if ( abs(arg1) < 1.0_prd .or. abs(arg2) < 1.0_prd ) then errz = abs(arg1-arg2) else x1 = atan2(aimag(arg1),real(arg1,prd)) x2 = atan2(aimag(arg2),real(arg2,prd)) x2 = mod(abs(x1-x2),2.0_prd*pid) errz = abs(cmplx(log(abs(arg1/arg2)),min(x2,2.0_prd*pid-x2),prd)) end if end function errz ! ! real single real(prs) function psis(arg) real(prs) arg, at at = abs(arg) if (at > 1.0_prs) then psis = sign(1.0_prs + log(at),arg) else psis = arg end if end function psis ! ! real double real(prd) function psid(arg) real(prd) arg, at at = abs(arg) if (at > 1.0_prd) then psid = sign(1.0_prd + log(at),arg) else psid = arg end if end function psid end module error ! ! ! program to test installation ! program test use airy_functions_real_single use airy_functions_real_double use airy_functions_complex_single use airy_functions_complex_double use error implicit none integer, parameter :: prs = kind(0.0e0) integer, parameter :: prd = kind(0.0d0) real(prd), parameter :: oned = 1.0_prd, thvd = 1.50_prd real(prs), parameter :: ones = 1.0_prs, thvs = 1.50_prs real(prs) aias, daias, aibs, daibs, rms, pis, g1, tols, ffacs, hs, rs real(prd) x, aix, daix, aiass, daiass, rmd, pid, f1, told, ffacd, hd, rd real (prd), dimension(5) :: aizv, daizv, aiassv, daiassv complex(prs) zs, ciunits, aiasz, daiasz, aibsz, daibsz complex(prd) z, aiz, daiz, aiassz, daiassz, ciunitd complex(prd), dimension(5) :: bizv, dbizv, biassv, dbiassv real(prd), dimension(20) :: ray complex(prd), dimension(20) :: airay, dairay character(3) ai, dai, bi, dbi integer ierr, itol, itol2 integer ierrcnt, i, iptcnt, ilocerrcnt, ilocptcnt, j, nt, na, np, nam, nap open(19,file="airyerr.log",status='unknown') pid = 3.14159265358979323846264338327950288419716939937510582097494459230782_prd pis = real(pid,prs) ciunits = cmplx(0.0_prs,1.0_prs,prs) ciunitd = cmplx(0.0_prd,1.0_prd,prd) ! ! factor for testing computed values agains stored values or for computing ! values at the same point in different precisions itol = 100 ! ! factor for extremely close, but different, argument values itol2 = 100 ! ! machine epsilons tols = exp((1-digits(tols))*log(real(radix(tols),prs))) told = exp((1-digits(told))*log(real(radix(told),prd))) ! ! set the counters and declare a few variables ierrcnt = 0 iptcnt = 0 ai = ' Ai' dai = 'dAi' bi = ' Bi' dbi = 'dBi' print*, ' ---------------------------------------------------' print*, ' ' print*, 'This is a program which tests the installation ' print*, ' the Airy function pacakge by B.R. Fabijonas.' print*, ' Please send all bug reports to bfabi@smu.edu' print*, ' ' ! ! write some preliminaries to the output file write(19,*) ' This is a log file for the installation of the Airy function' write(19,*) ' package by B.R. Fabijonas. Please send all bug reports' write(19,*) ' to bfabi@smu.edu. The reference for this package is ' write(19,*) ' B.R. Fabijonas 2004. Algorithm XXX: Airy functions. ' write(19,*) ' ACM Trans. Math. Softw., Vol.xxx, pp.xxx-xxx.' write(19,*) ' A detailed description of the computation method and parameter' write(19,*) ' determination can be found in B.R. Fabijonas, D.W. Lozier,' write(19,*) ' and F.W.J. Olver, Computation of complex Airy functions and ' write(19,*) ' their zeros using asymptotics and the differential equation. ' write(19,*) ' ACM Trans. Math. Softw., Vol.xxx, pp. xxx-xxx.' write(19,*) ' ' write(19,'("***Parameter values used in the code which are compiler specific.")') write(19,'(" We use the abbreviations no. for number of terms and ")') write(19,'(" a.e. for asymptotic expansion.")') write(19,'(" Single precision computations: standard Airy functions")') call airy_info(rms, hs, nt, na, np) write(19,'(A45,e24.16)') ' machine epsilon', epsilon(rms) write(19,'(A45,e24.16)') ' radius of the delineating circle rho', rms write(19,'(A45,e24.16)') ' integration step size h', hs write(19,'(A60,i4)') ' no. in the Taylor series N',nt write(19,'(A60,i4)') ' no. in the a.e.s S', na write(19,'(A60,i4/)') ' number of partitions of the integration ray P', np write(19,'(" Single precision computations: auxiliary functions")') call airy_aux_info ( rs, hs, nt, nam, nap, np) write(19,'(A45,e24.16)') ' radius of the delineating circle rho_m ',rs write(19,'(A45,e24.16)') ' integrationstep size h_m ',hs write(19,'(A60,i4)') ' no. of terms in the Taylor series N_m ',nt write(19,'(A60,i4)') ' no. of terms in the a.e.s of the modulus functions S_m ',nam write(19,'(A60,i4)') ' no. of terms in the a.e.s of the phase functions S_theta ',nap write(19,'(A60,i4/)') ' no. of partitions of the integration ray P_m ',np write(19,'(" Double precision computations: standard Airy functions")') call airy_info(rmd, hd, nt, na, np) write(19,'(A45,e24.16)') ' machine epsilon ', epsilon(rmd) write(19,'(A45,e24.16)') ' the parameter rho ', rmd write(19,'(A45,e24.16)') ' step size h ', hd write(19,'(A60,i4)') ' no. in the Taylor series N ',nt write(19,'(A60,i4)') ' no. used in the a.e.s S ', na write(19,'(A60,i4/)') ' number of partitions of the integration ray P ', np write(19,'(" Double precision computations: auxiliary functions")') call airy_aux_info ( rd, hd, nt, nam, nap, np) write(19,'(A45,e24.16)') ' the parameter rho_m ',rd write(19,'(A45,e24.16)') ' step size h_m ',hd write(19,'(A60,i4)') ' no. in the t.s. N_m ',nt write(19,'(A60,i4)') ' no. in the a.e. of the modulus functions S_m ',nam write(19,'(A60,i4)') ' no. in the a.e. of the phase functions S_theta ',nap write(19,'(A60,i4//)') ' number of partitions of the integration ray P_m ',np ! !***************** ! testing real routines !**************** ilocerrcnt = 0 ilocptcnt = 0 write(6,'(" Comparing the function values returned by the")') write(6,'(" real routines against stored values")',advance="NO") write(19,'( "***Comparing the function values returned by the")') write(19,'(" real routines against stored values")') x = 5.82340_prd ! integration routine ffacd = oned !10.0_prd**max(oned,thvd*log(abs(x))) call airy_ai(x, aix, daix, ierr) call airy_ai(real(x,prs),aias,daias,ierr) write(19,1000) ai, x write(19,1001) aias, aix, 0.1539272803479132e-04_prd f1 = errg(aix,0.1539272803479132e-04_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 f1 = errg(daix,-0.3777943926223622e-04_prd) write(19,1000) dai, x write(19,1001) daias, daix, -0.3777943926223622e-04_prd if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 call airy_bi(x, aix, daix, ierr) call airy_bi(real(x,prs),aias,daias,ierr) write(19,1000) bi, x write(19,1001) aias, aix, 0.4288114816683048e+04_prd f1 = errg(aix,0.4288114816683048e+04_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1000) dbi, x write(19,1001) daias, daix, 0.1015462058214280e+05_prd f1 = errg(daix,0.1015462058214280e+05_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 x = -17.345890_prd ! asymptotic expansion routine ffacd = oned !10.0_prd**max(oned,thvd*log(abs(x))) call airy_ai(x, aix, daix, ierr) call airy_ai(real(x,prs),aias,daias,ierr) write(19,1000) ai, x write(19,1001) aias, aix, -0.2677772589719993e+00_prd f1 = errg(aix,-0.2677772589719993e+00_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1000) dai, x write(19,1001) daias, daix, -0.2900295189410030e+00_prd f1 = errg(daix,-0.2900295189410030e+00_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 call airy_bi(x, aix, daix, ierr) call airy_bi(real(x,prs),aias,daias,ierr) write(19,1000) bi, x write(19,1001) aias, aix, 0.6870907209572491e-01_prd f1 = errg(aix,0.6870907209572491e-01_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1000) dbi, x write(19,1001) daias, daix, -0.1114292633371775e+01_prd f1 = errg(daix,-0.1114292633371775e+01_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...done with ",i3," errors")') ilocerrcnt ! !***************** ! testing complex routines !**************** ilocerrcnt = 0 ilocptcnt = 0 print*, ' ' write(6,'(" Comparing the function values returned by the")') write(6,'(" complex routines against stored values")',advance="NO") write(19,'( "***Comparing the function values returned by the complex routines")') write(19,'(" against stored values")') z = cmplx(1.5662_prd,3.681_prd,prd) ! integration routine call airy_ai(z, aiz, daiz, ierr) call airy_ai(cmplx(1.5662,3.681),aiasz,daiasz,ierr) ffacd = oned !10.0_prd**max(oned,thvd*log(abs(z))) write(19,2000) ai, real(z,prd),aimag(z) write(19,2001) real(aiasz), aimag(aiasz), real(aiz,prd), aimag(aiz), & 0.3807338825307106e+00_prd,0.3605171603516821e+00_prd f1 = errg(aiz,& cmplx(0.3807338825307106e+00_prd,0.3605171603516821e+00_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,2000) dai, real(z,prd),aimag(z) write(19,2001) real(daiasz), aimag(daiasz), real(daiz,prd), aimag(daiz), & -0.2685171251031442e+00_prd,-0.1010759493443000e+01_prd f1 = errg(daiz,& cmplx(-0.2685171251031442e+00_prd,-0.1010759493443000e+01_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 2 z = cmplx(-9.8342_prd,4.2811_prd,prd) ! asymptotic expansion routine call airy_ai(z, aiz, daiz, ierr) call airy_ai(cmplx(-9.8342,4.2811),aiasz,daiasz,ierr) write(19,2000) ai, real(z,prd),aimag(z) write(19,2001) real(aiasz), aimag(aiasz), real(aiz,prd), aimag(aiz), & 0.1069401631365655e+06_prd,-0.4772398758353830e+05_prd f1 = errg(aiz,& cmplx(0.1069401631365655e+06_prd,-0.4772398758353830e+05_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,2000) dai, real(z,prd),aimag(z) write(19,2001) real(daiasz), aimag(daiasz), real(daiz,prd), aimag(daiz), & -0.2216534220094920e+06_prd,-0.3110801974925891e+06_prd f1 = errg(daiz,& cmplx(-0.2216534220094920e+06_prd,-0.3110801974925891e+06_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 2 write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...done with ",i3," errors")') ilocerrcnt ! !***************** ! testing interface of the routines !**************** ! a perturbation analysis on the asymptotic expansion shows that ! |Ai(x(1+eps)) - Ai(x)| <= |Ai(x)|*(.25+x**(3/2))*eps. the ! same holds for the derivative. ! we will use this tolerance print*, ' ' write(6,'(" Comparing the function values on both side of the delineating circle")') write(6,'(" and verifying continuity within a tolerance")',advance="NO") write(19,'("***Testing function values on both side of the delineating circle")') write(19,'(" and verifying continuity within a tolerance")') write(19,'(" (no output printed unless an error is found)")') ilocerrcnt = 0 ilocptcnt = 0 ! for real variables call airy_info(rms) call airy_info(rmd) call airy_ai(rms+epsilon(rms)*2,aias,daias,ierr) call airy_ai(rms-epsilon(rms)*2,aibs,daibs,ierr) ffacs = abs(aias)*(.250_prs + rms*sqrt(rms))*4*epsilon(rms) g1 = errg(aias,aibs) if ( g1 > itol*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,3000) ai, rms write(19,3001) aias, aibs end if ffacs = abs(daias)*(.250_prs + rms*sqrt(rms))*4*epsilon(rms) g1 = errg(daias,daibs) if ( g1 > itol*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,3000) dai, rms write(19,3001) daias, daibs end if call airy_ai(rmd+epsilon(rmd)*2,aix,daix,ierr) call airy_ai(rmd-epsilon(rmd)*2,aiass,daiass,ierr) ffacd = abs(aix)*(.250_prd + rmd*sqrt(rmd))*4*epsilon(rmd) f1 = errg(aix,aiass) if ( f1 > itol*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,3000) ai, rmd write(19,3001) aix, aiass end if f1 = errg(daix,daiass) ffacd = abs(daix)*(.250_prd + rmd*sqrt(rmd))*4*epsilon(rmd) if ( f1 > itol*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,3000) dai, rmd write(19,3001) daix, daiass end if ilocptcnt = ilocptcnt + 4 ! for complex variables do i = -99,99 zs = (rms + epsilon(rms)*2)*exp(ciunits*i*0.01_prs*pis) call airy_ai(zs,aiasz,daiasz,ierr) zs = (rms - epsilon(rms)*2)*exp(ciunits*i*0.01_prs*pis) call airy_ai(zs,aibsz,daibsz,ierr) ffacs = abs(aiasz)*(0.250_prs+rms*sqrt(rms))*4*epsilon(rms) g1 = errg(aiasz,aibsz) if ( g1 > itol*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,4000) ai, rms*exp(ciunits*i*0.01_prs*pis) write(19,4001) aiasz,aibsz write(19,*) g1, itol*ffacs end if g1 = errg(daiasz,daibsz) ffacs = abs(daiasz)*(0.250_prs+rms*sqrt(rms))*4*epsilon(rms) if ( g1 > itol*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,*) g1, ffacs write(19,4000) dai, rms*exp(ciunits*i*0.01_prs*pis) write(19,4001) daiasz,daibsz end if z = (rmd + epsilon(rmd)*2)*exp(ciunitd*i*0.01_prd*pid) call airy_ai(z,aiz,daiz,ierr) z = (rmd - epsilon(rmd)*2)*exp(ciunitd*i*0.01_prd*pid) call airy_ai(z,aiassz,daiassz,ierr) f1 = errg(aiz,aiassz) ffacd = abs(aiz)*(0.250_prd+rmd*sqrt(rmd))*4*epsilon(rmd) if ( f1 > itol*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,4000) ai, rmd*exp(ciunitd*i*0.01_prd*pid) write(19,4002) aiasz,aiassz, aiassz end if f1 = errg(daiz,daiassz) ffacd = abs(daiz)*(0.250_prd+rmd*sqrt(rmd))*4*epsilon(rmd) if ( f1 > itol*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,4000) dai, rmd*exp(ciunitd*i*0.01_prd*pid) write(19,4002) daiz,daiassz end if ilocptcnt = ilocptcnt + 4 end do write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...done with ",i3," errors")') ilocerrcnt ! !***************** ! testing different precisions !**************** write(6,'(" ")') write(6,'(" Comparing the computation in different precisions at random points")') write(6,'(" on the real line and in the complex plane and verifying the")') write(6,'(" accuracy of the lower precision function values.")',advance='no') write(19,'("***Comparing the computation in different precisions at random points ")') write(19,'(" on the real line and in the complex plane and verifying the accuracy of")') write(19,'(" the lower precision function values (no output printed unless an error is found).")') ilocerrcnt = 0 ilocptcnt = 0 ! for real variables do i = 1,200 call random_number(rmd) rms = (rmd-0.50_prd)*20.0_prd rmd = rms call airy_ai(rmd, aix, daix, ierr) call airy_ai(rms, aias, daias, ierr) f1 = errg(aix,aias) ffacs = ones !10.0_prs**max(ones,thvs*log(abs(rms))) if ( f1 > itol2*tols*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,5000) ai, rmd write(19,5001) aias, aix end if f1 = errg(daix,daias) if ( f1 > itol2*tols*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,5000) dai, rmd write(19,5001) daias, daix end if ilocptcnt = ilocptcnt + 2 end do ! for complex variables do i = 1,200 call random_number(rmd) call random_number(x) zs = rmd*10*exp(ciunitd*(x-0.50_prd)*2.0_prd*pid) z = zs call airy_ai(zs,aiasz,daiasz,ierr) call airy_ai(z,aiz,daiz,ierr) ffacs = ones !10.0_prs**max(ones,thvs*log(abs(zs))) f1 = errg(aiz,aiasz) if ( f1 > itol*tols*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,6000) ai, z write(19,6001) aiasz,aiz end if f1 = errg(daiz,daiasz) if ( f1 > itol*tols*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,6000) dai, z write(19,6001) daiasz,daiz end if ilocptcnt = ilocptcnt + 2 end do write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...done with ",i3," errors")') ilocerrcnt ! !***************** ! testing vectorized airy_ai routine !**************** print*, ' ' write(6,'(" Testing the vectorized routine of the Airy functions")',advance='no') write(19,'("***Testing the vectorized routine of the Airy functions (no output printed")') write(19,'(" unless an error is found).")') ilocerrcnt = 0 ilocptcnt = 0 do j = -20,20 do i = 1,20 ray(i) = 0.5_prd*i ilocptcnt = ilocptcnt + 1 end do call airy_ai(ray,0.05_prd*j*pid,airay,dairay,ierr,.true.) do i = 1,20 call airy_ai(ray(i)*exp(ciunitd*j*0.05_prd*pid), aiz,daiz, ierr) g1 = errg (aiz,airay(i)) ffacd = oned !10.0_prd**max(oned,thvd*log(ray(i))) if ( g1 > itol2*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,2002) ai, ray(i), j*0.05_prd*pid write(19,2003) aiz, airay(i) end if g1 = errg (daiz,dairay(i)) if ( g1 > itol2*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,2002) dai, ray(i), j*0.05_prd*pid write(19,2003) daiz, dairay(i) end if end do end do write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...done with ",i3," errors")') ilocerrcnt ! !***************** ! testing zero routines !**************** print*, ' ' print*, 'Testing the routines for the zeros of the Airy functions' write(19,'("***Testing the routines for the zeros of the Airy functions")') write(19,'(" real and complex zeros--scalar version")') ilocerrcnt = 0 ilocptcnt = 0 write(6,'(" real zeros by table lookup")',advance='no') write(19,'(" **Zeros by table lookup")') call airy_ai_zero (18, aias, ierr, ai_assoc=daias, dai_zero=aibs, & dai_assoc=daibs) call airy_ai_zero (18, aix, ierr, ai_assoc=aiass, dai_zero=daix, & dai_assoc=daiass) ffacd = 10.0_prd**max(oned,thvd*log(abs(18.0_prd))) write(19,1002) 18, ai write(19,1001) aias, aix, -19.126380474246954_prd f1 = errg(aix, -19.126380474246954_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) ai, 18 write(19,1001) daias, aiass,-1.179880729870146_prd f1 = errg(aiass,-1.179880729870146_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1002) 18, dai write(19,1001) aibs, daix, -18.764798437665956_prd f1 = errg(daix,-18.764798437665956_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) dai, 18 write(19,1001) daibs, daiass,-0.2710702785769711_prd f1 = errg(daiass,-0.2710702785769711_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 call airy_bi_zero (18, aias, ierr, bi_assoc=daias, dbi_zero=aibs, & dbi_assoc=daibs) call airy_bi_zero (18, aix, ierr, bi_assoc=aiass, dbi_zero=daix, & dbi_assoc=daiass) write(19,1002) 18, bi write(19,1001) aias, aix, -18.765508284480081_prd f1 = errg(aix,-18.765508284480081_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) bi, 18 write(19,1001) daias, aiass, -1.174276253118059_prd f1 = errg(aiass,-1.174276253118059_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1002) 18, dbi write(19,1001) aibs, daix, -19.125697156412638_prd f1 = errg(daix,-19.125697156412638_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 18 write(19,1001) daibs, daiass,0.2697826139865426_prd f1 = errg(daiass,0.2697826139865426_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 call airy_bi_zero (18, aiasz, ierr, bi_assoc=daiasz, dbi_zero=aibsz, & dbi_assoc=daibsz) call airy_bi_zero (18, aiz, ierr, bi_assoc=aiassz, dbi_zero=daiz, & dbi_assoc=daiassz) write(19,1004) 18, bi write(19,2001) aiasz, aiz, 9.494603679187176_prd,16.603624606898762_prd f1 = errg(aiz,cmplx(9.494603679187176_prd,16.603624606898762_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1005) bi, 18 write(19,2001) daiasz, aiassz, 1.445920793639133_prd,-0.8328073286646578_prd f1 = errg(aiassz,cmplx(1.445920793639133_prd,-0.8328073286646578_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1004) 18, dbi write(19,2001) aibsz, daiz, 9.313152405593771_prd,16.290870307133684_prd f1 = errg(daiz,cmplx(9.313152405593771_prd,16.290870307133684_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1005) dbi, 18 write(19,2001) daibsz, daiassz, -0.3321948851433578_prd,-0.1913210621443901_prd f1 = errg(daiassz,& cmplx(-0.3321948851433578_prd,-0.1913210621443901_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...................done with ",i3," errors")') ilocerrcnt ! !***************** ! testing zero routines !**************** write(6,'(" real zeros by asymptotic expansion")',advance='no') write(19,'(" **Zeros by asymptotic expansion")') call airy_ai_zero (86, aias, ierr, ai_assoc=daias, dai_zero=aibs, & dai_assoc=daibs) call airy_ai_zero (86, aix, ierr, ai_assoc=aiass, dai_zero=daix, & dai_assoc=daiass) ffacd = oned !10.0_prd**max(oned,thvd*log(abs(86.0_prd))) write(19,1002) 86, ai write(19,1001) aias, aix, -54.657586491868699_prd f1 = errg(aix,-54.657586491868699_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) ai, 86 write(19,1001) daias, aiass,-1.534044239238300_prd f1 = errg(aiass,-1.534044239238300_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1002) 86, dai write(19,1001) aibs, daix, -54.444826792609810_prd f1 = errg(daix,-54.444826792609810_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) dai, 86 write(19,1001) daibs, daiass,-0.2076995780278862_prd f1 = errg(daiass,-0.2076995780278862_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 call airy_bi_zero (86, aias, ierr, bi_assoc=daias, dbi_zero=aibs, & dbi_assoc=daibs) call airy_bi_zero (86, aix, ierr, bi_assoc=aiass, dbi_zero=daix, & dbi_assoc=daiass) write(19,1002) 86, bi write(19,1001) aias, aix, -54.444911130586874_prd f1 = errg(aix,-54.444911130586874_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) bi, 86 write(19,1001) daias, aiass,-1.532549805062612_prd f1 = errg(aiass,-1.532549805062612_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1002) 86, dbi write(19,1001) aibs, daix,-54.657502808936158_prd f1 = errg(daix,-54.657502808936158_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 86 write(19,1001) daibs, daiass, 0.2074972409267743_prd f1 = errg(daiass, 0.2074972409267743_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 call airy_bi_zero (86, aiasz, ierr, bi_assoc=daiasz, dbi_zero=aibsz, & dbi_assoc=daibsz) call airy_bi_zero (86, aiz, ierr, bi_assoc=aiassz, dbi_zero=daiz, & dbi_assoc=daiassz) write(19,1004) 86, bi write(19,2001) aiasz, aiz, 27.288200667644368_prd,47.358306153862507_prd f1 = errg(aiz,cmplx(27.288200667644368_prd,47.358306153862507_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1005) bi, 86 write(19,2001) daiasz, aiassz, 1.879045614304762_prd,-1.084330361798879_prd f1 = errg(aiassz,cmplx(1.879045614304762_prd,-1.084330361798879_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1004) 86, dbi write(19,2001) aibsz, daiz, 27.181741517075782_prd,47.174096724925626_prd f1 = errg(daiz,cmplx(27.181741517075782_prd,47.174096724925626_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1005) dbi, 86 write(19,2001) daibsz, daiassz, -0.2544106266601735_prd,-0.1468108932911459_prd f1 = errg(daiassz,cmplx(-0.2544106266601735_prd,-0.1468108932911459_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...........done with ",i3," errors")') ilocerrcnt ! !***************** ! testing vectorized zero routines !**************** write(6,'(" real and complex zeros--vectorized version")',advance="no") write(19,'(" **Zeros--vectorized version (no output printed unless an error is found)")') ilocerrcnt = 0 ilocptcnt = 0 call airy_ai_zero (23, aizv, ierr, ai_assoc=aiassv) ffacd = oned !10.0_prd**max(oned,thvd*log(abs(27.0_prd))) f1 = errg(aizv(1),-22.567612917496504_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 23, ai write(19,1008) aizv(1), -22.567612917496504_prd end if f1 = errg(aizv(2),-23.224165001121680_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 24, ai write(19,1008) aizv(2), -23.224165001121680_prd end if f1 = errg(aizv(3),-23.871564455535918_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 25, ai write(19,1008) aizv(3), -23.871564455535918_prd end if f1 = errg(aizv(4),-24.510301236589672_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 26, ai write(19,1008) aizv(4), -24.510301236589672_prd end if f1 = errg(aizv(5),-25.140821166148957_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 27, ai write(19,1008) aizv(5), -25.140821166148957_prd end if f1 = errg(aiassv(1),1.229700701509681_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) ai, 23 write(19,1008) aiassv(1),1.229700701509681_prd end if f1 = errg(aiassv(2),-1.238547875329632_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) ai, 24 write(19,1008) aiassv(2),-1.238547875329632_prd end if f1 = errg(aiassv(3),1.247089945259408_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) ai, 25 write(19,1008) aiassv(3),1.247089945259408_prd end if f1 = errg(aiassv(4),-1.255349140475735_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) ai, 26 write(19,1008) aiassv(4),-1.255349140475735_prd end if f1 = errg(aiassv(5),1.263345282750799_prd) if ( f1 > itol*told*ffacd )then ilocerrcnt = ilocerrcnt+1 write(19,1003) ai, 27 write(19,1008) aiassv(5),1.263345282750799_prd end if ilocptcnt = ilocptcnt + 10 call airy_ai_zero (23, aizv, ierr, dai_zero=daizv, dai_assoc=daiassv) f1 = errg(daizv(1),-22.235232285348914_prd) if ( f1 > itol*told*ffacd )then ilocerrcnt = ilocerrcnt+1 write(19,1002) 23, dai write(19,1008) daizv(1),-22.235232285348914_prd end if f1 = errg(daizv(2),-22.896588738874620_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 24, dai write(19,1008) daizv(2),-22.896588738874620_prd end if f1 = errg(daizv(3),-23.548526295928802_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 25, dai write(19,1008) daizv(3),-23.548526295928802_prd end if f1 = errg(daizv(4),-24.191559709526349_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 26, dai write(19,1008) daizv(4),-24.191559709526349_prd end if f1 = errg(daizv(5),-24.826156425921152_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 27, dai write(19,1008) daizv(5),-24.826156425921152_prd end if f1 = errg(daiassv(1),0.259812670151466_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dai, 23 write(19,1008) daiassv(1),0.259812670151466_prd end if f1 = errg(daiassv(2),-0.257916075332572_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dai, 24 write(19,1008) daiassv(2),-0.257916075332572_prd end if f1 = errg(daiassv(3),0.256112333779654_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dai, 25 write(19,1008) daiassv(3),0.256112333779654_prd end if f1 = errg(daiassv(4),-0.254393342646825_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dai, 26 write(19,1008) daiassv(4),-0.254393342646825_prd end if f1 = errg(daiassv(5),0.252751992576574_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dai, 27 write(19,1008) daiassv(5),0.252751992576574_prd end if ilocptcnt = ilocptcnt + 10 ! !***************** ! testing vectorized zero routines !**************** call airy_bi_zero (23, aizv, ierr, bi_assoc=aiassv) f1 = errg(aizv(1),-22.235737881803384_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 23, bi write(19,1008) aizv(1),-22.235737881803384_prd end if f1 = errg(aizv(2),-22.897065554219793_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 24, bi write(19,1008) aizv(2),-22.897065554219793_prd end if f1 = errg(aizv(3),-23.548977079642448_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 25, bi write(19,1008) aizv(3),-23.548977079642448_prd end if f1 = errg(aizv(4),-24.191986850648995_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 26, bi write(19,1008) aizv(4),-24.191986850648995_prd end if f1 = errg(aizv(5),-24.826562012152888_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 27, bi write(19,1008) aizv(5),-24.826562012152888_prd end if f1 = errg(aiassv(1),1.225154995846960_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) bi, 23 write(19,1008) aiassv(1),1.225154995846960_prd end if f1 = errg(aiassv(2),-1.234163920487302_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) bi, 24 write(19,1008) aiassv(2),-1.234163920487302_prd end if f1 = errg(aiassv(3),1.242855598092317_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) bi, 25 write(19,1008) aiassv(3),1.242855598092317_prd end if f1 = errg(aiassv(4),-1.251253611221631_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) bi, 26 write(19,1008) aiassv(4),-1.251253611221631_prd end if f1 = errg(aiassv(5),1.259378938688538_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) bi, 27 write(19,1008) aiassv(5),1.259378938688538_prd end if ilocptcnt = ilocptcnt + 10 call airy_bi_zero (23, aizv, ierr, dbi_zero=daizv, dbi_assoc=daiassv) f1 = errg(daizv(1),-22.567122080497199_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 23, dbi write(19,1008) daizv(1),-22.567122080497199_prd end if f1 = errg(daizv(2),-23.223701521208962_prd) if ( f1 > itol*told*ffacd )then ilocerrcnt = ilocerrcnt+1 write(19,1002) 24, dbi write(19,1008) daizv(2),-23.223701521208962_prd end if f1 = errg(daizv(3),-23.871125771677974_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 25, dbi write(19,1008) daizv(3),-23.871125771677974_prd end if f1 = errg(daizv(4),-24.509885117016239_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 26, dbi write(19,1008) daizv(4),-24.509885117016239_prd end if f1 = errg(daizv(5),-25.140425655367874_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 27, dbi write(19,1008) daizv(5),-25.140425655367874_prd end if f1 = errg(daiassv(1),-0.258852215916922_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 23 write(19,1008) daiassv(1),-0.258852215916922_prd end if f1 = errg(daiassv(2),0.257003129647316_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 24 write(19,1008) daiassv(2),0.257003129647316_prd end if f1 = errg(daiassv(3),-0.255242710064815_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 25 write(19,1008) daiassv(3),-0.255242710064815_prd end if f1 = errg(daiassv(4),0.253563372437697_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 26 write(19,1008) daiassv(4),0.253563372437697_prd end if f1 = errg(daiassv(5),-0.251958444330784_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 27 write(19,1008) daiassv(5),-0.251958444330784_prd end if ilocptcnt = ilocptcnt + 10 ! !***************** ! testing vectorized zero routines !**************** call airy_bi_zero (23, bizv, ierr, bi_assoc=biassv) f1 = errg(bizv(1),& cmplx(11.220656371214879_prd,19.580653883038824_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1004) 23, bi write(19,2003) bizv(1),11.220656371214879_prd,19.580653883038824_prd end if f1 = errg(bizv(2),& cmplx(11.549830143889196_prd,20.148722567309367_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1004) 24, bi write(19,2003) bizv(2), 11.549830143889196_prd,20.148722567309367_prd end if f1 = errg(bizv(3),& cmplx(11.874378641347221_prd,20.708893464731311_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1004) 25, bi write(19,2003) bizv(3),11.874378641347221_prd,20.708893464731311_prd end if f1 = errg(bizv(4),& cmplx(12.194551330833875_prd,21.261588251953778_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1004) 26, bi write(19,2003) bizv(4),12.194551330833875_prd,21.261588251953778_prd end if f1 = errg(bizv(5),& cmplx(12.510575046777394_prd,21.807190718498749_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1004) 27, bi write(19,2003) bizv(5),12.510575046777394_prd,21.807190718498749_prd end if f1 = errg(biassv(1),& cmplx(-1.506774749698633_prd,0.868314074899650_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1005) bi, 23 write(19,2003) biassv(1),-1.506774749698633_prd,0.868314074899650_prd end if f1 = errg(biassv(2),& cmplx(1.517585357164310_prd,-0.874612709515018_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1005) bi, 24 write(19,2003) biassv(2),1.517585357164310_prd,-0.874612709515018_prd end if f1 = errg(biassv(3),& cmplx(-1.528024149210944_prd,0.880692431621191_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1005) bi, 25 write(19,2003) biassv(3),-1.528024149210944_prd,0.880692431621191_prd end if f1 = errg(biassv(4),& cmplx(1.538118145012280_prd,-0.886569309158123_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1005) bi, 26 write(19,2003) biassv(4),1.538118145012280_prd,-0.886569309158123_prd end if f1 = errg(biassv(5),& cmplx(-1.547891444975066_prd,0.892257657608977_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1005) bi, 27 write(19,2003) biassv(5),-1.547891444975066_prd,0.892257657608977_prd end if ilocptcnt = ilocptcnt + 10 call airy_bi_zero (23, bizv, ierr, dbi_zero=dbizv, dbi_assoc=dbiassv) f1 = errg(dbizv(1),& cmplx(11.053994360667899_prd,19.293078213959429_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 23, dbi write(19,2003) dbizv(1),11.053994360667899_prd,19.293078213959429_prd end if f1 = errg(dbizv(2),& cmplx(11.385596969302076_prd,19.865292017137630_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 24, dbi write(19,2003) dbizv(2),11.385596969302076_prd,19.865292017137630_prd end if f1 = errg(dbizv(3),& cmplx(11.712438641126221_prd,20.429378925421950_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 25, dbi write(19,2003) dbizv(3),11.712438641126221_prd,20.429378925421950_prd end if f1 = errg(dbizv(4),& cmplx(12.034781571133943_prd,20.985781895965161_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 26, dbi write(19,2003) dbizv(4),12.034781571133943_prd,20.985781895965161_prd end if f1 = errg(dbizv(5),& cmplx(12.352863681490561_prd,21.534903282918666_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 27, dbi write(19,2003) dbizv(5),12.352863681490561_prd,21.534903282918666_prd end if f1 = errg(dbiassv(1),& cmplx(0.318355274563740_prd,0.183451937317410_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 23 write(19,2003) dbiassv(1),0.318355274563740_prd,0.183451937317410_prd end if f1 = errg(dbiassv(2),& cmplx(-0.316024910474316_prd,-0.182124025592498_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 24 write(19,2003) dbiassv(2),-0.316024910474316_prd,-0.182124025592498_prd end if f1 = errg(dbiassv(3),& cmplx(0.313808934629479_prd,0.180860595992543_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 25 write(19,2003) dbiassv(3),0.313808934629479_prd,0.180860595992543_prd end if f1 = errg(dbiassv(4),& cmplx(-0.311697340122842_prd,-0.179656065952459_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 26 write(19,2003) dbiassv(4),-0.311697340122842_prd,-0.179656065952459_prd end if f1 = errg(dbiassv(5),& cmplx(0.309681349894832_prd,0.178505532129885_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 27 write(19,2003) dbiassv(5),0.309681349894832_prd,0.178505532129885_prd end if ilocptcnt = ilocptcnt + 10 write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...done with ",i3," errors")') ilocerrcnt ! !***************** ! testing modulus and phase routines !**************** print*, ' ' write(6,'(" Comparing the computed modulus, phase, and associated values")') write(6,'(" against stored values")',advance='no') write(19,'("***Comparing the computed modulus, phase, and associated values against stored values")') ilocerrcnt = 0 ilocptcnt = 0 x = -5.234_prd write(19,'(" **by ODE integration")') call airy_aux ( real(x,prs), aias, daias, ierr, aibs, daibs) call airy_aux ( x, aix, daix, ierr, aiass, daiass) ffacd = oned !10.0_prd**max(oned,thvd*log(abs(x))) write(19,1006) 'M(x)', x write(19,1001) aias, aix,0.3728081698402362_prd f1 = errg(aix,0.3728081698402362_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1007) 'theta(x)', x write(19,1001) daias, daix,8.759640547244738_prd f1 = errg(daix,8.759640547244738_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1006) 'M(x)', x write(19,1001) aibs, aiass,0.8540001773409781_prd f1 = errg(aiass,0.8540001773409781_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1007) 'phi(x)', x write(19,1001) daibs, daiass,7.209566736849307_prd f1 = errg(daiass,7.209566736849307_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 ! !***************** ! testing modulus and phase routines !**************** x = -25.392_prd write(19,'(" **by asymptotic expansion")') call airy_aux ( real(x,prs), aias, daias, ierr, aibs, daibs) call airy_aux ( x, aix, daix, ierr, aiass, daiass) ffacd = oned !10.0_prd**max(oned,thvd*log(abs(x))) write(19,1006) 'M(x)', x write(19,1001) aias, aix,0.2513325654640693_prd f1 = errg(aix,0.2513325654640693_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1007) 'theta(x)', x write(19,1001) daias, daix,86.085580681739742_prd f1 = errg(daix,86.085580681739742_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1006) 'N(x)', x write(19,1001) aibs, aiass,1.2664912447881540_prd f1 = errg(aiass,1.2664912447881540_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1007) 'phi(x)', x write(19,1001) daibs, daiass,84.516738087389200_prd f1 = errg(daiass,84.516738087389200_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 !***************** ! testing modulus and phase in different precisions !**************** write(19,'(" **in different precisions (no output printed unless an error is found)")') do i = 1,200 call random_number(rmd) rms = -rmd*10.0_prd rmd = rms call airy_aux ( rmd, aix, daix, ierr, aiass, daiass) call airy_aux ( rms, aias, daias, ierr, aibs, daibs) ffacd = ones !10.0_prs**max(ones,thvs*log(abs(rms))) f1 = errg(aix,aias) if ( f1 > itol*tols*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,5002) ai, rmd write(19,5001) aias,aix end if f1 = errg(daix,daias) if ( f1 > itol*tols*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,5002) ai, rmd write(19,5001) daias,daix end if f1 = errg(aiass,aibs) if ( f1 > itol*tols*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,5003) ai, rmd write(19,5001) aibs,aiass end if f1 = errg(daiass,daibs) if ( f1 > itol*tols*ffacs ) then ilocerrcnt = ilocerrcnt+11 write(19,5003) ai, rmd write(19,5001) daibs,daiass end if ilocptcnt = ilocptcnt + 4 end do write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...done with ",i3," errors")') ilocerrcnt write(19,'("Number of tests performed:",i6)') iptcnt write(19,'("Number of errors found:",i6)') ierrcnt print*, ' ' print*, 'Number of tests performed:',iptcnt print*, 'Number of errors found:',ierrcnt print*, ' ' print*, ' ' print*, 'End of run' print*, ' ' print*, ' ' close(21) 1500 format(' The number of points tested is ',i4) 1501 format(' The number of errors found is ',i4//) 1000 format(2x,'Computing',2x,A3,2x,'at x=',f9.5) 1001 format(17x,'We computed ',e15.8,' in single precision'/,17x,'We computed ', & e23.16,' in double precision',/,5x,'and the stored value is ',e23.16/) 1002 format(' Testing the ',i2,'th zero of ',A3) 1008 format(17x,'We computed ',e23.16,' in double precision',/,5x,'and the stored value is ',e23.16/) 1003 format(' Testing the associated function value of ',A3,' at the ',i2,' zero') 1004 format(' Testing the ',i2,'th complex zero of ',A3) 1005 format(' Testing the associated function value of ',A3,' at the ',i2,' complex zero') 1006 format(' Testing the modulus ',A4,' at x= ',f9.5) 1007 format(' Testing the phase ',A8,' at x= ',f9.5) 2000 format(2x,'Computing',2x,A3,2x,'at z=',f9.5,' +I*(',f9.5,')') 2001 format(2x,'In single precision, we computed ',e15.8,8x,' +I*(',e15.8,')',/,2x, & 'In double precision, we computed ',e23.16,' +I*(',e23.16,')',/,11x, & 'and the stored value is ',e23.16,' +I*(',e23.16,')',/) 2002 format('Error encountered testing'/' the vectorized version of ',2x,A3, & 2x,'at (r,theta)= (',f9.5,',',f9.5,')') 2003 format(2x,'In double precision, we computed ',e23.16,' +I*(',e23.16,')',/,11x, & 'and the stored value is ',e23.16,' +I*(',e23.16,')',/) 3000 format('Error encountered testing each side of ',2x,A3,2x,'at x=',f9.5) 3001 format(12x,'Below we got ',e23.16,/8x,'and above we got ',e23.16/) 4000 format('Error encountered testing each side of ',2x,A3,2x,'at z=',f9.5, & ' +I*(',f9.5,')') 4001 format(12x,'Below we got ',e15.8,' +I*(',e15.8,')',/& 8x,'and above we got ',e15.8,' +I*(',e15.8,')'/) 4002 format(12x,'Below we got ',e23.16,' +I*(',e23.16,')',/& 8x,'and above we got ',e23.16,' +I*(',e23.16,')'/) 5000 format('Error encountered testing',2x,A3,2x,'at x=',f9.5) 5001 format(9x,'Single returned ',e15.8,/5x,'and double returned ',e23.16/) 5002 format('Error encountered testing the modulus of ',2x,A3,2x,'at x=',f9.5) 5003 format('Error encountered testing the phase of ',2x,A3,2x,'at x=',f9.5) 6000 format('Error encountered testing',2x,A3,2x,'at z=',f9.5, & ' +I*(',f9.5,')') 6001 format(8x,'Single returned ',e15.8,9x,' +I*(',e15.8,')',/& 4x,'and double returned ',e23.16,' +I*(',e23.16,')'/) end program test SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'airy_complex' then echo shar: will not over-write existing file "'airy_complex'" else cat << "SHAR_EOF" > 'airy_complex' ! ! begin file airy_functions_complex ! version 1.0 ! edited 07/24/2002 ! !*** !************************************************************************ !*** subroutine airy_aic (z, ai, dai, ierr, argument_switch, modify_switch) complex(prd), intent(in) :: z complex(prd), intent(out) :: ai, dai integer, intent(out) :: ierr logical, intent(in), optional :: argument_switch, modify_switch !* ! The following subroutine is the driver routine that computes the ! Airy function Ai(z) and its derivative for complex arguments. ! The parameter PRD specifies the precision as determined by the ! module being used, e.g. single, double, quad. The subroutines ! called are: ! parameter_airy = determines the machine specific parameters ! flowsc = checks for over- and under- flow regions ! asympc = asymptotic expansion (direct evaluation), ! taylorc = taylor series expansion (integration). !* integer nn, iregion complex(prd) dG_z, dG_mz, G_z, G_mz, zeta_local, ai0s, ai1s complex(prd), dimension(2,2) :: tsm ! if (.not. is_init_airy) call parameter_airy iflag = 0 mod_local = .false. arg_local = .false. if (present(modify_switch)) mod_local = modify_switch if (present(argument_switch)) arg_local = argument_switch if (.not. arg_local) then r_global = abs(z) ! argument in rectangular form if ( r_global /= zero ) then theta_global = atan2(aimag(z),real(z,prd)) else theta_global = zero end if z_global = z else r_global = z ! argument in polar form theta_global = aimag(z) if (floor(abs(theta_global)/pi) /= 0 .or. r_global < zero) iflag = -10 z_global = r_global*exp(ciunit*theta_global) end if call flowsc if (iflag == 0 .or. iflag == 1) then zeta_global = two_third*z_global*sqrt(z_global) ! ! choose the region where z lies according to the paper iregion = 0 if (r_global >= r_min) then if (abs(theta_global) <= two_pi_third) then iregion = 1 else iregion = 2 end if else if (abs(theta_global) <= pi_thirds) then iregion = 3 if ( r_global <= r_min/n_parts ) iregion = 4 else iregion = 4 end if end if ! ! call the appropriate subroutine(s) to evaluate Ai(z) and/or Ai'(z) select case (iregion) case(1) call asympc (G_z, dG_z) ai0s = G_z ai1s = dG_z case(2) if (iflag == 1) then ! exp(2*zeta) < underflow limit call asympc (G_z, dG_z) ai0s = G_z ai1s = dG_z else call asympc (G_z, dG_z, G_mz, dG_mz) if (theta_global > zero) then ai0s = G_z + ciunit*G_mz*exp(two*zeta_global) ai1s = dG_z - ciunit*dG_mz*exp(two*zeta_global) else ai0s = G_z - ciunit*G_mz*exp(two*zeta_global) ai1s = dG_z + ciunit*dG_mz*exp(two*zeta_global) end if end if case(3) nn = ceiling(n_parts*(one-r_global/r_min)) call asympc (G_z, dG_z, r_in=r_min ) call taylorc (nn, r_min, r_global, tsm) zeta_local = two_third*r_min*sqrt(r_min) & *exp(ciunit*1.50_prd*theta_global) zeta_local = zeta_global - zeta_local ai0s = exp(zeta_local)*(tsm(1,1)*G_z + tsm(1,2)*dG_z) ai1s = exp(zeta_local)*(tsm(2,1)*G_z + tsm(2,2)*dG_z) case(4) nn = ceiling(n_parts*r_global/r_min) call taylorc (nn, zero, r_global, tsm) ai0s = tsm(1,1)*ai0zer + tsm(1,2)*ai1zer ai1s = tsm(2,1)*ai0zer + tsm(2,2)*ai1zer end select if ( iregion <= 3 .and. .not. mod_local ) then ai0s = exp(-zeta_global)*ai0s ai1s = exp(-zeta_global)*ai1s elseif ( iregion == 4 .and. mod_local) then ai0s = exp(zeta_global)*ai0s ai1s = exp(zeta_global)*ai1s end if end if ! ! error handling and return values select case (iflag) case(:-1) ai = czero dai = czero ierr = iflag case(0:1) ai = ai0s dai = ai1s ierr = 0 case(5) ai = cmplx(ai0zer,zero,prd) ! the value of Ai(0) dai = cmplx(ai1zer,zero,prd) ! the value of d(Ai(0))/dz ierr = 0 case default ai = czero dai = czero ierr = -2000 end select end subroutine airy_aic !*** !************************************************************************ !*** subroutine airy_ai_rayc (rin, zph, ai, dai, ierr, & argument_switch, modify_switch) real(prd), intent(in), dimension(:) :: rin real(prd), intent(in) :: zph complex(prd), intent(out), dimension(:) :: ai, dai integer, intent(out) :: ierr logical, intent(in) :: argument_switch logical, intent(in), optional :: modify_switch !* ! This subroutine is the driver routine that computes the ! Airy function Ai(z) and its derivative for complex arguments ! on rays directed from the origin. The subroutines called are: ! parameter_airy = determines the machine specific parameters ! flowsc = checks for over- and under- flow regions ! airy_aic = asymptotic expansion (direct evaluation), !* real(prd) h integer iplace, i, j, iregion, iflag2 complex(prd) ai0s, ai1s, G_z, dG_z, g_mz, dG_mz complex(prd), dimension(2,2) :: tsm ! if (.not. is_init_airy) call parameter_airy iflag = 0; iflag2 = 0 if ( argument_switch .eqv. .false. ) then ai = czero dai = czero ierr = -10 return end if arg_local = argument_switch mod_local = .false. if (present(modify_switch)) mod_local = modify_switch theta_global = zph ! ! run the routine for each point on the ray do j = 1,size(rin) r_global = rin(j) z_global = r_global*exp(ciunit*theta_global) if (r_global < zero) iflag = -10 call flowsc if (iflag == 0 .or. iflag == 1) then if (j == 1 .and. minval(rin) < r_min) call biggrid zeta_global = two_third*z_global*sqrt(z_global) ! ! choose the region where z lies according to the paper iregion = 0 if (r_global >= r_min) then if (abs(theta_global) <= two_pi_third) then iregion = 1 else iregion = 2 end if else if (abs(theta_global) <= pi_thirds) then iregion = 3 else iregion = 4 end if end if ! ! call the appropriate subroutine to evaluate Ai(z) and/or Ai'(z) select case (iregion) case(1) call asympc (G_z, dG_z) ai0s = G_z ai1s = dG_z case(2) if (iflag == 1) then ! exp(2*zeta) < underflow limit call asympc (G_z, dG_z) ai0s = G_z ai1s = dG_z else call asympc (G_z, dG_z, G_mz, dG_mz) if (theta_global > zero) then ai0s = G_z + ciunit*G_mz*exp(two*zeta_global) ai1s = dG_z - ciunit*dG_mz*exp(two*zeta_global) else ai0s = G_z - ciunit*G_mz*exp(two*zeta_global) ai1s = dG_z + ciunit*dG_mz*exp(two*zeta_global) end if end if case(3:4) h = r_min/n_parts if (iregion == 3) then do i = n_parts-1,0,-1 if ( r_global >= h*i ) then iplace = i+1 exit end if end do else do i = n_parts-1,0,-1 if ( r_global >= h*i ) then iplace = i exit end if end do end if call taylorc (1, h*iplace, r_global, tsm) ai0s = tsm(1,1)*aigrid(1,iplace) + tsm(1,2)*aigrid(2,iplace) ai1s = tsm(2,1)*aigrid(1,iplace) + tsm(2,2)*aigrid(2,iplace) case default ai0s = czero ai1s = czero iflag = -1000 end select if ( (iregion == 1 .or. iregion == 2) .and. .not. mod_local ) then ai0s = exp(-zeta_global)*ai0s ai1s = exp(-zeta_global)*ai1s elseif ( (iregion == 3 .or. iregion == 4) .and. mod_local) then ai0s = exp(zeta_global)*ai0s ai1s = exp(zeta_global)*ai1s end if end if ! ! error handling and return values select case (iflag) case(:-1) ai(j) = czero dai(j) = czero case(0:1) ai(j) = ai0s dai(j) = ai1s case(5) ai(j) = cmplx(ai0zer,prd) ! the value of Ai(0) dai(j) = cmplx(ai1zer,prd) ! the value of d(Ai(0))/dz case default ai(j) = czero dai(j) = czero end select iflag2 = min(iflag2,iflag) iflag = 0 end do ierr = 0 if ( iflag2 < 0 ) ierr = -20 end subroutine airy_ai_rayc !*** !************************************************************************ !*** subroutine flowsc !* ! PRIVATE subroutine to check overflow and underflow regions !* real(prd) tol, val ! if (iflag /= 0) return ! ! underflow for small |z| if (r_global <= r_lolimit) then iflag = 5 return end if ! ! test to see if r**(3/2) is out of the range of floating point if (r_global >= r_uplimit) then iflag = -6 return end if ! ! underflow in exp(2*xi) for large |z| and 2*pi/3 < |arg(z)| < pi. ! the value of the series is O(1), so once the difference between ! exp(2*xi) and G_z is larger than epsilon(tol), then the ! contribution is negligible. So we can make this tolerance ! significantly larger than underflow/overflow limit. tol = min(log(huge(tol)),-log(tiny(tol)))-15.0_prd val = abs(real(z_global*sqrt(z_global),prd)) if (abs(theta_global) > two_pi_third .and. val >= 0.75_prd*tol ) then iflag = 1 end if ! test to see if multiplication by exp(=/-xi) for unscaled option ! will over- or under-flow tol = min( -log(two*sqrt(pi)*abs(z_global**0.25_prd)) - log(tiny(tol)), & -log(two*sqrt(pi)/abs(z_global**0.25_prd)) - log(tiny(tol)), & log(two*sqrt(pi)/abs(z_global**0.25_prd)) + log(huge(tol)), & log(two*sqrt(pi)*abs(z_global**0.25_prd)) + log(huge(tol)) ) & - 15.0_prd if (.not. mod_local) then if ( val >= 1.50_prd*tol ) then iflag = -7 return end if end if end subroutine flowsc !*** !************************************************************************ !*** subroutine asympc (G_z, dG_z, G_mz, dG_mz, r_in) complex(prd), intent(out) :: G_z, dG_z complex(prd), intent(out), optional :: G_mz, dG_mz real(prd), intent(in), optional :: r_in !* ! PRIVATE subroutine to compute the asymptotic expansion for large |z| !* integer i complex(prd) zetar, z_local ! ! summation of asymptotic series is done using nested multiplication ! Compute G(zeta), G(-zeta), G'(zeta), G'(-zeta) if ( present(r_in) ) then z_local = r_in*exp(ciunit*theta_global) zetar = cunit/(two_third*z_local*sqrt(z_local)) else z_local = z_global zetar = cunit/zeta_global end if G_z = cmplx(ucoef(n_asymp),zero,prd) dG_z = cmplx(vcoef(n_asymp),zero,prd) do i = n_asymp-1, 1, -1 G_z = cmplx(ucoef(i),zero,prd) - zetar*G_z dG_z = cmplx(vcoef(i),zero,prd) - zetar*dG_z end do G_z = (cunit - zetar*G_z)/(two_sqrpi*z_local**fourth) dG_z = -(cunit - zetar*dG_z)/two_sqrpi*z_local**fourth if (present(G_mz)) then G_mz = cmplx(ucoef(n_asymp),zero,prd) dG_mz = cmplx(vcoef(n_asymp),zero,prd) do i = n_asymp-1, 1, -1 G_mz = cmplx(ucoef(i),zero,prd) + zetar*G_mz dG_mz = cmplx(vcoef(i),zero,prd) + zetar*dG_mz end do G_mz = (cunit + zetar*G_mz)/(two_sqrpi*z_local**fourth) dG_mz = -(cunit + zetar*dG_mz)/two_sqrpi*z_local**fourth end if end subroutine asympc !*** !************************************************************************ !*** subroutine taylorc (nn, r_start, r_stop, tsm, index) integer, intent(in) :: nn real(prd), intent(in) :: r_start, r_stop complex(prd), dimension(2,2), intent(out) :: tsm integer, optional, intent(in) :: index !* ! PRIVATE subroutine to integrate along a ray !* integer i, j complex(prd) h, xm, z_start, z_stop complex(prd), dimension(0:n_taylor) :: pterm, qterm complex(prd), dimension(nn,2,2) :: Phi ! j = 0 ! local error indicator--trap small step size if ( r_start == zero .and. abs(log(r_stop)) < epsilon(one) & .and. abs(r_stop-one) > two*epsilon(one) ) j = 1 if ( r_start /= zero .and. r_stop /= zero ) then if ( abs(log(r_start/r_stop)) < epsilon(one) ) j = 1 end if if ( j == 1 ) then tsm = reshape( (/ cunit, czero, czero, cunit /), (/ 2,2 /) ) return end if z_start = r_start*exp( ciunit*theta_global) z_stop = z_global if ( r_stop /= r_global ) z_stop = r_stop*exp( ciunit*theta_global) h = (z_stop-z_start)/nn ! ! compute the Taylor series in each partition do i = 1, nn xm = z_start+h*(i-1) ! ! compute the reduced derivatives: pterm(0) = cunit; pterm(1) = czero; pterm(2) = xm*pterm(0)*half qterm(0) = czero; qterm(1) = cunit; qterm(2) = czero do j = 3,n_taylor pterm(j) = (xm*pterm(j-2) + pterm(j-3))/real(j*j-j,prd) qterm(j) = (xm*qterm(j-2) + qterm(j-3))/real(j*j-j,prd) end do ! ! now sum the series Phi(i,1,1) = pterm(n_taylor) Phi(i,2,1) = pterm(n_taylor)*n_taylor Phi(i,1,2) = qterm(n_taylor) Phi(i,2,2) = qterm(n_taylor)*n_taylor do j = n_taylor-1,1,-1 Phi(i,1,1) = pterm(j) + h*Phi(i,1,1) Phi(i,2,1) = pterm(j)*j + h*Phi(i,2,1) Phi(i,1,2) = qterm(j) + h*Phi(i,1,2) Phi(i,2,2) = qterm(j)*j + h*Phi(i,2,2) end do Phi(i,1,1) = pterm(0) + h*Phi(i,1,1) Phi(i,1,2) = qterm(0) + h*Phi(i,1,2) end do ! ! multiply all the matrices together to obtain the fundamental solution matrix ! for the entire ray. do i = 1,nn-1 Phi(i+1,:,:) = matmul(Phi(i+1,:,:),Phi(i,:,:)) end do ! ! return values tsm(:,:) = Phi(nn,:,:) ! ! populuate a vector full of function and derivative values that ! lie along a given value of theta if (present(index)) then if (abs(theta_global) <= pi_thirds) then do i = 1,n_parts-1 aigrid(:,n_parts-i) = matmul(Phi(i,:,:),aigrid(:,n_parts)) end do else do i = 1,n_parts-1 aigrid(:,i) = matmul(Phi(i,:,:),aigrid(:,0)) end do end if end if end subroutine taylorc !*** !************************************************************************ !*** subroutine biggrid !* ! PRIVATE subroutine to generate the values of the Airy fucntions ! at the grid points along a ray. !* complex(prd), dimension(2,2) :: tsm complex(prd) G_z, dG_z, G_mz, dG_mz, zeta_local ! if (.not.allocated(aigrid)) allocate(aigrid(2,0:n_parts)) ! ! values at the origin are well known aigrid(1,0) = cmplx(ai0zer,zero,prd) aigrid(2,0) = cmplx(ai1zer,zero,prd) ! ! values on the circle are computed from the asymptotic expansions zeta_local = r_min*exp(ciunit*theta_global) ! dummy storage zeta_local = two_third*zeta_local*sqrt(zeta_local) if (abs(theta_global) <= two_pi_third) then ! region (I) call asympc (G_z, dG_z, r_in=r_min) aigrid(1,n_parts) = G_z aigrid(2,n_parts) = dG_z else ! region(II) call asympc (G_z, dG_z, G_mz, dG_mz, r_min) if (theta_global > zero) then aigrid(1,n_parts) = (G_z + ciunit*G_mz*exp(two*zeta_local)) aigrid(2,n_parts) = (dG_z - ciunit*dG_mz*exp(two*zeta_local)) else aigrid(1,n_parts) = (G_z - ciunit*G_mz*exp(two*zeta_local)) aigrid(2,n_parts) = (dG_z + ciunit*dG_mz*exp(two*zeta_local)) end if end if aigrid(:,n_parts) = exp(-zeta_local)*aigrid(:,n_parts) ! ! now integrate once and populate the vector of function and derivative ! values on the grid if (abs(theta_global) <= pi_thirds) then ! integrate towards the origin call taylorc (n_parts, r_min, zero, tsm, 1) else ! integrate away from the origin call taylorc (n_parts, zero, r_min, tsm, 1) end if end subroutine biggrid !*** !************************************************************************ !*** ! ! end file airy_functions_complex ! ! begin file airy_zeros ! version 1.0 ! edited 07/24/2002 ! !*** !************************************************************************ !*** subroutine ai_zeroc (n, ai_zero, ierr, ai_assoc, dai_zero, dai_assoc) integer, intent(in) :: n complex(prd), intent(out) :: ai_zero integer, intent(out) :: ierr complex(prd), intent(out), optional :: ai_assoc, dai_zero, dai_assoc !* ! PUBLIC subroutine which returns zeros with an error flag since ! Ai(z) and Ai'(z), do not have non-real zeros. !* ai_zero = czero if (present(ai_assoc)) ai_assoc = czero if (present(dai_zero)) dai_zero = czero if (present(dai_assoc)) dai_assoc = czero ierr = n; ierr = -3 ! the first is a dummy statement to use n ! end subroutine ai_zeroc !*** !************************************************************************ !*** subroutine bi_zeroc (n, bi_zero, ierr, bi_assoc, dbi_zero, dbi_assoc) integer, intent(in) :: n complex(prd), intent(out) :: bi_zero integer, intent(out) :: ierr complex(prd), intent(out), optional :: bi_assoc, dbi_zero, dbi_assoc !* ! PUBLIC subroutine which is the driver routine that computes the ! zeros of Bi(z) and optionally of Bi'(z), and the associated ! values Bi'(beta_n) and Bi(beta_n') for scalar arguments. The first 25 ! zeros are stored in an array and the rest are computed by summing ! approriate asymptotic expansions. The (PRIVATE) subroutines called are: ! parameters_airy = determines the machine specific parameters ! zero_parameters_airy = populates the vectors containing the ! stored zeros and the coefficients of the asymptotic expansions ! airy_bic = computes the associated function value for stored zeros ! ae_zero_c = sums the appropriate asymptotic expansion for the ! zeros of the Airy functions and their associated function values !* integer itemp complex(prd) xbi, xdbi, lam ! if (.not. is_zero_init_airy) call zero_parameter_airy iflag = 0 if (n <= 0) iflag = -3 itemp = floor(0.25*(huge(itemp)-1)) if (n >= itemp) iflag = -4 ! zero requested cannot be computed if (iflag < 0) then bi_zero = zero if (present(bi_assoc)) bi_assoc = czero if (present(dbi_zero)) dbi_zero = czero if (present(dbi_assoc)) dbi_assoc = czero ierr = iflag return end if if (n <= n_zeros) then bi_zero = bizc(n) if (present(bi_assoc)) then call airy_aic (bi_zero*exp(-ciunit*two_third*pi), lam, xdbi, ierr) call airy_aic (bi_zero*exp( ciunit*two_third*pi), lam, xbi, ierr) lam = 5.0_prd*pi*ciunit/six bi_assoc = exp(-lam)*xdbi + exp(lam)*xbi end if if (present(dbi_zero)) dbi_zero = dbizc(n) if (present(dbi_assoc)) then call airy_aic (dbizc(n)*exp(-ciunit*two_third*pi), xdbi, lam, ierr) call airy_aic (dbizc(n)*exp( ciunit*two_third*pi), xbi, lam, ierr) lam = pi*ciunit/six dbi_assoc = exp(-lam)*xdbi + exp(lam)*xbi end if else lam = cmplx(three_pi_ate*real(4*n-1,prd), & 0.75_prd*log(2.0_prd),prd) call ae_zero_c (xbi, lam, 1) bi_zero = exp(ciunit*pi/three)*xbi call ae_zero_c (xbi, lam, 3) if (present(bi_assoc)) then bi_assoc = sqrt(two)*exp(-ciunit*pi/six)*xbi if ( mod(n,2) == 1 ) bi_assoc = -bi_assoc end if if (present(dbi_zero) .or. present(dbi_assoc)) then lam = cmplx(three_pi_ate*real(4*n-3,prd),& 0.75_prd*log(2.0_prd),prd) call ae_zero_c (xbi, lam, 2) if (present(dbi_zero)) dbi_zero = exp(ciunit*pi/three)*xbi call ae_zero_c (xdbi, lam, 4) if (present(dbi_assoc)) then dbi_assoc = sqrt(two)*exp(ciunit*pi/six)*xdbi if ( mod(n-1,2) == 1 ) dbi_assoc = -dbi_assoc end if end if end if end subroutine bi_zeroc !*** !************************************************************************ !*** subroutine ae_zero_c (zr, lam, in) complex(prd), intent (out) :: zr complex(prd), intent (in) :: lam integer, intent (in) :: in !* ! PRIVATE subroutine to sum the asymptotic expansions for the zeros ! of the Airy functions and their associated functions !* complex(prd) zrsum, lamr integer i ! lamr = one/lam**2 select case (in) case(1) ! sum the T(x) expansion zrsum = cmplx(Tcoeff(n_asymp_zero),zero,prd) do i = n_asymp_zero-1,0,-1 zrsum = cmplx(Tcoeff(i),zero,prd) + zrsum*lamr end do zr = zrsum*lam**two_third case(2) ! sum the U(x) expansion zrsum = cmplx(Ucoeff(n_asymp_zero),zero,prd) do i = n_asymp_zero-1,0,-1 zrsum = cmplx(Ucoeff(i),zero,prd) + zrsum*lamr end do zr = zrsum*lam**two_third case(3) ! sum the V(x) expansion zrsum = cmplx(Vcoeff(n_asymp_asso),zero,prd) do i = n_asymp_asso-1,0,-1 zrsum = cmplx(Vcoeff(i),zero,prd) + zrsum*lamr end do zr = zrsum*lam**(half/three)/sqrpi case(4) ! sum the W(x) expansion zrsum = cmplx(Wcoeff(n_asymp_asso),zero,prd) do i = n_asymp_asso-1,0,-1 zrsum = cmplx(Wcoeff(i),zero,prd) + zrsum*lamr end do zr = zrsum/lam**(half/three)/sqrpi end select end subroutine ae_zero_c !*** !************************************************************************ !*** subroutine ai_zerocv (n, ai_zero, ierr, ai_assoc, dai_zero, dai_assoc) integer, intent(in) :: n complex(prd), intent(out), dimension(:) :: ai_zero integer, intent(out) :: ierr complex(prd), intent(out), dimension(:), optional :: & ai_assoc, dai_zero, dai_assoc !* !* ! PUBLIC subroutine which returns zeros with an error flag since ! Ai(z) and Ai'(z), do not have non-real zeros. !* ai_zero = czero if (present(ai_assoc)) ai_assoc = czero if (present(dai_zero)) dai_zero = czero if (present(dai_assoc)) dai_assoc = czero ierr = n; ierr = -3 ! the first is a dummy statement to use n ! ! end subroutine ai_zerocv !*** !************************************************************************ !*** subroutine bi_zerocv (n, bi_zero, ierr, bi_assoc, dbi_zero, dbi_assoc) integer, intent(in) :: n complex(prd), intent(out), dimension(:) :: bi_zero integer, intent(out) :: ierr complex(prd), intent(out), dimension(:), optional :: & bi_assoc, dbi_zero, dbi_assoc !* !* ! PUBLIC subroutine which is the driver routine that computes the ! zeros of Bi(z) and optionally of Bi'(z), and the associated ! values Bi'(beta_n) and Bi(beta_n') for vector arguments. The first 25 ! zeros are stored in an array and the rest are computed by summing ! approriate asymptotic expansions. The (PRIVATE) subroutines called are: ! parameters_airy = determines the machine specific parameters ! zero_parameters_airy = populates the vectors containing the ! stored zeros and the coefficients of the asymptotic expansions ! airy_bic = computes the associated function value for stored zeros ! ae_zero_cv = sums the appropriate asymptotic expansion for the ! zeros of the Airy functions and their associated function values !* integer itemp, K, i, iflag_zero complex(prd) xbi, xdbi, lam1 complex(prd), dimension(:), allocatable :: lam ! if (.not. is_zero_init_airy) call zero_parameter_airy iflag_zero = 0 if (n <= 0) iflag_zero = -3 itemp = floor(0.25*(huge(itemp)-1)) if (n >= itemp) iflag_zero = -4 ! zero requested cannot be computed K = size(bi_zero) ! ! check to see if the last zero requested is computable ! and change the last zero returned accordingly. if (n+K >= itemp .and. iflag_zero /= -4) then do i = n+K-1,n,-1 if ( i < itemp ) then K = i-n iflag_zero = 50 exit end if end do end if if (iflag_zero < 0) then bi_zero(:) = czero if (present(bi_assoc)) bi_assoc(:) = czero if (present(dbi_zero)) dbi_zero(:) = czero if (present(dbi_assoc)) dbi_assoc(:) = czero ierr = iflag_zero return end if if (n <= n_zeros) then itemp = min(n+K,n_zeros) bi_zero(1:itemp-n+1) = bizc(n:itemp) if (present(bi_assoc)) then do i = n, itemp call airy_aic (bizc(i)*exp(-ciunit*two_third*pi), lam1, xdbi, ierr) call airy_aic (bizc(i)*exp( ciunit*two_third*pi), lam1, xbi, ierr) lam1 = 5.0_prd*pi*ciunit/six bi_assoc(i-n+1) = exp(-lam1)*xdbi + exp(lam1)*xbi end do end if if (present(dbi_zero)) dbi_zero(1:itemp-n+1) = dbizc(n:itemp) if (present(dbi_assoc)) then do i = n, itemp call airy_aic (dbizc(i)*exp(-ciunit*two_third*pi), xdbi, lam1, ierr) call airy_aic (dbizc(i)*exp( ciunit*two_third*pi), xbi, lam1, ierr) lam1 = pi*ciunit/six dbi_assoc(i-n+1) = exp(-lam1)*xdbi + exp(lam1)*xbi end do end if end if if (n+K > n_zeros) then itemp = max(n,n_zeros+1) allocate (lam(itemp:n+K-1)) do i = itemp,n+K-1 lam(i) = cmplx(three_pi_ate*real(4*i-1,prd), & 0.75_prd*log(2.0_prd),prd) end do call ae_zero_cv (bi_zero(itemp-n+1:K), lam, 1) bi_zero(itemp-n+1:K) = exp(ciunit*pi/three)*bi_zero(itemp-n+1:K) if (present(bi_assoc)) then call ae_zero_cv (bi_assoc(itemp-n+1:K), lam, 3) bi_assoc(itemp-n+1:K) = & sqrt(two)*exp(-ciunit*pi/six)*bi_assoc(itemp-n+1:K) do i = itemp,n+K-1 if ( mod(i,2) == 1 ) bi_assoc(i-n+1) = -bi_assoc(i-n+1) end do end if deallocate (lam) if (present(dbi_zero) .or. present(dbi_assoc)) then allocate (lam(itemp:n+K-1)) do i = itemp,n+K-1 lam(i) = cmplx(three_pi_ate*real(4*i-3,prd), & 0.75_prd*log(2.0_prd),prd) end do if (present(dbi_zero)) then call ae_zero_cv (dbi_zero(itemp-n+1:K), lam, 2) dbi_zero(itemp-n+1:K) = exp(ciunit*pi/three)*dbi_zero(itemp-n+1:K) end if if (present(dbi_assoc)) then call ae_zero_cv (dbi_assoc(itemp-n+1:K), lam, 4) dbi_assoc(itemp-n+1:K) = & sqrt(two)*exp(ciunit*pi/six)*dbi_assoc(itemp-n+1:K) do i = itemp,n+K-1 if ( mod(i-1,2) == 1 ) dbi_assoc(i-n+1) = -dbi_assoc(i-n+1) end do end if deallocate (lam) end if end if ierr = iflag_zero end subroutine bi_zerocv !*** !************************************************************************ !*** subroutine ae_zero_cv (zr, lam, in) complex(prd), intent (out), dimension(1:) :: zr complex(prd), intent (in), dimension(1:) :: lam integer, intent (in) :: in complex(prd), dimension(size(zr)) :: zrsum, lamr !* ! PRIVATE subroutine to sum the asymptotic expansions for the zeros ! of the Airy functions and their associated functions for the ! vector case !* integer k, i, j ! k = size(zr) do i = 1,k lamr(i) = cunit/lam(i)**2 end do if (in == 1) then ! sum the T(x) expansion do j = 1,k zrsum(j) = cmplx(Tcoeff(n_asymp_zero),zero,prd) do i = n_asymp_zero-1,0,-1 zrsum(j) = cmplx(Tcoeff(i),zero,prd) + zrsum(j)*lamr(j) end do end do zr(:) = zrsum(:)*lam(:)**two_third elseif (in == 2) then ! sum the U(x) expansion do j = 1,k zrsum(j) = cmplx(Ucoeff(n_asymp_zero),zero,prd) do i = n_asymp_zero,0,-1 zrsum(j) = cmplx(Ucoeff(i),zero,prd) + zrsum(j)*lamr(j) end do end do zr(:) = zrsum(:)*lam(:)**two_third elseif (in == 3) then ! sum the V(x) expansion do j = 1,k zrsum(j) = cmplx(Vcoeff(n_asymp_asso),zero,prd) do i = n_asymp_asso-1,0,-1 zrsum(j) = cmplx(Vcoeff(i),zero,prd) + zrsum(j)*lamr(j) end do end do zr(:) = zrsum(:)*lam(:)**(half/three)/sqrpi else ! sum the W(x) expansion do j = 1,k zrsum(j) = cmplx(Wcoeff(n_asymp_asso),zero,prd) do i = n_asymp_asso-1,0,-1 zrsum(j) = cmplx(Wcoeff(i),zero,prd) + zrsum(j)*lamr(j) end do end do zr(:) = zrsum(:)/lam(:)**(half/three)/sqrpi end if end subroutine ae_zero_cv !*** !************************************************************************ !*** ! ! end file airy_zeros SHAR_EOF fi # end of overwriting check if test -f 'airy_functions.f90' then echo shar: will not over-write existing file "'airy_functions.f90'" else cat << "SHAR_EOF" > 'airy_functions.f90' ! ! version 1.0 ! built Fri Apr 16 18:03:05 CDT 2004 ! ! by B.R. Fabijonas ! Department of Mathematics ! Southern Methodist University ! bfabi@smu.edu ! ! see the file airy_README for an explanation !*** !************************************************************************ !*** module airy_functions_real_single implicit none private public :: airy_ai, airy_bi, airy_ai_zero, airy_bi_zero public :: airy_info, airy_aux, airy_aux_info integer, parameter :: prd = kind(0.0e0) include 'airy_head' !* ! interfaces !* interface airy_ai module procedure airy_air end interface interface airy_bi module procedure airy_bir end interface interface airy_ai_zero module procedure ai_zeror module procedure ai_zerorv end interface interface airy_bi_zero module procedure bi_zeror module procedure bi_zerorv end interface interface airy_info module procedure airy_paramsr end interface interface airy_aux module procedure airy_auxr end interface interface airy_aux_info module procedure airy_paramsr_aux end interface contains include 'airy_real' include 'airy_parameters' end module airy_functions_real_single !*** !************************************************************************ !*** module airy_functions_complex_single implicit none private public :: airy_ai, airy_ai_zero, airy_bi_zero integer, parameter :: prd = kind(0.0e0) include 'airy_head' !* ! interfaces !* interface airy_ai module procedure airy_aic module procedure airy_ai_rayc end interface interface airy_ai_zero module procedure ai_zeroc module procedure ai_zerocv end interface interface airy_bi_zero module procedure bi_zeroc module procedure bi_zerocv end interface contains include 'airy_complex' include 'airy_parameters' end module airy_functions_complex_single !*** !************************************************************************ !*** module airy_functions_real_double implicit none private public :: airy_ai, airy_bi, airy_ai_zero, airy_bi_zero public :: airy_info, airy_aux, airy_aux_info integer, parameter :: prd = kind(0.0d0) include 'airy_head' !* ! interfaces !* interface airy_ai module procedure airy_air end interface interface airy_bi module procedure airy_bir end interface interface airy_ai_zero module procedure ai_zeror module procedure ai_zerorv end interface interface airy_bi_zero module procedure bi_zeror module procedure bi_zerorv end interface interface airy_info module procedure airy_paramsr end interface interface airy_aux module procedure airy_auxr end interface interface airy_aux_info module procedure airy_paramsr_aux end interface contains include 'airy_real' include 'airy_parameters' end module airy_functions_real_double !*** !************************************************************************ !*** module airy_functions_complex_double implicit none private public :: airy_ai, airy_ai_zero, airy_bi_zero integer, parameter :: prd = kind(0.0d0) include 'airy_head' !* ! interfaces !* interface airy_ai module procedure airy_aic module procedure airy_ai_rayc end interface interface airy_ai_zero module procedure ai_zeroc module procedure ai_zerocv end interface interface airy_bi_zero module procedure bi_zeroc module procedure bi_zerocv end interface contains include 'airy_complex' include 'airy_parameters' end module airy_functions_complex_double !*** !************************************************************************ !*** ! module airy_functions_real_quad ! implicit none ! private ! public :: airy_ai, airy_bi, airy_ai_zero, airy_bi_zero ! public :: airy_info, airy_aux, airy_aux_info ! integer, parameter :: prd = kind(0.0q0) ! include 'airy_head' !!* !! interfaces !!* ! interface airy_ai ! module procedure airy_air ! end interface ! interface airy_bi ! module procedure airy_bir ! end interface ! interface airy_ai_zero ! module procedure ai_zeror ! module procedure ai_zerorv ! end interface ! interface airy_bi_zero ! module procedure bi_zeror ! module procedure bi_zerorv ! end interface ! interface airy_info ! module procedure airy_paramsr ! end interface ! interface airy_aux ! module procedure airy_auxr ! end interface ! interface airy_aux_info ! module procedure airy_paramsr_aux ! end interface ! contains ! include 'airy_real' ! include 'airy_parameters' ! end module airy_functions_real_quad !!*** !!************************************************************************ !!*** ! module airy_functions_complex_quad ! implicit none ! private ! public :: airy_ai, airy_ai_zero, airy_bi_zero ! integer, parameter :: prd = kind(0.0q0) ! include 'airy_head' !!* !! interfaces !!* ! interface airy_ai ! module procedure airy_aic ! module procedure airy_ai_rayc ! end interface ! interface airy_ai_zero ! module procedure ai_zeroc ! module procedure ai_zerocv ! end interface ! interface airy_bi_zero ! module procedure bi_zeroc ! module procedure bi_zerocv ! end interface ! contains ! include 'airy_complex' ! include 'airy_parameters' ! end module airy_functions_complex_quad !!*** !!************************************************************************ !!*** SHAR_EOF fi # end of overwriting check if test -f 'airy_head' then echo shar: will not over-write existing file "'airy_head'" else cat << "SHAR_EOF" > 'airy_head' ! ! begin file airy_head ! version 1.0 ! edited 12/27/2001 ! revised by DWL April 2001 ! revised by BRF January 2002 ! ! Global constants integer, parameter :: & n_coeff_asso = 63, & n_coeff_phase = 50, & n_coeff_zero = 56, & n_zeros = 25 real(prd), parameter :: & zero = 0.0_prd, & one = 1.0_prd, & two = 2.0_prd, & three = 3.0_prd, & four = 4.0_prd, & six = 6.0_prd, & fourth = 0.25_prd, & half = 0.5_prd, & two_third = & 0.66666666666666666666666666666666666666666666666666666666666666666667_prd, & ai0zer = & 0.35502805388781723926006318600418317639797917419917724058332651030081_prd, & ai1zer = & -0.25881940379280679840518356018920396347909113835493458221000181385610_prd, & bi0zer = & 0.61492662744600073515092236909361355359472818864859650504087875301430_prd, & bi1zer = & 0.44828835735382635791482371039882839086622679921226206108280877837233_prd, & pi = & 3.14159265358979323846264338327950288419716939937510582097494459230782_prd, & sqrpi = & 1.77245385090551602729816748334114518279754945612238712821380778985291_prd, & pi_thirds = pi/3, & two_pi_third = 2*pi/3, & two_sqrpi = 2*sqrpi, & three_pi_ate = 3*pi/8 complex(prd), parameter :: & czero = (0.0_prd,0.0_prd), & ! complex zero ciunit = (0.0_prd,1.0_prd), & ! complex imaginary number cunit = (1.0_prd,0.0_prd) ! complex one ! ! Global variables integer :: & iflag = 0, & n_asymp = 0, & n_asymp_asso = 0, & n_asymp_mod = 0, & n_asymp_phase = 0, & n_asymp_zero = 0, & n_parts = 0, & n_parts_aux = 0, & n_taylor = 0, & n_taylor_aux = 0 real(prd) :: & aux_min = 0.0_prd, & fstpsz = 0.0_prd, & mstpsz = 0.0_prd, & r_min = 0.0_prd, & r_uplimit = 0.0_prd, & r_lolimit = 0.0_prd, & theta_global, & xi_global, & x_global, & r_global complex(prd) :: & zeta_global, & z_global real(prd), dimension(:,:), allocatable :: & m2grid real(prd), dimension(:), allocatable :: & mcoef, & ncoef, & theta_grid, & ucoef, & vcoef complex(prd), dimension(n_zeros) :: & bizc, & ! complex zeros dbizc ! complex zeros complex(prd), dimension(:,:), allocatable :: & aigrid logical :: & big_integrate_aux = .false., & is_aux_init_airy = .false., & is_init_airy = .false., & is_zero_init_airy = .false., & arg_local = .false., & mod_local = .false. ! ! stored values of the first 25 zeros, the first 45 coefficients ! for the asymptotic expansions of the zeros, the first 50 ! coefficients for the asymptotic expansions for the associated ! functions, and the first 50 coefficients for the asymptotic ! expansions for the phase functions. real(prd), dimension(n_zeros) :: & aizr, & bizr, & daizr, & dbizr real(prd), dimension(0:n_coeff_zero-1) :: & Tcoeff, & Ucoeff real(prd), dimension(0:n_coeff_asso-1) :: & Vcoeff, & Wcoeff real(prd), dimension(0:n_coeff_phase-1) :: & thetacoeff, & phicoeff ! ! what follows is the hard coded coefficients of the ! asymptotic series of the zeros of the Airy functions. data Ucoeff(0:18) / & 1.0_prd, & -.14583333333333333333333333333333333333333333333333333333333333333333_prd,& .12152777777777777777777777777777777777777777777777777777777777777777_prd,& -.87395351080246913580246913580246913580246913580246913580246913580246_prd,& .15016855549125514403292181069958847736625514403292181069958847736625e2_prd,& -.47694644148817441754543895747599451303155006858710562414266117969821e3_prd,& .24241309444872289030901191917459921575147912596472267254160258275484e5_prd,& -.18036569459973376725878994930524510053402371646541742563690437490163e7_prd,& .18486906478642827292616911910089915509011302335490721407502386011759e9_prd,& -.24975761293716812237485135712074343571689388225751805089702927291450e11_prd,& .43010597766029133575957465515041682975116407364519869082404185791498e13_prd,& -.91966630646216224414462514444769376815822305958737503992528419520874e15_prd,& .23905780559290646307566900546728188233693016653754929837863494781017e18_prd,& -.74241332478952472290363916933866251290404112455691659500004579531135e20_prd,& .27148400107512125552064168030849313000106437718598339998802868110263e23_prd,& -.11546167442336934169526883463106999395192635312191589413884964066192e26_prd,& .56509197949692565122879784149625473179282442934941884667921655590141e28_prd,& -.31534468131244981806142918984101053489460643574011303209334451335403e31_prd,& .19903407540014516438638931066554027629203641348823621813687428803681e34_prd/ data Ucoeff(19:55) / 37*zero / ! ! the next set of coeffs are for T(x). data Tcoeff(0:18) / & 1.0_prd, & .10416666666666666666666666666666666666666666666666666666666666666667_prd,& -.13888888888888888888888888888888888888888888888888888888888888888889_prd,& .92984423225308641975308641975308641975308641975308641975308641975309_prd,& -.15509155201673647854203409758965314520870076425631981187536743092299e2_prd,& .48552909692595711928032529884381736233588085439937291789143640995492e3_prd,& -.24505300371559153682986228639520820590779438516064030467322648392607e5_prd,& .18166531815056219864955679723943719999961426573223555390907928630837e7_prd,& -.18581460679897715056836004159090076815148668306941218300549412966945e9_prd,& .25071304911268777333357715167823935456935326637881744159381973639923e11_prd,& -.43138527109842934644759345584413150635257003381911139923300864298731e13_prd,& .92185668129464141695987363212833552615161591992359885620965411318628e15_prd,& -.23952433840095928791094332571942871169658682388152814639232211028922e18_prd,& .74362270472474507345398633000288750016206828627391077277289564070702e20_prd,& -.27185882859962273722339835257780075234211548224517714800753203346127e23_prd,& .11559853167916871386933427721989223661195069934160109595184381087791e26_prd,& -.56567334133206565368102758661503886004417965273912945973100005362259e28_prd,& .31562894324838583270281305634852764403638314020995636590825013151460e31_prd,& -.19919258451508720424348231287353840437868650493382105342718540790297e34_prd / data Tcoeff(19:55) / 37*zero / ! ! the next set of data are the coefficients of the ! asymptotic series of the associated functions at the zeros of the ! Airy functions. data Wcoeff(0:18) / & 1.0_prd, & -.07291666666666666666666666666666666666666666666666666666666666666667_prd,& .27229817708333333333333333333333333333333333333333333333333333333333_prd,& -.31796569447458526234567901234567901234567901234567901234567901234568e1_prd,& .76556611494861021944524819958847736625514403292181069958847736625514e2_prd,& -.31342883391504333097748305081340020576131687242798353909465020576132e4_prd,& .19526467547467580356141263657907492119681468318299594019758628812127e6_prd,& -.17214881804894471760366662537203599796978117784592300195866725359181e8_prd,& .20404480428776782808583230846549783277836178166084920127306154826731e10_prd,& -.31299949629291202854512623939968839596125044797653615695230422712467e12_prd,& .60336461605090016658418262949946148185490712235984520011916410419179e14_prd,& -.14278069090774911120088839828530748584551695465524445535742732768897e17_prd,& .40694463078143798287130215270391777070050987066122969894115089940327e19_prd,& -.13750138163790283482202930118673116996392571046076704729266812849012e22_prd,& .54348926796785052801475752640893088472995254326308568359735532603746e24_prd,& -.24844792999331365214998636162748338476863223103457160185005725731021e27_prd,& .13006474905307843236575513358891917058326565765494198896278465517184e30_prd,& -.77308324290748324739057865230742197678773712235823741470353250780808e32_prd,& .51777850299807518173142784385934789871927480162907991284189078357799e35_prd/ data Wcoeff(19:62) / 44*zero / ! ! the second 50 coeffs are for V(x). data Vcoeff(0:18) / & 1.0_prd,& .10416666666666666666666666666666666666666666666666666666666666666667_prd,& -.33094618055555555555555555555555555555555555555555555555555555555556_prd,& .36136956862461419753086419753086419753086419753086419753086419753086e1_prd,& -.83984945326395275871729864785420340975896531452087007642563198118754e2_prd,& .33685630692930845753622830626592200666274740348814422888496962571037e4_prd,& -.20713231106607890252977377856845981926357440760732941802900650637276e6_prd,& .18097282066882778449431331373601761112628551971830435479269498473888e8_prd,& -.21309444395672364695106172827605375452442300629551305981428785008379e10_prd,& .32523869187619483120330380717331928218594181576469626622284133921492e12_prd,& -.62446698182714354178599585173443015731252810847722720974668646995208e14_prd,& .14729817433919221428748862989469425786509381441758339217882596234826e17_prd,& -.41870023262835387031831169307108942313981607413837696274312354570063e19_prd,& .14115582971009629114612866526991564864316743141214936159803917741306e22_prd,& -.55686486988532093135443555068430236099429127519474369927264196276483e24_prd,& .25414120415563205117290897011609704504117958895702715714297261718928e27_prd,& -.13285320976148994930374487403381382511677438741525714111410410227917e30_prd,& .78865443622886902650543525972805065726220987310540174043394169896799e32_prd,& -.52761245869440267730251126474749725416169924433173936574780403898729e35_prd/ data Vcoeff(19:62) / 44*zero / ! ! now follows the stored zeros. ! the real zeros of Ai(x) data aizr / & -2.33810741045976703848919725244673544063854014567238785248385443721367_prd, & -4.08794944413097061663698870145739106022476469910852975498416087602512_prd, & -5.52055982809555105912985551293129357379721428061752510483288757695750_prd, & -6.78670809007175899878024638449617696605388247739349361652352909355624_prd, & -7.94413358712085312313828055579826853214067439697221480864385428571645_prd, & -9.02265085334098038015819083988008925652467753515608251556068568311328_prd, & -10.04017434155808593059455673736251809404290256910583310437135385095343_prd, & -11.00852430373326289323543964959015101673082538150403750570001348624448_prd, & -11.93601556323626251700636490293058431557788623211982397092353864291238_prd, & -12.82877675286575720040672940724182447738641559957341994148696682189440_prd, & -13.69148903521071792829569677946692054166536980920076805254071472510959_prd, & -14.52782995177533498207398144299589337871416486983482524101054549228523_prd, & -15.34075513597799685714620851348148670511758332024803704969134367360914_prd, & -16.13268515694577143934598044720252179051827239707628246640238983564054_prd, & -16.90563399742994262703523877061147659909005109503171855603649023736814_prd, & -17.66130010569705750925365030401805595215321866811995330361255823961027_prd, & -18.40113259920711541586139792950433675459381460602005817693995569303975_prd, & -19.12638047424695214412414868973249468907545838475308007653292923976344_prd, & -19.83812989172149970094756361601140419833568249453885876463439190533102_prd, & -20.53733290767756635998268141130810174530421801473749208218771674302329_prd, & -21.22482994364209695519767672540520133128372022695142223538915255700301_prd, & -21.90136759558513070740823704210580626014944052125548559380507316477778_prd, & -22.56761291749650283145917485684307040940078566103896133980710912405921_prd, & -23.22416500112168106132095039028827471366842731266977699133091665595986_prd, & -23.87156445553591856711857620915978417026368251469855098111241320250607_prd/ ! ! the real zeros of Ai'(x) data daizr / & -1.01879297164747108901732478339974382421820544125443456387086141398280_prd, & -3.24819758217983653787542377077584338415362303651722056102813298651085_prd, & -4.82009921117873563940061626041637941446268823149559840636862039390690_prd, & -6.16330735563948654763784353309140759297448387720351451555576295688298_prd, & -7.37217725504777017709218227112798901002013456098406544685663881697862_prd, & -8.48848673401972213288099541573077295112964359668407460466365614498850_prd, & -9.53544905243354747070263342700063969538820217161052760208983694238694_prd, & -10.52766039695740728197814249142222864298654636972840675098825728664860_prd, & -11.47505663348024529493721652240563873475572084899288324564501046010262_prd, & -12.38478837184574732549333923538262410099237019478484289722429657302031_prd, & -13.26221896166521038243727842123812965366637123382426125928326035400567_prd, & -14.11150197046299528163274595367191475659596057504177550100311260702737_prd, & -14.93593719672051746650728205326279421117972255282091467187721012647215_prd, & -15.73820137369253830268645323011724506466987051724640138554427300411434_prd, & -16.52050382543379354220519308946062809031267401323065118962585800030048_prd, & -17.28469505021643735664373842332082107165009909874506202863636173764009_prd, & -18.03234462250439339526520862299153702097903336787137198226440397831148_prd, & -18.76479843766595474015165121077601065287558788488499214751351991943185_prd, & -19.48322165656723117752471804759502558594797583425851677835534147425693_prd, & -20.18863150946337315365453214288926902260240432967877744616116532852433_prd, & -20.88192275551673770086895619234965436951497044315845023182432979152615_prd, & -21.56388772319897495770118028517312743827608287640899482894099397157279_prd, & -22.23523228534891333075485328085016555368139534299890933471263856691166_prd, & -22.89658873887461900146658726220533486755445824339047035924534481515681_prd, & -23.54852629592880157396398563890195889083567420203896865761957496760415_prd/ ! ! the real zeros of Bi(x) data bizr / & -1.17371322270912792491997996247390210454364638917570309757920761150430_prd, & -3.27109330283635271568022824016641380630093596910028480148503239626113_prd, & -4.83073784166201593266770933990517817696614261732301026576850194149303_prd, & -6.16985212831025125983336452055593667996554943427563109765876409289291_prd, & -7.37676207936776371359995933044254122209152229939709572992390169564234_prd, & -8.49194884650938801344803949280977672860508755505545790582286779923401_prd, & -9.53819437934623888663298854515601962083907207638247038808917424104129_prd, & -10.52991350670535792440055559845314799952957759462139619609904123067107_prd, & -11.47695355127877943792346492473281967194825381488769633431302780799070_prd, & -12.38641713858273874556190150286328094825979838468560971551502681648386_prd, & -13.26363952294180555411074332439549077524115196098126631023507247982663_prd, & -14.11275680906865779158730978222401847168404282855088116185500171775926_prd, & -14.93705741215416404020321431049090463961217635177820922699751842280182_prd, & -15.73921035119048277089497847974818338071801627678405270593898207072179_prd, & -16.52141955063437905391794996521054571671103103705808118179905226958795_prd, & -17.28553162458124253293423669225353924252797536027103490380471228329517_prd, & -18.03311328722500157217111254333919200080872914164059812874417048491807_prd, & -18.76550828448008104134297892361051284402671895514206863403515927768504_prd, & -19.48388013298923401366599865924135751220629777936103551104169553178148_prd, & -20.18924478539620242022532322582753607646497835839342208483158243219733_prd, & -20.88249599419317576784241835353106077909763612165974937464360562852603_prd, & -21.56442528471297765275553097792876056885144913635808007380299667653380_prd, & -22.23573788180338516651471705596359912920512830316306238465359314928271_prd, & -22.89706555421979347392720167028728219215245475061019081243083088117395_prd, & -23.54897707964244826910809088370627200288704237097332714530518718245203_prd/ ! ! the real zeros of Bi'(x) data dbizr / & -2.29443968261412324662245867376910282907848946700969989589640271599098_prd, & -4.07315508907182821555236851509270995079318027043612840965664763255332_prd, & -5.51239572966359949625959252047918324308288262563228155846581341204411_prd, & -6.78129444599030539002866629612268711724774958242124874315576274421256_prd, & -7.94017868916857892668015506404469358236295026378631200761623213008077_prd, & -9.01958335879423906741485643920410500793021840458929511999436467689141_prd, & -10.03769633490854580175763552931088600719650708912364209574799165928857_prd, & -11.00646266771228994042178673901908858144008056513455219417559368827981_prd, & -11.93426164501484466280156966627884018620779786420344917372782764148357_prd, & -12.82725830917721764018582515874172632262398632000325301764491321749451_prd, & -13.69015582683504910123087530032843814837512065755814073750827235404234_prd, & -14.52664576348571140952971574483676917698490608428973245022236766674096_prd, & -15.33969308224240410861550731532421181298426399724957048054681874323756_prd, & -16.13172478238590057788823351351896352801125898819864759735671923382317_prd, & -16.90475941188964995823127191076849074551409618145881961325073707696760_prd, & -17.66049874311497610219306321601848781394657665773696300832840696681477_prd, & -18.40039436718170327981549177420453616844392469219826273367311029445693_prd, & -19.12569715641263806611390647390670371376640167393442425682883753076962_prd, & -19.83749471841591050293338761137287048759231913849216259169928462093709_prd, & -20.53674024145327398011335618110930387074716842210354620328691065812538_prd, & -21.22427504488926656905931114462097972114007365699962254630891208179793_prd, & -21.90084644513920828138642436184221845938113289053710021975475917896781_prd, & -22.56712208049720046972298553127734852587993691855079121202853774320013_prd, & -23.22370152120896211611253695016184095318132604851766457576200217045542_prd, & -23.87112577167797359456463262632343377384840108602012722903210953550931_prd/ ! ! the complex zeros of Bi(z) in upper half complex plane data bizc(1:10) / & ( 0.97754488673162068594699270603101291142795187683490068888647838406443_prd, & 2.14129070603874457574913922659732897573219215679089642486333363968462_prd ),& ( 1.89677501389533634662721727829884659615177393761103648449500365404271_prd, & 3.62729176435891941044049914398672742135138471558722329404568929501847_prd ),& ( 2.63315773935494659570801882716408027064344210717059360176219749413806_prd, & 4.85546817997984498317462803408683616985194839172584257833996797129166_prd ),& ( 3.27853123615674637108428003742990269304233294063643810889158089478166_prd, & 5.94450428117905212809591786731140536872984124248655143010029596308279_prd ),& ( 3.86585273173334614251511337081452708596976194632193282962879173253901_prd, & 6.94169220958211125352971091861420315862502036505329040289270784633450_prd ),& ( 4.41161187480932549691994423741928993680626265689626333136555188876233_prd, & 7.87183965948658196641675823684140806420708760005728700488878495074197_prd ),& ( 4.92552935386139721563708520907097320258698061300099737722695463508780_prd, & 8.74998254125672034599585620155415722642063525723209180381922538805753_prd ),& ( 5.41393680880765209985945724430184554308272025802972651933170020359738_prd, & 9.58609690055480830302647362820501805743586542389874387303025339162157_prd ),& ( 5.88124675398118307316367004091001520511362531706460173192832604227755_prd, & 10.38722739030382930057781919481935780886769028619723636729641393497882_prd ),& ( 6.33068856706306813563172144477503542412410562601369891961334059189483_prd, & 11.15858122676025326996998168923801510296971883724967951816234142038045_prd ) / data bizc(11:20) / & ( 6.76471521144954609083399562153578384959795677911006057480597886621212_prd, & 11.90414448594906222072272411594785648949550936704795271500663482128814_prd ),& ( 7.18524510851662117145512900549280668366932915483005252325615401580719_prd, & 12.62705408396200199225079630995220617488806584833153727277631175242449_prd ),& ( 7.59381439189641479433066378186480796937054378758561448826440279962457_prd, & 13.32983471285135517036416811676054372969934205799065028150167463568813_prd ),& ( 7.99167720544561504522079912457128011957657180816119083579125507265861_prd, & 14.01455643395353374169097997561315710183878207886644322747150624317981_prd ),& ( 8.37987427934857338983294059811360130455426039105626127581824088547107_prd, & 14.68294329580687780615451873066104265705853564957606800569369647237208_prd ),& ( 8.75928129548526042699156226353757020010251201662833861031173538511019_prd, & 15.33645046170969719696146138518801429832622098668132657334102587762053_prd ),& ( 9.13064390781380346430825266662829937041678699009428928404896199007373_prd, & 15.97632038201882754684149702154344284440620400452420536075074952023654_prd ),& ( 9.49460367918717573757647043022916164116961674089630284237584546194107_prd, & 16.60362460689876200031647099108460856282007525233639325310907992053669_prd ),& ( 9.85171767110785905414362066323798900771136815515095852131622183049713_prd, & 17.21929550678326657471985928759224748750728709435044879861903190323442_prd ),& (10.20247349615153585381393106348151306984924659387653537868799157421065_prd, & 17.82415074109532567681143632039537948506451838046048127526573110371406_prd ) / data bizc(21:25) / & (10.54730106104828820929581975439272426046752560015739777627890446398335_prd, & 18.41891241378788807771014119921099424202856244191869291450744037059967_prd ),& (10.88658185276360171020167400751349315297399389674160016524066373525175_prd, & 19.00422226814517965597246929483653836545986628099335717542098425442790_prd ),& (11.22065637121487978490697787477242008976642909239915124343495431079782_prd, & 19.58065388303882246532112096834251855905705006234152403408347378498074_prd ),& (11.54983014388919569873370603964904591527739018815080309763753398354921_prd, & 20.14872256730936655594126504241030771102838585784429195801964986472933_prd ),& (11.87437864134722050579925598219404690296309944637104360607741548188520_prd, & 20.70889346473131214175517609801909093363721604975867261189535114023168_prd ) / ! ! the complex zeros of Bi'(z) in upper half complex plane data dbizc(1:10) / & ( 0.21494707453743056760883287056250605187305224890138110720402076192263_prd, & 1.10060014330279788064719413742209192726369794831163048836404502078255_prd ),& ( 1.45816830922350739202821116056490933729951849260075434667616008817736_prd, & 2.91224936745844541923508257303691682899459230703577503871294374301562_prd ),& ( 2.27376076301348229979236177702212807123422671446087248798600277076873_prd, & 4.25452854921709786216701510112105958142147276518704610066771560830131_prd ),& ( 2.96105265783122604010312815478360930889882922797417247802418355700948_prd, & 5.40812896573005040310482366111553738384473168479697682273157331450683_prd ),& ( 3.57576947588239351634900591383787344574637462648477739545407286050269_prd, & 6.44882619664101667688766917305374968697692482658875449823132377110174_prd ),& ( 4.14140155009212801911513214050461357041551762168926679027184324265196_prd, & 7.41110342737289435155363610290071245471864460472101786854840838238931_prd ),& ( 4.67067019317851843587991658526182217779452060404735436259995520458139_prd, & 8.31435754027842112221325052343261560126706456395060377346876375645051_prd ),& ( 5.17144623057547702848684803654329091601603941065180958532011916973726_prd, & 9.17087255347546656390403896670854259095596858721013207891817710397333_prd ),& ( 5.64902779025260692772106205587541053902993530776471545880519832872898_prd, & 9.98904997303294846428882639497607244557317263335698186901433176286389_prd ),& ( 6.10719645778672314912565520434373606154087248897125587950087368586402_prd, & 10.77495654079679678728726857870537776563201114764383107437207137807667_prd ) / data dbizc(11:20) / & ( 6.54877066739736118676003806957370783043553142483036682383197139988049_prd, & 11.53315414598333046084793206006732808588154690786958448199022839129188_prd ),& ( 6.97592212792084487364430796442826835395571509025540741290291051849993_prd, & 12.26718263018804966430344819184555616181575729948779604216931474582405_prd ),& ( 7.39036907848466049228107321930774095776397818857316494262577487277162_prd, & 12.97985865379078415461882350476384435335059481487043374296045769473536_prd ),& ( 7.79350057516810578045134261640077023829274762550868956525380065099566_prd, & 13.67346999233605232897265431533977335746270196819866020912761759222910_prd ),& ( 8.18645982214873789501772438998924405082262072688599000055406755107824_prd, & 14.34990697064684323397347337125976171582065678522309064630023625752993_prd ),& ( 8.57020199560005738297958938286472624011678392077178524442698505383515_prd, & 15.01075434722730594102247347049141443546551439677981436949704284282558_prd ),& ( 8.94553554174992509923121042304493111410594720693234230010453783795201_prd, & 15.65735735534456505989719663535512163903202778433978793833497105463328_prd ),& ( 9.31315240559377120279006194179236218678382103286005699284034696552073_prd, & 16.29087030713368404673157848639547747437192335235981600587878627241399_prd ),& ( 9.67365063100573605685580820643318437333429764070748529694125424618393_prd, & 16.91229310481033620812543730559457811485843222338184756713409814473519_prd ),& (10.02755157233877751003651949722784369505767870153629086749814687666235_prd, & 17.52249916285219769188573833900744788846491409677167422613709212335609_prd ) / data dbizc(21:25) / & (10.37531321693404033189804291399040740722048977783034198713300877988027_prd, & 18.12225710100791573570080366782397788632586761581476795682060835313086_prd ),& (10.71734064688710698783676365059067829058738257679452688582483678994761_prd, & 18.71224783545241969047143044519901088119172969586176757125107451219527_prd ),& (11.05399436066789932335321274509712390948071324132825401515346039829534_prd, & 19.29307821395942947173251780131670821566064914355268076696448955332805_prd ),& (11.38559696930207684941841365020489192373274968894013891384043192965540_prd, & 19.86529201713763010184485681506434978210283080939980848299152471421751_prd ),& (11.71243864112622080163738620993012206688392259049239134593136924664076_prd, & 20.42937892542195221143659467613817502213630344239089786009258668146871_prd ) / ! ! what follows is the hard coded coefficients of the ! asymptotic series for the phase function theta(x) of the ! Airy functions. data thetacoeff(0:19) / & 1.0000000000000000000000000000000000000000000000000000000000000000000e+00_prd, & 1.5625000000000000000000000000000000000000000000000000000000000000000e-01_prd, & 1.7985026041666666666666666666666666666666666666666666666666666666670e-01_prd, & 1.2638092041015625000000000000000000000000000000000000000000000000000e+00_prd, & 2.1832866753850664411272321428571428571428571428571428571428571428570e+01_prd, & 6.9682822536884082688225640190972222222222222222222222222222222222220e+02_prd, & 3.5553981581391285131262107328935102982954545454545454545454545454550e+04_prd, & 2.6532584171176254655550180289607781630295973557692307692307692307690e+06_prd, & 2.7258321302415599342250374093055143021047115325927734375000000000000e+08_prd, & 3.6894027011098001901607477714297461341843967709471197689280790441180e+10_prd, & 6.3630204766133750347903292605201682576982217429567275470808932655740e+12_prd, & 1.3622439180500760352059703427764362988546306819701147192279763874550e+15_prd, & 3.5446828581078074793552458265811520091905839784775414983950826274890e+17_prd, & 1.1017990287911060069986723262330428844182307035279722983583949869060e+20_prd, & 4.0320823577747335533970347823492793114462319968284659862713417937250e+22_prd, & 1.7159646692716793386038283530387825075493716774249119449042189710410e+25_prd, & 8.4031035110193005527486652122292612450561213930864781736528036247480e+27_prd, & 4.6916652616593043608874014133910554079508141416060179915669243125790e+30_prd, & 2.9625489083114947131526274616960212993262723579957461622065841493850e+33_prd, & 2.1006649445213587744194779196040385627926101231259943031483704675450e+36_prd / data thetacoeff(20:39) / 20*zero / ! ! the second 50 coeffs are for phi(x). data phicoeff(0:19) / & 1.0000000000000000000000000000000000000000000000000000000000000000000e+00_prd, & -2.1875000000000000000000000000000000000000000000000000000000000000000e-01_prd, & -2.3811848958333333333333333333333333333333333333333333333333333333330e-01_prd, & -1.5114471435546875000000000000000000000000000000000000000000000000000e+00_prd, & -2.4620345711708068847656250000000000000000000000000000000000000000000e+01_prd, & -7.6086294204038050439622667100694444444444444444444444444444444444440e+02_prd, & -3.8092128868188640229742635380138050426136363636363636363636363636360e+04_prd, & -2.8082562244080195275903231679246975825383112980769230769230769230770e+06_prd, & -2.8607946453262385636628692964222864247858524322509765625000000000000e+08_prd, & -3.8480790905747179385766298433031973967017085455796297858743106617650e+10_prd, & -6.6050619870875374085141895250268697180548303405980971690855528178970e+12_prd, & -1.4087107203786266895474229097071387064564936449547379500775908430420e+15_prd, & -3.6542984737250173263822582147634056560257268026887776125454859973850e+17_prd, & -1.1329639450011994200743048058647230246736432348623124181236088148240e+20_prd, & -4.1371744404875883016541729567285839631712496380598048629271123618470e+22_prd, & -1.7574349642737663059031534378525691263268308654379228784525651291710e+25_prd, & -8.5924160506996514557379617019132398803302175426817653569682851557100e+27_prd, & -4.7906551025465751515669729015000637241952562513029744187686495732900e+30_prd, & -3.0213269421511401499687820977552075015245020527256257041378966249360e+33_prd, & -2.1399970937320017557354494737300168515390960226867407271664411223300e+36_prd / data phicoeff(20:39) / 20*zero / ! ! the following lines contain data statement for the coefficients of the asymptotic ! expansions for the zeros and associated function values needed to produce ! accurate answer for precisions greater than quad, i.e. 2^{-112}. enough ! coefficients are stored so to ensure octal precision, i.e. 2^{-224}. ! data Ucoeff(19:39) / & ! -.14107259800787276308427861823522682542378736666070937012298625588832e37_prd,& ! .11157559226376058388241181162252207780072675866049674911684381382474e40_prd,& ! -.97911168966451866992095894172730381052056195537079932657057595139361e42_prd,& ! .94842188231979773651687611788221221444304429864017500587580906466991e45_prd,& ! -.10093854348840004363421701540951189641241176805162426458768970023668e49_prd,& ! .11753309898206889209788907835323816047837770192802898892767943766799e52_prd,& ! -.14915259567389695154010536894836306258744098712656167509641682302319e55_prd,& ! .20555467603941926152452307773581023553674196606652536651923229197765e58_prd,& ! -.30664114828091931546036143434946049203640895201650759126568923937596e61_prd,& ! .49366102338787427090349132178175476785656098771026760487567560730389e64_prd,& ! -.85527872283563506905396193884078691515189421163063772903051195883576e67_prd,& ! .15905160493853017875103784933898084128023076207249471188546554365590e71_prd,& ! -.31671480020251898031362069511376069683729168870885141394131419726417e74_prd,& ! .67377747797473203058970344907089937116595231663364594832661747001111e77_prd,& ! -.15281403392070465160928233963603155027184060115379663405333917341838e81_prd,& ! .36876222351784858824732282512073352215729923468181255970009084742092e84_prd,& ! -.94505207024597841792499304832642666053323686278446400296025689886898e87_prd,& ! .25676049507927905036477828913103335013957477586456812093463169718887e91_prd,& ! -.73831918349991133688096513080694039594126308972102658977041222126047e94_prd,& ! .22434873650609465989916612597264280685508533191672309973604370953703e98_prd,& ! -.71932175167133289253941392314260564935469329383000953020665385234810e101_prd/ ! data Ucoeff(40:55) / & ! .24301499022274266015249236229059102519443259535429498372071300028184e105_prd,& ! -.86392200099184256673341753224207184715009085434527689795095151266352e108_prd,& ! .32277334408133278730034141082298966934072100188969252542528864933608e112_prd,& ! -.12658410453063860217634140848993855557208179826387248492220322139142e116_prd,& ! .52049980681802453665831970426777630104007284760271492623690409070701e119_prd,& ! -.22415396118038146187210801469410795705798174092757295000371600948203e123_prd,& ! .10099567175653350699936345320055405124203719932677954083892691964047e127_prd,& ! -.47561519985327666244537756889291936086598881976270611582042240452110e130_prd,& ! .23387844301594890965483902208040384641369089913587515033974449393995e134_prd,& ! -.11997993650041033217582177048624365958249639549393416213721216792541e138_prd,& ! .64154920040257821156504337755375966474958297697150253844482194020056e141_prd,& ! -.35726348720185093488728179482093361120107194523287613890849925654197e145_prd,& ! .20703014851226170160468376381877508456518993988984763894602747782221e149_prd,& ! -.12474625649281036421622524359179945285093096496161906265837264044813e153_prd,& ! .78099095291498984417321166641290057395110783507772223491883388065302e156_prd,& ! -.50766451666920499286704685800242398987056663490773056939748373429556e160_prd/ ! data Tcoeff(19:39) / & ! .14117258064737679737776630556275250408926431644493704674703400857145e37_prd,& ! -.11164642174072109561223350471691552374391433806310651183802082896951e40_prd,& ! .97967164158278846547438620160557824503041359438794841332630872374930e42_prd,& ! -.94891307580284755277548004980034901717550001626170852847573583227890e45_prd,& ! .10098610809248076249074036492916612325284053264374591096020664351942e49_prd,& ! -.11758370710196329175686452122768747739546685071872150784191369967919e52_prd,& ! .14921150984686502643918601323682722977762625135098393735932060011458e55_prd,& ! -.20562942407939128973788970107145954697239949579275764596681522368340e58_prd,& ! .30674414339703027984557756338766648623061371604751170077497853041336e61_prd,& ! -.49381464389065703514414376124090604029319195122151949195423702788440e64_prd,& ! .85552600050674514166037335169011816701693895322935495541549028524886e67_prd,& ! -.15909444082473299698105131270096108507658923102224119263418369560478e71_prd,& ! .31679445060129872330154447996576899972397808928310935503731008478804e74_prd,& ! -.67393606698205551356590737775952562838503676308476217226185171949385e77_prd,& ! .15284776889635128629958626667274826956672953066901360901017230535376e81_prd,& ! -.36883872857780404454641197328305961049114361355619869034021470060729e84_prd,& ! .94523667365207400774723728448041810267457304315127521834562030560915e87_prd,& ! -.25680780131046747934900503749660107326570832650658319996932836021029e91_prd,& ! .73844770114739079149892403045167309500572883020991309765436811625919e84_prd,& ! -.22438568978657876605574244940414470180612253866038210914591946322219e98_prd,& ! .71943403333451861202742470176857110279884665218133554269689930563214e101_prd / ! data Tcoeff(40:55) / & ! -.24305098890695512130164302767376074435501301582956470312823055101950e105_prd,& ! .86404361298136368246708141865214768397985490481404698406775721925269e108_prd,& ! -.32281657547773530130641638261260703492461229074644634987917760718557e112_prd,& ! .12660025574379689825390757955090645221822260387369027230864689425729e116_prd,& ! -.52056314573696814753098687759946586496965657569199355992702143018080e119_prd,& ! .22418000454826808504041968248387029081392043348043423531235015542335e123_prd,& ! -.10100688701029413642642584650224102377826199264338678362250408468416e127_prd,& ! .47566573022411084112382308739209248802018314978025349360661979655957e130_prd,& ! -.23390223843342824573557767032239217008469637672278461668515561594488e134_prd,& ! .11999163732445825923983675048550329573712483824248289866948173690662e138_prd,& ! -.64160922430691576929378390179082298613977749760199117434017886679318e141_prd,& ! .35729558206957913592353331316564747010295288658097259740502832900235e145_prd,& ! -.20704802099398561273399136516817367415999088767761277008018247421227e149_prd,& ! .12475661318287886645565405796047851203309187211403798809514130671614e153_prd,& ! -.78105335607384790780096263667849118744101236336911720704817696289821e156_prd,& ! .50770358431057167077679014158662551964839528596933932658427521736419e160_prd / ! data Wcoeffs(19:39) / & ! -.38814345591742149734668498349356519681026833217761869126091260492675e38_prd,& ! .32371427083820730592997422984643403177773522701009254386004811466743e41_prd,& ! -.29874986827080478884190074254990568828103267869569100228461781316986e44_prd,& ! .30360634740053951327896972827621158563037858664706363711081949353786e47_prd,& ! -.33825714675286665024140203614169399088635995626398761088900799825628e50_prd,& ! .41149164273551406585781481488408987782456057258383260258549416260778e53_prd,& ! -.54455986340331835374641221522823255837698570580066433295851004024933e56_prd,& ! .78130992767230544018414009790655958914552920815980627216875258213508e59_prd,& ! -.12115223587091684946329431505204521385450149548831782255139272257958e63_prd,& ! .20244592837085139181461609232296606192823980799932370000660702095961e66_prd,& ! -.36356845070938164932191419587043625426144667669878108642458711480165e69_prd,& ! .69996165929353579451877926547331080690884165497387913985052720353445e72_prd,& ! -.14413113297034043325245482652689576895508908029906780419542674771067e76_prd,& ! .31672871590345335266091698496695051388432534498125701861139662564097e79_prd,& ! -.74126513717880415662859852045921163300483589588651650409246140259593e82_prd,& ! .18440849873059519808034797163033300858530685710495395519130980839380e86_prd,& ! -.48676991919685415796644648701954794143568660492618142404503177084697e89_prd,& ! .13610103187692986437476716347372662183718353515873683134697916324211e93_prd,& ! -.40243417887979252089715153133466607795716506887314640636278752867259e96_prd,& ! .12565013777935480241848586316006283152879014921176384532124033031005e100_prd,& ! -.41365634264582272778566733965611481020861148242111761647744907714320e103_prd/ ! data Wcoeffs(40:59) / & ! .14339409364180527062863367710702909293538014155644831937447023847458e107_prd,& ! -.52272565695840736902310970199322716034896043916911897826711471370477e110_prd,& ! .20013873237902089546233112757343855340672136637249651667006730957441e114_prd,& ! -.80388278212226250229539050402383177921686252031963455625680997914909e117_prd,& ! .33835447680016009318131920319518924659697918603232614321420856900811e121_prd,& ! -.14907484080694089485057416124301611758212349344388791192461825578970e125_prd,& ! .68682543761177567812666841466566465804229617540851105036999084022986e128_prd,& ! -.33057783816698691465421284889296141763433202701355488119536813033424e132_prd,& ! .16606585680890010392955925191610451146159841700201344759512414035539e136_prd,& ! -.86991562451838247539863941467506400498257953847383079314747693710280e139_prd,& ! .47477840071035479070316437198311566899770847203711441099978992848440e143_prd,& ! -.26975139028563609163677868450212082245025416175786556168506457211935e147_prd,& ! .15942313127360607223331447173078019238938649855332049699892015189284e151_prd,& ! -.97931671254977175938624414519464954227804198086876729834671618925574e154_prd,& ! .62482875317012139806026314897328136965964374299906322451833626941862e158_prd,& ! -.41376954055676364528233672305143910272597687609280824732944653664345e162_prd,& ! .28419769324142769069351633035666203091678626692459848096036450446153e166_prd,& ! -.20233105203250752333992061204250782338499709868648694329341733953856e170_prd,& ! .14921413227682859768600979116852113769192082397308681022730906838764e174_prd,& ! -.11391942318647607278195227387147815982946299255584776067390546283266e178_prd/ ! data Wcoeffs(60:62) / & ! .89984973907118798973274600559964730505886269200285584307716001979151e181_prd,& ! -.73498618581016154263237896162942608870795994021129309881094563967910e185_prd,& ! .62042034816474031028845699860660079271705962703651243271640215253005e189_prd/ ! data Vcoeff(19:39) / & ! .39511749687738754978707186616441860115896523182485334297380695925561e38_prd,& ! -.32923290247609366681382018551721673465195880594473426214142086973247e41_prd,& ! .30359489943433278235144569435560837233966593957760177297566358869186e44_prd,& ! -.30830152332406270164377583091160563455014524379537218769617073966086e47_prd,& ! .34325609828559626619691428942715725232658096603901061837756302968756e50_prd,& ! -.41731456300646629062922610457970573676170714541136141310758967872890e53_prd,& ! .55195180893208079575182955709289525687994646475016335125669817883713e56_prd,& ! -.79150035232067397162841508858543292745122889767009874931157993833688e59_prd,& ! .12267286410014702060147710311493346807399372162743542592030229492395e63_prd,& ! -.20489465757544862177181168685362400780928284026545992418550424591155e66_prd,& ! .36781201995018744410586705497584162101937619239291503710300378853588e69_prd,& ! -.70785510224402391504410638973731615282432843867153838660600966190695e72_prd,& ! .14570328765817175201745344397239061049289639510667166652366926788623e76_prd,& ! -.32007402423228767995646231438096527047320321247086611792310326968654e79_prd,& ! .74885388950666402410596485880458556983033741343261318201018549562040e82_prd,& ! -.18624012221973572482406403572483861757301929301481060378790454939360e86_prd,& ! .49146479067201303283191057797168315897937697549653670834510914285792e89_prd,& ! -.13737679563804852903349620097905514534829657189061129536765282302984e93_prd,& ! .40610325317648937486635390453257288210644048252427882240862813659530e96_prd,& ! -.12676520941322158498322311515387022465292939801517893359670714044848e100_prd,& ! .41723208208662814474860048359114616247849879135056450584779419025583e103_prd/ ! data Vcoeff(40:59) / & ! -.14460228643992977278971715293625872465630777893677611595013226813385e107_prd,& ! .52702137141948184886975471002292506399568818205669740583669405862891e110_prd,& ! -.20174387491029824961205653853719597884391823201956355101548772468333e114_prd,& ! .81017853587955248434572097165582819569986494473379223283336066657050e117_prd,& ! -.34094352097537082769224279460429788581422930094799330429529116973719e121_prd,& ! .15018993996821870455724263440162625429659710408567120834504368444388e125_prd,& ! -.69185020550260382163148603962436395586032754664208939252674040930507e128_prd,& ! .33294437491009899625612070497418754367170356692121410097534665436023e132_prd,& ! -.16722968863179740044484559468409527225906823206213932532120401016256e136_prd,& ! .87588665788890625267255395566168013064348940579513872897521595626651e139_prd,& ! -.47797148292096096604328686586996148973995999656389939085032677629810e143_prd,& ! .27152969647873325149100043985011418074992909890944936124651916757659e147_prd,& ! -.16045372590074827247546573056862140802796751223989742959748936578369e151_prd,& ! .98552707357503715509221245702668173028253308928381456246148353572071e154_prd,& ! -.62871713691259881989197991638799674760450888458950498498758525415788e158_prd,& ! .41629728127968126225313221039832645027885966270218693567128257568167e162_prd,& ! -.28590262269286840162369262442491958517609091260359713489594756877926e166_prd,& ! .20352339449544412644852790819518468781704630082369339335234743302728e170_prd,& ! -.15007817824859664648921028573280959268038609339370901503470055457741e174_prd,& ! .11456782500469922643539779328058943363839320875262498410359460045474e178_prd/ ! data Vcoeff(60:62) / & ! -.90488547197595729727818007599162028667196587795950750687627446320423e181_prd,& ! .73903139126156436722227750117680141015843739487204686022681600834281e185_prd,& ! -.62377953850258650978985595488693003161283350731492126115310176459383e189_prd/ ! data thetacoeff(20:39) / & ! 1.6620414565199836046517241168226679575433914398847140159409431569400e+39_prd, & ! 1.4589768148863538676908338329952770803471443409116858702096237708490e+42_prd, & ! 1.4136706205926967227550391181089542632840298521113599879861436872780e+45_prd, & ! 1.5049528873233485453201090242148583430669331596967888813309286983220e+48_prd, & ! 1.7528120867435543022238293679384047783034582054789241354396665814190e+51_prd, & ! 2.2248798784215759381108095700088082121325362980595660350417563574900e+54_prd, & ! 3.0668742919686650506786294155949844120360694119732382794844730546530e+57_prd, & ! 4.5759897594438856162245163173704124179167067431675218696740068750570e+60_prd, & ! 7.3682324048945997036117127108954093722721949946103451048231900538500e+63_prd, & ! 1.2767812907658639911193679647694549374414268530237897473653855408160e+67_prd, & ! 2.3747413389528063711999564032375036732632838473204649286532194299650e+70_prd, & ! 4.7294599973899606439620494750909194373558101517093261954177214594550e+73_prd, & ! 1.0062840158864202353776785670460809939692662665374140432947800580520e+77_prd, & ! 2.2825722023880581733954080677391959965131449243240891757158807086650e+80_prd, & ! 5.5088575969201769653275088819145994402997486982906377742850251255490e+83_prd, & ! 1.4119574904108397754635319483178352586648916671672757537207984396580e+87_prd, & ! 3.8365603338945820651443721086714857583386907373667502273506193817610e+90_prd, & ! 1.1033246832978464491034681560090717113928219873397757024970797143410e+94_prd, & ! 3.3529399615471696053962311884989125864086206752781566570555462826900e+97_prd, & ! 1.0751426905865614366516033970490417988024797993593296124425204481830e+101_prd / ! data thetacoeff(40:49) / & ! 3.6325759548575941054331531628741592937638711725735797561598681234640e+104_prd, & ! 1.2914958421278365429073565165388490025468863490549642689956743754280e+108_prd, & ! 4.8255989348832298734268326774116915672615870738562487397862875989870e+111_prd, & ! 1.8926320413720882742332339252335781415823127647556584051535791102430e+115_prd, & ! 7.7828653867243241612219242316942961428397330907284697950596299672640e+118_prd, & ! 3.3519372264089624493413195232679066342813817418642368984048705399450e+122_prd, & ! 1.5103632391716812602137633280398852067548262843999179493462357088140e+126_prd, & ! 7.1131556794721283742180920205235953698098271487901179565894695165920e+129_prd, & ! 3.4980302564753714551386271025304835533159465936045523336850630596900e+133_prd, & ! 1.7946000747753342325010185977181141233793656634866192538659890317680e+137_prd / ! data phicoeff(20:39) / & !-1.6915031381769804103850733555848209612224387842573724614493225104370e+39_prd, & !-1.4835311087410397492083286676216470036366696291606158301927652269280e+42_prd, & !-1.4363173833472507948203636105983791470914532060599594653238608453020e+45_prd, & !-1.5279550958145288037798152480016932257564787055599234667936240458810e+48_prd, & !-1.7784267946552010527555830684846510215218772178165602333517334638230e+51_prd, & !-2.2560262434454646323482238905021187368735965582739309243375596398990e+54_prd, & !-3.1080758268011686600204807563047338098476704125999060010703485766130e+57_prd, & !-4.6350816361828754091843816015741307407256970145816581103438065469270e+60_prd, & !-7.4598303934675995061007145208437781401846859095965526853293166387510e+63_prd, & !-1.2920824883839973738038478102778590132033803794125290485499782279180e+67_prd, & !-2.4022123670554103535271126405427971447224944735511581081191319646570e+70_prd, & !-4.7823342996246230871145477576426171365849607028316240761103920136300e+73_prd, & !-1.0171687597738396277948309273735238244126224786530205309841241554100e+77_prd, & !-2.3064858246307166204738030138103766616236470199051511300921739815190e+80_prd, & !-5.5648121928700061582464391625164386860081156002157737088990604346560e+83_prd, & !-1.4258747425358444009570638240678952414106824061428143662684567461640e+87_prd, & !-3.8732895967824595508759346648192032120618162030750584346035466479860e+90_prd, & !-1.1135923286968508441850331229673794325234269780229843580847122796070e+94_prd, & !-3.3832949359638733345296896827515103441106509123369820515812386269470e+97_prd, & !-1.0846187498103346877861073639932215030464972161028749862302896323670e+101_prd / ! data phicoeff(40:49) / & !-3.6637675896351569613989314774169363221116993261202308931493480926110e+104_prd, & !-1.3023068683764358850217284313843312241882676687996376761813939376260e+108_prd, & !-4.8650038107772889723933269035595710492136466820336895421677762354480e+111_prd, & !-1.9077172322290849802626455124991995446193522466030919121989161940300e+115_prd, & !-7.8434494940022838542204878691135615680589924478607909220870915929650e+118_prd, & !-3.3774340943788767024750346018657171672375785744278296573692058467690e+122_prd, & !-1.5215955884734982790544475883850926084837078168255872870631652278060e+126_prd, & !-7.1649004190192480023524589812450546706749489342612772092982680431290e+129_prd, & !-3.5229331179784842655076831037707769835979262101176203206020619820250e+133_prd, & !-1.8071088131997102426584141687123445370794488782743819670027768650070e+137_prd / ! ! ! end file airy_head SHAR_EOF fi # end of overwriting check if test -f 'airy_parameters' then echo shar: will not over-write existing file "'airy_parameters'" else cat << "SHAR_EOF" > 'airy_parameters' ! ! begin file airy_functions_parameters ! version 1.0 ! edited 12/28/2001 ! subroutine airy_paramsr (Radius, max_step_size, & n_terms_taylor, n_terms_asymp, n_partitions) real(prd), intent (out) :: Radius real(prd), intent (out), optional :: max_step_size integer, intent (out), optional :: n_terms_taylor, & n_terms_asymp, n_partitions !**************************************************************************! ! This subroutine returns to the user the parameters used specific to the ! machine and compiler used !**************************************************************************! if (.not. is_init_airy) call parameter_airy Radius = r_min if (present(max_step_size)) max_step_size = fstpsz if (present(n_terms_taylor)) n_terms_taylor = n_taylor if (present(n_terms_asymp)) n_terms_asymp = n_asymp if (present(n_partitions)) n_partitions = n_parts end subroutine airy_paramsr ! subroutine parameter_airy !***********************************************************************! ! This subroutine computes and stores the global variables ! r_lolimit ! r_uplimit ! r_min (rho) ! n_asymp (S) ! n_taylor (N) ! n_parts (P) ! fstpsz (h) ! and the asymptotic expansion coefficients ! ucoef (Eq. 2.6) ! vcoef (Eq. 2.6) ! Parenthetical notes refer to the ACM TOMS paper (submitted?). !***********************************************************************! integer M, i real(prd) eta, Sfloat, xinew, xiold, chi, fun, dfun, basep real(prd), dimension(0:2) :: lambda ! ! required machine constants M = digits(one) basep = real(radix(one),prd) eta = epsilon(one)/basep ! ! three global variables: the upper and lower limits on the ! computation region and the number of terms in the ! asymptotic expansions for large |z| r_lolimit = epsilon(eta)/abs(ai1zer) r_uplimit = exp(two_third*log(huge(eta)*0.95_prd)) n_asymp = M*log(basep) + half ! ! generate the coefficients for the asymptotic expansions for large |z| allocate ( ucoef(n_asymp) ) allocate ( vcoef(n_asymp) ) ucoef(1) = 5.0_prd/72.0_prd vcoef(1) = -7.0_prd/72.0_prd do i = 1, n_asymp-1 ucoef(i+1) = real((6*i+5)*(6*i+1),prd)/real((i+1)*72,prd)*ucoef(i) vcoef(i+1) = -real(6*i+7,prd)/real(6*i+5,prd)*ucoef(i+1) end do ! ! determine the parameter r_min (rho) using Newton's method chi = one if ( mod(n_asymp,2) == 0 ) then do i = 1,n_asymp/2 chi = chi*i/(i-half) end do else chi = pi*half do i = 1,(n_asymp-1)/2 chi = chi*(i+half)/i end do end if Sfloat = real(n_asymp,prd) xinew = half*Sfloat chi = M*log(basep) + log(abs(vcoef(n_asymp))) + log(two*chi) ! dummy fstpsz = 7.0_prd*pi/72.0_prd ! dummy do xiold = xinew fun = Sfloat*log(xiold) - fstpsz/xiold - chi dfun = Sfloat/xiold + fstpsz/xiold**2 xinew = xiold - fun/dfun if ( abs((xinew-xiold)/xinew) < four*eta ) exit end do ! ! three of the global variables can now be computed r_min = (1.50_prd*xinew)**two_third n_parts = ceiling(r_min) fstpsz = r_min/n_parts ! ! the number of terms in the taylor series lambda(0) = one/fstpsz lambda(1) = one lambda(2) = r_min*half*fstpsz i = 2 do i = i+1 lambda(mod(i,3)) = (r_min*lambda(mod(i-2,3)) + & fstpsz*lambda(mod(i-3,3)))/real(i*(i-1),prd)*fstpsz**2 if ( i*lambda(mod(i,3)) < half*eta ) exit end do ! ! the final global variable n_taylor = i ! ! reset the logical and exit is_init_airy = .true. ! ! the following are dummy statements to use variables in the ! header file that pertain to the complex routines. arg_local = .false. allocate(aigrid(1,1)) deallocate(aigrid) z_global = dbizc(1) zeta_global = bizc(1) theta_global = zero ! ! the following are dummy statements to use variables in the ! header file that pertain to the real routines. big_integrate_aux = .false. allocate(m2grid(1,1)) deallocate(m2grid) allocate(theta_grid(1)) deallocate(theta_grid) x_global = dbizr(1) xi_global = bizr(1) x_global = daizr(1) xi_global = aizr(1) ! end subroutine parameter_airy ! ! end file airy_functions_parameters ! ! begin file airy_aux_parameters ! version 1.0 ! edited 01/07/2002 ! !*** !************************************************************************ !*** subroutine airy_paramsr_aux (Radius, max_step_size, & n_terms_taylor, n_terms_asymp_mod, n_terms_asymp_phase, n_partitions) real(prd), intent (out) :: Radius real(prd), intent (out), optional :: max_step_size integer, intent (out), optional :: n_terms_taylor, & n_terms_asymp_mod, n_terms_asymp_phase, n_partitions !* ! This PUBLIC subroutine returns to the user the parameters used for the ! modulus and phase functions specific to the ! machine and compiler used. !* if (.not. is_aux_init_airy) call aux_parameter_airy Radius = aux_min if (present(max_step_size)) max_step_size = mstpsz if (present(n_terms_taylor)) n_terms_taylor = n_taylor_aux if (present(n_terms_asymp_mod)) n_terms_asymp_mod = n_asymp_mod if (present(n_terms_asymp_phase)) n_terms_asymp_phase = n_asymp_phase if (present(n_partitions)) n_partitions = n_parts_aux end subroutine airy_paramsr_aux !*** !************************************************************************ !*** subroutine aux_parameter_airy !* ! This PRIVATE subroutine computes and stores the global variables ! aux_min (hat rho) ! n_asymp_phase (hat T) ! n_asymp_mod (hat S) ! n_taylor_aux (hat N) ! n_parts_aux (hat P) ! mstpsz (hat h) ! and the asymptotic expansion coefficients ! mcoef (Eq. 7.7) ! ncoef (Eq. 7.8) ! Parenthetical notes refer to the ACM TOMS paper (submitted?). !* real(prd), dimension(0:2) :: lambda real(prd) h, eta integer i, M ! ! determine the compiler dependent constants if (.not. is_init_airy) call parameter_airy eta = epsilon(one)/real(radix(one),prd) M = digits(one) ! ! determine the number of terms needed in the asymptotic ! expansions for the modulus functions n_asymp_mod = half*M*log(two) - fourth*log(one*M) + half allocate ( mcoef(n_asymp_mod) ) allocate ( ncoef(n_asymp_mod) ) mcoef(1) = 5.0_prd/32.0_prd ncoef(1) = -1.40_prd*mcoef(1) do i = 1, n_asymp_mod-1 mcoef(i+1) = real((6*i+5)*(6*i+3)*(6*i+1),prd)/real((i+1)*96,prd)& *mcoef(i) ncoef(i+1) = -real(6*i+7,prd)/real(6*i+5,prd)*mcoef(i+1) end do ! ! we now determine the parameter hat rho aux_min = abs(ncoef(n_asymp_mod-1))/eta aux_min = aux_min**(one/real(3*(n_asymp_mod-1),prd)) ! ! determine the number of terms needed in asymptotic expansion ! of the phase functions i = 1 do if ( thetacoeff(i) /= zero .and. i /= n_coeff_phase ) then if (two*max(abs(thetacoeff(i)),abs(phicoeff(i)))/ & aux_min**(3*i)/eta < one) exit i = i+1 else aux_min = aux_min + one i = 1 end if end do n_asymp_phase = i ! ! the number of terms in the taylor series h = aux_min/ceiling(aux_min) lambda(0) = one/h**2 lambda(1) = one/h lambda(2) = half i = 3 do lambda(mod(i,3)) = h**4*(four*aux_min*(mod(i-2,3))*lambda(mod(i-2,3)) & + lambda(mod(i-3,3))*(4*i-10)*h)/real(i*(i-1)*(i-2),prd) if (two*i*(i-1)*lambda(mod(i,3))/eta < one) exit i = i+1 end do ! ! the machine specific constants are determined. n_taylor_aux = i n_parts_aux = ceiling(aux_min) mstpsz = h ! ! reset the logical and exit the routine is_aux_init_airy = .true. end subroutine aux_parameter_airy !*** !************************************************************************ !*** ! ! end file airy_aux_parameters ! ! begin file airy_zeros_parameters ! version 1.0 ! edited 01/03/2002 ! !*** !************************************************************************ !*** subroutine zero_parameter_airy !* ! PRIVATE subroutine which determines the number of terms to use ! in the asymptotic expansion for the zeros !* integer NN real(prd) tol, term1, term2 ! if (.not. is_init_airy) call parameter_airy tol = epsilon(tol)/two do NN = 1,n_coeff_zero term1 = Ucoeff(NN)/(three_pi_ate*real(4*n_zeros-1,prd))**(2*NN) term2 = Tcoeff(NN)/(three_pi_ate*real(4*n_zeros-3,prd))**(2*NN) if ( abs(two*max(term1,term2)/tol) < one ) then n_asymp_zero = NN exit end if end do do NN = 1,n_coeff_asso term1 = Vcoeff(NN)/(three_pi_ate*real(4*n_zeros-1,prd))**(2*NN) term2 = Wcoeff(NN)/(three_pi_ate*real(4*n_zeros-3,prd))**(2*NN) if ( abs(two*max(term1,term2)/tol) < one ) then n_asymp_asso = NN exit end if end do ! ! reset the logical and exit the routine is_zero_init_airy = .true. end subroutine zero_parameter_airy !*** !************************************************************************ !*** ! ! end file airy_zeros_parameters SHAR_EOF fi # end of overwriting check if test -f 'airy_real' then echo shar: will not over-write existing file "'airy_real'" else cat << "SHAR_EOF" > 'airy_real' ! ! begin file airy_functions_real ! version 1.0 ! edited 07/24/2003 ! !*** !************************************************************************ !*** subroutine airy_air (x, ai, dai, ierr, modify_switch) real(prd), intent(out) :: ai, dai integer, intent(out) :: ierr real(prd), intent(in) :: x logical, intent(in), optional :: modify_switch !* ! PUBLIC subroutine which is the driver routine that computes Ai(x) ! and its derivative. The (PRIVATE) subroutines called are: ! parameter_airy = determines the machine specific parameters ! flows = checks for over- and under- flow regions ! asymp1r = asymptotic expansion for x > 0 ! asymp2r = asymptotic expansion for x < 0 ! taylorr = taylor series expansion (integration) !* integer nn, iregion real(prd) Pxi, Qxi, Rxi, Sxi, ai0s, ai1s, a, b, xi0 real(prd), dimension(2,2) :: tsm ! ! check for some quick under- and over- flow regions if (.not. is_init_airy) call parameter_airy iflag = 0 mod_local = .false. if (present(modify_switch)) mod_local = modify_switch r_global = abs(x) x_global = x if (x_global < zero .and. mod_local) ierr = -1 call flows if (iflag == 0) then xi_global = two_third*r_global*sqrt(r_global) ! ! choose the region where x lies and call the appropriate subroutine ! to evaluate Ai(x) and/or Ai'(x) iregion = 0 if (r_global >= r_min) then if (x_global > zero) then iregion = 1 else iregion = 2 end if else if (x_global > zero) then iregion = 3 else iregion = 4 end if end if ! ! call the appropriate subroutine(s) to evaluate Ai(x) and/or Ai'(x) select case (iregion) case(1) ai1s = two*sqrpi ! dummy storage of number ai0s = x_global**fourth ! dummy storage of number a = one/(ai1s*ai0s) b = -ai0s/ai1s iflag = 10 ! indicates Ai functions to asymp1r call asymp1r (ai0s, ai1s) ai0s = a*ai0s ai1s = b*ai1s if (.not. mod_local) then ai0s = exp(-xi_global)*ai0s ai1s = exp(-xi_global)*ai1s end if case(2) ai1s = r_global**fourth ! dummy storage of number a = one/(sqrpi*ai1s) b = ai1s/sqrpi call asymp2r (Pxi, Qxi, Rxi, Sxi) ai1s = xi_global-pi*fourth ! dummy storage of number ai0s = a*(Pxi*cos(ai1s) + Qxi*sin(ai1s)) ai1s = b*(Rxi*sin(ai1s) - Sxi*cos(ai1s)) case(3) nn = ceiling(n_parts*(one-r_global/r_min)) ai1s = two*sqrpi ! dummy storage of number ai0s = r_min**fourth ! dummy storage of number a = one/(ai1s*ai0s) b = -ai0s/ai1s xi0 = two_third*r_min*sqrt(r_min) iflag = 10 ! indicates Ai function to subroutine asymp1r call asymp1r (ai0s, ai1s, xi0) call taylorr (nn, r_min, tsm) a = a*exp(-xi0)*ai0s b = b*exp(-xi0)*ai1s ai0s = tsm(1,1)*a + tsm(1,2)*b ai1s = tsm(2,1)*a + tsm(2,2)*b if (mod_local) then ai0s = exp(xi_global)*ai0s ai1s = exp(xi_global)*ai1s end if case(4) nn = ceiling(n_parts*r_global/r_min) call taylorr (nn, zero, tsm) ai0s = tsm(1,1)*ai0zer + tsm(1,2)*ai1zer ai1s = tsm(2,1)*ai0zer + tsm(2,2)*ai1zer end select end if ! ! error handling and return values select case (iflag) case(:-1) ai = zero dai = zero ierr = iflag case(0:2) ai = ai0s dai = ai1s ierr = 0 case(5) ai = ai0zer ! the value of Ai(0) dai = ai1zer ! the value of d(Ai(0))/dx ierr = 0 end select end subroutine airy_air !*** !************************************************************************ !*** subroutine airy_bir (x, bi, dbi, ierr, modify_switch) real(prd), intent(out) :: bi, dbi integer, intent(out) :: ierr real(prd), intent(in) :: x logical, intent(in), optional :: modify_switch !* ! PUBLIC subroutine which is the driver routine that computes Bi(x) and ! its derivative. The (PRIVATE) subroutines called are: ! parameter_airy = determines the machine specific parameters ! flows = checks for over- and under- flow regions ! asymp1r = asymptotic expansion for x > 0 ! asymp2r = asymptotic expansion for x < 0 ! taylorr = taylor series expansion (integration) !* integer nn, iregion real(prd) Pxi, Qxi, Rxi, Sxi, bi0s, bi1s, a, b real(prd), dimension(2,2) :: tsm ! ! check for some quick under- and over- flow regions if (.not. is_init_airy) call parameter_airy iflag = 0 mod_local = .false. if (present(modify_switch)) mod_local = modify_switch r_global = abs(x) x_global = x if (x < zero .and. mod_local) iflag = -1 call flows if (iflag == 0) then xi_global = two_third*r_global*sqrt(r_global) ! ! choose the region where x lies and call the appropriate subroutine ! to evaluate Bi(x) and/or Bi'(x) iregion = 0 if (r_global >= r_min) then if (x_global > zero) then iregion = 1 else iregion = 2 end if else iregion = 3 end if select case (iregion) case(1) bi0s = x_global**fourth ! dummy storage of number a = one/(sqrpi*bi0s) b = bi0s/sqrpi call asymp1r (bi0s, bi1s) bi0s = a*bi0s bi1s = b*bi1s if (.not. mod_local) then bi0s = exp(xi_global)*bi0s bi1s = exp(xi_global)*bi1s end if case(2) bi1s = r_global**fourth ! dummy storage of number a = one/(sqrpi*bi1s) b = bi1s/sqrpi call asymp2r (Pxi, Qxi, Rxi, Sxi) bi1s = xi_global-pi*fourth ! dummy storage of number bi0s = a*(-Pxi*sin(bi1s) + Qxi*cos(bi1s)) bi1s = b*( Rxi*cos(bi1s) + Sxi*sin(bi1s)) case(3) nn = ceiling(n_parts*r_global/r_min) call taylorr (nn, zero, tsm) bi0s = tsm(1,1)*bi0zer + tsm(1,2)*bi1zer bi1s = tsm(2,1)*bi0zer + tsm(2,2)*bi1zer if (mod_local) then bi0s = exp(-xi_global)*bi0s bi1s = exp(-xi_global)*bi1s end if end select end if ! ! error handling and return values select case (iflag) case(:-1) bi = zero dbi = zero ierr = iflag case(0:2) bi = bi0s dbi = bi1s ierr = 0 case(5) bi = bi0zer ! the value of Bi(0) dbi = bi1zer ! the value of d(Bi(0))/dx ierr = 0 end select end subroutine airy_bir !*** !************************************************************************ !*** subroutine flows !* ! PRIVATE subroutine to check overflow and underflow regions !* real(prd) tol ! if (iflag /= 0) return ! test for underflow for small |z| if (r_global <= r_lolimit) then iflag = 5 return end if ! ! test for overflow for large |z| if (r_global >= r_uplimit) then iflag = -6 return end if ! ! test the unscaled functions if (.not.mod_local .and. x_global > zero) then ! ! test to see if multiplication by exp(-xi) will underflow ! or if multiplication by exp(xi) for x>0 will over-flow tol = min( -log(two_sqrpi*r_global**fourth) - log(tiny(tol)), & -log(two_sqrpi/r_global**fourth) - log(tiny(tol)), & log(two_sqrpi/r_global**fourth) + log(huge(tol)), & log(two_sqrpi*r_global**fourth) + log(huge(tol)) ) & - 15.0_prd if (r_global >= (1.50_prd*tol)**two_third) then iflag = -7 return end if end if end subroutine flows !*** !************************************************************************ !*** subroutine asymp1r (F1, F2, xi_in) real(prd), intent(out) :: F1, F2 real(prd), intent(in), optional :: xi_in !* ! PRIVATE subroutine to compute the asymptotic expansions as x->infty !* integer i real(prd) xir ! ! summation of asymptotic series is done using nested multiplication ! Compute G(xi), G'(xi) F1 = ucoef(n_asymp) F2 = vcoef(n_asymp) if ( present(xi_in) ) then xir = one/xi_in else xir = one/xi_global end if if (iflag == 10) xir = -xir ! Ai(x) is an alternating series do i = n_asymp-1, 1, -1 F1 = ucoef(i) + xir*F1 F2 = vcoef(i) + xir*F2 end do F1 = (one + xir*F1) F2 = (one + xir*F2) iflag = 0 end subroutine asymp1r !*** !************************************************************************ !*** subroutine asymp2r (Pxi, Qxi, Rxi, Sxi) real(prd), intent(out) :: Pxi, Qxi, Rxi, Sxi !* ! PRIVATE subroutine to compute the asymptotic expansions as x->-infty !* integer i, itemp real(prd) xir ! ! summation of asymptotic series is done using nested multiplication ! Compute G(xi), G'(xi) itemp = floor(half*n_asymp - half) if (mod(itemp,2) /= 0) itemp = itemp - 1 Pxi = ucoef(itemp) Qxi = ucoef(itemp-1) Rxi = vcoef(itemp) Sxi = vcoef(itemp-1) xir = one/xi_global**2 do i = itemp-2, 2, -2 Pxi = ucoef(i ) - xir*Pxi Rxi = vcoef(i ) - xir*Rxi Qxi = ucoef(i-1) - xir*Qxi Sxi = vcoef(i-1) - xir*Sxi end do Pxi = one - xir*Pxi Rxi = one - xir*Rxi Qxi = Qxi/xi_global Sxi = Sxi/xi_global end subroutine asymp2r !*** !************************************************************************ !*** subroutine taylorr (nn, x_start, tsm) integer, intent(in) :: nn real(prd), intent(in) :: x_start real(prd), dimension(2,2), intent(out) :: tsm !* ! PRIVATE subroutine to integrate along a ray !* integer i, j real(prd) h, xm real (prd), dimension(0:n_taylor) :: pterm, qterm real (prd), dimension(nn,2,2) :: Phi ! h = (x_global - x_start)/nn ! step size and integration direction ! ! compute the Taylor series in each partition do i = 1, nn xm = (x_start+h*(i-1)) ! xm always appears with h**2 ! ! compute the reduced derivatives: pterm(0) = one pterm(1) = zero pterm(2) = xm*pterm(0)*half qterm(0) = zero qterm(1) = one qterm(2) = zero do j = 3, n_taylor ! ! compute the next Taylor series term for each partition. pterm(j) = (xm*pterm(j-2) + pterm(j-3))/real(j*j-j,prd) qterm(j) = (xm*qterm(j-2) + qterm(j-3))/real(j*j-j,prd) end do ! ! now sum the series Phi(i,1,1) = pterm(n_taylor) Phi(i,2,1) = pterm(n_taylor)*n_taylor Phi(i,1,2) = qterm(n_taylor) Phi(i,2,2) = qterm(n_taylor)*n_taylor do j = n_taylor-1,1,-1 Phi(i,1,1) = pterm(j) + h*Phi(i,1,1) Phi(i,2,1) = pterm(j)*j + h*Phi(i,2,1) Phi(i,1,2) = qterm(j) + h*Phi(i,1,2) Phi(i,2,2) = qterm(j)*j + h*Phi(i,2,2) end do Phi(i,1,1) = pterm(0) + h*Phi(i,1,1) Phi(i,1,2) = qterm(0) + h*Phi(i,1,2) end do ! ! multiply all the matrices together do i = 1,nn-1 Phi(i+1,:,:) = matmul(Phi(i+1,:,:),Phi(i,:,:)) end do ! ! return values tsm(:,:) = Phi(nn,:,:) end subroutine taylorr !*** !************************************************************************ !*** ! ! end file airy_functions_real ! ! begin file airy_aux ! version 1.0 ! edited 07/25/2003 ! !*** !************************************************************************ !*** subroutine airy_auxr (x, airy_mod, airy_phase, ierr, & dairy_mod, dairy_phase) real(prd), intent(in) :: x real(prd), intent(out) :: airy_mod, airy_phase real(prd), intent(out), optional :: dairy_mod, dairy_phase integer, intent(out) :: ierr !* ! PUBLIC subroutine is the driver routine that computes the modulus ! and phase functions for Ai(x) and Bi(x). The (PRIVATE) subroutines ! called are: ! parameter_airy = determines the machine specific parameters ! aux_parameter_airy = determines the machine specific parameters ! specific to the auxiliary functions ! asymp_aux = asymptotic expansion for auxiliary functions ! taylor_modulus = taylor series integration for modulus functions M(x) ! and N(x) ! taylor_phase = taylor series integration for phase functions theta(x) ! and phi(x) !* real(prd) m_inter, n_inter, theta_inter, phi_inter real(prd), dimension(3) :: tsm ! ! initialize all parameters if (.not. is_aux_init_airy) call aux_parameter_airy ! ! check to see if the main grid has been allocated and populated. if ( .not.big_integrate_aux ) then call taylor_modulus (n_parts_aux, tsm) call taylor_phase (n_parts_aux, theta_inter) end if iflag = 0 x_global = x r_global = abs(x) call flows if (x > zero) iflag = -2 if (iflag == 0) then ! ! choose the region where x lies and call the appropriate subroutine ! to evaluate Ai(x) and/or Ai'(x) rtest: if (x_global <= -aux_min) then call asymp_aux (m_inter, n_inter, theta_inter, phi_inter) else ! ! integrate the odes to obtain the functions and phases values ! using one step from the grid point to the desired point call taylor_modulus (1, tsm) call taylor_phase (1, theta_inter) m_inter = sqrt(tsm(1)) if ( present(dairy_mod) ) & n_inter = sqrt(tsm(2)**2*fourth/tsm(1) + one/(pi**2*tsm(1))) if ( present(dairy_phase) ) & phi_inter = theta_inter - atan(two/tsm(2)/pi) end if rtest end if ! ! error handling and return values select case (iflag) case(:-1) airy_mod = zero airy_phase = zero if ( present(dairy_mod) ) dairy_mod = zero if ( present(dairy_phase) ) dairy_phase = zero ierr = iflag case(0) airy_mod = m_inter airy_phase = theta_inter if ( present(dairy_mod) ) dairy_mod = n_inter if ( present(dairy_phase) ) dairy_phase = phi_inter ierr = 0 case(5) airy_mod = sqrt(m2grid(0,1)) airy_phase = pi/six if ( present(dairy_mod) ) dairy_mod = sqrt(ai1zer**2 + bi1zer**2) if ( present(dairy_phase) ) dairy_phase = -pi/six ierr = 0 end select end subroutine airy_auxr !*** !************************************************************************ !*** subroutine asymp_aux (m_inter, n_inter, theta_inter, phi_inter) real(prd), intent(out) :: m_inter, n_inter, theta_inter, phi_inter !* ! PRIVATE subroutine to compute the asymptotic expansions as x->-infty !* integer i real(prd) darg1, darg2, darg3 ! ! summation of asymptotic series is done using nested multiplication ! Compute the modulus functions darg1 = mcoef(n_asymp_mod) darg2 = ncoef(n_asymp_mod) darg3 = one/(x_global**3) do i = n_asymp_mod-1, 1, -1 darg1 = mcoef(i) + darg3*darg1 darg2 = ncoef(i) + darg3*darg2 end do darg1 = (one + darg3*darg1)/pi/sqrt(-x_global) darg2 = (one + darg3*darg2)*sqrt(-x_global)/pi m_inter = sqrt(darg1) n_inter = sqrt(darg2) ! ! Compute the phase functions darg1 = thetacoeff(n_asymp_phase) darg2 = phicoeff(n_asymp_phase) do i = n_asymp_phase-1, 1, -1 darg1 = thetacoeff(i) + darg3*darg1 darg2 = phicoeff(i) + darg3*darg2 end do darg1 = (thetacoeff(0) + darg3*darg1)*two_third*sqrt(-x_global)**3 darg2 = (phicoeff(0) + darg3*darg2)*two_third*sqrt(-x_global)**3 theta_inter = pi*fourth + darg1 phi_inter = -pi*fourth + darg2 ! end subroutine asymp_aux !*** !************************************************************************ !*** subroutine taylor_modulus (nn, tsm) integer, intent(in) :: nn real(prd), dimension(3), intent(out) :: tsm !* ! PRIVATE subroutine to integrate for the modulus function M^2(x) ! and the corresponding phase function theta(x). ! we check to see if the main grid has been computed and ! integrate from the closest partition point to the one desired. !* integer iplace, j, i real(prd) h, xright, xm real(prd), dimension(0:n_taylor_aux) :: pterm, qterm, rterm real(prd), dimension(nn,3,3) :: Phi ! h = -aux_min/n_parts_aux if ( nn == 1 ) then ! ! find the grid points between which x lies do i = n_parts_aux-1,0,-1 if ( x_global <= h*i ) then xright = h*i iplace = i exit end if end do h = x_global-xright end if ! ! compute the Taylor series in each partition do i = 1, nn xm = h*(i-1) if ( nn == 1 ) xm = xright ! ! compute the reduced derivatives: pterm(0) = one; pterm(1) = zero; pterm(2) = zero qterm(0) = zero; qterm(1) = one; qterm(2) = zero rterm(0) = zero; rterm(1) = zero; rterm(2) = half do j = 3, n_taylor_aux pterm(j) = (pterm(j-3)*(4*j-10) & + four*xm*(j-2)*pterm(j-2))/real(j*(j-1)*(j-2),prd) qterm(j) = (qterm(j-3)*(4*j-10) & + four*xm*(j-2)*qterm(j-2))/real(j*(j-1)*(j-2),prd) rterm(j) = (rterm(j-3)*(4*j-10) & + four*xm*(j-2)*rterm(j-2))/real(j*(j-1)*(j-2),prd) end do ! ! sum the series Phi(i,:,1) = pterm(n_taylor_aux) Phi(i,:,2) = qterm(n_taylor_aux) Phi(i,:,3) = rterm(n_taylor_aux) Phi(i,2,:) = Phi(i,2,:)*n_taylor_aux Phi(i,3,:) = Phi(i,3,:)*n_taylor_aux*(n_taylor_aux-1) do j = n_taylor_aux-1,2,-1 Phi(i,1,1) = h*Phi(i,1,1) + pterm(j) Phi(i,2,1) = h*Phi(i,2,1) + pterm(j)*j Phi(i,3,1) = h*Phi(i,3,1) + pterm(j)*j*(j-1) Phi(i,1,2) = h*Phi(i,1,2) + qterm(j) Phi(i,2,2) = h*Phi(i,2,2) + qterm(j)*j Phi(i,3,2) = h*Phi(i,3,2) + qterm(j)*j*(j-1) Phi(i,1,3) = h*Phi(i,1,3) + rterm(j) Phi(i,2,3) = h*Phi(i,2,3) + rterm(j)*j Phi(i,3,3) = h*Phi(i,3,3) + rterm(j)*j*(j-1) end do ! j = 1 term Phi(i,1,1) = h*Phi(i,1,1) + pterm(1) Phi(i,2,1) = h*Phi(i,2,1) + pterm(1) Phi(i,1,2) = h*Phi(i,1,2) + qterm(1) Phi(i,2,2) = h*Phi(i,2,2) + qterm(1) Phi(i,1,3) = h*Phi(i,1,3) + rterm(1) Phi(i,2,3) = h*Phi(i,2,3) + rterm(1) ! j = 0 term Phi(i,1,1) = h*Phi(i,1,1) + pterm(0) Phi(i,1,2) = h*Phi(i,1,2) + qterm(0) Phi(i,1,3) = h*Phi(i,1,3) + rterm(0) end do ! ! return values if ( nn == 1 ) then tsm(:) = Phi(1,:,1)*m2grid(iplace,1) + Phi(1,:,2)*m2grid(iplace,2) & + Phi(1,:,3)*m2grid(iplace,3) else ! ! obtain the function, first, and second derivative ! values at each partition point and store them ! into a global module variable allocate ( m2grid(0:n_parts_aux,3) ) ! ! initial conditions for M^2(0) m2grid(0,1) = ai0zer**2 + bi0zer**2 m2grid(0,2) = two*(ai0zer*ai1zer + bi0zer*bi1zer) m2grid(0,3) = two*(ai1zer**2 + bi1zer**2) ! ! build the rest of the two grids do i = 1,n_parts_aux m2grid(i,:) = Phi(i,:,1)*m2grid(i-1,1) + Phi(i,:,2)*m2grid(i-1,2) & + Phi(i,:,3)*m2grid(i-1,3) end do end if end subroutine taylor_modulus !*** !************************************************************************ !*** subroutine taylor_phase (nn, theta_inter) integer, intent(in) :: nn real(prd), intent(out) :: theta_inter !* ! PRIVATE subroutine to integrate for the phase functions M^2(x) ! and the corresponding phase function theta(x). ! we check to see if the main grid has been computed and ! integrate from the closest partition point to the one desired. !* integer iplace, j, i, k real(prd) h, xright, xm, xloc, darg real(prd), dimension(0:n_taylor_aux) :: M2taylor, theta ! h = -aux_min/n_parts_aux if ( nn == 1 ) then ! ! find the partition point closest to the desire value do i = n_parts_aux-1,0,-1 if ( x_global <= h*i ) then xright = h*i iplace = i exit end if end do h = x_global-xright else ! ! compute and populate a vector of function values for ! theta(x) on the same grid points. ! build the Taylor series expansion for M^2(x) ! at each grid point, then calculate theta(x) from the ! first order ODE allocate ( theta_grid(0:n_parts_aux) ) ! ! initial conditions for theta(0) theta_grid(0) = pi/six end if ! ! compute the Taylor series for M^2 in each partition do i = 1,nn if (nn == 1) then xm = xright xloc = x_global else iplace = i-1 xm = h*(i-1) xloc = h*i end if M2taylor(0) = m2grid(iplace,1) M2taylor(1) = m2grid(iplace,2) M2taylor(2) = m2grid(iplace,3)*half do j = 3,n_taylor_aux M2taylor(j) = (four*(j-2)*xm*M2taylor(j-2) & + M2taylor(j-3)*(4*j-10))/real(j*(j-1)*(j-2),prd) end do ! ! now build the Taylor series for the phase function theta(x) theta(0) = theta_grid(iplace) theta(1) = -one/pi/M2taylor(0) theta(2) = -theta(1)*M2taylor(1)*half/M2taylor(0) do k = 3,n_taylor_aux darg = zero do j = 0,k-2 darg = darg + theta(j+1)*(j+1)*M2taylor(k-j-1) end do theta(k) = -darg/M2taylor(0)/real(k,prd) end do ! ! now sum the Taylor series to obtain the value of theta. darg = theta(n_taylor_aux) do j = n_taylor_aux-1,1,-1 darg = h*darg + theta(j) end do theta_inter = h*darg + theta(0) if ( nn == n_parts_aux ) theta_grid(i) = theta_inter end do ! ! reset the logical and exit if ( nn == n_parts_aux ) big_integrate_aux = .true. end subroutine taylor_phase !*** !************************************************************************ !*** ! ! end file airy_aux ! ! begin file airy_zeros ! version 1.0 ! edited 07/24/2002 ! !*** !************************************************************************ !*** subroutine ai_zeror (n, ai_zero, ierr, ai_assoc, dai_zero, dai_assoc) integer, intent(in) :: n real(prd), intent(out) :: ai_zero integer, intent(out) :: ierr real(prd), intent(out), optional :: ai_assoc, dai_zero, dai_assoc !* ! PUBLIC subroutine which is the driver routine that computes the ! zeros of Ai(x) and optionally of Ai'(x), and the associated ! values Ai'(a_n) and Ai(a_n') for scalar arguments. The first 25 ! zeros are stored in an array and the rest are computed by summing ! approriate asymptotic expansions. The (PRIVATE) subroutines called are: ! parameters_airy = determines the machine specific parameters ! zero_parameters_airy = populates the vectors containing the ! stored zeros and the coefficients of the asymptotic expansions ! airy_air = computes the associated function value for stored zeros ! ae_zero_r = sums the appropriate asymptotic expansion for the ! zeros of the Airy functions and their associated function values !* integer itemp real(prd) xai, xdai, lam ! if (.not. is_zero_init_airy) call zero_parameter_airy iflag = 0 if (n <= 0) iflag = - 3 itemp = floor(0.25*(huge(itemp)-1)) if (n >= itemp) iflag = -4 ! zero requested cannot be computed if (iflag < 0) then ai_zero = zero if (present(ai_assoc)) ai_assoc = zero if (present(dai_zero)) dai_zero = zero if (present(dai_assoc)) dai_assoc = zero ierr = iflag return end if if (n <= n_zeros) then ai_zero = aizr(n) if (present(ai_assoc)) then call airy_air (aizr(n), xai, xdai, ierr) ai_assoc = xdai end if if (present(dai_zero)) dai_zero = daizr(n) if (present(dai_assoc)) then call airy_air (daizr(n), xai, xdai, ierr) dai_assoc = xai end if else lam = three_pi_ate*real(4*n-1,prd) call ae_zero_r (xai, lam, 1) ai_zero = -xai call ae_zero_r (xai, lam, 3) if (present(ai_assoc)) then ai_assoc = xai if ( mod(n-1,2) == 1 ) ai_assoc = -ai_assoc end if if (present(dai_zero) .or. present(dai_assoc)) then lam = three_pi_ate*(4*n-3) call ae_zero_r (xai, lam, 2) if (present(dai_zero)) dai_zero = -xai call ae_zero_r (xdai, lam, 4) if (present(dai_assoc)) then dai_assoc = xdai if ( mod(n-1,2) == 1 ) dai_assoc = -dai_assoc end if end if end if end subroutine ai_zeror !*** !************************************************************************ !*** subroutine bi_zeror (n, bi_zero, ierr, bi_assoc, dbi_zero, dbi_assoc) integer, intent(in) :: n real(prd), intent(out) :: bi_zero integer, intent(out) :: ierr real(prd), intent(out), optional :: bi_assoc, dbi_zero, dbi_assoc !* ! PUBLIC subroutine which is the driver routine that computes the ! zeros of Bi(x) and optionally of Bi'(x), and the associated ! values Bi'(b_n) and Bi(b_n') for scalar arguments. The first 25 ! zeros are stored in an array and the rest are computed by summing ! approriate asymptotic expansions. The (PRIVATE) subroutines called are: ! parameters_airy = determines the machine specific parameters ! zero_parameters_airy = populates the vectors containing the ! stored zeros and the coefficients of the asymptotic expansions ! airy_bir = computes the associated function value for stored zeros ! ae_zero_r = sums the appropriate asymptotic expansion for the ! zeros of the Airy functions and their associated function values !* integer itemp real(prd) xbi, xdbi, lam ! if (.not. is_zero_init_airy) call zero_parameter_airy iflag = 0 if (n <= 0) iflag = -3 itemp = floor(0.25*(huge(itemp)-1)) if (n >= itemp) iflag = -4 ! zero requested cannot be computed if (iflag < 0) then bi_zero = zero if (present(bi_assoc)) bi_assoc = zero if (present(dbi_zero)) dbi_zero = zero if (present(dbi_assoc)) dbi_assoc = zero ierr = iflag return end if if (n <= n_zeros) then bi_zero = bizr(n) if (present(bi_assoc)) then call airy_bir (bi_zero, xbi, xdbi, ierr) bi_assoc = xdbi end if if (present(dbi_zero)) dbi_zero = dbizr(n) if (present(dbi_assoc)) then call airy_bir (dbizr(n), xbi, xdbi, ierr) dbi_assoc = xbi end if else lam = three_pi_ate*real(4*n-3,prd) call ae_zero_r (xbi, lam, 1) bi_zero = -xbi call ae_zero_r (xbi, lam, 3) if (present(bi_assoc)) then bi_assoc = xbi if ( mod(n-1,2) == 1 ) bi_assoc = -bi_assoc end if if (present(dbi_zero) .or. present(dbi_assoc)) then lam = three_pi_ate*(4*n-1) call ae_zero_r (xbi, lam, 2) if (present(dbi_zero)) dbi_zero = -xbi call ae_zero_r (xdbi, lam, 4) if (present(dbi_assoc)) then dbi_assoc = xdbi if ( mod(n,2) == 1 ) dbi_assoc = -dbi_assoc end if end if end if end subroutine bi_zeror !*** !************************************************************************ !*** subroutine ae_zero_r (zr, lam, in) real(prd), intent (out) :: zr real(prd), intent (in) :: lam integer, intent (in) :: in !* ! PRIVATE subroutine to sum the asymptotic expansions for the zeros ! of the Airy functions and their associated functions !* real(prd) zrsum, lamr integer i ! lamr = one/lam**2 select case (in) case(1) ! sum the T(x) expansion zrsum = Tcoeff(n_asymp_zero) do i = n_asymp_zero-1,0,-1 zrsum = Tcoeff(i) + zrsum*lamr end do zr = zrsum*lam**two_third case(2) ! sum the U(x) expansion zrsum = Ucoeff(n_asymp_zero) do i = n_asymp_zero-1,0,-1 zrsum = Ucoeff(i) + zrsum*lamr end do zr = zrsum*lam**two_third case(3) ! sum the V(x) expansion zrsum = Vcoeff(n_asymp_asso) do i = n_asymp_asso-1,0,-1 zrsum = Vcoeff(i) + zrsum*lamr end do zr = zrsum*lam**(half/three)/sqrpi case(4) ! sum the W(x) expansion zrsum = Wcoeff(n_asymp_asso) do i = n_asymp_asso-1,0,-1 zrsum = Wcoeff(i) + zrsum*lamr end do zr = zrsum/lam**(half/three)/sqrpi end select end subroutine ae_zero_r !*** !************************************************************************ !*** subroutine ai_zerorv (n, ai_zero, ierr, ai_assoc, dai_zero, dai_assoc) integer, intent(in) :: n real(prd), intent(out), dimension(:) :: ai_zero integer, intent(out) :: ierr real(prd), intent(out), dimension(:), optional :: & ai_assoc, dai_zero, dai_assoc !* !* ! PUBLIC subroutine which is the driver routine that computes the ! zeros of Ai(x) and optionally of Ai'(x), and the associated ! values Ai'(a_n) and Ai(a_n') for vector arguments. The first 25 ! zeros are stored in an array and the rest are computed by summing ! approriate asymptotic expansions. The (PRIVATE) subroutines called are: ! parameters_airy = determines the machine specific parameters ! zero_parameters_airy = populates the vectors containing the ! stored zeros and the coefficients of the asymptotic expansions ! airy_air = computes the associated function value for stored zeros ! ae_zero_rv = sums the appropriate asymptotic expansion for the ! zeros of the Airy functions and their associated function values !* integer itemp, K, i, iflag_zero real(prd) xai, xdai real(prd), dimension(:), allocatable :: lam ! ! if (.not. is_zero_init_airy) call zero_parameter_airy iflag_zero = 0 if (n <= 0) iflag_zero = -3 itemp = floor(0.25*(huge(itemp)-1)) if (n >= itemp) iflag_zero = -4 ! zero requested cannot be computed K = size(ai_zero) ! ! check to see if the last zero requested is computable ! and change the last zero returned accordingly. if (n+K >= itemp .and. iflag_zero /= -4) then do i = n+K-1,n,-1 if ( i < itemp ) then K = i-n iflag_zero = 50 exit end if end do end if if ( iflag_zero < 0) then ai_zero(:) = zero if (present(ai_assoc)) ai_assoc(:) = zero if (present(dai_zero)) dai_zero(:) = zero if (present(dai_assoc)) dai_assoc(:) = zero ierr = iflag_zero return end if if (n <= n_zeros) then itemp = min(n+K,n_zeros) ai_zero(1:itemp-n+1) = aizr(n:itemp) if (present(ai_assoc)) then do i = n, itemp call airy_air (aizr(i), xai, xdai, ierr) ai_assoc(i-n+1) = xdai end do end if if (present(dai_zero)) dai_zero(1:itemp-n+1) = daizr(n:itemp) if (present(dai_assoc)) then do i = n, itemp call airy_air (daizr(i), xai, xdai, ierr) dai_assoc(i-n+1) = xai end do end if end if if (n+K > n_zeros) then itemp = max(n,n_zeros+1) allocate (lam(itemp:n+K-1)) do i = itemp,n+K-1 lam(i) = three_pi_ate*real(4*i-1,prd) end do call ae_zero_rv (ai_zero(itemp-n+1:K), lam, 1) ai_zero(itemp-n+1:K) = -ai_zero(itemp-n+1:K) if (present(ai_assoc)) then call ae_zero_rv (ai_assoc(itemp-n+1:K), lam, 3) do i = itemp,n+K-1 if ( mod(i-1,2) == 1 ) ai_assoc(i-n+1) = -ai_assoc(i-n+1) end do end if deallocate (lam) if (present(dai_zero) .or. present(dai_assoc)) then allocate (lam(itemp:n+K-1)) do i = itemp,n+K-1 lam(i) = three_pi_ate*real(4*i-3,prd) end do if (present(dai_zero)) then call ae_zero_rv (dai_zero(itemp-n+1:K), lam, 2) dai_zero(itemp-n+1:K) = -dai_zero(itemp-n+1:K) end if if (present(dai_assoc)) then call ae_zero_rv (dai_assoc(itemp-n+1:K), lam, 4) do i = itemp,n+K-1 if ( mod(i-1,2) == 1 ) dai_assoc(i-n+1) = -dai_assoc(i-n+1) end do end if deallocate (lam) end if end if ierr = iflag_zero end subroutine ai_zerorv !*** !************************************************************************ !*** subroutine bi_zerorv (n, bi_zero, ierr, bi_assoc, dbi_zero, dbi_assoc) integer, intent(in) :: n real(prd), intent(out), dimension(:) :: bi_zero integer, intent(out) :: ierr real(prd), intent(out), dimension(:), optional :: & bi_assoc, dbi_zero, dbi_assoc !* !* ! PUBLIC subroutine which is the driver routine that computes the ! zeros of Bi(x) and optionally of Bi'(x), and the associated ! values Bi'(b_n) and Bi(b_n') for vector arguments. The first 25 ! zeros are stored in an array and the rest are computed by summing ! approriate asymptotic expansions. The (PRIVATE) subroutines called are: ! parameters_airy = determines the machine specific parameters ! zero_parameters_airy = populates the vectors containing the ! stored zeros and the coefficients of the asymptotic expansions ! airy_bir = computes the associated function value for stored zeros ! ae_zero_rv = sums the appropriate asymptotic expansion for the ! zeros of the Airy functions and their associated function values !* integer itemp, K, i, iflag_zero real(prd) xbi, xdbi real(prd), dimension(:), allocatable :: lam ! if (.not. is_zero_init_airy) call zero_parameter_airy iflag_zero = 0 if (n <= 0) iflag_zero = -3 itemp = floor(0.25*(huge(itemp)-1)) if (n >= itemp) iflag_zero = -4 ! zero requested cannot be computed K = size(bi_zero) ! ! check to see if the last zero requested is computable ! and change the last zero returned accordingly. if (n+K >= itemp .and. iflag_zero /= -4) then do i = n+K-1,n,-1 if ( i < itemp ) then K = i-n iflag_zero = 50 exit end if end do end if if (iflag_zero < 0) then bi_zero(:) = zero if (present(bi_assoc)) bi_assoc(:) = zero if (present(dbi_zero)) dbi_zero(:) = zero if (present(dbi_assoc)) dbi_assoc(:) = zero ierr = iflag_zero return end if if (n <= n_zeros) then itemp = min(n+K,n_zeros) bi_zero(1:itemp-n+1) = bizr(n:itemp) if (present(bi_assoc)) then do i = n, itemp call airy_bir (bizr(i), xbi, xdbi, ierr) bi_assoc(i-n+1) = xdbi end do end if if (present(dbi_zero)) dbi_zero(1:itemp-n+1) = dbizr(n:itemp) if (present(dbi_assoc)) then do i = n, itemp call airy_bir (dbizr(i), xbi, xdbi, ierr) dbi_assoc(i-n+1) = xbi end do end if end if if (n+K > n_zeros) then itemp = max(n,n_zeros+1) allocate (lam(itemp:n+K-1)) do i = itemp,n+K-1 lam(i) = three_pi_ate*real(4*i-3,prd) end do call ae_zero_rv (bi_zero(itemp-n+1:K), lam, 1) bi_zero(itemp-n+1:K) = -bi_zero(itemp-n+1:K) if (present(bi_assoc)) then call ae_zero_rv (bi_assoc(itemp-n+1:K), lam, 3) do i = itemp,n+K-1 if ( mod(i-1,2) == 1 ) bi_assoc(i-n+1) = -bi_assoc(i-n+1) end do end if deallocate (lam) if (present(dbi_zero) .or. present(dbi_assoc)) then allocate (lam(itemp:n+K-1)) do i = itemp,n+K-1 lam(i) = three_pi_ate*real(4*i-1,prd) end do if (present(dbi_zero)) then call ae_zero_rv (dbi_zero(itemp-n+1:K), lam, 2) dbi_zero(itemp-n+1:K) = -dbi_zero(itemp-n+1:K) end if if (present(dbi_assoc)) then call ae_zero_rv (dbi_assoc(itemp-n+1:K), lam, 4) do i = itemp,n+K-1 if ( mod(i,2) == 1 ) dbi_assoc(i-n+1) = -dbi_assoc(i-n+1) end do end if deallocate (lam) end if end if ierr = iflag_zero end subroutine bi_zerorv !*** !************************************************************************ !*** subroutine ae_zero_rv (zr, lam, in) real(prd), intent (out), dimension(1:) :: zr real(prd), intent (in), dimension(1:) :: lam integer, intent (in) :: in real(prd), dimension(size(zr)) :: zrsum, lamr !* ! PRIVATE subroutine to sum the asymptotic expansions for the zeros ! of the Airy functions and their associated functions for the ! vector case !* integer k, i, j ! k = size(zr) do i = 1,k lamr(i) = one/lam(i)**2 end do if (in == 1) then ! sum the T(x) expansion do j = 1,k zrsum(j) = Tcoeff(n_asymp_zero) do i = n_asymp_zero-1,0,-1 zrsum(j) = Tcoeff(i) + zrsum(j)*lamr(j) end do end do zr(:) = zrsum(:)*lam(:)**two_third elseif (in == 2) then ! sum the U(x) expansion do j = 1,k zrsum(j) = Ucoeff(n_asymp_zero) do i = n_asymp_zero,0,-1 zrsum(j) = Ucoeff(i) + zrsum(j)*lamr(j) end do end do zr(:) = zrsum(:)*lam(:)**two_third elseif (in == 3) then ! sum the V(x) expansion do j = 1,k zrsum(j) = Vcoeff(n_asymp_asso) do i = n_asymp_asso-1,0,-1 zrsum(j) = Vcoeff(i) + zrsum(j)*lamr(j) end do end do zr(:) = zrsum(:)*lam(:)**(half/three)/sqrpi else ! sum the W(x) expansion do j = 1,k zrsum(j) = Wcoeff(n_asymp_asso) do i = n_asymp_asso-1,0,-1 zrsum(j) = Wcoeff(i) + zrsum(j)*lamr(j) end do end do zr(:) = zrsum(:)/lam(:)**(half/three)/sqrpi end if end subroutine ae_zero_rv !*** !************************************************************************ !*** ! ! end file airy_zeros SHAR_EOF fi # end of overwriting check if test -f 'airyerr.log' then echo shar: will not over-write existing file "'airyerr.log'" else cat << "SHAR_EOF" > 'airyerr.log' This is a log file for the installation of the Airy function package by B.R. Fabijonas. Please send all bug reports to bfabi@smu.edu. The reference for this package is B.R. Fabijonas 2004. Algorithm XXX: Airy functions. ACM Trans. Math. Softw., Vol.xxx, pp.xxx-xxx. A detailed description of the computation method and parameter determination can be found in B.R. Fabijonas, D.W. Lozier, and F.W.J. Olver, Computation of complex Airy functions and their zeros using asymptotics and the differential equation. ACM Trans. Math. Softw., Vol.xxx, pp. xxx-xxx. ***Parameter values used in the code which are compiler specific. We use the abbreviations no. for number of terms and a.e. for asymptotic expansion. Single precision computations: standard Airy functions machine epsilon 0.1192092895507812E-06 radius of the delineating circle rho 0.5394496440887451E+01 integration step size h 0.8990827202796936E+00 no. in the Taylor series N 19 no. in the a.e.s S 17 number of partitions of the integration ray P 6 Single precision computations: auxiliary functions radius of the delineating circle rho_m 0.5058679103851318E+01 integrationstep size h_m 0.8431131839752197E+00 no. of terms in the Taylor series N_m 20 no. of terms in the a.e.s of the modulus functions S_m 8 no. of terms in the a.e.s of the phase functions S_theta 5 no. of partitions of the integration ray P_m 6 Double precision computations: standard Airy functions machine epsilon 0.2220446049250313E-15 the parameter rho 0.9127616308750468E+01 step size h 0.9127616308750468E+00 no. in the Taylor series N 31 no. used in the a.e.s S 37 number of partitions of the integration ray P 10 Double precision computations: auxiliary functions the parameter rho_m 0.8818819021035965E+01 step size h_m 0.9798687801151073E+00 no. in the t.s. N_m 35 no. in the a.e. of the modulus functions S_m 17 no. in the a.e. of the phase functions S_theta 12 number of partitions of the integration ray P_m 9 ***Comparing the function values returned by the real routines against stored values Computing Ai at x= 5.82340 We computed 0.15392729E-04 in single precision We computed 0.1539272803479132E-04 in double precision and the stored value is 0.1539272803479132E-04 Computing dAi at x= 5.82340 We computed -0.37779439E-04 in single precision We computed -0.3777943926223622E-04 in double precision and the stored value is -0.3777943926223622E-04 Computing Bi at x= 5.82340 We computed 0.42881152E+04 in single precision We computed 0.4288114816683048E+04 in double precision and the stored value is 0.4288114816683048E+04 Computing dBi at x= 5.82340 We computed 0.10154620E+05 in single precision We computed 0.1015462058214280E+05 in double precision and the stored value is 0.1015462058214280E+05 Computing Ai at x=-17.34589 We computed -0.26777703E+00 in single precision We computed -0.2677772589719993E+00 in double precision and the stored value is -0.2677772589719993E+00 Computing dAi at x=-17.34589 We computed -0.29003367E+00 in single precision We computed -0.2900295189410030E+00 in double precision and the stored value is -0.2900295189410030E+00 Computing Bi at x=-17.34589 We computed 0.68710074E-01 in single precision We computed 0.6870907209572491E-01 in double precision and the stored value is 0.6870907209572491E-01 Computing dBi at x=-17.34589 We computed -0.11142915E+01 in single precision We computed -0.1114292633371775E+01 in double precision and the stored value is -0.1114292633371775E+01 The number of points tested is 8 The number of errors found is 0 ***Comparing the function values returned by the complex routines against stored values Computing Ai at z= 1.56620 +I*( 3.68100) In single precision, we computed 0.38073376E+00 +I*( 0.36051720E+00) In double precision, we computed 0.3807338825307106E+00 +I*( 0.3605171603516821E+00) and the stored value is 0.3807338825307106E+00 +I*( 0.3605171603516821E+00) Computing dAi at z= 1.56620 +I*( 3.68100) In single precision, we computed -0.26851702E+00 +I*(-0.10107594E+01) In double precision, we computed -0.2685171251031442E+00 +I*(-0.1010759493443000E+01) and the stored value is -0.2685171251031442E+00 +I*(-0.1010759493443000E+01) Computing Ai at z= -9.83420 +I*( 4.28110) In single precision, we computed 0.10694007E+06 +I*(-0.47724039E+05) In double precision, we computed 0.1069401631365655E+06 +I*(-0.4772398758353833E+05) and the stored value is 0.1069401631365655E+06 +I*(-0.4772398758353830E+05) Computing dAi at z= -9.83420 +I*( 4.28110) In single precision, we computed -0.22165348E+06 +I*(-0.31107988E+06) In double precision, we computed -0.2216534220094920E+06 +I*(-0.3110801974925891E+06) and the stored value is -0.2216534220094920E+06 +I*(-0.3110801974925891E+06) The number of points tested is 4 The number of errors found is 0 ***Testing function values on both side of the delineating circle and verifying continuity within a tolerance (no output printed unless an error is found) The number of points tested is 800 The number of errors found is 0 ***Comparing the computation in different precisions at random points on the real line and in the complex plane and verifying the accuracy of the lower precision function values (no output printed unless an error is found). The number of points tested is 800 The number of errors found is 0 ***Testing the vectorized routine of the Airy functions (no output printed unless an error is found). The number of points tested is 820 The number of errors found is 0 ***Testing the routines for the zeros of the Airy functions real and complex zeros--scalar version **Zeros by table lookup Testing the 18th zero of Ai We computed -0.19126381E+02 in single precision We computed -0.1912638047424695E+02 in double precision and the stored value is -0.1912638047424695E+02 Testing the associated function value of Ai at the 18 zero We computed -0.11798807E+01 in single precision We computed -0.1179880729870146E+01 in double precision and the stored value is -0.1179880729870146E+01 Testing the 18th zero of dAi We computed -0.18764799E+02 in single precision We computed -0.1876479843766596E+02 in double precision and the stored value is -0.1876479843766596E+02 Testing the associated function value of dAi at the 18 zero We computed -0.27107027E+00 in single precision We computed -0.2710702785769711E+00 in double precision and the stored value is -0.2710702785769711E+00 Testing the 18th zero of Bi We computed -0.18765509E+02 in single precision We computed -0.1876550828448008E+02 in double precision and the stored value is -0.1876550828448008E+02 Testing the associated function value of Bi at the 18 zero We computed -0.11742761E+01 in single precision We computed -0.1174276253118059E+01 in double precision and the stored value is -0.1174276253118059E+01 Testing the 18th zero of dBi We computed -0.19125698E+02 in single precision We computed -0.1912569715641264E+02 in double precision and the stored value is -0.1912569715641264E+02 Testing the associated function value of dBi at the 18 zero We computed 0.26978263E+00 in single precision We computed 0.2697826139865426E+00 in double precision and the stored value is 0.2697826139865426E+00 Testing the 18th complex zero of Bi In single precision, we computed 0.94946041E+01 +I*( 0.16603624E+02) In double precision, we computed 0.9494603679187176E+01 +I*( 0.1660362460689876E+02) and the stored value is 0.9494603679187176E+01 +I*( 0.1660362460689876E+02) Testing the associated function value of Bi at the 18 complex zero In single precision, we computed 0.14459233E+01 +I*(-0.83281225E+00) In double precision, we computed 0.1445920793639127E+01 +I*(-0.8328073286646603E+00) and the stored value is 0.1445920793639133E+01 +I*(-0.8328073286646578E+00) Testing the 18th complex zero of dBi In single precision, we computed 0.93131523E+01 +I*( 0.16290871E+02) In double precision, we computed 0.9313152405593771E+01 +I*( 0.1629087030713368E+02) and the stored value is 0.9313152405593771E+01 +I*( 0.1629087030713368E+02) Testing the associated function value of dBi at the 18 complex zero In single precision, we computed -0.33219612E+00 +I*(-0.19132091E+00) In double precision, we computed -0.3321948851433578E+00 +I*(-0.1913210621443876E+00) and the stored value is -0.3321948851433578E+00 +I*(-0.1913210621443901E+00) The number of points tested is 12 The number of errors found is 0 **Zeros by asymptotic expansion Testing the 86th zero of Ai We computed -0.54657593E+02 in single precision We computed -0.5465758649186870E+02 in double precision and the stored value is -0.5465758649186870E+02 Testing the associated function value of Ai at the 86 zero We computed -0.15340441E+01 in single precision We computed -0.1534044239238300E+01 in double precision and the stored value is -0.1534044239238300E+01 Testing the 86th zero of dAi We computed -0.54444836E+02 in single precision We computed -0.5444482679260981E+02 in double precision and the stored value is -0.5444482679260981E+02 Testing the associated function value of dAi at the 86 zero We computed -0.20769954E+00 in single precision We computed -0.2076995780278862E+00 in double precision and the stored value is -0.2076995780278862E+00 Testing the 86th zero of Bi We computed -0.54444920E+02 in single precision We computed -0.5444491113058687E+02 in double precision and the stored value is -0.5444491113058687E+02 Testing the associated function value of Bi at the 86 zero We computed -0.15325499E+01 in single precision We computed -0.1532549805062612E+01 in double precision and the stored value is -0.1532549805062612E+01 Testing the 86th zero of dBi We computed -0.54657509E+02 in single precision We computed -0.5465750280893616E+02 in double precision and the stored value is -0.5465750280893616E+02 Testing the associated function value of dBi at the 86 zero We computed 0.20749725E+00 in single precision We computed 0.2074972409267743E+00 in double precision and the stored value is 0.2074972409267743E+00 Testing the 86th complex zero of Bi In single precision, we computed 0.27288200E+02 +I*( 0.47358311E+02) In double precision, we computed 0.2728820066764437E+02 +I*( 0.4735830615386251E+02) and the stored value is 0.2728820066764437E+02 +I*( 0.4735830615386251E+02) Testing the associated function value of Bi at the 86 complex zero In single precision, we computed 0.18790455E+01 +I*(-0.10843303E+01) In double precision, we computed 0.1879045614304762E+01 +I*(-0.1084330361798879E+01) and the stored value is 0.1879045614304762E+01 +I*(-0.1084330361798879E+01) Testing the 86th complex zero of dBi In single precision, we computed 0.27181744E+02 +I*( 0.47174103E+02) In double precision, we computed 0.2718174151707578E+02 +I*( 0.4717409672492563E+02) and the stored value is 0.2718174151707578E+02 +I*( 0.4717409672492563E+02) Testing the associated function value of dBi at the 86 complex zero In single precision, we computed -0.25441056E+00 +I*(-0.14681086E+00) In double precision, we computed -0.2544106266601734E+00 +I*(-0.1468108932911459E+00) and the stored value is -0.2544106266601735E+00 +I*(-0.1468108932911459E+00) The number of points tested is 24 The number of errors found is 0 **Zeros--vectorized version (no output printed unless an error is found) The number of points tested is 60 The number of errors found is 0 ***Comparing the computed modulus, phase, and associated values against stored values **by ODE integration Testing the modulus M(x) at x= -5.23400 We computed 0.37280816E+00 in single precision We computed 0.3728081698402362E+00 in double precision and the stored value is 0.3728081698402362E+00 Testing the phase theta(x) at x= -5.23400 We computed 0.87596397E+01 in single precision We computed 0.8759640547244738E+01 in double precision and the stored value is 0.8759640547244738E+01 Testing the modulus M(x) at x= -5.23400 We computed 0.85400015E+00 in single precision We computed 0.8540001773409781E+00 in double precision and the stored value is 0.8540001773409781E+00 Testing the phase phi(x) at x= -5.23400 We computed 0.72095666E+01 in single precision We computed 0.7209566736849307E+01 in double precision and the stored value is 0.7209566736849307E+01 **by asymptotic expansion Testing the modulus M(x) at x= -25.39200 We computed 0.25133255E+00 in single precision We computed 0.2513325654640693E+00 in double precision and the stored value is 0.2513325654640693E+00 Testing the phase theta(x) at x= -25.39200 We computed 0.86085594E+02 in single precision We computed 0.8608558068173974E+02 in double precision and the stored value is 0.8608558068173974E+02 Testing the modulus N(x) at x= -25.39200 We computed 0.12664913E+01 in single precision We computed 0.1266491244788154E+01 in double precision and the stored value is 0.1266491244788154E+01 Testing the phase phi(x) at x= -25.39200 We computed 0.84516747E+02 in single precision We computed 0.8451673808738920E+02 in double precision and the stored value is 0.8451673808738920E+02 **in different precisions (no output printed unless an error is found) The number of points tested is 808 The number of errors found is 0 Number of tests performed: 3336 Number of errors found: 0 SHAR_EOF fi # end of overwriting check if test -f 'res' then echo shar: will not over-write existing file "'res'" else cat << "SHAR_EOF" > 'res' --------------------------------------------------- This is a program which tests the installation the Airy function pacakge by B.R. Fabijonas. Please send all bug reports to bfabi@smu.edu Comparing the function values returned by the real routines against stored values...done with 0 errors Comparing the function values returned by the complex routines against stored values...done with 0 errors Comparing the function values on both side of the delineating circle and verifying continuity within a tolerance...done with 0 errors Comparing the computation in different precisions at random points on the real line and in the complex plane and verifying the accuracy of the lower precision function values....done with 0 errors Testing the vectorized routine of the Airy functions...done with 0 errors Testing the routines for the zeros of the Airy functions real zeros by table lookup...................done with 0 errors real zeros by asymptotic expansion...........done with 0 errors real and complex zeros--vectorized version...done with 0 errors Comparing the computed modulus, phase, and associated values against stored values...done with 0 errors Number of tests performed: 3336 Number of errors found: 0 End of run SHAR_EOF fi # end of overwriting check if test -f 'test_install.f90' then echo shar: will not over-write existing file "'test_install.f90'" else cat << "SHAR_EOF" > 'test_install.f90' ! updated 04/26/2004 ! This is a program which tests the installation ! The Airy function pacakge by B.R.Fabijonas. ! reference: B.R. Fabijonas, Algorithm XXX: Airy functions, ! ACM Trans. Math. Softw., Vol. xxx, pp. xxx-xxx. ! Please report all bug reports to bfabi@mail.smu.edu ! ! This programs tests all featuers of the package: ! Airy functions of real argument computed by integration of ODE ! Airy functions or real argument computed by summation of ! asymptotic expansions ! Airy functions of complex argument computed by integration of ODE ! Airy functions or complex argument computed by summation of ! asymptotic expansions ! Real and complex zeros of the Airy functions computed by table lookup ! Real and complex zeros of the Airy functions computed by ! summation of asymptotic expansions ! Vectorized version of zeros of Ai(x) and its derivative with ! associated values ! Vectorized version of zeros of Bi(x) and its derivative with ! associated values ! Vectorized version of zeros of Bi(z) and its derivative with ! associated values ! Auxiliary functions for negative real arguments computed by ! integration of ODE ! Auxiliary functions for negative real arguments computed by ! summation of asymptotic expansions ! ! A module to interface the error computation and switch ! between relative precision and absolute error. ! the endings: ! s = real single ! d = real double ! c = complex single ! z = complex double ! module error implicit none private public errg integer, parameter :: prs = kind(0.0e0) integer, parameter :: prd = kind(0.0d0) real(prd), parameter :: zerod = 0.0_prd, & pid = 3.1415926535897932384626433832795028841971_prd real(prs), parameter :: zeros = 0.0_prs, pis = pid interface errg module procedure errds module procedure errs module procedure errd module procedure errzc module procedure errz module procedure errc end interface interface psi module procedure psis module procedure psid end interface contains ! ! real double and single real(prd) function errds(arg1,arg2) real(prd) arg1 real(prs) arg2 errds = abs(psi(abs(arg1)) - psi(abs(arg2))) end function errds ! ! real single real(prs) function errs(arg1,arg2) real(prs) arg1, arg2 errs = abs(psi(arg1) - psi(arg2)) end function errs ! ! real double real(prd) function errd(arg1,arg2) real(prd) arg1, arg2 errd = abs(psi(arg1) - psi(arg2)) end function errd ! ! complex double and single real(prd) function errzc(arg1,arg2) complex(prd) arg1, arg3 complex(prs) arg2 real(prd) x1,x2 arg3 = arg2 if ( abs(arg1) < 1.0_prd .or. abs(arg3) < 1.0_prd ) then errzc = abs(arg1-arg2) else x1 = atan2(aimag(arg1),real(arg1,prd)) x2 = atan2(aimag(arg3),real(arg3,prd)) x2 = mod(abs(x1-x2),2.0_prd*pid) errzc = abs(cmplx(log(abs(arg1/arg3)),min(x2,2.0_prd*pid-x2),prd)) end if end function errzc ! ! complex single real(prs) function errc(arg1,arg2) complex(prs) arg1, arg2 real(prs) x1,x2 if ( abs(arg1) < 1.0_prs .or. abs(arg2) < 1.0_prs ) then errc = abs(arg1-arg2) else x1 = atan2(aimag(arg1),real(arg1)) x2 = atan2(aimag(arg2),real(arg2)) x2 = mod(abs(x1-x2),2.0_prs*pis) errc = abs(cmplx(log(abs(arg1/arg2)),min(x2,2.0_prs*pis-x2))) end if end function errc ! ! complex double real(prd) function errz(arg1,arg2) complex(prd) arg1, arg2 real(prd) x1,x2 if ( abs(arg1) < 1.0_prd .or. abs(arg2) < 1.0_prd ) then errz = abs(arg1-arg2) else x1 = atan2(aimag(arg1),real(arg1,prd)) x2 = atan2(aimag(arg2),real(arg2,prd)) x2 = mod(abs(x1-x2),2.0_prd*pid) errz = abs(cmplx(log(abs(arg1/arg2)),min(x2,2.0_prd*pid-x2),prd)) end if end function errz ! ! real single real(prs) function psis(arg) real(prs) arg, at at = abs(arg) if (at > 1.0_prs) then psis = sign(1.0_prs + log(at),arg) else psis = arg end if end function psis ! ! real double real(prd) function psid(arg) real(prd) arg, at at = abs(arg) if (at > 1.0_prd) then psid = sign(1.0_prd + log(at),arg) else psid = arg end if end function psid end module error ! ! ! program to test installation ! program test use airy_functions_real_single use airy_functions_real_double use airy_functions_complex_single use airy_functions_complex_double use error implicit none integer, parameter :: prs = kind(0.0e0) integer, parameter :: prd = kind(0.0d0) real(prd), parameter :: oned = 1.0_prd, thvd = 1.50_prd real(prs), parameter :: ones = 1.0_prs, thvs = 1.50_prs real(prs) aias, daias, aibs, daibs, rms, pis, g1, tols, ffacs, hs, rs real(prd) x, aix, daix, aiass, daiass, rmd, pid, f1, told, ffacd, hd, rd real (prd), dimension(5) :: aizv, daizv, aiassv, daiassv complex(prs) zs, ciunits, aiasz, daiasz, aibsz, daibsz complex(prd) z, aiz, daiz, aiassz, daiassz, ciunitd complex(prd), dimension(5) :: bizv, dbizv, biassv, dbiassv real(prd), dimension(20) :: ray complex(prd), dimension(20) :: airay, dairay character(3) ai, dai, bi, dbi integer ierr, itol, itol2 integer ierrcnt, i, iptcnt, ilocerrcnt, ilocptcnt, j, nt, na, np, nam, nap open(19,file="airyerr.log",status='unknown') pid = 3.14159265358979323846264338327950288419716939937510582097494459230782_prd pis = real(pid,prs) ciunits = cmplx(0.0_prs,1.0_prs,prs) ciunitd = cmplx(0.0_prd,1.0_prd,prd) ! ! factor for testing computed values agains stored values or for computing ! values at the same point in different precisions itol = 100 ! ! factor for extremely close, but different, argument values itol2 = 100 ! ! machine epsilons tols = exp((1-digits(tols))*log(real(radix(tols),prs))) told = exp((1-digits(told))*log(real(radix(told),prd))) ! ! set the counters and declare a few variables ierrcnt = 0 iptcnt = 0 ai = ' Ai' dai = 'dAi' bi = ' Bi' dbi = 'dBi' print*, ' ---------------------------------------------------' print*, ' ' print*, 'This is a program which tests the installation ' print*, ' the Airy function pacakge by B.R. Fabijonas.' print*, ' Please send all bug reports to bfabi@smu.edu' print*, ' ' ! ! write some preliminaries to the output file write(19,*) ' This is a log file for the installation of the Airy function' write(19,*) ' package by B.R. Fabijonas. Please send all bug reports' write(19,*) ' to bfabi@smu.edu. The reference for this package is ' write(19,*) ' B.R. Fabijonas 2004. Algorithm XXX: Airy functions. ' write(19,*) ' ACM Trans. Math. Softw., Vol.xxx, pp.xxx-xxx.' write(19,*) ' A detailed description of the computation method and parameter' write(19,*) ' determination can be found in B.R. Fabijonas, D.W. Lozier,' write(19,*) ' and F.W.J. Olver, Computation of complex Airy functions and ' write(19,*) ' their zeros using asymptotics and the differential equation. ' write(19,*) ' ACM Trans. Math. Softw., Vol.xxx, pp. xxx-xxx.' write(19,*) ' ' write(19,'("***Parameter values used in the code which are compiler specific.")') write(19,'(" We use the abbreviations no. for number of terms and ")') write(19,'(" a.e. for asymptotic expansion.")') write(19,'(" Single precision computations: standard Airy functions")') call airy_info(rms, hs, nt, na, np) write(19,'(A45,e24.16)') ' machine epsilon', epsilon(rms) write(19,'(A45,e24.16)') ' radius of the delineating circle rho', rms write(19,'(A45,e24.16)') ' integration step size h', hs write(19,'(A60,i4)') ' no. in the Taylor series N',nt write(19,'(A60,i4)') ' no. in the a.e.s S', na write(19,'(A60,i4/)') ' number of partitions of the integration ray P', np write(19,'(" Single precision computations: auxiliary functions")') call airy_aux_info ( rs, hs, nt, nam, nap, np) write(19,'(A45,e24.16)') ' radius of the delineating circle rho_m ',rs write(19,'(A45,e24.16)') ' integrationstep size h_m ',hs write(19,'(A60,i4)') ' no. of terms in the Taylor series N_m ',nt write(19,'(A60,i4)') ' no. of terms in the a.e.s of the modulus functions S_m ',nam write(19,'(A60,i4)') ' no. of terms in the a.e.s of the phase functions S_theta ',nap write(19,'(A60,i4/)') ' no. of partitions of the integration ray P_m ',np write(19,'(" Double precision computations: standard Airy functions")') call airy_info(rmd, hd, nt, na, np) write(19,'(A45,e24.16)') ' machine epsilon ', epsilon(rmd) write(19,'(A45,e24.16)') ' the parameter rho ', rmd write(19,'(A45,e24.16)') ' step size h ', hd write(19,'(A60,i4)') ' no. in the Taylor series N ',nt write(19,'(A60,i4)') ' no. used in the a.e.s S ', na write(19,'(A60,i4/)') ' number of partitions of the integration ray P ', np write(19,'(" Double precision computations: auxiliary functions")') call airy_aux_info ( rd, hd, nt, nam, nap, np) write(19,'(A45,e24.16)') ' the parameter rho_m ',rd write(19,'(A45,e24.16)') ' step size h_m ',hd write(19,'(A60,i4)') ' no. in the t.s. N_m ',nt write(19,'(A60,i4)') ' no. in the a.e. of the modulus functions S_m ',nam write(19,'(A60,i4)') ' no. in the a.e. of the phase functions S_theta ',nap write(19,'(A60,i4//)') ' number of partitions of the integration ray P_m ',np ! !***************** ! testing real routines !**************** ilocerrcnt = 0 ilocptcnt = 0 write(6,'(" Comparing the function values returned by the")') write(6,'(" real routines against stored values")',advance="NO") write(19,'( "***Comparing the function values returned by the")') write(19,'(" real routines against stored values")') x = 5.82340_prd ! integration routine ffacd = oned !10.0_prd**max(oned,thvd*log(abs(x))) call airy_ai(x, aix, daix, ierr) call airy_ai(real(x,prs),aias,daias,ierr) write(19,1000) ai, x write(19,1001) aias, aix, 0.1539272803479132e-04_prd f1 = errg(aix,0.1539272803479132e-04_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 f1 = errg(daix,-0.3777943926223622e-04_prd) write(19,1000) dai, x write(19,1001) daias, daix, -0.3777943926223622e-04_prd if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 call airy_bi(x, aix, daix, ierr) call airy_bi(real(x,prs),aias,daias,ierr) write(19,1000) bi, x write(19,1001) aias, aix, 0.4288114816683048e+04_prd f1 = errg(aix,0.4288114816683048e+04_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1000) dbi, x write(19,1001) daias, daix, 0.1015462058214280e+05_prd f1 = errg(daix,0.1015462058214280e+05_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 x = -17.345890_prd ! asymptotic expansion routine ffacd = oned !10.0_prd**max(oned,thvd*log(abs(x))) call airy_ai(x, aix, daix, ierr) call airy_ai(real(x,prs),aias,daias,ierr) write(19,1000) ai, x write(19,1001) aias, aix, -0.2677772589719993e+00_prd f1 = errg(aix,-0.2677772589719993e+00_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1000) dai, x write(19,1001) daias, daix, -0.2900295189410030e+00_prd f1 = errg(daix,-0.2900295189410030e+00_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 call airy_bi(x, aix, daix, ierr) call airy_bi(real(x,prs),aias,daias,ierr) write(19,1000) bi, x write(19,1001) aias, aix, 0.6870907209572491e-01_prd f1 = errg(aix,0.6870907209572491e-01_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1000) dbi, x write(19,1001) daias, daix, -0.1114292633371775e+01_prd f1 = errg(daix,-0.1114292633371775e+01_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...done with ",i3," errors")') ilocerrcnt ! !***************** ! testing complex routines !**************** ilocerrcnt = 0 ilocptcnt = 0 print*, ' ' write(6,'(" Comparing the function values returned by the")') write(6,'(" complex routines against stored values")',advance="NO") write(19,'( "***Comparing the function values returned by the complex routines")') write(19,'(" against stored values")') z = cmplx(1.5662_prd,3.681_prd,prd) ! integration routine call airy_ai(z, aiz, daiz, ierr) call airy_ai(cmplx(1.5662,3.681),aiasz,daiasz,ierr) ffacd = oned !10.0_prd**max(oned,thvd*log(abs(z))) write(19,2000) ai, real(z,prd),aimag(z) write(19,2001) real(aiasz), aimag(aiasz), real(aiz,prd), aimag(aiz), & 0.3807338825307106e+00_prd,0.3605171603516821e+00_prd f1 = errg(aiz,& cmplx(0.3807338825307106e+00_prd,0.3605171603516821e+00_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,2000) dai, real(z,prd),aimag(z) write(19,2001) real(daiasz), aimag(daiasz), real(daiz,prd), aimag(daiz), & -0.2685171251031442e+00_prd,-0.1010759493443000e+01_prd f1 = errg(daiz,& cmplx(-0.2685171251031442e+00_prd,-0.1010759493443000e+01_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 2 z = cmplx(-9.8342_prd,4.2811_prd,prd) ! asymptotic expansion routine call airy_ai(z, aiz, daiz, ierr) call airy_ai(cmplx(-9.8342,4.2811),aiasz,daiasz,ierr) write(19,2000) ai, real(z,prd),aimag(z) write(19,2001) real(aiasz), aimag(aiasz), real(aiz,prd), aimag(aiz), & 0.1069401631365655e+06_prd,-0.4772398758353830e+05_prd f1 = errg(aiz,& cmplx(0.1069401631365655e+06_prd,-0.4772398758353830e+05_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,2000) dai, real(z,prd),aimag(z) write(19,2001) real(daiasz), aimag(daiasz), real(daiz,prd), aimag(daiz), & -0.2216534220094920e+06_prd,-0.3110801974925891e+06_prd f1 = errg(daiz,& cmplx(-0.2216534220094920e+06_prd,-0.3110801974925891e+06_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 2 write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...done with ",i3," errors")') ilocerrcnt ! !***************** ! testing interface of the routines !**************** ! a perturbation analysis on the asymptotic expansion shows that ! |Ai(x(1+eps)) - Ai(x)| <= |Ai(x)|*(.25+x**(3/2))*eps. the ! same holds for the derivative. ! we will use this tolerance print*, ' ' write(6,'(" Comparing the function values on both side of the delineating circle")') write(6,'(" and verifying continuity within a tolerance")',advance="NO") write(19,'("***Testing function values on both side of the delineating circle")') write(19,'(" and verifying continuity within a tolerance")') write(19,'(" (no output printed unless an error is found)")') ilocerrcnt = 0 ilocptcnt = 0 ! for real variables call airy_info(rms) call airy_info(rmd) call airy_ai(rms+epsilon(rms)*2,aias,daias,ierr) call airy_ai(rms-epsilon(rms)*2,aibs,daibs,ierr) ffacs = abs(aias)*(.250_prs + rms*sqrt(rms))*4*epsilon(rms) g1 = errg(aias,aibs) if ( g1 > itol*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,3000) ai, rms write(19,3001) aias, aibs end if ffacs = abs(daias)*(.250_prs + rms*sqrt(rms))*4*epsilon(rms) g1 = errg(daias,daibs) if ( g1 > itol*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,3000) dai, rms write(19,3001) daias, daibs end if call airy_ai(rmd+epsilon(rmd)*2,aix,daix,ierr) call airy_ai(rmd-epsilon(rmd)*2,aiass,daiass,ierr) ffacd = abs(aix)*(.250_prd + rmd*sqrt(rmd))*4*epsilon(rmd) f1 = errg(aix,aiass) if ( f1 > itol*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,3000) ai, rmd write(19,3001) aix, aiass end if f1 = errg(daix,daiass) ffacd = abs(daix)*(.250_prd + rmd*sqrt(rmd))*4*epsilon(rmd) if ( f1 > itol*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,3000) dai, rmd write(19,3001) daix, daiass end if ilocptcnt = ilocptcnt + 4 ! for complex variables do i = -99,99 zs = (rms + epsilon(rms)*2)*exp(ciunits*i*0.01_prs*pis) call airy_ai(zs,aiasz,daiasz,ierr) zs = (rms - epsilon(rms)*2)*exp(ciunits*i*0.01_prs*pis) call airy_ai(zs,aibsz,daibsz,ierr) ffacs = abs(aiasz)*(0.250_prs+rms*sqrt(rms))*4*epsilon(rms) g1 = errg(aiasz,aibsz) if ( g1 > itol*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,4000) ai, rms*exp(ciunits*i*0.01_prs*pis) write(19,4001) aiasz,aibsz write(19,*) g1, itol*ffacs end if g1 = errg(daiasz,daibsz) ffacs = abs(daiasz)*(0.250_prs+rms*sqrt(rms))*4*epsilon(rms) if ( g1 > itol*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,*) g1, ffacs write(19,4000) dai, rms*exp(ciunits*i*0.01_prs*pis) write(19,4001) daiasz,daibsz end if z = (rmd + epsilon(rmd)*2)*exp(ciunitd*i*0.01_prd*pid) call airy_ai(z,aiz,daiz,ierr) z = (rmd - epsilon(rmd)*2)*exp(ciunitd*i*0.01_prd*pid) call airy_ai(z,aiassz,daiassz,ierr) f1 = errg(aiz,aiassz) ffacd = abs(aiz)*(0.250_prd+rmd*sqrt(rmd))*4*epsilon(rmd) if ( f1 > itol*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,4000) ai, rmd*exp(ciunitd*i*0.01_prd*pid) write(19,4002) aiasz,aiassz, aiassz end if f1 = errg(daiz,daiassz) ffacd = abs(daiz)*(0.250_prd+rmd*sqrt(rmd))*4*epsilon(rmd) if ( f1 > itol*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,4000) dai, rmd*exp(ciunitd*i*0.01_prd*pid) write(19,4002) daiz,daiassz end if ilocptcnt = ilocptcnt + 4 end do write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...done with ",i3," errors")') ilocerrcnt ! !***************** ! testing different precisions !**************** write(6,'(" ")') write(6,'(" Comparing the computation in different precisions at random points")') write(6,'(" on the real line and in the complex plane and verifying the")') write(6,'(" accuracy of the lower precision function values.")',advance='no') write(19,'("***Comparing the computation in different precisions at random points ")') write(19,'(" on the real line and in the complex plane and verifying the accuracy of")') write(19,'(" the lower precision function values (no output printed unless an error is found).")') ilocerrcnt = 0 ilocptcnt = 0 ! for real variables do i = 1,200 call random_number(rmd) rms = (rmd-0.50_prd)*20.0_prd rmd = rms call airy_ai(rmd, aix, daix, ierr) call airy_ai(rms, aias, daias, ierr) f1 = errg(aix,aias) ffacs = ones !10.0_prs**max(ones,thvs*log(abs(rms))) if ( f1 > itol2*tols*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,5000) ai, rmd write(19,5001) aias, aix end if f1 = errg(daix,daias) if ( f1 > itol2*tols*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,5000) dai, rmd write(19,5001) daias, daix end if ilocptcnt = ilocptcnt + 2 end do ! for complex variables do i = 1,200 call random_number(rmd) call random_number(x) zs = rmd*10*exp(ciunitd*(x-0.50_prd)*2.0_prd*pid) z = zs call airy_ai(zs,aiasz,daiasz,ierr) call airy_ai(z,aiz,daiz,ierr) ffacs = ones !10.0_prs**max(ones,thvs*log(abs(zs))) f1 = errg(aiz,aiasz) if ( f1 > itol*tols*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,6000) ai, z write(19,6001) aiasz,aiz end if f1 = errg(daiz,daiasz) if ( f1 > itol*tols*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,6000) dai, z write(19,6001) daiasz,daiz end if ilocptcnt = ilocptcnt + 2 end do write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...done with ",i3," errors")') ilocerrcnt ! !***************** ! testing vectorized airy_ai routine !**************** print*, ' ' write(6,'(" Testing the vectorized routine of the Airy functions")',advance='no') write(19,'("***Testing the vectorized routine of the Airy functions (no output printed")') write(19,'(" unless an error is found).")') ilocerrcnt = 0 ilocptcnt = 0 do j = -20,20 do i = 1,20 ray(i) = 0.5_prd*i ilocptcnt = ilocptcnt + 1 end do call airy_ai(ray,0.05_prd*j*pid,airay,dairay,ierr,.true.) do i = 1,20 call airy_ai(ray(i)*exp(ciunitd*j*0.05_prd*pid), aiz,daiz, ierr) g1 = errg (aiz,airay(i)) ffacd = oned !10.0_prd**max(oned,thvd*log(ray(i))) if ( g1 > itol2*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,2002) ai, ray(i), j*0.05_prd*pid write(19,2003) aiz, airay(i) end if g1 = errg (daiz,dairay(i)) if ( g1 > itol2*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,2002) dai, ray(i), j*0.05_prd*pid write(19,2003) daiz, dairay(i) end if end do end do write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...done with ",i3," errors")') ilocerrcnt ! !***************** ! testing zero routines !**************** print*, ' ' print*, 'Testing the routines for the zeros of the Airy functions' write(19,'("***Testing the routines for the zeros of the Airy functions")') write(19,'(" real and complex zeros--scalar version")') ilocerrcnt = 0 ilocptcnt = 0 write(6,'(" real zeros by table lookup")',advance='no') write(19,'(" **Zeros by table lookup")') call airy_ai_zero (18, aias, ierr, ai_assoc=daias, dai_zero=aibs, & dai_assoc=daibs) call airy_ai_zero (18, aix, ierr, ai_assoc=aiass, dai_zero=daix, & dai_assoc=daiass) ffacd = 10.0_prd**max(oned,thvd*log(abs(18.0_prd))) write(19,1002) 18, ai write(19,1001) aias, aix, -19.126380474246954_prd f1 = errg(aix, -19.126380474246954_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) ai, 18 write(19,1001) daias, aiass,-1.179880729870146_prd f1 = errg(aiass,-1.179880729870146_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1002) 18, dai write(19,1001) aibs, daix, -18.764798437665956_prd f1 = errg(daix,-18.764798437665956_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) dai, 18 write(19,1001) daibs, daiass,-0.2710702785769711_prd f1 = errg(daiass,-0.2710702785769711_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 call airy_bi_zero (18, aias, ierr, bi_assoc=daias, dbi_zero=aibs, & dbi_assoc=daibs) call airy_bi_zero (18, aix, ierr, bi_assoc=aiass, dbi_zero=daix, & dbi_assoc=daiass) write(19,1002) 18, bi write(19,1001) aias, aix, -18.765508284480081_prd f1 = errg(aix,-18.765508284480081_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) bi, 18 write(19,1001) daias, aiass, -1.174276253118059_prd f1 = errg(aiass,-1.174276253118059_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1002) 18, dbi write(19,1001) aibs, daix, -19.125697156412638_prd f1 = errg(daix,-19.125697156412638_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 18 write(19,1001) daibs, daiass,0.2697826139865426_prd f1 = errg(daiass,0.2697826139865426_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 call airy_bi_zero (18, aiasz, ierr, bi_assoc=daiasz, dbi_zero=aibsz, & dbi_assoc=daibsz) call airy_bi_zero (18, aiz, ierr, bi_assoc=aiassz, dbi_zero=daiz, & dbi_assoc=daiassz) write(19,1004) 18, bi write(19,2001) aiasz, aiz, 9.494603679187176_prd,16.603624606898762_prd f1 = errg(aiz,cmplx(9.494603679187176_prd,16.603624606898762_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1005) bi, 18 write(19,2001) daiasz, aiassz, 1.445920793639133_prd,-0.8328073286646578_prd f1 = errg(aiassz,cmplx(1.445920793639133_prd,-0.8328073286646578_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1004) 18, dbi write(19,2001) aibsz, daiz, 9.313152405593771_prd,16.290870307133684_prd f1 = errg(daiz,cmplx(9.313152405593771_prd,16.290870307133684_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1005) dbi, 18 write(19,2001) daibsz, daiassz, -0.3321948851433578_prd,-0.1913210621443901_prd f1 = errg(daiassz,& cmplx(-0.3321948851433578_prd,-0.1913210621443901_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...................done with ",i3," errors")') ilocerrcnt ! !***************** ! testing zero routines !**************** write(6,'(" real zeros by asymptotic expansion")',advance='no') write(19,'(" **Zeros by asymptotic expansion")') call airy_ai_zero (86, aias, ierr, ai_assoc=daias, dai_zero=aibs, & dai_assoc=daibs) call airy_ai_zero (86, aix, ierr, ai_assoc=aiass, dai_zero=daix, & dai_assoc=daiass) ffacd = oned !10.0_prd**max(oned,thvd*log(abs(86.0_prd))) write(19,1002) 86, ai write(19,1001) aias, aix, -54.657586491868699_prd f1 = errg(aix,-54.657586491868699_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) ai, 86 write(19,1001) daias, aiass,-1.534044239238300_prd f1 = errg(aiass,-1.534044239238300_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1002) 86, dai write(19,1001) aibs, daix, -54.444826792609810_prd f1 = errg(daix,-54.444826792609810_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) dai, 86 write(19,1001) daibs, daiass,-0.2076995780278862_prd f1 = errg(daiass,-0.2076995780278862_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 call airy_bi_zero (86, aias, ierr, bi_assoc=daias, dbi_zero=aibs, & dbi_assoc=daibs) call airy_bi_zero (86, aix, ierr, bi_assoc=aiass, dbi_zero=daix, & dbi_assoc=daiass) write(19,1002) 86, bi write(19,1001) aias, aix, -54.444911130586874_prd f1 = errg(aix,-54.444911130586874_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) bi, 86 write(19,1001) daias, aiass,-1.532549805062612_prd f1 = errg(aiass,-1.532549805062612_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1002) 86, dbi write(19,1001) aibs, daix,-54.657502808936158_prd f1 = errg(daix,-54.657502808936158_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 86 write(19,1001) daibs, daiass, 0.2074972409267743_prd f1 = errg(daiass, 0.2074972409267743_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 call airy_bi_zero (86, aiasz, ierr, bi_assoc=daiasz, dbi_zero=aibsz, & dbi_assoc=daibsz) call airy_bi_zero (86, aiz, ierr, bi_assoc=aiassz, dbi_zero=daiz, & dbi_assoc=daiassz) write(19,1004) 86, bi write(19,2001) aiasz, aiz, 27.288200667644368_prd,47.358306153862507_prd f1 = errg(aiz,cmplx(27.288200667644368_prd,47.358306153862507_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1005) bi, 86 write(19,2001) daiasz, aiassz, 1.879045614304762_prd,-1.084330361798879_prd f1 = errg(aiassz,cmplx(1.879045614304762_prd,-1.084330361798879_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1004) 86, dbi write(19,2001) aibsz, daiz, 27.181741517075782_prd,47.174096724925626_prd f1 = errg(daiz,cmplx(27.181741517075782_prd,47.174096724925626_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1005) dbi, 86 write(19,2001) daibsz, daiassz, -0.2544106266601735_prd,-0.1468108932911459_prd f1 = errg(daiassz,cmplx(-0.2544106266601735_prd,-0.1468108932911459_prd,prd)) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...........done with ",i3," errors")') ilocerrcnt ! !***************** ! testing vectorized zero routines !**************** write(6,'(" real and complex zeros--vectorized version")',advance="no") write(19,'(" **Zeros--vectorized version (no output printed unless an error is found)")') ilocerrcnt = 0 ilocptcnt = 0 call airy_ai_zero (23, aizv, ierr, ai_assoc=aiassv) ffacd = oned !10.0_prd**max(oned,thvd*log(abs(27.0_prd))) f1 = errg(aizv(1),-22.567612917496504_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 23, ai write(19,1008) aizv(1), -22.567612917496504_prd end if f1 = errg(aizv(2),-23.224165001121680_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 24, ai write(19,1008) aizv(2), -23.224165001121680_prd end if f1 = errg(aizv(3),-23.871564455535918_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 25, ai write(19,1008) aizv(3), -23.871564455535918_prd end if f1 = errg(aizv(4),-24.510301236589672_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 26, ai write(19,1008) aizv(4), -24.510301236589672_prd end if f1 = errg(aizv(5),-25.140821166148957_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 27, ai write(19,1008) aizv(5), -25.140821166148957_prd end if f1 = errg(aiassv(1),1.229700701509681_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) ai, 23 write(19,1008) aiassv(1),1.229700701509681_prd end if f1 = errg(aiassv(2),-1.238547875329632_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) ai, 24 write(19,1008) aiassv(2),-1.238547875329632_prd end if f1 = errg(aiassv(3),1.247089945259408_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) ai, 25 write(19,1008) aiassv(3),1.247089945259408_prd end if f1 = errg(aiassv(4),-1.255349140475735_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) ai, 26 write(19,1008) aiassv(4),-1.255349140475735_prd end if f1 = errg(aiassv(5),1.263345282750799_prd) if ( f1 > itol*told*ffacd )then ilocerrcnt = ilocerrcnt+1 write(19,1003) ai, 27 write(19,1008) aiassv(5),1.263345282750799_prd end if ilocptcnt = ilocptcnt + 10 call airy_ai_zero (23, aizv, ierr, dai_zero=daizv, dai_assoc=daiassv) f1 = errg(daizv(1),-22.235232285348914_prd) if ( f1 > itol*told*ffacd )then ilocerrcnt = ilocerrcnt+1 write(19,1002) 23, dai write(19,1008) daizv(1),-22.235232285348914_prd end if f1 = errg(daizv(2),-22.896588738874620_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 24, dai write(19,1008) daizv(2),-22.896588738874620_prd end if f1 = errg(daizv(3),-23.548526295928802_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 25, dai write(19,1008) daizv(3),-23.548526295928802_prd end if f1 = errg(daizv(4),-24.191559709526349_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 26, dai write(19,1008) daizv(4),-24.191559709526349_prd end if f1 = errg(daizv(5),-24.826156425921152_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 27, dai write(19,1008) daizv(5),-24.826156425921152_prd end if f1 = errg(daiassv(1),0.259812670151466_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dai, 23 write(19,1008) daiassv(1),0.259812670151466_prd end if f1 = errg(daiassv(2),-0.257916075332572_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dai, 24 write(19,1008) daiassv(2),-0.257916075332572_prd end if f1 = errg(daiassv(3),0.256112333779654_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dai, 25 write(19,1008) daiassv(3),0.256112333779654_prd end if f1 = errg(daiassv(4),-0.254393342646825_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dai, 26 write(19,1008) daiassv(4),-0.254393342646825_prd end if f1 = errg(daiassv(5),0.252751992576574_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dai, 27 write(19,1008) daiassv(5),0.252751992576574_prd end if ilocptcnt = ilocptcnt + 10 ! !***************** ! testing vectorized zero routines !**************** call airy_bi_zero (23, aizv, ierr, bi_assoc=aiassv) f1 = errg(aizv(1),-22.235737881803384_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 23, bi write(19,1008) aizv(1),-22.235737881803384_prd end if f1 = errg(aizv(2),-22.897065554219793_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 24, bi write(19,1008) aizv(2),-22.897065554219793_prd end if f1 = errg(aizv(3),-23.548977079642448_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 25, bi write(19,1008) aizv(3),-23.548977079642448_prd end if f1 = errg(aizv(4),-24.191986850648995_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 26, bi write(19,1008) aizv(4),-24.191986850648995_prd end if f1 = errg(aizv(5),-24.826562012152888_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 27, bi write(19,1008) aizv(5),-24.826562012152888_prd end if f1 = errg(aiassv(1),1.225154995846960_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) bi, 23 write(19,1008) aiassv(1),1.225154995846960_prd end if f1 = errg(aiassv(2),-1.234163920487302_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) bi, 24 write(19,1008) aiassv(2),-1.234163920487302_prd end if f1 = errg(aiassv(3),1.242855598092317_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) bi, 25 write(19,1008) aiassv(3),1.242855598092317_prd end if f1 = errg(aiassv(4),-1.251253611221631_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) bi, 26 write(19,1008) aiassv(4),-1.251253611221631_prd end if f1 = errg(aiassv(5),1.259378938688538_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) bi, 27 write(19,1008) aiassv(5),1.259378938688538_prd end if ilocptcnt = ilocptcnt + 10 call airy_bi_zero (23, aizv, ierr, dbi_zero=daizv, dbi_assoc=daiassv) f1 = errg(daizv(1),-22.567122080497199_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 23, dbi write(19,1008) daizv(1),-22.567122080497199_prd end if f1 = errg(daizv(2),-23.223701521208962_prd) if ( f1 > itol*told*ffacd )then ilocerrcnt = ilocerrcnt+1 write(19,1002) 24, dbi write(19,1008) daizv(2),-23.223701521208962_prd end if f1 = errg(daizv(3),-23.871125771677974_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 25, dbi write(19,1008) daizv(3),-23.871125771677974_prd end if f1 = errg(daizv(4),-24.509885117016239_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 26, dbi write(19,1008) daizv(4),-24.509885117016239_prd end if f1 = errg(daizv(5),-25.140425655367874_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 27, dbi write(19,1008) daizv(5),-25.140425655367874_prd end if f1 = errg(daiassv(1),-0.258852215916922_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 23 write(19,1008) daiassv(1),-0.258852215916922_prd end if f1 = errg(daiassv(2),0.257003129647316_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 24 write(19,1008) daiassv(2),0.257003129647316_prd end if f1 = errg(daiassv(3),-0.255242710064815_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 25 write(19,1008) daiassv(3),-0.255242710064815_prd end if f1 = errg(daiassv(4),0.253563372437697_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 26 write(19,1008) daiassv(4),0.253563372437697_prd end if f1 = errg(daiassv(5),-0.251958444330784_prd) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 27 write(19,1008) daiassv(5),-0.251958444330784_prd end if ilocptcnt = ilocptcnt + 10 ! !***************** ! testing vectorized zero routines !**************** call airy_bi_zero (23, bizv, ierr, bi_assoc=biassv) f1 = errg(bizv(1),& cmplx(11.220656371214879_prd,19.580653883038824_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1004) 23, bi write(19,2003) bizv(1),11.220656371214879_prd,19.580653883038824_prd end if f1 = errg(bizv(2),& cmplx(11.549830143889196_prd,20.148722567309367_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1004) 24, bi write(19,2003) bizv(2), 11.549830143889196_prd,20.148722567309367_prd end if f1 = errg(bizv(3),& cmplx(11.874378641347221_prd,20.708893464731311_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1004) 25, bi write(19,2003) bizv(3),11.874378641347221_prd,20.708893464731311_prd end if f1 = errg(bizv(4),& cmplx(12.194551330833875_prd,21.261588251953778_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1004) 26, bi write(19,2003) bizv(4),12.194551330833875_prd,21.261588251953778_prd end if f1 = errg(bizv(5),& cmplx(12.510575046777394_prd,21.807190718498749_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1004) 27, bi write(19,2003) bizv(5),12.510575046777394_prd,21.807190718498749_prd end if f1 = errg(biassv(1),& cmplx(-1.506774749698633_prd,0.868314074899650_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1005) bi, 23 write(19,2003) biassv(1),-1.506774749698633_prd,0.868314074899650_prd end if f1 = errg(biassv(2),& cmplx(1.517585357164310_prd,-0.874612709515018_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1005) bi, 24 write(19,2003) biassv(2),1.517585357164310_prd,-0.874612709515018_prd end if f1 = errg(biassv(3),& cmplx(-1.528024149210944_prd,0.880692431621191_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1005) bi, 25 write(19,2003) biassv(3),-1.528024149210944_prd,0.880692431621191_prd end if f1 = errg(biassv(4),& cmplx(1.538118145012280_prd,-0.886569309158123_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1005) bi, 26 write(19,2003) biassv(4),1.538118145012280_prd,-0.886569309158123_prd end if f1 = errg(biassv(5),& cmplx(-1.547891444975066_prd,0.892257657608977_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1005) bi, 27 write(19,2003) biassv(5),-1.547891444975066_prd,0.892257657608977_prd end if ilocptcnt = ilocptcnt + 10 call airy_bi_zero (23, bizv, ierr, dbi_zero=dbizv, dbi_assoc=dbiassv) f1 = errg(dbizv(1),& cmplx(11.053994360667899_prd,19.293078213959429_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 23, dbi write(19,2003) dbizv(1),11.053994360667899_prd,19.293078213959429_prd end if f1 = errg(dbizv(2),& cmplx(11.385596969302076_prd,19.865292017137630_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 24, dbi write(19,2003) dbizv(2),11.385596969302076_prd,19.865292017137630_prd end if f1 = errg(dbizv(3),& cmplx(11.712438641126221_prd,20.429378925421950_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 25, dbi write(19,2003) dbizv(3),11.712438641126221_prd,20.429378925421950_prd end if f1 = errg(dbizv(4),& cmplx(12.034781571133943_prd,20.985781895965161_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 26, dbi write(19,2003) dbizv(4),12.034781571133943_prd,20.985781895965161_prd end if f1 = errg(dbizv(5),& cmplx(12.352863681490561_prd,21.534903282918666_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1002) 27, dbi write(19,2003) dbizv(5),12.352863681490561_prd,21.534903282918666_prd end if f1 = errg(dbiassv(1),& cmplx(0.318355274563740_prd,0.183451937317410_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 23 write(19,2003) dbiassv(1),0.318355274563740_prd,0.183451937317410_prd end if f1 = errg(dbiassv(2),& cmplx(-0.316024910474316_prd,-0.182124025592498_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 24 write(19,2003) dbiassv(2),-0.316024910474316_prd,-0.182124025592498_prd end if f1 = errg(dbiassv(3),& cmplx(0.313808934629479_prd,0.180860595992543_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 25 write(19,2003) dbiassv(3),0.313808934629479_prd,0.180860595992543_prd end if f1 = errg(dbiassv(4),& cmplx(-0.311697340122842_prd,-0.179656065952459_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 26 write(19,2003) dbiassv(4),-0.311697340122842_prd,-0.179656065952459_prd end if f1 = errg(dbiassv(5),& cmplx(0.309681349894832_prd,0.178505532129885_prd,prd)) if ( f1 > itol*told*ffacd ) then ilocerrcnt = ilocerrcnt+1 write(19,1003) dbi, 27 write(19,2003) dbiassv(5),0.309681349894832_prd,0.178505532129885_prd end if ilocptcnt = ilocptcnt + 10 write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...done with ",i3," errors")') ilocerrcnt ! !***************** ! testing modulus and phase routines !**************** print*, ' ' write(6,'(" Comparing the computed modulus, phase, and associated values")') write(6,'(" against stored values")',advance='no') write(19,'("***Comparing the computed modulus, phase, and associated values against stored values")') ilocerrcnt = 0 ilocptcnt = 0 x = -5.234_prd write(19,'(" **by ODE integration")') call airy_aux ( real(x,prs), aias, daias, ierr, aibs, daibs) call airy_aux ( x, aix, daix, ierr, aiass, daiass) ffacd = oned !10.0_prd**max(oned,thvd*log(abs(x))) write(19,1006) 'M(x)', x write(19,1001) aias, aix,0.3728081698402362_prd f1 = errg(aix,0.3728081698402362_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1007) 'theta(x)', x write(19,1001) daias, daix,8.759640547244738_prd f1 = errg(daix,8.759640547244738_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1006) 'M(x)', x write(19,1001) aibs, aiass,0.8540001773409781_prd f1 = errg(aiass,0.8540001773409781_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1007) 'phi(x)', x write(19,1001) daibs, daiass,7.209566736849307_prd f1 = errg(daiass,7.209566736849307_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 ! !***************** ! testing modulus and phase routines !**************** x = -25.392_prd write(19,'(" **by asymptotic expansion")') call airy_aux ( real(x,prs), aias, daias, ierr, aibs, daibs) call airy_aux ( x, aix, daix, ierr, aiass, daiass) ffacd = oned !10.0_prd**max(oned,thvd*log(abs(x))) write(19,1006) 'M(x)', x write(19,1001) aias, aix,0.2513325654640693_prd f1 = errg(aix,0.2513325654640693_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1007) 'theta(x)', x write(19,1001) daias, daix,86.085580681739742_prd f1 = errg(daix,86.085580681739742_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1006) 'N(x)', x write(19,1001) aibs, aiass,1.2664912447881540_prd f1 = errg(aiass,1.2664912447881540_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 write(19,1007) 'phi(x)', x write(19,1001) daibs, daiass,84.516738087389200_prd f1 = errg(daiass,84.516738087389200_prd) if ( f1 > itol*told*ffacd ) ilocerrcnt = ilocerrcnt+1 ilocptcnt = ilocptcnt + 4 !***************** ! testing modulus and phase in different precisions !**************** write(19,'(" **in different precisions (no output printed unless an error is found)")') do i = 1,200 call random_number(rmd) rms = -rmd*10.0_prd rmd = rms call airy_aux ( rmd, aix, daix, ierr, aiass, daiass) call airy_aux ( rms, aias, daias, ierr, aibs, daibs) ffacd = ones !10.0_prs**max(ones,thvs*log(abs(rms))) f1 = errg(aix,aias) if ( f1 > itol*tols*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,5002) ai, rmd write(19,5001) aias,aix end if f1 = errg(daix,daias) if ( f1 > itol*tols*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,5002) ai, rmd write(19,5001) daias,daix end if f1 = errg(aiass,aibs) if ( f1 > itol*tols*ffacs ) then ilocerrcnt = ilocerrcnt+1 write(19,5003) ai, rmd write(19,5001) aibs,aiass end if f1 = errg(daiass,daibs) if ( f1 > itol*tols*ffacs ) then ilocerrcnt = ilocerrcnt+11 write(19,5003) ai, rmd write(19,5001) daibs,daiass end if ilocptcnt = ilocptcnt + 4 end do write(19,1500) ilocptcnt write(19,1501) ilocerrcnt iptcnt = iptcnt + ilocptcnt ierrcnt = ierrcnt + ilocerrcnt write(6,'("...done with ",i3," errors")') ilocerrcnt write(19,'("Number of tests performed:",i6)') iptcnt write(19,'("Number of errors found:",i6)') ierrcnt print*, ' ' print*, 'Number of tests performed:',iptcnt print*, 'Number of errors found:',ierrcnt print*, ' ' print*, ' ' print*, 'End of run' print*, ' ' print*, ' ' close(21) 1500 format(' The number of points tested is ',i4) 1501 format(' The number of errors found is ',i4//) 1000 format(2x,'Computing',2x,A3,2x,'at x=',f9.5) 1001 format(17x,'We computed ',e15.8,' in single precision'/,17x,'We computed ', & e23.16,' in double precision',/,5x,'and the stored value is ',e23.16/) 1002 format(' Testing the ',i2,'th zero of ',A3) 1008 format(17x,'We computed ',e23.16,' in double precision',/,5x,'and the stored value is ',e23.16/) 1003 format(' Testing the associated function value of ',A3,' at the ',i2,' zero') 1004 format(' Testing the ',i2,'th complex zero of ',A3) 1005 format(' Testing the associated function value of ',A3,' at the ',i2,' complex zero') 1006 format(' Testing the modulus ',A4,' at x= ',f9.5) 1007 format(' Testing the phase ',A8,' at x= ',f9.5) 2000 format(2x,'Computing',2x,A3,2x,'at z=',f9.5,' +I*(',f9.5,')') 2001 format(2x,'In single precision, we computed ',e15.8,8x,' +I*(',e15.8,')',/,2x, & 'In double precision, we computed ',e23.16,' +I*(',e23.16,')',/,11x, & 'and the stored value is ',e23.16,' +I*(',e23.16,')',/) 2002 format('Error encountered testing'/' the vectorized version of ',2x,A3, & 2x,'at (r,theta)= (',f9.5,',',f9.5,')') 2003 format(2x,'In double precision, we computed ',e23.16,' +I*(',e23.16,')',/,11x, & 'and the stored value is ',e23.16,' +I*(',e23.16,')',/) 3000 format('Error encountered testing each side of ',2x,A3,2x,'at x=',f9.5) 3001 format(12x,'Below we got ',e23.16,/8x,'and above we got ',e23.16/) 4000 format('Error encountered testing each side of ',2x,A3,2x,'at z=',f9.5, & ' +I*(',f9.5,')') 4001 format(12x,'Below we got ',e15.8,' +I*(',e15.8,')',/& 8x,'and above we got ',e15.8,' +I*(',e15.8,')'/) 4002 format(12x,'Below we got ',e23.16,' +I*(',e23.16,')',/& 8x,'and above we got ',e23.16,' +I*(',e23.16,')'/) 5000 format('Error encountered testing',2x,A3,2x,'at x=',f9.5) 5001 format(9x,'Single returned ',e15.8,/5x,'and double returned ',e23.16/) 5002 format('Error encountered testing the modulus of ',2x,A3,2x,'at x=',f9.5) 5003 format('Error encountered testing the phase of ',2x,A3,2x,'at x=',f9.5) 6000 format('Error encountered testing',2x,A3,2x,'at z=',f9.5, & ' +I*(',f9.5,')') 6001 format(8x,'Single returned ',e15.8,9x,' +I*(',e15.8,')',/& 4x,'and double returned ',e23.16,' +I*(',e23.16,')'/) end program test SHAR_EOF fi # end of overwriting check cd .. cd .. cd .. # End of shell archive exit 0