#! /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.2737607630134822997923617770221