subroutine airy_ai_values ( n_data, x, ai ) !*****************************************************************************80 ! !! AIRY_AI_VALUES returns some values of the Airy Ai(x) function. ! ! Discussion: ! ! The Airy functions Ai(X) and Bi(X) are a pair of linearly independent ! solutions of the differential equation: ! ! W'' - X * W = 0 ! ! In Mathematica, the function can be evaluated by: ! ! AiryAi[x] ! ! Modified: ! ! 11 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Stephen Wolfram, ! The Mathematica Book, ! Fourth Edition, ! Wolfram Media / Cambridge University Press, 1999. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) AI, the value of the Airy AI function. ! implicit none integer, parameter :: n_max = 11 real ( kind = 8 ) ai real ( kind = 8 ), save, dimension ( n_max ) :: ai_vec = (/ & 0.3550280538878172D+00, & 0.3292031299435381D+00, & 0.3037031542863820D+00, & 0.2788064819550049D+00, & 0.2547423542956763D+00, & 0.2316936064808335D+00, & 0.2098000616663795D+00, & 0.1891624003981501D+00, & 0.1698463174443649D+00, & 0.1518868036405444D+00, & 0.1352924163128814D+00 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.0D+00, & 0.1D+00, & 0.2D+00, & 0.3D+00, & 0.4D+00, & 0.5D+00, & 0.6D+00, & 0.7D+00, & 0.8D+00, & 0.9D+00, & 1.0D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 ai = 0.0D+00 else x = x_vec(n_data) ai = ai_vec(n_data) end if return end subroutine airy_ai_prime_values ( n_data, x, aip ) !*****************************************************************************80 ! !! AIRY_AI_PRIME_VALUES returns some values of the Airy function Ai'(x). ! ! Discussion: ! ! The Airy functions Ai(X) and Bi(X) are a pair of linearly independent ! solutions of the differential equation: ! ! W'' - X * W = 0 ! ! In Mathematica, the function can be evaluated by: ! ! AiryAiPrime[x] ! ! Modified: ! ! 11 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Stephen Wolfram, ! The Mathematica Book, ! Fourth Edition, ! Wolfram Media / Cambridge University Press, 1999. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) AIP, the derivative of the Airy AI function. ! implicit none integer, parameter :: n_max = 11 real ( kind = 8 ) aip real ( kind = 8 ), save, dimension ( n_max ) :: aip_vec = (/ & -0.2588194037928068D+00, & -0.2571304219075862D+00, & -0.2524054702856195D+00, & -0.2451463642190548D+00, & -0.2358320344192082D+00, & -0.2249105326646839D+00, & -0.2127932593891585D+00, & -0.1998511915822805D+00, & -0.1864128638072717D+00, & -0.1727638434616347D+00, & -0.1591474412967932D+00 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.0D+00, & 0.1D+00, & 0.2D+00, & 0.3D+00, & 0.4D+00, & 0.5D+00, & 0.6D+00, & 0.7D+00, & 0.8D+00, & 0.9D+00, & 1.0D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 aip = 0.0D+00 else x = x_vec(n_data) aip = aip_vec(n_data) end if return end subroutine airy_bi_values ( n_data, x, bi ) !*****************************************************************************80 ! !! AIRY_BI_VALUES returns some values of the Airy Bi(x) function. ! ! Discussion: ! ! The Airy functions Ai(X) and Bi(X) are a pair of linearly independent ! solutions of the differential equation: ! ! W'' - X * W = 0 ! ! In Mathematica, the function can be evaluated by: ! ! AiryBi[x] ! ! Modified: ! ! 11 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Stephen Wolfram, ! The Mathematica Book, ! Fourth Edition, ! Wolfram Media / Cambridge University Press, 1999. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) BI, the value of the Airy BI function. ! implicit none integer, parameter :: n_max = 11 real ( kind = 8 ) bi real ( kind = 8 ), save, dimension ( n_max ) :: bi_vec = (/ & 0.6149266274460007D+00, & 0.6598616901941892D+00, & 0.7054642029186612D+00, & 0.7524855850873156D+00, & 0.8017730000135972D+00, & 0.8542770431031555D+00, & 0.9110633416949405D+00, & 0.9733286558781659D+00, & 0.1042422171231561D+01, & 0.1119872813134447D+01, & 0.1207423594952871D+01 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.0D+00, & 0.1D+00, & 0.2D+00, & 0.3D+00, & 0.4D+00, & 0.5D+00, & 0.6D+00, & 0.7D+00, & 0.8D+00, & 0.9D+00, & 1.0D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 bi = 0.0D+00 else x = x_vec(n_data) bi = bi_vec(n_data) end if return end subroutine airy_bi_prime_values ( n_data, x, bip ) !*****************************************************************************80 ! !! AIRY_BI_PRIME_VALUES returns some values of the Airy function Bi'(x). ! ! Discussion: ! ! The Airy functions Ai(X) and Bi(X) are a pair of linearly independent ! solutions of the differential equation: ! ! W'' - X * W = 0 ! ! In Mathematica, the function can be evaluated by: ! ! AiryBiPrime[x] ! ! Modified: ! ! 11 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Stephen Wolfram, ! The Mathematica Book, ! Fourth Edition, ! Wolfram Media / Cambridge University Press, 1999. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) BIP, the derivative of the Airy BI function. ! implicit none integer, parameter :: n_max = 11 real ( kind = 8 ) bip real ( kind = 8 ), save, dimension ( n_max ) :: bip_vec = (/ & 0.4482883573538264D+00, & 0.4515126311496465D+00, & 0.4617892843621509D+00, & 0.4800490287524480D+00, & 0.5072816760506224D+00, & 0.5445725641405923D+00, & 0.5931444786342857D+00, & 0.6544059191721400D+00, & 0.7300069016152518D+00, & 0.8219038903072090D+00, & 0.9324359333927756D+00 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.0D+00, & 0.1D+00, & 0.2D+00, & 0.3D+00, & 0.4D+00, & 0.5D+00, & 0.6D+00, & 0.7D+00, & 0.8D+00, & 0.9D+00, & 1.0D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 bip = 0.0D+00 else x = x_vec(n_data) bip = bip_vec(n_data) end if return end function alngam ( x ) !*****************************************************************************80 ! !! ALNGAM computes the log of the absolute value of the Gamma function. ! ! Definition: ! ! The Gamma function is defined as ! ! GAMMA(Z) = INTEGRAL ( 0 <= T < Infinity) T**(Z-1) EXP ( -T ) DT ! ! If Z is a positive integer, GAMMA(Z) = (Z-1)!, the factorial. ! ! There is a special value: ! ! GAMMA(0.5) = SQRT ( PI ). ! ! Modified: ! ! 31 May 2000 ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the argument of the gamma function. ! ! Output, real ( kind = 8 ) ALNGAM, the logarithm of the absolute ! value of GAMMA(X). ! implicit none real ( kind = 8 ) alngam real ( kind = 8 ) d9lgmc real ( kind = 8 ), save :: dxrel = 0.0D+00 real ( kind = 8 ) gamma real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) sinpiy real ( kind = 8 ), parameter :: sq2pil = 0.91893853320467274D+00 real ( kind = 8 ), parameter :: sqpi2l = 0.22579135264472743D+00 real ( kind = 8 ) x real ( kind = 8 ), save :: xmax = 0.0D+00 real ( kind = 8 ) y if ( xmax == 0.0D+00 ) then xmax = huge ( xmax ) / log ( huge ( xmax ) ) dxrel = sqrt ( epsilon ( dxrel ) ) end if y = abs ( x ) if ( y <= 10.0D+00 ) then alngam = log ( abs ( gamma ( x ) ) ) return end if if ( xmax < y ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ALNGAM - Fatal error!' write ( *, '(a)' ) ' |X| is so big that ALNGAM will overflow.' stop end if if ( 0.0D+00 < x ) then alngam = sq2pil + ( x - 0.5D+00 ) * log ( x ) - x + d9lgmc ( y ) return end if sinpiy = abs ( sin ( pi * y ) ) if ( sinpiy == 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ALNGAM - Fatal error!' write ( *, '(a)' ) ' X is a negative integer.' stop end if if ( abs ( ( x - real ( int ( x - 0.5D+00 ), kind = 8 ) ) / x ) < dxrel ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ALNGAM - Warning:' write ( *, '(a)' ) ' Answer has less than half usual precision.' write ( *, '(a)' ) ' X is very near a negative integer.' end if alngam = sqpi2l + ( x - 0.5D+00 ) * log ( y ) - x - log ( sinpiy ) & - d9lgmc ( y ) return end subroutine asyjy ( funjy, x, fnu, flgjy, in, y, wk, iflw ) !*****************************************************************************80 ! !! ASYJY computes high order Bessel functions J and Y. ! ! Description: ! ! ASYJY implements the uniform asymptotic expansion of ! the J and Y Bessel functions for 35 <= FNU and 0.0 < X. ! ! The forms are identical except for a change ! in sign of some of the terms. This change in sign is ! accomplished by means of the flag FLGJY = 1 or -1. ! ! On FLGJY = 1 the Airy functions AI(X) and DAI(X) are ! supplied by the external function JAIRY, and on ! FLGJY = -1 the Airy functions BI(X) and DBI(X) are ! supplied by the external funtion YAIRY. ! ! Modified: ! ! 25 August 2001 ! ! Author: ! ! Donald Amos ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, external FUNJY, is the function JAIRY or YAIRY. ! ! Input, real ( kind = 8 ) X, the argument, which must be greater than 0. ! ! Input, real ( kind = 8 ) FNU, the order of the first Bessel function. ! FNU is generally at least 35. ! ! Input, real ( kind = 8 ) FLGJY, a selection flag ! 1.0D+00 gives the J function ! -1.0D+00 gives the Y function ! ! Input, integer IN, the number of functions desired, which should be ! 1 or 2. ! ! Output, real ( kind = 8 ) Y(IN), contains the desired function values. ! ! Output, integer IFLW, a flag indicating underflow or overflow ! return variables for BESJ only. ! ! Output, real ( kind = 8 ) WK(7), contains the following values: ! ! wk(1) = 1 - (x/fnu)**2 = w**2 ! wk(2) = sqrt ( abs ( wk(1) ) ) ! wk(3) = abs ( wk(2) - atan ( wk(2) ) ) or ! abs ( ln((1 + wk(2) )/ ( x / fnu ) ) - wk(2)) ! = abs ( (2/3)*zeta**(3/2)) ! wk(4) = fnu*wk(3) ! wk(5) = (1.5*wk(3) * fnu)**(1/3) = sqrt ( zeta ) * fnu**(1/3) ! wk(6) = sign ( 1.0, w**2 ) * wk(5)**2 ! = sign ( 1.0, w**2 ) * zeta * fnu**(2/3) ! wk(7) = fnu**(1/3) ! implicit none real ( kind = 8 ) abw2 real ( kind = 8 ) akm real ( kind = 8 ) alfa(26,4) real ( kind = 8 ) alfa1(26,2) real ( kind = 8 ) alfa2(26,2) real ( kind = 8 ) ap real ( kind = 8 ), parameter, dimension ( 8 ) :: ar = (/ & 8.35503472222222D-02, 1.28226574556327D-01, & 2.91849026464140D-01, 8.81627267443758D-01, 3.32140828186277D+00, & 1.49957629868626D+01, 7.89230130115865D+01, 4.74451538868264D+02 /) real ( kind = 8 ) asum real ( kind = 8 ) az real ( kind = 8 ) beta(26,5) real ( kind = 8 ) beta1(26,2) real ( kind = 8 ) beta2(26,2) real ( kind = 8 ) beta3(26,1) real ( kind = 8 ) br(10) real ( kind = 8 ) bsum real ( kind = 8 ) c(65) real ( kind = 8 ), parameter :: con1 = 6.66666666666667D-01 real ( kind = 8 ), parameter :: con2 = 3.33333333333333D-01 real ( kind = 8 ), parameter :: con548 = 1.04166666666667D-01 real ( kind = 8 ) cr(10) real ( kind = 8 ) crz32 real ( kind = 8 ) d1mach real ( kind = 8 ) dfi real ( kind = 8 ) elim real ( kind = 8 ) dr(10) real ( kind = 8 ) fi real ( kind = 8 ) flgjy real ( kind = 8 ) fn real ( kind = 8 ) fnu real ( kind = 8 ) fn2 external funjy real ( kind = 8 ) gama(26) integer i integer i1mach integer iflw integer in integer j integer jn integer jr integer ju integer k integer kb integer klast integer kmax(5) integer kp1 integer ks integer ksp1 integer kstemp integer l integer lr integer lrp1 real ( kind = 8 ) phi real ( kind = 8 ) rcz real ( kind = 8 ) rden real ( kind = 8 ) relb real ( kind = 8 ) rfn2 real ( kind = 8 ) rtz real ( kind = 8 ) rzden real ( kind = 8 ) sa real ( kind = 8 ) sb real ( kind = 8 ) suma real ( kind = 8 ) sumb real ( kind = 8 ) s1 real ( kind = 8 ) ta real ( kind = 8 ) tau real ( kind = 8 ) tb real ( kind = 8 ) tfn real ( kind = 8 ) tol real ( kind = 8 ), save :: tols = -6.90775527898214D+00 real ( kind = 8 ) t2 real ( kind = 8 ) upol(10) real ( kind = 8 ) wk(*) real ( kind = 8 ) x real ( kind = 8 ) xx real ( kind = 8 ) y(*) real ( kind = 8 ) z real ( kind = 8 ) z32 equivalence (alfa(1,1),alfa1(1,1)) equivalence (alfa(1,3),alfa2(1,1)) equivalence (beta(1,1),beta1(1,1)) equivalence (beta(1,3),beta2(1,1)) equivalence (beta(1,5),beta3(1,1)) data br(1), br(2), br(3), br(4), br(5), br(6), br(7), br(8), br(9), br(10) & /-1.45833333333333D-01,-9.87413194444444D-02, & -1.43312053915895D-01,-3.17227202678414D-01,-9.42429147957120D-01, & -3.51120304082635D+00,-1.57272636203680D+01,-8.22814390971859D+01, & -4.92355370523671D+02,-3.31621856854797D+03/ data c(1), c(2), c(3), c(4), c(5), c(6), c(7), c(8), c(9), c(10), & c(11), c(12), c(13), c(14), c(15), c(16), c(17), c(18), & c(19), c(20), c(21), c(22), c(23), c(24)/ & -2.08333333333333D-01, 1.25000000000000D-01, & 3.34201388888889D-01, -4.01041666666667D-01, & 7.03125000000000D-02, -1.02581259645062D+00, & 1.84646267361111D+00, -8.91210937500000D-01, & 7.32421875000000D-02, 4.66958442342625D+00, & -1.12070026162230D+01, 8.78912353515625D+00, & -2.36408691406250D+00, 1.12152099609375D-01, & -2.82120725582002D+01, 8.46362176746007D+01, & -9.18182415432400D+01, 4.25349987453885D+01, & -7.36879435947963D+00, 2.27108001708984D-01, & 2.12570130039217D+02, -7.65252468141182D+02, & 1.05999045252800D+03, -6.99579627376133D+02/ data c(25), c(26), c(27), c(28), c(29), c(30), c(31), c(32), & c(33), c(34), c(35), c(36), c(37), c(38), c(39), c(40), & c(41), c(42), c(43), c(44), c(45), c(46), c(47), c(48)/ & 2.18190511744212D+02, -2.64914304869516D+01, & 5.72501420974731D-01, -1.91945766231841D+03, & 8.06172218173731D+03, -1.35865500064341D+04, & 1.16553933368645D+04, -5.30564697861340D+03, & 1.20090291321635D+03, -1.08090919788395D+02, & 1.72772750258446D+00, 2.02042913309661D+04, & -9.69805983886375D+04, 1.92547001232532D+05, & -2.03400177280416D+05, 1.22200464983017D+05, & -4.11926549688976D+04, 7.10951430248936D+03, & -4.93915304773088D+02, 6.07404200127348D+00, & -2.42919187900551D+05, 1.31176361466298D+06, & -2.99801591853811D+06, 3.76327129765640D+06/ data c(49), c(50), c(51), c(52), c(53), c(54), c(55), c(56), & c(57), c(58), c(59), c(60), c(61), c(62), c(63), c(64), & c(65)/ & -2.81356322658653D+06, 1.26836527332162D+06, & -3.31645172484564D+05, 4.52187689813627D+04, & -2.49983048181121D+03, 2.43805296995561D+01, & 3.28446985307204D+06, -1.97068191184322D+07, & 5.09526024926646D+07, -7.41051482115327D+07, & 6.63445122747290D+07, -3.75671766607634D+07, & 1.32887671664218D+07, -2.78561812808645D+06, & 3.08186404612662D+05, -1.38860897537170D+04, & 1.10017140269247D+02/ data alfa1(1,1), alfa1(2,1), alfa1(3,1), alfa1(4,1), alfa1(5,1), & alfa1(6,1), alfa1(7,1), alfa1(8,1), alfa1(9,1), alfa1(10,1), & alfa1(11,1),alfa1(12,1),alfa1(13,1),alfa1(14,1),alfa1(15,1), & alfa1(16,1),alfa1(17,1),alfa1(18,1),alfa1(19,1),alfa1(20,1), & alfa1(21,1),alfa1(22,1),alfa1(23,1),alfa1(24,1),alfa1(25,1), & alfa1(26,1) /-4.44444444444444D-03,-9.22077922077922D-04, & -8.84892884892885D-05, 1.65927687832450D-04, 2.46691372741793D-04, & 2.65995589346255D-04, 2.61824297061501D-04, 2.48730437344656D-04, & 2.32721040083232D-04, 2.16362485712365D-04, 2.00738858762752D-04, & 1.86267636637545D-04, 1.73060775917876D-04, 1.61091705929016D-04, & 1.50274774160908D-04, 1.40503497391270D-04, 1.31668816545923D-04, & 1.23667445598253D-04, 1.16405271474738D-04, 1.09798298372713D-04, & 1.03772410422993D-04, 9.82626078369363D-05, 9.32120517249503D-05, & 8.85710852478712D-05, 8.42963105715700D-05, 8.03497548407791D-05/ data alfa1(1,2), alfa1(2,2), alfa1(3,2), alfa1(4,2), alfa1(5,2), & alfa1(6,2), alfa1(7,2), alfa1(8,2), alfa1(9,2), alfa1(10,2), & alfa1(11,2),alfa1(12,2),alfa1(13,2),alfa1(14,2),alfa1(15,2), & alfa1(16,2),alfa1(17,2),alfa1(18,2),alfa1(19,2),alfa1(20,2), & alfa1(21,2),alfa1(22,2),alfa1(23,2),alfa1(24,2),alfa1(25,2), & alfa1(26,2) / 6.93735541354589D-04, 2.32241745182922D-04, & -1.41986273556691D-05,-1.16444931672049D-04,-1.50803558053049D-04,& -1.55121924918096D-04,-1.46809756646466D-04,-1.33815503867491D-04, & -1.19744975684254D-04,-1.06184319207974D-04,-9.37699549891194D-05, & -8.26923045588193D-05,-7.29374348155221D-05,-6.44042357721016D-05, & -5.69611566009369D-05,-5.04731044303562D-05,-4.48134868008883D-05, & -3.98688727717599D-05,-3.55400532972042D-05,-3.17414256609022D-05, & -2.83996793904175D-05,-2.54522720634871D-05,-2.28459297164725D-05, & -2.05352753106481D-05,-1.84816217627666D-05,-1.66519330021394D-05/ data alfa2(1,1), alfa2(2,1), alfa2(3,1), alfa2(4,1), alfa2(5,1), & alfa2(6,1), alfa2(7,1), alfa2(8,1), alfa2(9,1), alfa2(10,1), & alfa2(11,1),alfa2(12,1),alfa2(13,1),alfa2(14,1),alfa2(15,1), & alfa2(16,1),alfa2(17,1),alfa2(18,1),alfa2(19,1),alfa2(20,1), & alfa2(21,1),alfa2(22,1),alfa2(23,1),alfa2(24,1),alfa2(25,1), & alfa2(26,1) /-3.54211971457744D-04,-1.56161263945159D-04, & 3.04465503594936D-05, 1.30198655773243D-04, 1.67471106699712D-04, & 1.70222587683593D-04, 1.56501427608595D-04, 1.36339170977445D-04, & 1.14886692029825D-04, 9.45869093034688D-05, 7.64498419250898D-05, & 6.07570334965197D-05, 4.74394299290509D-05, 3.62757512005344D-05, & 2.69939714979225D-05, 1.93210938247939D-05, 1.30056674793963D-05, & 7.82620866744497D-06, 3.59257485819352D-06, 1.44040049814252D-07, & -2.65396769697939D-06,-4.91346867098486D-06,-6.72739296091248D-06, & -8.17269379678658D-06,-9.31304715093561D-06,-1.02011418798016D-05/ data alfa2(1,2), alfa2(2,2), alfa2(3,2), alfa2(4,2), alfa2(5,2), & alfa2(6,2), alfa2(7,2), alfa2(8,2), alfa2(9,2), alfa2(10,2), & alfa2(11,2),alfa2(12,2),alfa2(13,2),alfa2(14,2),alfa2(15,2), & alfa2(16,2),alfa2(17,2),alfa2(18,2),alfa2(19,2),alfa2(20,2), & alfa2(21,2),alfa2(22,2),alfa2(23,2),alfa2(24,2),alfa2(25,2), & alfa2(26,2) / 3.78194199201773D-04, 2.02471952761816D-04, & -6.37938506318862D-05,-2.38598230603006D-04,-3.10916256027362D-04, & -3.13680115247576D-04,-2.78950273791323D-04,-2.28564082619141D-04, & -1.75245280340847D-04,-1.25544063060690D-04,-8.22982872820208D-05, & -4.62860730588116D-05,-1.72334302366962D-05, 5.60690482304602D-06, & 2.31395443148287D-05, 3.62642745856794D-05, 4.58006124490189D-05, & 5.24595294959114D-05, 5.68396208545815D-05, 5.94349820393104D-05, & 6.06478527578422D-05, 6.08023907788436D-05, 6.01577894539460D-05, & 5.89199657344698D-05, 5.72515823777593D-05, 5.52804375585853D-05/ data beta1(1,1), beta1(2,1), beta1(3,1), beta1(4,1), beta1(5,1), & beta1(6,1), beta1(7,1), beta1(8,1), beta1(9,1), beta1(10,1), & beta1(11,1),beta1(12,1),beta1(13,1),beta1(14,1),beta1(15,1), & beta1(16,1),beta1(17,1),beta1(18,1),beta1(19,1),beta1(20,1), & beta1(21,1),beta1(22,1),beta1(23,1),beta1(24,1),beta1(25,1), & beta1(26,1) / 1.79988721413553D-02, 5.59964911064388D-03, & 2.88501402231133D-03, 1.80096606761054D-03, 1.24753110589199D-03, & 9.22878876572938D-04, 7.14430421727287D-04, 5.71787281789705D-04, & 4.69431007606482D-04, 3.93232835462917D-04, 3.34818889318298D-04, & 2.88952148495752D-04, 2.52211615549573D-04, 2.22280580798883D-04, & 1.97541838033063D-04, 1.76836855019718D-04, 1.59316899661821D-04, & 1.44347930197334D-04, 1.31448068119965D-04, 1.20245444949303D-04, & 1.10449144504599D-04, 1.01828770740567D-04, 9.41998224204238D-05, & 8.74130545753834D-05, 8.13466262162801D-05, 7.59002269646219D-05/ data beta1(1,2), beta1(2,2), beta1(3,2), beta1(4,2), beta1(5,2), & beta1(6,2), beta1(7,2), beta1(8,2), beta1(9,2), beta1(10,2), & beta1(11,2),beta1(12,2),beta1(13,2),beta1(14,2),beta1(15,2), & beta1(16,2),beta1(17,2),beta1(18,2),beta1(19,2),beta1(20,2), & beta1(21,2),beta1(22,2),beta1(23,2),beta1(24,2),beta1(25,2), & beta1(26,2) /-1.49282953213429D-03,-8.78204709546389D-04, & -5.02916549572035D-04,-2.94822138512746D-04,-1.75463996970783D-04, & -1.04008550460816D-04,-5.96141953046458D-05,-3.12038929076098D-05, & -1.26089735980230D-05,-2.42892608575730D-07, 8.05996165414274D-06, & 1.36507009262147D-05, 1.73964125472926D-05, 1.98672978842134D-05, & 2.14463263790823D-05, 2.23954659232457D-05, 2.28967783814713D-05, & 2.30785389811178D-05, 2.30321976080909D-05, 2.28236073720349D-05, & 2.25005881105292D-05, 2.20981015361991D-05, 2.16418427448104D-05, & 2.11507649256221D-05, 2.06388749782171D-05, 2.01165241997082D-05/ data beta2(1,1), beta2(2,1), beta2(3,1), beta2(4,1), beta2(5,1), & beta2(6,1), beta2(7,1), beta2(8,1), beta2(9,1), beta2(10,1), & beta2(11,1),beta2(12,1),beta2(13,1),beta2(14,1),beta2(15,1), & beta2(16,1),beta2(17,1),beta2(18,1),beta2(19,1),beta2(20,1), & beta2(21,1),beta2(22,1),beta2(23,1),beta2(24,1),beta2(25,1), & beta2(26,1) / 5.52213076721293D-04, 4.47932581552385D-04, & 2.79520653992021D-04, 1.52468156198447D-04, 6.93271105657044D-05, & 1.76258683069991D-05,-1.35744996343269D-05,-3.17972413350427D-05, & -4.18861861696693D-05,-4.69004889379141D-05,-4.87665447413787D-05, & -4.87010031186735D-05,-4.74755620890087D-05,-4.55813058138628D-05, & -4.33309644511266D-05,-4.09230193157750D-05,-3.84822638603221D-05, & -3.60857167535411D-05,-3.37793306123367D-05,-3.15888560772110D-05, & -2.95269561750807D-05,-2.75978914828336D-05,-2.58006174666884D-05, & -2.41308356761280D-05,-2.25823509518346D-05,-2.11479656768913D-05/ data beta2(1,2), beta2(2,2), beta2(3,2), beta2(4,2), beta2(5,2), & beta2(6,2), beta2(7,2), beta2(8,2), beta2(9,2), beta2(10,2), & beta2(11,2),beta2(12,2),beta2(13,2),beta2(14,2),beta2(15,2), & beta2(16,2),beta2(17,2),beta2(18,2),beta2(19,2),beta2(20,2), & beta2(21,2),beta2(22,2),beta2(23,2),beta2(24,2),beta2(25,2), & beta2(26,2) /-4.74617796559960D-04,-4.77864567147321D-04, & -3.20390228067038D-04,-1.61105016119962D-04,-4.25778101285435D-05, & 3.44571294294968D-05, 7.97092684075675D-05, 1.03138236708272D-04, & 1.12466775262204D-04, 1.13103642108481D-04, 1.08651634848774D-04, & 1.01437951597662D-04, 9.29298396593364D-05, 8.40293133016090D-05, & 7.52727991349134D-05, 6.69632521975731D-05, 5.92564547323195D-05, & 5.22169308826976D-05, 4.58539485165361D-05, 4.01445513891487D-05, & 3.50481730031328D-05, 3.05157995034347D-05, 2.64956119950516D-05, & 2.29363633690998D-05, 1.97893056664022D-05, 1.70091984636413D-05/ data beta3(1,1), beta3(2,1), beta3(3,1), beta3(4,1), beta3(5,1), & beta3(6,1), beta3(7,1), beta3(8,1), beta3(9,1), beta3(10,1), & beta3(11,1),beta3(12,1),beta3(13,1),beta3(14,1),beta3(15,1), & beta3(16,1),beta3(17,1),beta3(18,1),beta3(19,1),beta3(20,1), & beta3(21,1),beta3(22,1),beta3(23,1),beta3(24,1),beta3(25,1), & beta3(26,1) / 7.36465810572578D-04, 8.72790805146194D-04, & 6.22614862573135D-04, 2.85998154194304D-04, 3.84737672879366D-06, & -1.87906003636972D-04,-2.97603646594555D-04,-3.45998126832656D-04, & -3.53382470916038D-04,-3.35715635775049D-04,-3.04321124789040D-04, & -2.66722723047613D-04,-2.27654214122820D-04,-1.89922611854562D-04, & -1.55058918599094D-04,-1.23778240761874D-04,-9.62926147717644D-05, & -7.25178327714425D-05,-5.22070028895634D-05,-3.50347750511901D-05, & -2.06489761035552D-05,-8.70106096849767D-06, 1.13698686675100D-06, & 9.16426474122779D-06, 1.56477785428873D-05, 2.08223629482467D-05/ data gama(1), gama(2), gama(3), gama(4), gama(5), & gama(6), gama(7), gama(8), gama(9), gama(10), & gama(11), gama(12), gama(13), gama(14), gama(15), & gama(16), gama(17), gama(18), gama(19), gama(20), & gama(21), gama(22), gama(23), gama(24), gama(25), & gama(26) / 6.29960524947437D-01, 2.51984209978975D-01, & 1.54790300415656D-01, 1.10713062416159D-01, 8.57309395527395D-02, & 6.97161316958684D-02, 5.86085671893714D-02, 5.04698873536311D-02, & 4.42600580689155D-02, 3.93720661543510D-02, 3.54283195924455D-02, & 3.21818857502098D-02, 2.94646240791158D-02, 2.71581677112934D-02, & 2.51768272973862D-02, 2.34570755306079D-02, 2.19508390134907D-02, & 2.06210828235646D-02, 1.94388240897881D-02, 1.83810633800683D-02, & 1.74293213231963D-02, 1.65685837786612D-02, 1.57865285987918D-02, & 1.50729501494096D-02, 1.44193250839955D-02, 1.38184805735342D-02/ ! ! I1MACH(14) replaces I1MACH(11) in a double precision code ! I1MACH(15) replaces I1MACH(12) in a double precision code ! ta = epsilon ( ta ) tol = max ( ta, 1.0D-15 ) tb = d1mach(5) ju = i1mach(15) if ( flgjy /= 1.0D+00 ) then jr = i1mach(14) elim = 2.303D+00 * tb * ( real ( - ju ) - real ( jr ) ) else elim = 2.303D+00 * ( tb * real ( - ju ) - 3.0D+00 ) end if fn = fnu iflw = 0 do jn = 1, in xx = x / fn wk(1) = 1.0D+00 - xx * xx abw2 = abs ( wk(1) ) wk(2) = sqrt ( abw2 ) wk(7) = fn**con2 if ( 0.27750D+00 < abw2 ) then go to 80 end if ! ! Asymptotic expansion. ! ! Cases near x=fn, abs ( 1-(x/fn)**2 ) <= 0.2775 ! coefficients of asymptotic expansion by series ! ! ZETA and truncation for a(zeta) and b(zeta) series ! ! KMAX is the truncation index for a(zeta) and b(zeta) series = max ( 2, sa ) ! if ( abw2 == 0.0D+00 ) then sa = 0.0D+00 else sa = tols / log ( abw2 ) end if sb = sa do i = 1, 5 akm = max ( sa, 2.0D+00 ) kmax(i) = int ( akm ) sa = sa + sb end do kb = kmax(5) klast = kb - 1 sa = gama(kb) do k = 1, klast kb = kb - 1 sa = sa * wk(1) + gama(kb) end do z = wk(1) * sa az = abs ( z ) rtz = sqrt ( az ) wk(3) = con1 * az * rtz wk(4) = wk(3) * fn wk(5) = rtz * wk(7) wk(6) = - wk(5) * wk(5) if ( 0.0D+00 < z ) then if ( elim < wk(4) ) then iflw = 1 return end if wk(6) = -wk(6) end if phi = sqrt ( sqrt ( sa + sa + sa + sa ) ) ! ! b(zeta) for s=0 ! kb = kmax(5) klast = kb - 1 sb = beta(kb,1) do k = 1, klast kb = kb - 1 sb = sb * wk(1) + beta(kb,1) end do ksp1 = 1 fn2 = fn * fn rfn2 = 1.0D+00 / fn2 rden = 1.0D+00 asum = 1.0D+00 relb = tol * abs ( sb ) bsum = sb do ks = 1, 4 ksp1 = ksp1 + 1 rden = rden * rfn2 ! ! a(zeta) and b(zeta) for s=1,2,3,4 ! kstemp = 5 - ks kb = kmax(kstemp) klast = kb - 1 sa = alfa(kb,ks) sb = beta(kb,ksp1) do k = 1, klast kb = kb - 1 sa = sa * wk(1) + alfa(kb,ks) sb = sb * wk(1) + beta(kb,ksp1) end do ta = sa * rden tb = sb * rden asum = asum + ta bsum = bsum + tb if ( abs ( ta ) <= tol .and. abs ( tb ) <= relb ) then exit end if end do bsum = bsum / ( fn * wk(7) ) go to 160 80 continue upol(1) = 1.0D+00 tau = 1.0D+00 / wk(2) t2 = 1.0D+00 / wk(1) ! ! Cases for sqrt ( 1.2775 ) < (x/fn). ! if ( wk(1) < 0.0D+00 ) then wk(3) = abs ( wk(2) - atan ( wk(2) ) ) wk(4) = wk(3) * fn rcz = -con1 / wk(4) z32 = 1.5D+00 * wk(3) rtz = z32**con2 wk(5) = rtz * wk(7) wk(6) = -wk(5) * wk(5) ! ! Cases for (x/fn) < sqrt ( 0.7225 ) ! else wk(3) = abs ( log ( ( 1.0D+00 + wk(2) ) / xx ) - wk(2) ) wk(4) = wk(3) * fn rcz = con1 / wk(4) if ( elim < wk(4) ) then iflw = 1 return end if z32 = 1.5D+00 * wk(3) rtz = z32**con2 wk(7) = fn**con2 wk(5) = rtz * wk(7) wk(6) = wk(5) * wk(5) end if phi = sqrt ( ( rtz + rtz ) * tau ) tb = 1.0D+00 asum = 1.0D+00 tfn = tau / fn upol(2) = ( c(1) * t2 + c(2) ) * tfn crz32 = con548 * rcz bsum = upol(2) + crz32 relb = tol * abs ( bsum ) ap = tfn ks = 0 kp1 = 2 rzden = rcz l = 2 do lr = 2, 8, 2 ! ! Compute two U polynomials for next a(zeta) and b(zeta) ! lrp1 = lr + 1 do k = lr, lrp1 ks = ks + 1 kp1 = kp1 + 1 l = l + 1 s1 = c(l) do j = 2, kp1 l = l + 1 s1 = s1 * t2 + c(l) end do ap = ap * tfn upol(kp1) = ap * s1 cr(ks) = br(ks) * rzden rzden = rzden * rcz dr(ks) = ar(ks) * rzden end do suma = upol(lrp1) sumb = upol(lr+2) + upol(lrp1) * crz32 ju = lrp1 do jr = 1, lr ju = ju - 1 suma = suma + cr(jr) * upol(ju) sumb = sumb + dr(jr) * upol(ju) end do tb = -tb if ( 0.0D+00 < wk(1) ) then tb = abs ( tb ) end if asum = asum + suma * tb bsum = bsum + sumb * tb if ( abs ( suma ) <= tol .and. abs ( sumb ) <= relb ) then exit end if end do tb = wk(5) if ( 0.0D+00 < wk(1) ) then tb = -tb end if bsum = bsum / tb 160 continue call funjy ( wk(6), wk(5), wk(4), fi, dfi ) y(jn) = flgjy * phi * ( fi * asum + dfi * bsum ) / wk(7) fn = fn - flgjy end do return end subroutine bakslv ( nr, n, a, x, b ) !*****************************************************************************80 ! !! BAKSLV solves A'*x=b where A is a lower triangular matrix. ! ! Discussion: ! ! BAKSLV solves the linear equations A'*X=B, where A is a ! lower triangular matrix and A' is the transpose of A. ! ! This routine is required by the UNCMIN minimization program. ! ! If B is no longer required by calling routine, then vectors B and ! X may share the same storage, and the output value of X will ! overwrite the input value of B. ! ! Reference: ! ! John Dennis, Robert Schnabel, ! Numerical Methods for Unconstrained Optimization ! and Nonlinear Equations, ! SIAM, 1996, ! ISBN13: 978-0-898713-64-0, ! LC: QA402.5.D44. ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, integer NR, the leading dimension of A. ! ! Input, integer N, the number of rows and columns in A. ! ! Input, real ( kind = 8 ) A(NR,N), the N by N matrix, containing the lower ! triangular matrix. A is not altered by this routine. ! ! Output, real ( kind = 8 ) X(N), the solution vector. ! ! Input, real ( kind = 8 ) B(N), the right hand side vector. ! implicit none integer n integer nr real ( kind = 8 ) a(nr,n) real ( kind = 8 ) b(n) integer i integer ip1 real ( kind = 8 ) x(n) ! ! Solve L' * x = b. ! i = n x(i) = b(i) / a(i,i) if ( n == 1 ) then return end if do ip1 = i i = i - 1 x(i) = ( b(i) - dot_product ( x(ip1:n), a(ip1:n,i) ) ) / a(i,i) if ( i == 1 ) then exit end if end do return end subroutine bernstein_poly_values ( n_data, n, k, x, b ) !*****************************************************************************80 ! !! BERNSTEIN_POLY_VALUES returns some values of the Bernstein polynomials. ! ! Discussion: ! ! The Bernstein polynomials are assumed to be based on [0,1]. ! ! The formula for the Bernstein polynomials is ! ! B(N,I)(X) = [N!/(I!*(N-I)!)] * (1-X)**(N-I) * X**I ! ! In Mathematica, the function can be evaluated by: ! ! Binomial[n,i] * (1-x)^(n-i) * x^i ! ! First values: ! ! B(0,0)(X) = 1 ! ! B(1,0)(X) = 1-X ! B(1,1)(X) = X ! ! B(2,0)(X) = (1-X)**2 ! B(2,1)(X) = 2 * (1-X) * X ! B(2,2)(X) = X**2 ! ! B(3,0)(X) = (1-X)**3 ! B(3,1)(X) = 3 * (1-X)**2 * X ! B(3,2)(X) = 3 * (1-X) * X**2 ! B(3,3)(X) = X**3 ! ! B(4,0)(X) = (1-X)**4 ! B(4,1)(X) = 4 * (1-X)**3 * X ! B(4,2)(X) = 6 * (1-X)**2 * X**2 ! B(4,3)(X) = 4 * (1-X) * X**3 ! B(4,4)(X) = X**4 ! ! Special values: ! ! B(N,I)(X) has a unique maximum value at X = I/N. ! ! B(N,I)(X) has an I-fold zero at 0 and and N-I fold zero at 1. ! ! B(N,I)(1/2) = C(N,K) / 2**N ! ! For a fixed X and N, the polynomials add up to 1: ! ! Sum ( 0 <= I <= N ) B(N,I)(X) = 1 ! ! Modified: ! ! 19 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Stephen Wolfram, ! The Mathematica Book, ! Fourth Edition, ! Wolfram Media / Cambridge University Press, 1999. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, integer N, the degree of the polynomial. ! ! Output, integer K, the index of the polynomial. ! ! Output, real ( kind = 8 ) X, the argument of the polynomial. ! ! Output, real ( kind = 8 ) B, the value of the polynomial B(N,K)(X). ! implicit none integer, parameter :: n_max = 15 real ( kind = 8 ) b real ( kind = 8 ), save, dimension ( n_max ) :: b_vec = (/ & 0.1000000000000000D+01, & 0.7500000000000000D+00, & 0.2500000000000000D+00, & 0.5625000000000000D+00, & 0.3750000000000000D+00, & 0.6250000000000000D-01, & 0.4218750000000000D+00, & 0.4218750000000000D+00, & 0.1406250000000000D+00, & 0.1562500000000000D-01, & 0.3164062500000000D+00, & 0.4218750000000000D+00, & 0.2109375000000000D+00, & 0.4687500000000000D-01, & 0.3906250000000000D-02 /) integer k integer, save, dimension ( n_max ) :: k_vec = (/ & 0, & 0, 1, & 0, 1, 2, & 0, 1, 2, 3, & 0, 1, 2, 3, 4 /) integer n integer n_data integer, save, dimension ( n_max ) :: n_vec = (/ & 0, & 1, 1, & 2, 2, 2, & 3, 3, 3, 3, & 4, 4, 4, 4, 4 /) real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 n = 0 k = 0 x = 0.0D+00 b = 0.0D+00 else n = n_vec(n_data) k = k_vec(n_data) x = x_vec(n_data) b = b_vec(n_data) end if return end function besi0 ( x ) !*****************************************************************************80 ! !! BESI0 computes the hyperbolic Bessel function of the first kind, order zero. ! ! Modified: ! ! 25 August 2001 ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the argument of the Bessel function. ! ! Output, real ( kind = 8 ) BESI0, the value of the Bessel function at X. ! implicit none real ( kind = 8 ) besi0 real ( kind = 8 ) besi0e real ( kind = 8 ), parameter, dimension ( 12 ) :: bi0cs = (/ & -0.07660547252839144951D+00, 1.927337953993808270D+00, & 0.2282644586920301339D+00, 0.01304891466707290428D+00, & 0.00043442709008164874D+00, 0.00000942265768600193D+00, & 0.00000014340062895106D+00, 0.00000000161384906966D+00, & 0.00000000001396650044D+00, 0.00000000000009579451D+00, & 0.00000000000000053339D+00, 0.00000000000000000245D+00 /) real ( kind = 8 ) csevl integer inits integer, save :: nti0 = 0 real ( kind = 8 ) x real ( kind = 8 ), save :: xmax = 0.0D+00 real ( kind = 8 ), save :: xsml = 0.0D+00 real ( kind = 8 ) y if ( nti0 == 0 ) then nti0 = inits ( bi0cs, 12, 0.1D+00 * epsilon ( bi0cs ) ) xsml = 2.0D+00 * sqrt ( epsilon ( xsml ) ) xmax = log ( huge ( xmax ) ) end if y = abs ( x ) if ( y <= 3.0D+00 ) then if ( xsml < y ) then besi0 = 2.75D+00 + csevl ( y * y / 4.5D+00 - 1.0D+00, bi0cs, nti0 ) else besi0 = 1.0D+00 end if return end if if ( xmax < y ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BESI0 - Fatal error!' write ( *, '(a)' ) ' |X| is so big that BESI0 will overflow.' stop end if besi0 = exp ( y ) * besi0e ( x ) return end function besi0e ( x ) !*****************************************************************************80 ! !! BESI0E computes the scaled hyperbolic Bessel function I0(X). ! ! Discussion: ! ! BESI0E calculates the exponentially scaled modified hyperbolic ! Bessel function of the first kind of order zero for real argument X. ! ! besi0e(x) = exp ( - abs ( x ) ) * i0 ( x ). ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the argument of the Bessel function. ! ! Output, real ( kind = 8 ) BESI0E, the value of the Bessel function at X. ! implicit none real ( kind = 8 ) ai02cs(22) real ( kind = 8 ) ai0cs(21) real ( kind = 8 ) besi0e real ( kind = 8 ) bi0cs(12) real ( kind = 8 ) csevl integer inits integer, save :: ntai0 = 0 integer, save :: ntai02 = 0 integer, save :: nti0 = 0 real ( kind = 8 ) x real ( kind = 8 ), save :: xsml = 0.0D+00 real ( kind = 8 ) y data bi0cs( 1) / -0.07660547252839144951D+00 / data bi0cs( 2) / 1.927337953993808270D+00 / data bi0cs( 3) / 0.2282644586920301339D+00 / data bi0cs( 4) / 0.01304891466707290428D+00 / data bi0cs( 5) / 0.00043442709008164874D+00 / data bi0cs( 6) / 0.00000942265768600193D+00 / data bi0cs( 7) / 0.00000014340062895106D+00 / data bi0cs( 8) / 0.00000000161384906966D+00 / data bi0cs( 9) / 0.00000000001396650044D+00 / data bi0cs(10) / 0.00000000000009579451D+00 / data bi0cs(11) / 0.00000000000000053339D+00 / data bi0cs(12) / 0.00000000000000000245D+00 / data ai0cs( 1) / 0.07575994494023796D+00 / data ai0cs( 2) / 0.00759138081082334D+00 / data ai0cs( 3) / 0.00041531313389237D+00 / data ai0cs( 4) / 0.00001070076463439D+00 / data ai0cs( 5) / -0.00000790117997921D+00 / data ai0cs( 6) / -0.00000078261435014D+00 / data ai0cs( 7) / 0.00000027838499429D+00 / data ai0cs( 8) / 0.00000000825247260D+00 / data ai0cs( 9) / -0.00000001204463945D+00 / data ai0cs(10) / 0.00000000155964859D+00 / data ai0cs(11) / 0.00000000022925563D+00 / data ai0cs(12) / -0.00000000011916228D+00 / data ai0cs(13) / 0.00000000001757854D+00 / data ai0cs(14) / 0.00000000000112822D+00 / data ai0cs(15) / -0.00000000000114684D+00 / data ai0cs(16) / 0.00000000000027155D+00 / data ai0cs(17) / -0.00000000000002415D+00 / data ai0cs(18) / -0.00000000000000608D+00 / data ai0cs(19) / 0.00000000000000314D+00 / data ai0cs(20) / -0.00000000000000071D+00 / data ai0cs(21) / 0.00000000000000007D+00 / data ai02cs( 1) / 0.05449041101410882D+00 / data ai02cs( 2) / 0.00336911647825569D+00 / data ai02cs( 3) / 0.00006889758346918D+00 / data ai02cs( 4) / 0.00000289137052082D+00 / data ai02cs( 5) / 0.00000020489185893D+00 / data ai02cs( 6) / 0.00000002266668991D+00 / data ai02cs( 7) / 0.00000000339623203D+00 / data ai02cs( 8) / 0.00000000049406022D+00 / data ai02cs( 9) / 0.00000000001188914D+00 / data ai02cs(10) / -0.00000000003149915D+00 / data ai02cs(11) / -0.00000000001321580D+00 / data ai02cs(12) / -0.00000000000179419D+00 / data ai02cs(13) / 0.00000000000071801D+00 / data ai02cs(14) / 0.00000000000038529D+00 / data ai02cs(15) / 0.00000000000001539D+00 / data ai02cs(16) / -0.00000000000004151D+00 / data ai02cs(17) / -0.00000000000000954D+00 / data ai02cs(18) / 0.00000000000000382D+00 / data ai02cs(19) / 0.00000000000000176D+00 / data ai02cs(20) / -0.00000000000000034D+00 / data ai02cs(21) / -0.00000000000000027D+00 / data ai02cs(22) / 0.00000000000000003D+00 / if ( nti0 == 0 ) then nti0 = inits ( bi0cs, 12, 0.1D+00 * epsilon ( bi0cs ) ) ntai0 = inits ( ai0cs, 21, 0.1D+00 * epsilon ( ai0cs ) ) ntai02 = inits ( ai02cs, 22, 0.1D+00 * epsilon ( ai02cs ) ) xsml = 2.0D+00 * sqrt ( epsilon ( xsml ) ) end if y = abs ( x ) if ( y <= xsml ) then besi0e = 1.0D+00 else if ( y <= 3.0D+00 ) then besi0e = exp ( -y ) * & ( 2.75D+00 + csevl ( y*y/4.5D+00 - 1.0D+00, bi0cs, nti0 ) ) else if ( y <= 8.0D+00 ) then besi0e = ( 0.375D+00 + & csevl ( ( 48.0D+00 / y - 11.0D+00 ) / 5.0D+00, ai0cs, ntai0 ) ) & / sqrt ( y ) else if ( 8.0D+00 < y ) then besi0e = ( 0.375D+00 + csevl ( 16.0D+00 / y - 1.0D+00, ai02cs, ntai02 ) ) & / sqrt ( y ) end if return end subroutine besj ( x, alpha, n, y, nz ) !*****************************************************************************80 ! !! BESJ computes a sequence of J Bessel functions of increasing order. ! ! Discussion: ! ! BESJ computes an N member sequence of J Bessel functions ! ! J(ALPHA+K-1) (X) ! ! for K=1,..,N for non-negative order ALPHA and X. ! ! A combination of the power series, the asymptotic expansion for X ! to infinity and the uniform asymptotic expansion for NU to infinity ! are applied over subdivisions of the (NU,X) plane. For values of ! (NU,X) not covered by one of these formulae, the order is ! incremented or decremented by integer values into a region where ! one of the formulas apply. ! ! Backward recursion is applied to reduce orders by integer values ! except where the entire sequence lies in the oscillatory region. ! In this case forward recursion is stable and values from the ! asymptotic expansion for X to infinity start the recursion when it ! is efficient to do so. ! ! Leading terms of the series and uniform expansion are tested for ! underflow. If a sequence is requested and the last member would ! underflow, the result is set to zero and the next lower order ! tried, until a member comes on scale or all members are set ! to zero. ! ! Overflow cannot occur. ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Donald Amos, SL Daniel, MK Weston, ! CDC 6600 subroutines IBESS and JBESS for Bessel functions ! I(NU,X) and J(NU,X), X .GE. 0, NU .GE. 0 ! ACM Transactions on Mathematical Software, ! Volume 3, pages 76-92, 1977. ! ! Frank Olver, ! Tables of Bessel Functions of Moderate or Large Orders, ! NPL Mathematical Tables, Volume 6, ! Her Majesty's Stationery Office, London, 1962. ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the argument of the Bessel function. ! X must be nonnegative. ! ! Input, real ( kind = 8 ) ALPHA, the order of the first member of ! the sequence. ALPHA must be at least 0.0. ! ! Input, integer N, the number of members in the sequence, ! N must be at least 1. ! ! Output, real ( kind = 8 ) Y(N), a vector whose first N components contain ! values for J(ALPHA+K-1)(X), K=1,...,N ! ! Output, integer NZ, the number of components of Y set to zero ! due to underflow. ! ! NZ=0, normal return, computation completed ! ! NZ /= 0, Y(N-NZ+1) through Y(N) were set to 0. ! implicit none integer n real ( kind = 8 ) ak real ( kind = 8 ) akm real ( kind = 8 ) alngam real ( kind = 8 ) alpha real ( kind = 8 ) ans real ( kind = 8 ) ap real ( kind = 8 ) arg real ( kind = 8 ) coef real ( kind = 8 ) d1mach real ( kind = 8 ) dalpha real ( kind = 8 ) dfn real ( kind = 8 ) dtm real ( kind = 8 ) earg real ( kind = 8 ) elim1 real ( kind = 8 ) etx real ( kind = 8 ) fidal real ( kind = 8 ) flgjy real ( kind = 8 ) fn real ( kind = 8 ) fnf real ( kind = 8 ) fni real ( kind = 8 ) fnp1 real ( kind = 8 ) fnu real ( kind = 8 ), parameter, dimension ( 2 ) :: fnulim = (/ & 100.0D+00, 60.0D+00 /) real ( kind = 8 ) gln integer i integer i1 integer i1mach integer i2 integer ialp integer idalp integer iflw integer in integer, parameter :: inlim = 150 integer is external jairy integer k integer kk integer km integer kt integer nn integer ns integer nz real ( kind = 8 ), parameter :: pdf = 0.785398163397448D+00 real ( kind = 8 ), parameter :: pidt = 1.57079632679490D+00 real ( kind = 8 ), parameter, dimension ( 4 ) :: pp = (/ & 8.72909153935547D+00, 2.65693932265030D-01, & 1.24578576865586D-01, 7.70133747430388D-04 /) real ( kind = 8 ) rden real ( kind = 8 ) relb real ( kind = 8 ), parameter :: rttp = 7.97884560802865D-01 real ( kind = 8 ), parameter :: rtwo = 1.34839972492648D+00 real ( kind = 8 ) rtx real ( kind = 8 ) rzden real ( kind = 8 ) s real ( kind = 8 ) sa real ( kind = 8 ) sb real ( kind = 8 ) sxo2 real ( kind = 8 ) s1 real ( kind = 8 ) s2 real ( kind = 8 ) t real ( kind = 8 ) ta real ( kind = 8 ) tau real ( kind = 8 ) tb real ( kind = 8 ) temp(3) real ( kind = 8 ) tfn real ( kind = 8 ) tm real ( kind = 8 ) tol real ( kind = 8 ) tolln real ( kind = 8 ) trx real ( kind = 8 ) tx real ( kind = 8 ) t1 real ( kind = 8 ) t2 real ( kind = 8 ) wk(7) real ( kind = 8 ) x real ( kind = 8 ) xo2 real ( kind = 8 ) xo2l real ( kind = 8 ) y(n) nz = 0 kt = 1 ! ! I1MACH(14) replaces I1MACH(11) in a double precision code ! I1MACH(15) replaces I1MACH(12) in a double precision code ! ta = epsilon ( ta ) tol = max ( ta, 1.0D-15 ) i1 = i1mach(14) + 1 i2 = i1mach(15) tb = d1mach(5) elim1 = 2.303D+00 * ( real ( -i2, kind = 8 ) * tb - 3.0D+00 ) ! ! TOLLN = -ln(tol) ! tolln = 2.303D+00 * tb * real ( i1, kind = 8 ) tolln = min ( tolln, 34.5388D+00 ) if ( n < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BESJ - Fatal error!' write ( *, '(a)' ) ' N is less than 1.' return end if if ( n == 1 ) then kt = 2 end if nn = n if ( x < 0.0D+00 ) then call xerror ( 'BESJ - X less than zero.', 2, 1 ) return end if if ( x == 0.0D+00 ) then if ( alpha < 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BESJ - Fatal error!' write ( *, '(a)' ) ' ALPHA less than zero.' return end if if ( alpha == 0.0D+00 ) then y(1) = 1.0D+00 if ( n == 1 ) then return end if i1 = 2 else i1 = 1 end if y(i1:n) = 0.0D+00 return end if if ( alpha < 0.0D+00 ) then call xerror ( 'BESJ - order, alpha, less than zero.', 2, 1) return end if ialp = int ( alpha ) fni = real ( ialp + n - 1, kind = 8 ) fnf = alpha - real ( ialp, kind = 8 ) dfn = fni + fnf fnu = dfn xo2 = x * 0.5D+00 sxo2 = xo2 * xo2 ! ! Decision tree for region where series, asymptotic expansion for x ! to infinity and asymptotic expansion for nu to infinity are applied. ! if ( sxo2 <= ( fnu+1.0D+00 ) ) then go to 90 end if ta = max ( 20.0D+00, fnu ) if ( ta < x ) then go to 120 end if if ( 12.0D+00 < x ) then go to 110 end if xo2l = log ( xo2 ) ns = int ( sxo2 - fnu ) + 1 go to 100 90 continue fn = fnu fnp1 = fn + 1.0D+00 xo2l = log ( xo2 ) is = kt if ( x <= 0.50D+00 ) then go to 330 end if ns = 0 100 continue fni = fni + real ( ns, kind = 8 ) dfn = fni + fnf fn = dfn fnp1 = fn + 1.0D+00 is = kt if ( 0 < n - 1 + ns ) then is = 3 end if go to 330 110 continue ans = max ( 36.0D+00 - fnu, 0.0D+00 ) ns = int ( ans ) fni = fni + real ( ns, kind = 8 ) dfn = fni + fnf fn = dfn is = kt if ( 0 < n - 1 + ns ) then is = 3 end if go to 130 120 continue rtx = sqrt ( x ) tau = rtwo * rtx ta = tau + fnulim(kt) if ( fnu <= ta ) then go to 480 end if fn = fnu is = kt ! ! Uniform asymptotic expansion for NU to infinity. ! 130 continue i1 = abs ( 3 - is ) i1 = max ( i1, 1 ) flgjy = 1.0D+00 call asyjy ( jairy, x, fn, flgjy, i1, temp(is), wk, iflw ) if ( iflw /= 0 ) then go to 380 end if go to (320, 450, 620), is 310 continue temp(1) = temp(3) kt = 1 320 continue is = 2 fni = fni - 1.0D+00 dfn = fni + fnf fn = dfn if ( i1 == 2 ) then go to 450 end if go to 130 ! ! Series for (x/2)**2<=nu+1 ! 330 continue gln = alngam ( fnp1 ) arg = fn * xo2l - gln if ( arg < (-elim1) ) then go to 400 end if earg = exp ( arg ) 340 continue s = 1.0D+00 if ( x < tol ) then go to 360 end if ak = 3.0D+00 t2 = 1.0D+00 t = 1.0D+00 s1 = fn do k = 1, 17 s2 = t2 + s1 t = - t * sxo2 / s2 s = s + t if ( abs ( t ) < tol ) then exit end if t2 = t2 + ak ak = ak + 2.0D+00 s1 = s1 + fn end do 360 continue temp(is) = s * earg go to (370, 450, 610), is 370 continue earg = earg * fn / xo2 fni = fni - 1.0D+00 dfn = fni + fnf fn = dfn is = 2 go to 340 ! ! Set underflow value and update parameters ! 380 continue y(nn) = 0.0D+00 nn = nn - 1 fni = fni - 1.0D+00 dfn = fni + fnf fn = dfn if ( nn-1 ) 440, 390, 130 390 continue kt = 2 is = 2 go to 130 400 continue y(nn) = 0.0D+00 nn = nn - 1 fnp1 = fn fni = fni - 1.0D+00 dfn = fni + fnf fn = dfn if ( nn-1 ) 440, 410, 420 410 continue kt = 2 is = 2 420 continue if ( sxo2 <= fnp1 ) then go to 430 end if go to 130 430 continue arg = arg - xo2l + log ( fnp1 ) if ( arg < (-elim1) ) then go to 400 end if go to 330 440 nz = n - nn return ! ! Backward recursion section ! 450 continue nz = n - nn if ( kt == 2 ) then go to 470 end if ! ! Backward recur from index ALPHA+NN-1 to ALPHA. ! y(nn) = temp(1) y(nn-1) = temp(2) if ( nn == 2 ) then return end if trx = 2.0D+00 / x dtm = fni tm = ( dtm + fnf ) * trx k = nn + 1 do i = 3, nn k = k - 1 y(k-2) = tm * y(k-1) - y(k) dtm = dtm - 1.0D+00 tm = ( dtm + fnf ) * trx end do return 470 continue y(1) = temp(2) return ! ! Asymptotic expansion for X to infinity with forward recursion in ! oscillatory region max ( 20, NU ) < X, provided the last member ! of the sequence is also in the region. ! 480 continue in = int ( alpha - tau + 2.0D+00 ) if ( in <= 0 ) then go to 490 end if idalp = ialp - in - 1 kt = 1 go to 500 490 continue idalp = ialp in = 0 500 continue is = kt fidal = real ( idalp, kind = 8 ) dalpha = fidal + fnf arg = x - pidt * dalpha - pdf sa = sin ( arg ) sb = cos ( arg ) coef = rttp / rtx etx = 8.0D+00 * x 510 continue dtm = fidal + fidal dtm = dtm * dtm tm = 0.0D+00 if ( fidal == 0.0D+00 .and. abs ( fnf ) < tol ) then go to 520 end if tm = 4.0D+00 * fnf * ( fidal + fidal + fnf ) 520 continue trx = dtm - 1.0D+00 t2 = ( trx + tm ) / etx s2 = t2 relb = tol * abs ( t2 ) t1 = etx s1 = 1.0D+00 fn = 1.0D+00 ak = 8.0D+00 do k = 1, 13 t1 = t1 + etx fn = fn + ak trx = dtm - fn ap = trx + tm t2 = -t2 * ap / t1 s1 = s1 + t2 t1 = t1 + etx ak = ak + 8.0D+00 fn = fn + ak trx = dtm - fn ap = trx + tm t2 = t2 * ap / t1 s2 = s2 + t2 if ( abs ( t2 ) <= relb ) then exit end if ak = ak + 8.0D+00 end do 540 continue temp(is) = coef * ( s1 * sb - s2 * sa ) if ( is == 2 ) then go to 560 end if 550 continue fidal = fidal + 1.0D+00 dalpha = fidal + fnf is = 2 tb = sa sa = -sb sb = tb go to 510 ! ! Forward recursion section ! 560 continue if ( kt == 2 ) then go to 470 end if s1 = temp(1) s2 = temp(2) tx = 2.0D+00 / x tm = dalpha * tx if ( in == 0 ) then go to 580 end if ! ! Forward recursion to index alpha ! do i = 1, in s = s2 s2 = tm * s2 - s1 tm = tm + tx s1 = s end do if ( nn == 1 ) then go to 600 end if s = s2 s2 = tm * s2 - s1 tm = tm + tx s1 = s 580 continue ! ! Forward recursion from index ALPHA to ALPHA+N-1. ! y(1) = s1 y(2) = s2 do i = 3, nn y(i) = tm * y(i-1) - y(i-2) tm = tm + tx end do return 600 continue y(1) = s2 return ! ! Backward recursion with normalization by ! asymptotic expansion for nu to infinity or power series. ! 610 continue ! ! Computation of last order for series normalization ! akm = max ( 3.0D+00 - fn, 0.0D+00 ) km = int ( akm ) tfn = fn + real ( km, kind = 8 ) ta = ( gln + tfn - 0.9189385332D+00 - 0.0833333333D+00 / tfn ) & / ( tfn + 0.5D+00 ) ta = xo2l - ta tb = - ( 1.0D+00 -1.5D+00 / tfn ) / tfn akm = tolln / ( - ta + sqrt ( ta * ta - tolln * tb ) ) + 1.5D+00 in = km + int ( akm ) go to 660 620 continue ! ! Computation of last order for asymptotic expansion normalization ! gln = wk(3) + wk(2) if ( 30.0D+00 < wk(6) ) then go to 640 end if rden = ( pp(4) * wk(6) + pp(3) ) * wk(6) + 1.0D+00 rzden = pp(1) + pp(2) * wk(6) ta = rzden / rden if ( wk(1) < 0.10D+00 ) then go to 630 end if tb = gln / wk(5) go to 650 630 continue tb = ( 1.259921049D+00 + ( 0.1679894730D+00 + 0.0887944358D+00 & * wk(1) ) * wk(1) ) / wk(7) go to 650 640 continue ta = 0.5D+00 * tolln / wk(4) ta=( ( 0.0493827160D+00 * ta - 0.1111111111D+00 ) * ta & + 0.6666666667D+00 ) * ta * wk(6) if ( wk(1) < 0.10D+00 ) then go to 630 end if tb = gln / wk(5) 650 continue in = int ( ta / tb + 1.5D+00 ) if ( inlim < in ) then go to 310 end if 660 continue dtm = fni + real ( in, kind = 8 ) trx = 2.0D+00 / x tm = ( dtm + fnf ) * trx ta = 0.0D+00 tb = tol kk = 1 670 continue ! ! Backward recur unindexed ! do i = 1, in s = tb tb = tm * tb - ta ta = s dtm = dtm - 1.0D+00 tm = ( dtm + fnf ) * trx end do ! ! Normalization. ! if ( kk == 1 ) then ta = ( ta / tb ) * temp(3) tb = temp(3) kk = 2 in = ns if ( ns /= 0 ) then go to 670 end if end if y(nn) = tb nz = n - nn if ( nn == 1 ) then return end if k = nn - 1 y(k) = tm * tb - ta if ( nn == 2 ) then return end if dtm = dtm - 1.0D+00 tm = ( dtm + fnf ) * trx km = k - 1 ! ! Backward recur indexed ! do i = 1, km y(k-1) = tm * y(k) - y(k+1) dtm = dtm - 1.0D+00 tm = ( dtm + fnf ) * trx k = k - 1 end do return end subroutine bessel_i0_values ( n_data, x, fx ) !*****************************************************************************80 ! !! BESSEL_I0_VALUES returns some values of the I0 Bessel function. ! ! Discussion: ! ! The modified Bessel functions In(Z) and Kn(Z) are solutions of ! the differential equation ! ! Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. ! ! The modified Bessel function I0(Z) corresponds to N = 0. ! ! In Mathematica, the function can be evaluated by: ! ! BesselI[0,x] ! ! Modified: ! ! 20 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Stephen Wolfram, ! The Mathematica Book, ! Fourth Edition, ! Wolfram Media / Cambridge University Press, 1999. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.1000000000000000D+01, & 0.1010025027795146D+01, & 0.1040401782229341D+01, & 0.1092045364317340D+01, & 0.1166514922869803D+01, & 0.1266065877752008D+01, & 0.1393725584134064D+01, & 0.1553395099731217D+01, & 0.1749980639738909D+01, & 0.1989559356618051D+01, & 0.2279585302336067D+01, & 0.3289839144050123D+01, & 0.4880792585865024D+01, & 0.7378203432225480D+01, & 0.1130192195213633D+02, & 0.1748117185560928D+02, & 0.2723987182360445D+02, & 0.6723440697647798D+02, & 0.4275641157218048D+03, & 0.2815716628466254D+04 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.00D+00, & 0.20D+00, & 0.40D+00, & 0.60D+00, & 0.80D+00, & 0.10D+01, & 0.12D+01, & 0.14D+01, & 0.16D+01, & 0.18D+01, & 0.20D+01, & 0.25D+01, & 0.30D+01, & 0.35D+01, & 0.40D+01, & 0.45D+01, & 0.50D+01, & 0.60D+01, & 0.80D+01, & 0.10D+02 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_j0_values ( n_data, x, fx ) !*****************************************************************************80 ! !! BESSEL_J0_VALUES returns some values of the J0 Bessel function. ! ! Discussion: ! ! In Mathematica, the function can be evaluated by: ! ! BesselJ[0,x] ! ! Modified: ! ! 10 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Stephen Wolfram, ! The Mathematica Book, ! Fourth Edition, ! Wolfram Media / Cambridge University Press, 1999. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 21 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & -0.1775967713143383D+00, & -0.3971498098638474D+00, & -0.2600519549019334D+00, & 0.2238907791412357D+00, & 0.7651976865579666D+00, & 0.1000000000000000D+01, & 0.7651976865579666D+00, & 0.2238907791412357D+00, & -0.2600519549019334D+00, & -0.3971498098638474D+00, & -0.1775967713143383D+00, & 0.1506452572509969D+00, & 0.3000792705195556D+00, & 0.1716508071375539D+00, & -0.9033361118287613D-01, & -0.2459357644513483D+00, & -0.1711903004071961D+00, & 0.4768931079683354D-01, & 0.2069261023770678D+00, & 0.1710734761104587D+00, & -0.1422447282678077D-01 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & -5.0D+00, & -4.0D+00, & -3.0D+00, & -2.0D+00, & -1.0D+00, & 0.0D+00, & 1.0D+00, & 2.0D+00, & 3.0D+00, & 4.0D+00, & 5.0D+00, & 6.0D+00, & 7.0D+00, & 8.0D+00, & 9.0D+00, & 10.0D+00, & 11.0D+00, & 12.0D+00, & 13.0D+00, & 14.0D+00, & 15.0D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_j1_values ( n_data, x, fx ) !*****************************************************************************80 ! !! BESSEL_J1_VALUES returns some values of the J1 Bessel function. ! ! Discussion: ! ! In Mathematica, the function can be evaluated by: ! ! BesselJ[1,x] ! ! Modified: ! ! 12 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Stephen Wolfram, ! The Mathematica Book, ! Fourth Edition, ! Wolfram Media / Cambridge University Press, 1999. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 21 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.3275791375914652D+00, & 0.6604332802354914D-01, & -0.3390589585259365D+00, & -0.5767248077568734D+00, & -0.4400505857449335D+00, & 0.0000000000000000D+00, & 0.4400505857449335D+00, & 0.5767248077568734D+00, & 0.3390589585259365D+00, & -0.6604332802354914D-01, & -0.3275791375914652D+00, & -0.2766838581275656D+00, & -0.4682823482345833D-02, & 0.2346363468539146D+00, & 0.2453117865733253D+00, & 0.4347274616886144D-01, & -0.1767852989567215D+00, & -0.2234471044906276D+00, & -0.7031805212177837D-01, & 0.1333751546987933D+00, & 0.2051040386135228D+00 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & -5.0D+00, & -4.0D+00, & -3.0D+00, & -2.0D+00, & -1.0D+00, & 0.0D+00, & 1.0D+00, & 2.0D+00, & 3.0D+00, & 4.0D+00, & 5.0D+00, & 6.0D+00, & 7.0D+00, & 8.0D+00, & 9.0D+00, & 10.0D+00, & 11.0D+00, & 12.0D+00, & 13.0D+00, & 14.0D+00, & 15.0D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_jn_values ( n_data, nu, x, fx ) !*****************************************************************************80 ! !! BESSEL_JN_VALUES returns some values of the Jn Bessel function. ! ! Discussion: ! ! In Mathematica, the function can be evaluated by: ! ! BesselJ[n,x] ! ! Modified: ! ! 29 April 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Stephen Wolfram, ! The Mathematica Book, ! Fourth Edition, ! Wolfram Media / Cambridge University Press, 1999. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, integer NU, the order of the function. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.1149034849319005D+00, & 0.3528340286156377D+00, & 0.4656511627775222D-01, & 0.2546303136851206D+00, & -0.5971280079425882D-01, & 0.2497577302112344D-03, & 0.7039629755871685D-02, & 0.2611405461201701D+00, & -0.2340615281867936D+00, & -0.8140024769656964D-01, & 0.2630615123687453D-09, & 0.2515386282716737D-06, & 0.1467802647310474D-02, & 0.2074861066333589D+00, & -0.1138478491494694D+00, & 0.3873503008524658D-24, & 0.3918972805090754D-18, & 0.2770330052128942D-10, & 0.1151336924781340D-04, & -0.1167043527595797D+00 /) integer n_data integer nu integer, save, dimension ( n_max ) :: nu_vec = (/ & 2, 2, 2, 2, & 2, 5, 5, 5, & 5, 5, 10, 10, & 10, 10, 10, 20, & 20, 20, 20, 20 /) real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 nu = 0 x = 0.0D+00 fx = 0.0D+00 else nu = nu_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bp01 ( n, x, b ) !*****************************************************************************80 ! !! BP01 evaluates the N+1 Bernstein basis functions of degree N on [0,1]. ! ! Definition: ! ! The I-th Bernstein basis polynomial of degree N is defined as: ! ! B(N,I,X)= N!/(I!*(N-I)!) * (1-X)**(N-I) * X**I ! ! although this is not how the values are computed. ! ! Modified: ! ! 08 February 2003 ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, integer N, should be 0 or greater. ! ! Input, real ( kind = 8 ) X, the point where the functions should be ! evaluated. ! ! Output, real ( kind = 8 ) B(0:N), the values of the Bernstein polynomials ! at the point X. ! implicit none integer n real ( kind = 8 ) b(0:n) integer i integer j real ( kind = 8 ) x if ( n == 0 ) then b(0) = 1.0D+00 else if ( 0 < n ) then do i = 1, n if ( i == 1 ) then b(1) = x else b(i) = x * b(i-1) end if do j = i-1, 1, -1 b(j) = x * b(j-1) + ( 1.0D+00 - x ) * b(j) end do if ( i == 1 ) then b(0) = 1.0D+00 - x else b(0) = ( 1.0D+00 - x ) * b(0) end if end do end if return end subroutine c8vec_print_some ( n, x, i_lo, i_hi, title ) !*****************************************************************************80 ! !! C8VEC_PRINT_SOME prints some of a C8VEC. ! ! Discussion: ! ! A C8VEC is a vector of complex ( kind = 8 ) values. ! ! Modified: ! ! 18 October 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of entries of the vector. ! ! Input, complex ( kind = 8 ) X(N), the vector to be printed. ! ! Input, integer I_LO, I_HI, the first and last entries to print. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer n integer i integer i_hi integer i_lo character ( len = * ) title complex ( kind = 8 ) x(n) if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if write ( *, '(a)' ) ' ' do i = max ( 1, i_lo ), min ( n, i_hi ) write ( *, '(2x,i8,2x,2g14.6)' ) i, x(i) end do return end subroutine c8vec_uniform_01 ( n, seed, c ) !*****************************************************************************80 ! !! C8VEC_UNIFORM_01 returns a unit pseudorandom C8VEC. ! ! Discussion: ! ! A C8VEC is a vector of complex ( kind = 8 ) values. ! ! For now, the input quantity SEED is an integer ( kind = 4 ) variable. ! ! The angles should be uniformly distributed between 0 and 2 * PI, ! the square roots of the radius uniformly distributed between 0 and 1. ! ! This results in a uniform distribution of values in the unit circle. ! ! Modified: ! ! 15 March 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of values to compute. ! ! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which should NOT be 0. ! On output, SEED has been updated. ! ! Output, complex ( kind = 8 ) C(N), the pseudorandom complex vector. ! implicit none integer ( kind = 4 ) n complex ( kind = 8 ) c(n) integer ( kind = 4 ) i real ( kind = 8 ) r integer ( kind = 4 ) k real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 integer ( kind = 4 ) seed real ( kind = 8 ) theta if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'C8VEC_UNIFORM_01 - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop end if do i = 1, n k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + 2147483647 end if r = sqrt ( real ( seed, kind = 8 ) * 4.656612875D-10 ) k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + 2147483647 end if theta = 2.0D+00 * pi * ( real ( seed, kind = 8 ) * 4.656612875D-10 ) c(i) = r * cmplx ( cos ( theta ), sin ( theta ), kind = 8 ) end do return end subroutine chfdv ( x1, x2, f1, f2, d1, d2, ne, xe, fe, de, next, ierr ) !*****************************************************************************80 ! !! CHFDV evaluates a cubic polynomial and its derivative given in Hermite form. ! ! Discussion: ! ! CHFDV evaluates a cubic polynomial and its first derivative. ! The cubic polynomial is given in Hermite form. The evaluation ! is carried out at an array of points. ! ! This routine was designed for use by PCHFD, but it may also be ! useful directly as an evaluator for a piecewise cubic Hermite ! function in applications, such as graphing, where the interval ! is known in advance. ! ! If only function values are required, use CHFEV instead. ! ! Author: ! ! Fred Fritsch, ! Mathematics and Statistics Division, ! Lawrence Livermore National Laboratory. ! ! FORTRAN90 translation by John Burkardt. ! ! Reference: ! ! Fred Fritsch, Ralph Carlson, ! Monotone Piecewise Cubic Interpolation, ! SIAM Journal on Numerical Analysis, ! Volume 17, Number 2, April 1980, pages 238-246. ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, real ( kind = 8 ) X1, X2, the endpoints of the interval of ! definition of the cubic. X1 and X2 must be distinct. ! ! Input, real ( kind = 8 ) F1, F2, the values of the function at X1 and ! X2, respectively. ! ! Input, real ( kind = 8 ) D1, D2, the derivative values at the ends ! of the interval. ! ! Input, integer NE, the number of evaluation points. ! ! Input, real ( kind = 8 ) XE(NE), the points at which the functions are to ! be evaluated. If any of the XE are outside the interval ! [X1,X2], a warning error is returned in next. ! ! Output, real ( kind = 8 ) FE(NE), DE(NE), the values of the cubic ! function and its derivative at the points XE(*). ! ! Output, integer NEXT(2), indicates the number of extrapolation points: ! NEXT(1) = number of evaluation points to left of interval. ! NEXT(2) = number of evaluation points to right of interval. ! ! Output, integer IERR, error flag. ! 0, no errors. ! -1, NE < 1. ! -2, X1 == X2. ! implicit none integer ne real ( kind = 8 ) c2 real ( kind = 8 ) c2t2 real ( kind = 8 ) c3 real ( kind = 8 ) c3t3 real ( kind = 8 ) d1 real ( kind = 8 ) d2 real ( kind = 8 ) de(ne) real ( kind = 8 ) del1 real ( kind = 8 ) del2 real ( kind = 8 ) delta real ( kind = 8 ) f1 real ( kind = 8 ) f2 real ( kind = 8 ) fe(ne) real ( kind = 8 ) h integer i integer ierr integer next(2) real ( kind = 8 ) x real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) xe(ne) real ( kind = 8 ) xma real ( kind = 8 ) xmi ! ! Check arguments. ! if ( ne < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHFDV - Fatal error!' write ( *, '(a)' ) ' The number of evaluation points was less than 1.' stop end if h = x2 - x1 if ( h == 0.0D+00 ) then ierr = -2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHFDV - Fatal error!' write ( *, '(a)' ) ' The interval endpoints are equal.' return end if ! ! Initialize. ! ierr = 0 next(1) = 0 next(2) = 0 xmi = min ( 0.0D+00, h ) xma = max ( 0.0D+00, h ) ! ! Compute cubic coefficients expanded about X1. ! delta = ( f2 - f1 ) / h del1 = ( d1 - delta ) / h del2 = ( d2 - delta ) / h c2 = -( del1 + del1 + del2 ) c2t2 = c2 + c2 c3 = ( del1 + del2 ) / h c3t3 = c3 + c3 + c3 ! ! Evaluation loop. ! do i = 1, ne x = xe(i) - x1 fe(i) = f1 + x * ( d1 + x * ( c2 + x * c3 ) ) de(i) = d1 + x * ( c2t2 + x * c3t3 ) ! ! Count extrapolation points. ! if ( x < xmi ) then next(1) = next(1) + 1 end if if ( xma < x ) then next(2) = next(2) + 1 end if end do return end subroutine chfev ( x1, x2, f1, f2, d1, d2, ne, xe, fe, next, ierr ) !*****************************************************************************80 ! !! CHFEV evaluates a cubic polynomial given in Hermite form. ! ! Discussion: ! ! This routine evaluates a cubic polynomial given in Hermite form at an ! array of points. While designed for use by PCHFE, it may ! be useful directly as an evaluator for a piecewise cubic ! Hermite function in applications, such as graphing, where ! the interval is known in advance. ! ! The cubic polynomial is determined by function values ! F1, F2 and derivatives D1, D2 on the interval [X1,X2]. ! ! Author: ! ! Fred Fritsch, ! Mathematics and Statistics Division, ! Lawrence Livermore National Laboratory. ! ! FORTRAN90 translation by John Burkardt. ! ! Reference: ! ! Fred Fritsch, Ralph Carlson, ! Monotone Piecewise Cubic Interpolation, ! SIAM Journal on Numerical Analysis, ! Volume 17, Number 2, April 1980, pages 238-246. ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, real ( kind = 8 ) X1, X2, the endpoints of the interval of ! definition of the cubic. X1 and X2 must be distinct. ! ! Input, real ( kind = 8 ) F1, F2, the values of the function at X1 and ! X2, respectively. ! ! Input, real ( kind = 8 ) D1, D2, the derivative values at X1 and ! X2, respectively. ! ! Input, integer NE, the number of evaluation points. ! ! Input, real ( kind = 8 ) XE(NE), the points at which the function is to ! be evaluated. If any of the XE are outside the interval ! [X1,X2], a warning error is returned in NEXT. ! ! Output, real ( kind = 8 ) FE(NE), the value of the cubic function ! at the points XE. ! ! Output, integer NEXT(2), indicates the number of extrapolation points: ! NEXT(1) = number of evaluation points to the left of interval. ! NEXT(2) = number of evaluation points to the right of interval. ! ! Output, integer IERR, error flag. ! 0, no errors. ! -1, NE < 1. ! -2, X1 == X2. ! implicit none integer ne real ( kind = 8 ) c2 real ( kind = 8 ) c3 real ( kind = 8 ) d1 real ( kind = 8 ) d2 real ( kind = 8 ) del1 real ( kind = 8 ) del2 real ( kind = 8 ) delta real ( kind = 8 ) f1 real ( kind = 8 ) f2 real ( kind = 8 ) fe(ne) real ( kind = 8 ) h integer i integer ierr integer next(2) real ( kind = 8 ) x real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) xe(ne) real ( kind = 8 ) xma real ( kind = 8 ) xmi if ( ne < 1 ) then ierr = -1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHFEV - Fatal error!' write ( *, '(a)' ) ' Number of evaluation points is less than 1.' write ( *, '(a,i6)' ) ' NE = ', ne stop end if h = x2 - x1 if ( h == 0.0D+00 ) then ierr = -2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHFEV - Fatal error!' write ( *, '(a)' ) ' The interval [X1,X2] is of zero length.' stop end if ! ! Initialize. ! ierr = 0 next(1) = 0 next(2) = 0 xmi = min ( 0.0D+00, h ) xma = max ( 0.0D+00, h ) ! ! Compute cubic coefficients expanded about X1. ! delta = ( f2 - f1 ) / h del1 = ( d1 - delta ) / h del2 = ( d2 - delta ) / h c2 = -( del1 + del1 + del2 ) c3 = ( del1 + del2 ) / h ! ! Evaluation loop. ! do i = 1, ne x = xe(i) - x1 fe(i) = f1 + x * ( d1 + x * ( c2 + x * c3 ) ) ! ! Count the extrapolation points. ! if ( x < xmi ) then next(1) = next(1) + 1 end if if ( xma < x ) then next(2) = next(2) + 1 end if end do return end function chfiv ( x1, x2, f1, f2, d1, d2, a, b, ierr ) !*****************************************************************************80 ! !! CHFIV evaluates the integral of a cubic polynomial in Hermite form. ! ! Discussion: ! ! CHFIV is called by PCHIA to evaluate the integral of a single cubic (in ! Hermite form) over an arbitrary interval (A,B). ! ! Author: ! ! Fred Fritsch, ! Mathematics and Statistics Division, ! Lawrence Livermore National Laboratory. ! ! FORTRAN90 translation by John Burkardt. ! ! Reference: ! ! Fred Fritsch, Ralph Carlson, ! Monotone Piecewise Cubic Interpolation, ! SIAM Journal on Numerical Analysis, ! Volume 17, Number 2, April 1980, pages 238-246. ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Output, real ( kind = 8 ) VALUE, the value of the requested integral. ! ! Input, real ( kind = 8 ) X1, X2, the endpoints of the interval of ! definition of the cubic. X1 and X2 must be distinct. ! ! Input, real ( kind = 8 ) F1, F2, the values of the function at X1 ! and X2, respectively. ! ! Input, real ( kind = 8 ) D1, D2, the derivative values at the ends ! of the interval. ! ! Input, real ( kind = 8 ) A, B, the endpoints of interval of integration. ! ! Output, integer IERR, error flag. ! 0, no errors. ! -1, X1 == X2. ! implicit none real ( kind = 8 ) a real ( kind = 8 ) b real ( kind = 8 ) chfiv real ( kind = 8 ) d1 real ( kind = 8 ) d2 real ( kind = 8 ) dterm real ( kind = 8 ) f1 real ( kind = 8 ) f2 real ( kind = 8 ) fterm real ( kind = 8 ) h integer ierr real ( kind = 8 ) phia1 real ( kind = 8 ) phia2 real ( kind = 8 ) phib1 real ( kind = 8 ) phib2 real ( kind = 8 ) psia1 real ( kind = 8 ) psia2 real ( kind = 8 ) psib1 real ( kind = 8 ) psib2 real ( kind = 8 ) ta1 real ( kind = 8 ) ta2 real ( kind = 8 ) tb1 real ( kind = 8 ) tb2 real ( kind = 8 ) ua1 real ( kind = 8 ) ua2 real ( kind = 8 ) ub1 real ( kind = 8 ) ub2 real ( kind = 8 ) x1 real ( kind = 8 ) x2 ! ! Check input. ! if ( x1 == x2 ) then ierr = -1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHFIV - Fatal error!' write ( *, '(a)' ) ' X1 = X2.' stop end if ierr = 0 ! ! Compute integral. ! h = x2 - x1 ta1 = ( a - x1 ) / h ta2 = ( x2 - a ) / h tb1 = ( b - x1 ) / h tb2 = ( x2 - b ) / h ua1 = ta1**3 phia1 = ua1 * ( 2.0D+00 - ta1 ) psia1 = ua1 * ( 3.0D+00 * ta1 - 4.0D+00 ) ua2 = ta2**3 phia2 = ua2 * ( 2.0D+00 - ta2) psia2 = -ua2 * ( 3.0D+00 * ta2 - 4.0D+00 ) ub1 = tb1**3 phib1 = ub1 * ( 2.0D+00 - tb1 ) psib1 = ub1 * ( 3.0D+00 * tb1 - 4.0D+00 ) ub2 = tb2**3 phib2 = ub2 * ( 2.0D+00 - tb2 ) psib2 = -ub2 * ( 3.0D+00 * tb2 - 4.0D+00 ) fterm = f1 * ( phia2 - phib2 ) + f2 * ( phib1 - phia1 ) dterm = ( d1 * ( psia2 - psib2 ) + d2 * ( psib1 - psia1 ) ) * ( h / 6.0D+00 ) chfiv = 0.5D+00 * h * ( fterm + dterm ) return end function chfmc ( d1, d2, delta ) !*****************************************************************************80 ! !! CHFMC determines the monotonicity properties of a cubic polynomial. ! ! Discussion: ! ! CHFMC is called by PCHMC to determine the monotonicity properties ! of the cubic with boundary derivative values D1, D2 and chord ! slope DELTA. ! ! Author: ! ! Fred Fritsch, ! Mathematics and Statistics Division, ! Lawrence Livermore National Laboratory. ! ! FORTRAN90 translation by John Burkardt. ! ! Reference: ! ! Fred Fritsch, Ralph Carlson, ! Monotone Piecewise Cubic Interpolation, ! SIAM Journal on Numerical Analysis, ! Volume 17, Number 2, April 1980, pages 238-246. ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, real ( kind = 8 ) D1, D2, the derivative values at the ends ! of the interval. ! ! Input, real ( kind = 8 ) DELTA, the data slope over that interval. ! ! Output, integer CHFMC, indicates the monotonicity of the cubic segment: ! -1, if function is strictly decreasing; ! 0, if function is constant; ! 1, if function is strictly increasing; ! 2, if function is non-monotonic; ! 3, if unable to determine. ! implicit none real ( kind = 8 ) a real ( kind = 8 ) b integer chfmc real ( kind = 8 ) d1 real ( kind = 8 ) d2 real ( kind = 8 ) delta real ( kind = 8 ) eps integer ismon integer itrue real ( kind = 8 ) phi eps = 10.0D+00 * epsilon ( eps ) ! ! Make the check. ! if ( delta == 0.0D+00 ) then if ( d1 == 0.0D+00 .and. d2 == 0.0D+00 ) then ismon = 0 else ismon = 2 end if else itrue = sign ( 1.0D+00, delta) a = d1 / delta b = d2 / delta if ( a < 0.0D+00 .or. b < 0.0D+00 ) then ismon = 2 else if ( a <= 3.0D+00 - eps .and. b <= 3.0D+00 -eps ) then ! ! Inside square (0,3)x(0,3) implies OK. ! ismon = itrue else if ( 4.0D+00 + eps < a .and. 4.0D+00 + eps < b ) then ! ! Outside square (0,4)x(0,4) implies nonmonotonic. ! ismon = 2 else ! ! Must check against boundary of ellipse. ! a = a - 2.0D+00 b = b - 2.0D+00 phi = ( ( a * a + b * b ) + a * b ) - 3.0D+00 if ( phi < -eps ) then ismon = itrue else if ( eps < phi ) then ismon = 2 else ! ! Too close to boundary to tell, ! in the presence of round-off errors. ! ismon = 3 end if end if end if chfmc = ismon return end subroutine chkder ( m, n, x, fvec, fjac, ldfjac, xp, fvecp, mode, err ) !*****************************************************************************80 ! !! CHKDER checks the gradients of M functions of N variables. ! ! Discussion: ! ! CHKDER checks the gradients of M nonlinear functions in N variables, ! evaluated at a point X, for consistency with the functions themselves. ! ! The user calls CHKDER twice, first with MODE = 1 and then with MODE = 2. ! ! MODE = 1. ! On input, ! X contains the point of evaluation. ! On output, ! XP is set to a neighboring point. ! ! Now the user must evaluate the function and gradients at X, and the ! function at XP. Then the subroutine is called again: ! ! MODE = 2. ! On input, ! FVEC contains the function values at X, ! FJAC contains the function gradients at X. ! FVECP contains the functions evaluated at XP. ! On output, ! ERR contains measures of correctness of the respective gradients. ! ! The subroutine does not perform reliably if cancellation or ! rounding errors cause a severe loss of significance in the ! evaluation of a function. Therefore, none of the components ! of X should be unusually small (in particular, zero) or any ! other value which may cause loss of significance. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Argonne National Laboratory, ! Argonne, Illinois. ! ! Parameters: ! ! Input, integer M, is the number of functions. ! ! Input, integer N, is the number of variables. ! ! Input, real ( kind = 8 ) X(N), the point at which the jacobian is ! to be evaluated. ! ! Input, real ( kind = 8 ) FVEC(M), is used only when MODE = 2. ! In that case, it should contain the function values at X. ! ! Input, real ( kind = 8 ) FJAC(LDFJAC,N), an M by N array. When MODE = 2, ! FJAC(I,J) should contain the value of dF(I)/dX(J). ! ! Input, integer LDFJAC, the leading dimension of the array FJAC. ! LDFJAC must be at least M. ! ! Output, real ( kind = 8 ) XP(N), on output with MODE = 1, is a ! neighboring point of X, at which the function is to be evaluated. ! ! Input, real ( kind = 8 ) FVECP(M), on input with MODE = 2, is the ! function value at XP. ! ! Input, integer MODE, should be set to 1 on the first call, and ! 2 on the second. ! ! Output, real ( kind = 8 ) ERR(M). On output when MODE = 2, ERR ! contains measures of correctness of the respective gradients. If ! there is no severe loss of significance, then if ERR(I): ! = 1.0D+00, the I-th gradient is correct, ! = 0.0D+00, the I-th gradient is incorrect. ! > 0.5D+00, the I-th gradient is probably correct. ! < 0.5D+00, the I-th gradient is probably incorrect. ! implicit none integer ldfjac integer m integer n real ( kind = 8 ) eps real ( kind = 8 ) epsf real ( kind = 8 ) epslog real ( kind = 8 ) epsmch real ( kind = 8 ) err(m) real ( kind = 8 ) fjac(ldfjac,n) real ( kind = 8 ) fvec(m) real ( kind = 8 ) fvecp(m) integer i integer j integer mode real ( kind = 8 ) temp real ( kind = 8 ) x(n) real ( kind = 8 ) xp(n) epsmch = epsilon ( epsmch ) eps = sqrt ( epsmch ) ! ! MODE = 1. ! if ( mode == 1 ) then do j = 1, n temp = eps * abs ( x(j) ) if ( temp == 0.0D+00 ) then temp = eps end if xp(j) = x(j) + temp end do ! ! MODE = 2. ! else if ( mode == 2 ) then epsf = 100.0D+00 * epsmch epslog = log10 ( eps ) err = 0.0D+00 do j = 1, n temp = abs ( x(j) ) if ( temp == 0.0D+00 ) then temp = 1.0D+00 end if err(1:m) = err(1:m) + temp * fjac(1:m,j) end do do i = 1, m temp = 1.0D+00 if ( fvec(i) /= 0.0D+00 .and. fvecp(i) /= 0.0D+00 .and. & epsf * abs ( fvec(i) ) <= abs ( fvecp(i) - fvec(i) ) ) then temp = eps * abs ( ( fvecp(i) - fvec(i) ) / eps - err(i) ) & / ( abs ( fvec(i) ) + abs ( fvecp(i) ) ) end if err(i) = 1.0D+00 if ( epsmch < temp .and. temp < eps ) then err(i) = ( log10 ( temp ) - epslog ) / epslog end if if ( eps <= temp ) then err(i) = 0.0D+00 end if end do end if return end subroutine chlhsn ( nr, n, a, epsm, sx, udiag ) !*****************************************************************************80 ! !! CHLHSN finds the L*L' decomposition of the perturbed model hessian matrix. ! ! Discussion: ! ! The perturbed model Hessian matrix has the form ! ! A + MU * I ! ! (where 0 <= MU and I is the identity matrix) which is safely ! positive definite. ! ! If A is safely positive definite upon entry, then MU=0. ! ! 1. If A has any negative diagonal elements, then choose 0 < MU ! such that the diagonal of A:=A+MU*I is all positive ! with the ratio of its smallest to largest element on the ! order of sqrt ( EPSM ). ! ! 2. A undergoes a perturbed Cholesky decomposition which ! results in an LL+ decomposition of A+D, where D is a ! non-negative diagonal matrix which is implicitly added to ! A during the decomposition if A is not positive definite. ! A is retained and not changed during this process by ! copying L into the upper triangular part of A and the ! diagonal into UDIAG. Then the Cholesky decomposition routine ! is called. On return, ADDMAX contains the maximum element of D. ! ! 3. If ADDMAX=0, A was positive definite going into step 2 ! and return is made to calling program. Otherwise, ! the minimum number SDD which must be added to the ! diagonal of A to make it safely strictly diagonally dominant ! is calculated. Since A + ADDMAX * I and A + SDD * I are safely ! positive definite, choose MU = min ( ADDMAX, SDD ) and decompose ! A + MU * I to obtain L. ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, integer NR, the row dimension of the matrix. ! ! Input, integer N, the dimension of the problem. ! ! Input/output, real A(NR,N), contains an N by N matrix. ! On input, A is the model hessian. Only the lower triangular part and ! diagonal are stored. On output, A contains the factor L of the ! LL+ decomposition of the perturbed model hessian in the lower triangular ! part and diagonal, and contains the hessian in the upper triangular part ! and UDIAG. ! ! Input, real ( kind = 8 ) EPSM, the machine epsilon. ! ! Input, real ( kind = 8 ) SX(N), the diagonal scaling matrix for X. ! ! Output, real ( kind = 8 ) UDIAG(N), the diagonal of the hessian. ! ! Local variables: ! ! tol tolerance ! diagmn minimum element on diagonal of a ! diagmx maximum element on diagonal of a ! offmax maximum off-diagonal element of a ! offrow sum of off-diagonal elements in a row of a ! evmin minimum eigenvalue of a ! evmax maximum eigenvalue of a ! implicit none integer n integer nr real ( kind = 8 ) a(nr,n) real ( kind = 8 ) addmax real ( kind = 8 ) amu real ( kind = 8 ) diagmx real ( kind = 8 ) diagmn real ( kind = 8 ) epsm real ( kind = 8 ) evmax real ( kind = 8 ) evmin integer i integer j real ( kind = 8 ) offmax real ( kind = 8 ) offrow real ( kind = 8 ) posmax real ( kind = 8 ) sdd real ( kind = 8 ) sx(n) real ( kind = 8 ) tol real ( kind = 8 ) udiag(n) ! ! Scale the hessian. ! do j = 1, n do i = j, n a(i,j) = a(i,j) / ( sx(i) * sx(j) ) end do end do ! ! Step1 ! tol = sqrt ( epsm ) diagmx = a(1,1) diagmn = a(1,1) do i = 2, n if ( a(i,i) < diagmn ) then diagmn = a(i,i) end if if ( diagmx < a(i,i) ) then diagmx = a(i,i) end if end do posmax = max ( diagmx, 0.0D+00 ) if ( diagmn <= posmax * tol ) then amu = tol * ( posmax - diagmn ) - diagmn ! ! Find the largest off-diagonal element of A. ! if ( amu == 0.0D+00 ) then offmax = 0.0D+00 do i = 2, n do j = 1, i-1 if ( offmax < abs ( a(i,j) ) ) then offmax = abs ( a(i,j) ) end if end do end do amu = offmax if ( amu == 0.0D+00 ) then amu = 1.0D+00 else amu = amu * ( 1.0D+00 + tol ) end if end if ! ! A = A + MU*I ! do i = 1, n a(i,i) = a(i,i) + amu end do diagmx = diagmx + amu end if ! ! Step2 ! ! Copy lower triangular part of A to upper triangular part ! and diagonal of A to udiag ! do j = 1, n udiag(j) = a(j,j) do i = j+1, n a(j,i) = a(i,j) end do end do call choldc ( nr, n, a, diagmx, tol, addmax ) ! ! Step3 ! ! If ADDMAX=0, A was positive definite going into step 2, ! the ll+ decomposition has been done, and we return. ! ! Otherwise, 0 < ADDMAX. perturb A so that it is safely ! diagonally dominant and find ll+ decomposition ! if ( 0.0D+00 < addmax ) then ! ! Restore original A (lower triangular part and diagonal) ! do j = 1, n a(j,j) = udiag(j) do i = j+1, n a(i,j) = a(j,i) end do end do ! ! Find SDD such that A+sdd*i is safely positive definite ! note: evmin<0 since A is not positive definite; ! evmin = 0.0D+00 evmax = a(1,1) do i = 1, n offrow = sum ( abs ( a(i,1:i-1) ) ) + sum ( abs ( a(i+1:n,i) ) ) evmin = min ( evmin, a(i,i)-offrow ) evmax = max ( evmax, a(i,i)+offrow ) end do sdd = tol * ( evmax - evmin ) - evmin ! ! Perturb A and decompose again. ! amu = min ( sdd, addmax ) do i = 1, n a(i,i) = a(i,i) + amu udiag(i) = a(i,i) end do ! ! A is now guaranteed safely positive definite ! call choldc ( nr, n, a, 0.0D+00, tol, addmax ) end if ! ! Unscale the hessian and Cholesky decomposition matrix. ! do j = 1, n a(j:n,j) = sx(j:n) * a(j:n,j) do i = 1, j-1 a(i,j) = sx(i) * sx(j) * a(i,j) end do udiag(j) = udiag(j) * sx(j) * sx(j) end do return end subroutine choldc ( nr, n, a, diagmx, tol, addmax ) !*****************************************************************************80 ! !! CHOLDC finds the perturbed L*L' decomposition of A+D. ! ! Discussion: ! ! D is a non-negative diagonal matrix added to A if ! necessary to allow the Cholesky decomposition to continue. ! ! The normal Cholesky decomposition is performed. However, if at any ! point the algorithm would attempt to set ! L(I,I) = sqrt ( TEMP ) ! with ! TEMP < TOL * DIAGMX, ! then L(I,I) is set to sqrt ( TOL * DIAGMX ) ! instead. This is equivalent to adding TOL * DIAGMX-TEMP to A(I,I) ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, integer NR, the row dimension of the matrix. ! ! Input, integer N, the dimension of the problem. ! ! Input/output, real ( kind = 8 ) A(NR,N), the N by N matrix. ! On input, the matrix for which to find the perturbed ! Cholesky decomposition. ! On output, the lower triangular part contains the L factor, ! and the diagonal of A. ! ! Input, real ( kind = 8 ) DIAGMX, the maximum diagonal element of A. ! ! Input, real ( kind = 8 ) TOL, a tolerance. ! ! Output, real ( kind = 8 ) ADDMAX, the maximum amount implicitly added to ! the diagonal of A in forming the Cholesky decomposition of A+D. ! ! Local variables: ! ! aminl smallest element allowed on diagonal of L. ! ! amnlsq =aminl**2 ! ! offmax maximum off-diagonal element in column of a ! implicit none integer n integer nr real ( kind = 8 ) a(nr,n) real ( kind = 8 ) addmax real ( kind = 8 ) aminl real ( kind = 8 ) amnlsq real ( kind = 8 ) diagmx integer i integer j integer k real ( kind = 8 ) offmax real ( kind = 8 ) sum2 real ( kind = 8 ) temp real ( kind = 8 ) tol addmax = 0.0D+00 aminl = sqrt ( diagmx * tol ) amnlsq = aminl**2 ! ! Form column J of L. ! do j = 1, n ! ! Find diagonal elements of L. ! sum2 = sum ( a(j,1:j-1)**2 ) temp = a(j,j) - sum2 if ( amnlsq <= temp ) then a(j,j) = sqrt ( temp ) ! ! Find maximum off-diagonal element in column. ! else offmax = 0.0D+00 do i = j+1, n if ( offmax < abs ( a(i,j) ) ) then offmax = abs ( a(i,j) ) end if end do if ( offmax <= amnlsq ) then offmax = amnlsq end if ! ! Add to diagonal element to allow Cholesky decomposition to continue ! a(j,j) = sqrt ( offmax ) addmax = max ( addmax, offmax - temp ) end if ! ! Find (I,J) element of lower triangular matrix. ! do i = j+1, n sum2 = 0.0D+00 do k = 1, j-1 sum2 = sum2 + a(i,k) * a(j,k) end do a(i,j) = ( a(i,j) - sum2 ) / a(j,j) end do end do return end subroutine cosqb ( n, x, wsave ) !*****************************************************************************80 ! !! COSQB computes the fast cosine transform of quarter wave data. ! ! Discussion: ! ! COSQB computes a sequence from its representation in terms of a cosine ! series with odd wave numbers. ! ! The transform is defined by: ! ! X_out(I) = sum ( 1 <= K <= N ) ! ! 4 * X_in(K) * cos ( ( 2 * K - 1 ) * ( I - 1 ) * PI / ( 2 * N ) ) ! ! COSQB is the unnormalized inverse of COSQF since a call of COSQB ! followed by a call of COSQF will multiply the input sequence X by 4*N. ! ! The array WSAVE must be initialized by calling COSQI. ! ! Modified: ! ! 09 March 2001 ! ! Author: ! ! Paul Swarztrauber, ! National Center for Atmospheric Research ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Paul Swarztrauber, ! Vectorizing the FFT's, ! in Parallel Computations (G. Rodrigue, editor), ! Academic Press, 1982, pages 51-83. ! ! Bill Buzbee, ! The SLATEC Common Math Library, ! in Sources and Development of Mathematical Software, ! W. Cowell, editor, ! Prentice Hall, 1984, pages 302-318. ! ! Parameters: ! ! Input, integer N, the length of the array X. The method is ! more efficient when N is the product of small primes. ! ! Input/output, real ( kind = 8 ) X(N). ! On input, the cosine series coefficients. ! On output, the corresponding data vector. ! ! Input, real WSAVE(3*N+15), contains data, depending on N, and ! required by the algorithm. The WSAVE array must be initialized by ! calling COSQI. A different WSAVE array must be used for each different ! value of N. ! implicit none integer n real ( kind = 8 ), parameter :: tsqrt2 = 2.82842712474619D+00 real ( kind = 8 ) wsave(3*n+15) real ( kind = 8 ) x(n) real ( kind = 8 ) x1 if ( n < 2 ) then x(1) = 4.0D+00 * x(1) else if ( n == 2 ) then x1 = 4.0D+00 * ( x(1) + x(2) ) x(2) = tsqrt2 * ( x(1) - x(2) ) x(1) = x1 else call cosqb1 ( n, x, wsave(1), wsave(n+1) ) end if return end subroutine cosqb1 ( n, x, w, xh ) !*****************************************************************************80 ! !! COSQB1 is a lower level routine used by COSQB. ! ! Modified: ! ! 09 March 2001 ! ! Author: ! ! Paul Swarztrauber, ! National Center for Atmospheric Research ! ! Parameters: ! ! Input, integer N, the length of the array. ! ! Input/output, real ( kind = 8 ) X(N). ! On input, the cosine series coefficients. ! On output, the corresponding data vector. ! ! Input, real ( kind = 8 ) W(N). ! ! Input, real ( kind = 8 ) XH(2*N+15). ! implicit none integer n integer i integer k integer kc integer ns2 real ( kind = 8 ) w(n) real ( kind = 8 ) x(n) real ( kind = 8 ) xh(2*n+15) real ( kind = 8 ) xim1 ns2 = ( n + 1 ) / 2 do i = 3, n, 2 xim1 = x(i-1) + x(i) x(i) = x(i) - x(i-1) x(i-1) = xim1 end do x(1) = x(1) + x(1) if ( mod ( n, 2 ) == 0 ) then x(n) = 2.0D+00 * x(n) end if call dfftb ( n, x, xh ) do k = 2, ns2 kc = n + 2 - k xh(k) = w(k-1) * x(kc) + w(kc-1) * x(k) xh(kc) = w(k-1) * x(k) - w(kc-1) * x(kc) end do if ( mod ( n, 2 ) == 0 ) then x(ns2+1) = w(ns2) * ( x(ns2+1) + x(ns2+1) ) end if do k = 2, ns2 kc = n + 2 - k x(k) = xh(k) + xh(kc) x(kc) = xh(k) - xh(kc) end do x(1) = 2.0D+00 * x(1) return end subroutine cosqf ( n, x, wsave ) !*****************************************************************************80 ! !! COSQF computes the fast cosine transform of quarter wave data. ! ! Discussion: ! ! COSQF computes the coefficients in a cosine series representation ! with only odd wave numbers. ! ! COSQF is the unnormalized inverse of COSQB since a call of COSQF ! followed by a call of COSQB will multiply the input sequence X ! by 4*N. ! ! The array WSAVE must be initialized by calling COSQI. ! ! The transform is defined by: ! ! X_out(I) = X_in(1) + sum ( 2 <= K <= N ) ! ! 2 * X_in(K) * cos ( ( 2 * I - 1 ) * ( K - 1 ) * PI / ( 2 * N ) ) ! ! Modified: ! ! 09 March 2001 ! ! Author: ! ! Paul Swarztrauber, ! National Center for Atmospheric Research ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Paul Swarztrauber, ! Vectorizing the FFT's, ! in Parallel Computations (G. Rodrigue, editor), ! Academic Press, 1982, pages 51-83. ! ! Bill Buzbee, ! The SLATEC Common Math Library, ! in Sources and Development of Mathematical Software, ! W. Cowell, editor, ! Prentice Hall, 1984, pages 302-318. ! ! Parameters: ! ! Input, integer N, the length of the array X. The method is ! more efficient when N is the product of small primes. ! ! Input/output, real ( kind = 8 ) X(N). ! On input, the data to be transformed. ! On output, the transformed data. ! ! Input, real ( kind = 8 ) WSAVE(3*N+15), contains data, depending on N, and ! required by the algorithm. The WSAVE array must be initialized by ! calling COSQI. A different WSAVE array must be used for each different ! value of N. ! implicit none integer n real ( kind = 8 ), parameter :: sqrt2 = 1.4142135623731D+00 real ( kind = 8 ) tsqx real ( kind = 8 ) wsave(3*n+15) real ( kind = 8 ) x(n) if ( n < 2 ) then else if ( n == 2 ) then tsqx = sqrt2 * x(2) x(2) = x(1) - tsqx x(1) = x(1) + tsqx else call cosqf1 ( n, x, wsave(1), wsave(n+1) ) end if return end subroutine cosqf1 ( n, x, w, xh ) !*****************************************************************************80 ! !! COSQF1 is a lower level routine used by COSQF. ! ! Modified: ! ! 09 March 2001 ! ! Author: ! ! Paul Swarztrauber, ! National Center for Atmospheric Research ! ! Parameters: ! ! Input, integer N, the length of the array to be transformed. ! ! Input/output, real ( kind = 8 ) X(N). ! On input, the data to be transformed. ! On output, the transformed data. ! ! Input, real ( kind = 8 ) W(N). ! ! Input, real ( kind = 8 ) XH(2*N+15). ! implicit none integer n integer i integer k integer kc integer ns2 real ( kind = 8 ) w(n) real ( kind = 8 ) x(n) real ( kind = 8 ) xh(2*n+15) real ( kind = 8 ) xim1 ns2 = ( n + 1 ) / 2 do k = 2, ns2 kc = n + 2 - k xh(k) = x(k) + x(kc) xh(kc) = x(k) - x(kc) end do if ( mod ( n, 2 ) == 0 ) then xh(ns2+1) = x(ns2+1) + x(ns2+1) end if do k = 2, ns2 kc = n+2-k x(k) = w(k-1) * xh(kc) + w(kc-1) * xh(k) x(kc) = w(k-1) * xh(k) - w(kc-1) * xh(kc) end do if ( mod ( n, 2 ) == 0 ) then x(ns2+1) = w(ns2) * xh(ns2+1) end if call dfftf ( n, x, xh ) do i = 3, n, 2 xim1 = x(i-1) - x(i) x(i) = x(i-1) + x(i) x(i-1) = xim1 end do return end subroutine cosqi ( n, wsave ) !*****************************************************************************80 ! !! COSQI initializes WSAVE, used in COSQF and COSQB. ! ! Discussion: ! ! The prime factorization of N together with a tabulation of the ! trigonometric functions are computed and stored in WSAVE. ! ! Modified: ! ! 09 March 2001 ! ! Author: ! ! Paul Swarztrauber, ! National Center for Atmospheric Research ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Paul Swarztrauber, ! Vectorizing the FFT's, ! in Parallel Computations (G. Rodrigue, editor), ! Academic Press, 1982, pages 51-83. ! ! Bill Buzbee, ! The SLATEC Common Math Library, ! in Sources and Development of Mathematical Software, ! W. Cowell, editor, ! Prentice Hall, 1984, pages 302-318. ! ! Parameters: ! ! Input, integer N, the length of the array to be transformed. The method ! is more efficient when N is the product of small primes. ! ! Output, real ( kind = 8 ) WSAVE(3*N+15), contains data, depending on N, ! and required by the COSQB and COSQF algorithms. ! implicit none integer n real ( kind = 8 ) :: pi = 3.141592653589793D+00 real ( kind = 8 ) dt integer k real ( kind = 8 ) wsave(3*n+15) dt = 0.5D+00 * pi / real ( n, kind = 8 ) do k = 1, n wsave(k) = cos ( real ( k, kind = 8 ) * dt ) end do call dffti ( n, wsave(n+1) ) return end subroutine cost ( n, x, wsave ) !*****************************************************************************80 ! !! COST computes the discrete Fourier cosine transform of an even sequence. ! ! Discussion: ! ! COST is the unnormalized inverse of itself since a call of COST ! followed by another call of COST will multiply the input sequence ! X by 2*(N-1). ! ! The array WSAVE must be initialized by calling COSTI. ! ! The transform is defined by: ! ! X_out(I) = X_in(1) + (-1) **(I-1) * X_in(N) + sum ( 2 <= K <= N-1 ) ! ! 2 * X_in(K) * cos ( ( K - 1 ) * ( I - 1 ) * PI / ( N - 1 ) ) ! ! Modified: ! ! 09 March 2001 ! ! Author: ! ! Paul Swarztrauber, ! National Center for Atmospheric Research ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Paul Swarztrauber, ! Vectorizing the FFT's, ! in Parallel Computations (G. Rodrigue, editor), ! Academic Press, 1982, pages 51-83. ! ! Bill Buzbee, ! The SLATEC Common Math Library, ! in Sources and Development of Mathematical Software, ! W. Cowell, editor, ! Prentice Hall, 1984, pages 302-318. ! ! Parameters: ! ! Input, integer N, the length of the sequence to be transformed. The ! method is more efficient when N-1 is the product of small primes. ! ! Input/output, real ( kind = 8 ) X(N). ! On input, the sequence to be transformed. ! On output, the transformed sequence. ! ! Input, real ( kind = 8 ) WSAVE(3*N+15). ! The WSAVE array must be initialized by calling COSTI. A different ! array must be used for each different value of N. ! implicit none integer n real ( kind = 8 ) c1 integer i integer k integer kc integer ns2 real ( kind = 8 ) t1 real ( kind = 8 ) t2 real ( kind = 8 ) tx2 real ( kind = 8 ) wsave(3*n+15) real ( kind = 8 ) x(n) real ( kind = 8 ) x1h real ( kind = 8 ) x1p3 real ( kind = 8 ) xi real ( kind = 8 ) xim2 ns2 = n / 2 if ( n <= 1 ) then return end if if ( n == 2 ) then x1h = x(1) + x(2) x(2) = x(1) - x(2) x(1) = x1h return end if if ( n == 3 ) then x1p3 = x(1) + x(3) tx2 = x(2) + x(2) x(2) = x(1) - x(3) x(1) = x1p3 + tx2 x(3) = x1p3 - tx2 return end if c1 = x(1) - x(n) x(1) = x(1) + x(n) do k = 2, ns2 kc = n + 1 - k t1 = x(k) + x(kc) t2 = x(k) - x(kc) c1 = c1 + wsave(kc) * t2 t2 = wsave(k) * t2 x(k) = t1 - t2 x(kc) = t1 + t2 end do if ( mod ( n, 2 ) /= 0 ) then x(ns2+1) = x(ns2+1) + x(ns2+1) end if call dfftf ( n-1, x, wsave(n+1) ) xim2 = x(2) x(2) = c1 do i = 4, n, 2 xi = x(i) x(i) = x(i-2) - x(i-1) x(i-1) = xim2 xim2 = xi end do if ( mod ( n, 2 ) /= 0 ) then x(n) = xim2 end if return end subroutine costi ( n, wsave ) !*****************************************************************************80 ! !! COSTI initializes WSAVE, used in COST. ! ! Discussion: ! ! The prime factorization of N together with a tabulation of the ! trigonometric functions are computed and stored in WSAVE. ! ! Modified: ! ! 09 March 2001 ! ! Author: ! ! Paul Swarztrauber, ! National Center for Atmospheric Research ! ! Reference: ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Paul Swarztrauber, ! Vectorizing the FFT's, ! in Parallel Computations (G. Rodrigue, editor), ! Academic Press, 1982, pages 51-83. ! ! Bill Buzbee, ! The SLATEC Common Math Library, ! in Sources and Development of Mathematical Software, ! W. Cowell, editor, ! Prentice Hall, 1984, pages 302-318. ! ! Parameters: ! ! Input, integer N, the length of the sequence to be transformed. The ! method is more efficient when N-1 is the product of small primes. ! ! Output, real ( kind = 8 ) WSAVE(3*N+15), contains data, depending on N, and ! required by the COST algorithm. ! implicit none integer n real ( kind = 8 ) :: pi = 3.141592653589793D+00 real ( kind = 8 ) dt integer k real ( kind = 8 ) wsave(3*n+15) if ( n <= 3 ) then return end if dt = pi / real ( n - 1, kind = 8 ) do k = 2, ( n / 2 ) wsave(k) = 2.0D+00 * sin ( real ( k - 1, kind = 8 ) * dt ) wsave(n+1-k) = 2.0D+00 * cos ( real ( k - 1, kind = 8 ) * dt ) end do call dffti ( n-1, wsave(n+1) ) return end function csevl ( x, cs, n ) !*****************************************************************************80 ! !! CSEVL evaluates an N term Chebyshev series. ! ! Modified: ! ! 15 April 2003 ! ! Reference: ! ! R Broucke, ! Algorithm 446, ! Communications of the ACM, ! Volume 16, page 254, 1973. ! ! Leslie Fox, Ian Parker, ! Chebyshev Polynomials in Numerical Analysis, ! Oxford Press, page 56. ! ! David Kahaner, Cleve Moler, Steven Nash, ! Numerical Methods and Software, ! Prentice Hall, 1989, ! ISBN: 0-13-627258-4, ! LC: TA345.K34. ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the argument at which the series is to be ! evaluated. X must satisfy -1.0 <= X <= 1.0. ! ! Input, real ( kind = 8 ) CS(N), the array of N terms of a Chebyshev series. ! In evaluating CS, only half the first coefficient is summed. ! ! Input, integer N, the number of terms in array CS. ! N must be at least 1, and no more than 1000. ! ! Output, real ( kind = 8 ) CSEVL, the value of the Chebyshev series. ! implicit none integer n real ( kind = 8 ) b0 real ( kind = 8 ) b1 real ( kind = 8 ) b2 real ( kind = 8 ) cs(n) real ( kind = 8 ) csevl integer i real ( kind = 8 ) x if ( n < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CSEVL - Fatal error!' write ( *, '(a)' ) ' Number of terms N is less than 1.' stop end if if ( 1000 < n ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CSEVL - Fatal error!' write ( *, '(a)' ) ' The number of terms is more than 1000.' stop end if if ( x < -1.0D+00 .or. 1.0D+00 < x ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CSEVL - Fatal error!' write ( *, '(a)' ) ' The input argument X is outside the interval [-1,1].' stop end if b1 = 0.0D+00 b0 = 0.0D+00 do i = n, 1, -1 b2 = b1 b1 = b0 b0 = 2.0D+00 * x * b1 - b2 + cs(i) end do csevl = 0.5D+00 * ( b0 - b2 ) return end subroutine d1fcn ( n, x, g ) !*****************************************************************************80 ! !! D1FCN is a dummy routine for evaluating the gradient vector. ! ! Discussion: ! ! We assume that F is a scalar function of N variables X. The routine ! is to compute the vector G where G(I) = d F/d X(I). ! ! Modified: ! ! 16 April 2001 ! ! Parameters: ! ! Input, integer N, the dimension of X, and order of A. ! ! Input, real ( kind = 8 ) X(N), the point at which the gradient ! is to be evaluated. ! ! Output, real ( kind = 8 ) G(N), the gradient vector.. ! implicit none integer n real ( kind = 8 ) g(n) real ( kind = 8 ) x(n) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'D1FCN - Fatal error!' write ( *, '(a)' ) ' This is a dummy routine.' write ( *, '(a)' ) ' The user is required to replace it with a' write ( *, '(a)' ) ' routine that computes the gradient of F.' stop end function d1mach ( i ) !*****************************************************************************80 ! !! D1MACH returns double precision real machine constants. ! ! Discussion: ! ! Assuming that the internal representation of a double precision real ! number is in base B, with T the number of base-B digits in the mantissa, ! and EMIN the smallest possible exponent and EMAX the largest possible ! exponent, then ! ! D1MACH(1) = B**(EMIN-1), the smallest positive magnitude. ! D1MACH(2) = B**EMAX*(1-B**(-T)), the largest magnitude. ! D1MACH(3) = B**(-T), the smallest relative spacing. ! D1MACH(4) = B**(1-T), the largest relative spacing. ! D1MACH(5) = log10(B). ! ! Modified: ! ! 24 April 2007 ! ! Author: ! ! Phyllis Fox, Andrew Hall, Norman Schryer ! ! Reference: ! ! Phyllis Fox, Andrew Hall, Norman Schryer, ! Algorithm 528: ! Framework for a Portable Library, ! ACM Transactions on Mathematical Software, ! Volume 4, Number 2, June 1978, page 176-188. ! ! Parameters: ! ! Input, integer I, chooses the parameter to be returned. ! 1 <= I <= 5. ! ! Output, real ( kind = 8 ) D1MACH, the value of the chosen parameter. ! implicit none real ( kind = 8 ) d1mach integer i if ( i < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'D1MACH - Fatal error!' write ( *, '(a)' ) ' The input argument I is out of bounds.' write ( *, '(a)' ) ' Legal values satisfy 1 <= I <= 5.' write ( *, '(a,i12)' ) ' I = ', i d1mach = 0.0D+00 stop else if ( i == 1 ) then d1mach = 4.450147717014403D-308 else if ( i == 2 ) then d1mach = 8.988465674311579D+307 else if ( i == 3 ) then d1mach = 1.110223024625157D-016 else if ( i == 4 ) then d1mach = 2.220446049250313D-016 else if ( i == 5 ) then d1mach = 0.301029995663981D+000 else if ( 5 < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'D1MACH - Fatal error!' write ( *, '(a)' ) ' The input argument I is out of bounds.' write ( *, '(a)' ) ' Legal values satisfy 1 <= I <= 5.' write ( *, '(a,i12)' ) ' I = ', i d1mach = 0.0D+00 stop end if return end subroutine d1mpyq ( m, n, a, lda, v, w ) !*****************************************************************************80 ! !! D1MPYQ computes A*Q, where Q is the product of Householder transformations. ! ! Discussion: ! ! Given an M by N matrix A, this subroutine computes A * Q where ! Q is the product of 2 * (N - 1) transformations ! ! GV(N-1) * ... * GV(1) * GW(1) * ... * GW(N-1) ! ! and GV(I), GW(I) are Givens rotations in the (I,N) plane which ! eliminate elements in the I-th and N-th planes, respectively. ! Q itself is not given, rather the information to recover the ! GV, GW rotations is supplied. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1 ! Argonne National Laboratory, ! Argonne, Illinois. ! ! Parameters: ! ! Input, integer M, the number of rows of A. ! ! Input, integer N, the number of columns of A. ! ! Input/output, real ( kind = 8 ) A(LDA,N), the M by N array. ! On input, the matrix A to be postmultiplied by the orthogonal matrix Q. ! On output, the value of A*Q. ! ! Input, integer LDA, the leading dimension of A, which must not ! be less than M. ! ! Input, real ( kind = 8 ) V(N), W(N), contain the information necessary ! to recover the Givens rotations GV and GW. ! implicit none integer lda integer m integer n real ( kind = 8 ) a(lda,n) real ( kind = 8 ) c integer i integer j real ( kind = 8 ) s real ( kind = 8 ) temp real ( kind = 8 ) v(n) real ( kind = 8 ) w(n) ! ! Apply the first set of Givens rotations to A. ! do j = n-1, 1, -1 if ( 1.0D+00 < abs ( v(j) ) ) then c = 1.0D+00 / v(j) s = sqrt ( 1.0D+00 - c * c ) else s = v(j) c = sqrt ( 1.0D+00 - s * s ) end if do i = 1, m temp = c * a(i,j) - s * a(i,n) a(i,n) = s * a(i,j) + c * a(i,n) a(i,j) = temp end do end do ! ! Apply the second set of Givens rotations to A. ! do j = 1, n-1 if ( 1.0D+00 < abs ( w(j) ) ) then c = 1.0D+00 / w(j) s = sqrt ( 1.0D+00 - c * c ) else s = w(j) c = sqrt ( 1.0D+00 - s * s ) end if do i = 1, m temp = c * a(i,j) + s * a(i,n) a(i,n) = - s * a(i,j) + c * a(i,n) a(i,j) = temp end do end do return end subroutine d2fcn ( nr, n, x, a ) !*****************************************************************************80 ! !! D2FCN is a dummy version of a routine that computes the second derivative. ! ! Discussion: ! ! We assume that F is a scalar function of N variables X. The routine ! is to compute the matrix H where H(I,J) = d d F / d X(I) d X(J). ! ! Modified: ! ! 16 April 2001 ! ! Parameters: ! ! Input, integer NR, the leading dimension of A, which must be at ! least N. ! ! Input, integer N, the dimension of X, and order of A. ! ! Input, real ( kind = 8 ) X(N), the point at which the Hessian matrix ! is to be evaluated. ! ! Output, real ( kind = 8 ) A(NR,N), the N by N Hessian matrix. ! implicit none integer n integer nr real ( kind = 8 ) a(nr,n) real ( kind = 8 ) x(n) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'D2FCN - Fatal error!' write ( *, '(a)' ) ' This is a dummy routine.' write ( *, '(a)' ) ' The user is required to replace it with a' write ( *, '(a)' ) ' routine that computes the Hessian matrix of F.' stop end function d9lgmc ( x ) !*****************************************************************************80 ! !! D9LGMC computes the log gamma correction factor. ! ! Discussion: ! ! The routine computes the log gamma correction factor for 10 <= X ! so that ! ! log ( gamma ( x ) ) = ! log ( sqrt ( 2 * pi ) ) + ( x - 0.5 ) * log ( x ) - x + d9lgmc ( x ) ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the argument of the log gamma function. ! X must be at least 10. ! ! Output, real ( kind = 8 ) D9LGMC, the correction. ! implicit none real ( kind = 8 ), save, dimension ( 6 ) :: algmcs = (/ & 0.166638948045186D+00, -0.0000138494817606D+00, 0.0000000098108256D+00, & -0.0000000000180912D+00, 0.0000000000000622D+00, -0.0000000000000003D+00 /) real ( kind = 8 ) arg real ( kind = 8 ) csevl real ( kind = 8 ) d9lgmc integer ( kind = 8 ) inits integer ( kind = 8 ), save :: nalgm = 0 real ( kind = 8 ) x real ( kind = 8 ), save :: xbig = 0.0D+00 real ( kind = 8 ), save :: xmax = 0.0D+00 if ( nalgm == 0 ) then nalgm = inits ( algmcs, 6, epsilon ( algmcs ) ) xbig = 1.0D+00 / sqrt ( epsilon ( xbig ) ) xmax = exp ( min ( log ( huge ( xmax ) / 12.0D+00 ), & -log ( 12.0D+00 * tiny ( xmax ) ) ) ) end if if ( x < 10.0D+00 ) then call xerror ( 'D9LGMC - 10 <= x required', 1, 2 ) else if ( x < xbig ) then arg = 2.0D+00 * ( 10.0D+00 / x )**2 - 1.0D+00 d9lgmc = csevl ( arg, algmcs, nalgm ) / x else if ( x < xmax ) then d9lgmc = 1.0D+00 / ( 12.0D+00 * x ) else d9lgmc = 0.0D+00 call xerror ( 'D9LGMC - X so big d9lgmc underflows', 2, 1) end if return end function damax ( n, x, incx ) !*****************************************************************************80 ! !! DAMAX returns the maximum absolute value of the entries in a vector. ! ! Modified: ! ! 08 April 1999 ! ! Parameters: ! ! Input, integer N, the number of entries in the vector. ! ! Input, real ( kind = 8 ) X(*), the vector to be examined. ! ! Input, integer INCX, the increment between successive entries of X. ! ! Output, real ( kind = 8 ) DAMAX, the maximum absolute value of ! an element of X. ! implicit none integer i integer incx integer ix integer n real ( kind = 8 ) damax real ( kind = 8 ) x(*) if ( n <= 0 ) then damax = 0.0D+00 else if ( n == 1 ) then damax = abs ( x(1) ) else if ( incx == 1 ) then damax = abs ( x(1) ) do i = 2, n if ( damax < abs ( x(i) ) ) then damax = abs ( x(i) ) end if end do else if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if damax = abs ( x(ix) ) ix = ix + incx do i = 2, n if ( damax < abs ( x(ix) ) ) then damax = abs ( x(ix) ) end if ix = ix + incx end do end if return end function dasum ( n, x, incx ) !*****************************************************************************80 ! !! DASUM takes the sum of the absolute values of a vector. ! ! Modified: ! ! 15 February 2001 ! ! Author: ! ! Jack Dongarra ! ! Reference: ! ! Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979. ! ! Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, ! Basic Linear Algebra Subprograms for Fortran Usage, ! Algorithm 539, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Parameters: ! ! Input, integer N, the number of entries in the vector. ! ! Input, real X(*), the vector to be examined. ! ! Input, integer INCX, the increment between successive entries of X. ! INCX must not be negative. ! ! Output, real ( kind = 8 ) DASUM, the sum of the absolute values of X. ! implicit none real ( kind = 8 ) dasum integer incx integer n real x(*) dasum = sum ( abs ( x(1:1+(n-1)*incx:incx) ) ) return end subroutine daxpy ( n, da, dx, incx, dy, incy ) !*****************************************************************************80 ! !! DAXPY computes constant times a vector plus a vector. ! ! Discussion: ! ! Uses unrolled loops for increments equal to one. ! ! Author: ! ! Jack Dongarra ! ! Reference: ! ! Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979. ! ! Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, ! Basic Linear Algebra Subprograms for Fortran Usage, ! Algorithm 539, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Parameters: ! ! Input, integer N, the number of elements in DX and DY. ! ! Input, real ( kind = 8 ) DA, the multiplier of DX. ! ! Input, real ( kind = 8 ) DX(*), the first vector. ! ! Input, integer INCX, the increment between successive entries of DX. ! ! Input/output, real ( kind = 8 ) DY(*), the second vector. ! On output, DY(*) has been replaced by DY(*) + DA * DX(*). ! ! Input, integer INCY, the increment between successive entries of DY. ! implicit none real ( kind = 8 ) da real ( kind = 8 ) dx(*) real ( kind = 8 ) dy(*) integer i integer incx integer incy integer ix integer iy integer m integer n if ( n <= 0 ) then return end if if ( da == 0.0D+00 ) then return end if ! ! Code for unequal increments or equal increments ! not equal to 1. ! if ( incx /= 1 .or. incy /= 1 ) then if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( 0 <= incy ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if do i = 1, n dy(iy) = dy(iy) + da * dx(ix) ix = ix + incx iy = iy + incy end do ! ! Code for both increments equal to 1. ! else m = mod ( n, 4 ) do i = 1, m dy(i) = dy(i) + da * dx(i) end do do i = m+1, n, 4 dy(i ) = dy(i ) + da * dx(i ) dy(i+1) = dy(i+1) + da * dx(i+1) dy(i+2) = dy(i+2) + da * dx(i+2) dy(i+3) = dy(i+3) + da * dx(i+3) end do end if return end subroutine ddcor ( dfdy, el, fa, h, impl, ipvt, matdim, miter, ml, mu, n, & nde, nq, t, users, y, yh, ywt, evalfa, save1, save2, a, d, jstate ) !*****************************************************************************80 ! !! DDCOR computes corrections to the Y array of DDRIV3. ! ! Discussion: ! ! In the case of functional iteration, update Y directly from the ! result of the last call to F. ! ! In the case of the chord method, compute the corrector error and ! solve the linear system with that as right hand side and DFDY as ! coefficient matrix, using the lu decomposition if miter is 1, 2, 4, ! or 5. ! ! Parameters: ! implicit none integer matdim integer n real ( kind = 8 ) a(matdim,*) real ( kind = 8 ) d real ( kind = 8 ) dfdy(matdim,*) real ( kind = 8 ) dnrm2 real ( kind = 8 ) el(13,12) logical evalfa external fa real ( kind = 8 ) h integer i integer i1 integer i2 integer i3 integer iflag integer impl integer ipvt(*) integer j integer jstate integer miter integer ml integer mu integer mw integer nde integer nq real ( kind = 8 ) save1(*) real ( kind = 8 ) save2(*) real ( kind = 8 ) t external users real ( kind = 8 ) y(*) real ( kind = 8 ) yh(n,*) real ( kind = 8 ) ywt(*) if ( miter == 0 ) then save1(1:n) = ( h * save2(1:n) - yh(1:n,2) - save1(1:n) ) / ywt(1:n) d = dnrm2 ( n, save1, 1 ) / sqrt ( real ( n, kind = 8 ) ) save1(1:n) = h * save2(1:n) - yh(1:n,2) else if ( miter == 1 .or. miter == 2 ) then if ( impl == 0 ) then save2(1:n) = h * save2(1:n) - yh(1:n,2) - save1(1:n) else if ( impl == 1 ) then if ( evalfa ) then call fa ( n, t, y, a, matdim, ml, mu, nde ) if ( n == 0 ) then jstate = 9 return end if else evalfa = .true. end if save2(1:n) = h * save2(1:n) do j = 1,n save2(1:n) = save2(1:n) - a(1:n,j) * ( yh(j,2) + save1(j) ) end do else if ( impl == 2 ) then if ( evalfa ) then call fa ( n, t, y, a, matdim, ml, mu, nde ) if ( n == 0 ) then jstate = 9 return end if else evalfa = .true. end if save2(1:n) = h * save2(1:n) - a(1:n,1) * (yh(i,2) + save1(1:n)) end if call dgesl ( dfdy, matdim, n, ipvt, save2, 0 ) save1(1:n) = save1(1:n) + save2(1:n) save2(1:n) = save2(1:n) / ywt(1:n) d = dnrm2 ( n, save2, 1 ) / sqrt ( real ( n, kind = 8 ) ) else if ( miter == 4 .or. miter == 5 ) then if ( impl == 0 ) then save2(1:n) = h * save2(1:n) - yh(1:n,2) - save1(1:n) else if ( impl == 1 ) then if ( evalfa ) then call fa ( n, t, y, a(ml+1,1), matdim, ml, mu, nde ) if ( n == 0 ) then jstate = 9 return end if