slatec-bessel-cpp/fortran/dyairy.f

*DECK DYAIRY
      SUBROUTINE DYAIRY (X, RX, C, BI, DBI)
C***BEGIN PROLOGUE  DYAIRY
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DBESJ and DBESY
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (YAIRY-S, DYAIRY-D)
C***AUTHOR  Amos, D. E., (SNLA)
C           Daniel, S. L., (SNLA)
C***DESCRIPTION
C
C                  DYAIRY computes the Airy function BI(X)
C                   and its derivative DBI(X) for DASYJY
C
C                                     INPUT
C
C         X  - Argument, computed by DASYJY, X unrestricted
C        RX  - RX=SQRT(ABS(X)), computed by DASYJY
C         C  - C=2.*(ABS(X)**1.5)/3., computed by DASYJY
C
C                                    OUTPUT
C        BI  - Value of function BI(X)
C       DBI  - Value of the derivative DBI(X)
C
C***SEE ALSO  DBESJ, DBESY
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   750101  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C   910408  Updated the AUTHOR section.  (WRB)
C***END PROLOGUE  DYAIRY
C
      INTEGER I, J, M1, M1D, M2, M2D, M3, M3D, M4D, N1, N1D, N2, N2D,
     1 N3, N3D, N4D
      DOUBLE PRECISION AA,AX,BB,BI,BJN,BJP,BK1,BK2,BK3,BK4,C,CON1,CON2,
     1 CON3, CV, DAA, DBB, DBI, DBJN, DBJP, DBK1, DBK2, DBK3, DBK4, D1,
     2 D2, EX, E1, E2, FPI12, F1, F2, RTRX, RX, SPI12, S1, S2, T, TC,
     3 TEMP1, TEMP2, TT, X
      DIMENSION BK1(20), BK2(20), BK3(20), BK4(14)
      DIMENSION BJP(19), BJN(19), AA(14), BB(14)
      DIMENSION DBK1(21), DBK2(20), DBK3(20), DBK4(14)
      DIMENSION DBJP(19), DBJN(19), DAA(14), DBB(14)
      SAVE N1, N2, N3, M1, M2, M3, N1D, N2D, N3D, N4D,
     1 M1D, M2D, M3D, M4D, FPI12, SPI12, CON1, CON2, CON3,
     2 BK1, BK2, BK3, BK4, BJN, BJP, AA, BB, DBK1, DBK2, DBK3, DBK4,
     3 DBJP, DBJN, DAA, DBB
      DATA N1,N2,N3/20,19,14/
      DATA M1,M2,M3/18,17,12/
      DATA N1D,N2D,N3D,N4D/21,20,19,14/
      DATA M1D,M2D,M3D,M4D/19,18,17,12/
      DATA FPI12,SPI12,CON1,CON2,CON3/
     1 1.30899693899575D+00, 1.83259571459405D+00, 6.66666666666667D-01,
     2 7.74148278841779D+00, 3.64766105490356D-01/
      DATA BK1(1),  BK1(2),  BK1(3),  BK1(4),  BK1(5),  BK1(6),
     1     BK1(7),  BK1(8),  BK1(9),  BK1(10), BK1(11), BK1(12),
     2     BK1(13), BK1(14), BK1(15), BK1(16), BK1(17), BK1(18),
     3     BK1(19), BK1(20)/ 2.43202846447449D+00, 2.57132009754685D+00,
     4 1.02802341258616D+00, 3.41958178205872D-01, 8.41978629889284D-02,
     5 1.93877282587962D-02, 3.92687837130335D-03, 6.83302689948043D-04,
     6 1.14611403991141D-04, 1.74195138337086D-05, 2.41223620956355D-06,
     7 3.24525591983273D-07, 4.03509798540183D-08, 4.70875059642296D-09,
     8 5.35367432585889D-10, 5.70606721846334D-11, 5.80526363709933D-12,
     9 5.76338988616388D-13, 5.42103834518071D-14, 4.91857330301677D-15/
      DATA BK2(1),  BK2(2),  BK2(3),  BK2(4),  BK2(5),  BK2(6),
     1     BK2(7),  BK2(8),  BK2(9),  BK2(10), BK2(11), BK2(12),
     2     BK2(13), BK2(14), BK2(15), BK2(16), BK2(17), BK2(18),
     3     BK2(19), BK2(20)/ 5.74830555784088D-01,-6.91648648376891D-03,
     4 1.97460263052093D-03,-5.24043043868823D-04, 1.22965147239661D-04,
     5-2.27059514462173D-05, 2.23575555008526D-06, 4.15174955023899D-07,
     6-2.84985752198231D-07, 8.50187174775435D-08,-1.70400826891326D-08,
     7 2.25479746746889D-09,-1.09524166577443D-10,-3.41063845099711D-11,
     8 1.11262893886662D-11,-1.75542944241734D-12, 1.36298600401767D-13,
     9 8.76342105755664D-15,-4.64063099157041D-15, 7.78772758732960D-16/
      DATA BK3(1),  BK3(2),  BK3(3),  BK3(4),  BK3(5),  BK3(6),
     1     BK3(7),  BK3(8),  BK3(9),  BK3(10), BK3(11), BK3(12),
     2     BK3(13), BK3(14), BK3(15), BK3(16), BK3(17), BK3(18),
     3     BK3(19), BK3(20)/ 5.66777053506912D-01, 2.63672828349579D-03,
     4 5.12303351473130D-05, 2.10229231564492D-06, 1.42217095113890D-07,
     5 1.28534295891264D-08, 7.28556219407507D-10,-3.45236157301011D-10,
     6-2.11919115912724D-10,-6.56803892922376D-11,-8.14873160315074D-12,
     7 3.03177845632183D-12, 1.73447220554115D-12, 1.67935548701554D-13,
     8-1.49622868806719D-13,-5.15470458953407D-14, 8.75741841857830D-15,
     9 7.96735553525720D-15,-1.29566137861742D-16,-1.11878794417520D-15/
      DATA BK4(1),  BK4(2),  BK4(3),  BK4(4),  BK4(5),  BK4(6),
     1     BK4(7),  BK4(8),  BK4(9),  BK4(10), BK4(11), BK4(12),
     2     BK4(13), BK4(14)/ 4.85444386705114D-01,-3.08525088408463D-03,
     3 6.98748404837928D-05,-2.82757234179768D-06, 1.59553313064138D-07,
     4-1.12980692144601D-08, 9.47671515498754D-10,-9.08301736026423D-11,
     5 9.70776206450724D-12,-1.13687527254574D-12, 1.43982917533415D-13,
     6-1.95211019558815D-14, 2.81056379909357D-15,-4.26916444775176D-16/
      DATA BJP(1),  BJP(2),  BJP(3),  BJP(4),  BJP(5),  BJP(6),
     1     BJP(7),  BJP(8),  BJP(9),  BJP(10), BJP(11), BJP(12),
     2     BJP(13), BJP(14), BJP(15), BJP(16), BJP(17), BJP(18),
     3     BJP(19)         / 1.34918611457638D-01,-3.19314588205813D-01,
     4 5.22061946276114D-02, 5.28869112170312D-02,-8.58100756077350D-03,
     5-2.99211002025555D-03, 4.21126741969759D-04, 8.73931830369273D-05,
     6-1.06749163477533D-05,-1.56575097259349D-06, 1.68051151983999D-07,
     7 1.89901103638691D-08,-1.81374004961922D-09,-1.66339134593739D-10,
     8 1.42956335780810D-11, 1.10179811626595D-12,-8.60187724192263D-14,
     9-5.71248177285064D-15, 4.08414552853803D-16/
      DATA BJN(1),  BJN(2),  BJN(3),  BJN(4),  BJN(5),  BJN(6),
     1     BJN(7),  BJN(8),  BJN(9),  BJN(10), BJN(11), BJN(12),
     2     BJN(13), BJN(14), BJN(15), BJN(16), BJN(17), BJN(18),
     3     BJN(19)         / 6.59041673525697D-02,-4.24905910566004D-01,
     4 2.87209745195830D-01, 1.29787771099606D-01,-4.56354317590358D-02,
     5-1.02630175982540D-02, 2.50704671521101D-03, 3.78127183743483D-04,
     6-7.11287583284084D-05,-8.08651210688923D-06, 1.23879531273285D-06,
     7 1.13096815867279D-07,-1.46234283176310D-08,-1.11576315688077D-09,
     8 1.24846618243897D-10, 8.18334132555274D-12,-8.07174877048484D-13,
     9-4.63778618766425D-14, 4.09043399081631D-15/
      DATA AA(1),   AA(2),   AA(3),   AA(4),   AA(5),   AA(6),
     1     AA(7),   AA(8),   AA(9),   AA(10),  AA(11),  AA(12),
     2     AA(13),  AA(14) /-2.78593552803079D-01, 3.52915691882584D-03,
     3 2.31149677384994D-05,-4.71317842263560D-06, 1.12415907931333D-07,
     4 2.00100301184339D-08,-2.60948075302193D-09, 3.55098136101216D-11,
     5 3.50849978423875D-11,-5.83007187954202D-12, 2.04644828753326D-13,
     6 1.10529179476742D-13,-2.87724778038775D-14, 2.88205111009939D-15/
      DATA BB(1),   BB(2),   BB(3),   BB(4),   BB(5),   BB(6),
     1     BB(7),   BB(8),   BB(9),   BB(10),  BB(11),  BB(12),
     2     BB(13),  BB(14) /-4.90275424742791D-01,-1.57647277946204D-03,
     3 9.66195963140306D-05,-1.35916080268815D-07,-2.98157342654859D-07,
     4 1.86824767559979D-08, 1.03685737667141D-09,-3.28660818434328D-10,
     5 2.57091410632780D-11, 2.32357655300677D-12,-9.57523279048255D-13,
     6 1.20340828049719D-13, 2.90907716770715D-15,-4.55656454580149D-15/
      DATA DBK1(1), DBK1(2), DBK1(3), DBK1(4), DBK1(5), DBK1(6),
     1     DBK1(7), DBK1(8), DBK1(9), DBK1(10),DBK1(11),DBK1(12),
     2     DBK1(13),DBK1(14),DBK1(15),DBK1(16),DBK1(17),DBK1(18),
     3     DBK1(19),DBK1(20),
     4     DBK1(21)        / 2.95926143981893D+00, 3.86774568440103D+00,
     5 1.80441072356289D+00, 5.78070764125328D-01, 1.63011468174708D-01,
     6 3.92044409961855D-02, 7.90964210433812D-03, 1.50640863167338D-03,
     7 2.56651976920042D-04, 3.93826605867715D-05, 5.81097771463818D-06,
     8 7.86881233754659D-07, 9.93272957325739D-08, 1.21424205575107D-08,
     9 1.38528332697707D-09, 1.50190067586758D-10, 1.58271945457594D-11,
     1 1.57531847699042D-12, 1.50774055398181D-13, 1.40594335806564D-14,
     2 1.24942698777218D-15/
      DATA DBK2(1), DBK2(2), DBK2(3), DBK2(4), DBK2(5), DBK2(6),
     1     DBK2(7), DBK2(8), DBK2(9), DBK2(10),DBK2(11),DBK2(12),
     2     DBK2(13),DBK2(14),DBK2(15),DBK2(16),DBK2(17),DBK2(18),
     3    DBK2(19),DBK2(20)/ 5.49756809432471D-01, 9.13556983276901D-03,
     4-2.53635048605507D-03, 6.60423795342054D-04,-1.55217243135416D-04,
     5 3.00090325448633D-05,-3.76454339467348D-06,-1.33291331611616D-07,
     6 2.42587371049013D-07,-8.07861075240228D-08, 1.71092818861193D-08,
     7-2.41087357570599D-09, 1.53910848162371D-10, 2.56465373190630D-11,
     8-9.88581911653212D-12, 1.60877986412631D-12,-1.20952524741739D-13,
     9-1.06978278410820D-14, 5.02478557067561D-15,-8.68986130935886D-16/
      DATA DBK3(1), DBK3(2), DBK3(3), DBK3(4), DBK3(5), DBK3(6),
     1     DBK3(7), DBK3(8), DBK3(9), DBK3(10),DBK3(11),DBK3(12),
     2     DBK3(13),DBK3(14),DBK3(15),DBK3(16),DBK3(17),DBK3(18),
     3    DBK3(19),DBK3(20)/ 5.60598509354302D-01,-3.64870013248135D-03,
     4-5.98147152307417D-05,-2.33611595253625D-06,-1.64571516521436D-07,
     5-2.06333012920569D-08,-4.27745431573110D-09,-1.08494137799276D-09,
     6-2.37207188872763D-10,-2.22132920864966D-11, 1.07238008032138D-11,
     7 5.71954845245808D-12, 7.51102737777835D-13,-3.81912369483793D-13,
     8-1.75870057119257D-13, 6.69641694419084D-15, 2.26866724792055D-14,
     9 2.69898141356743D-15,-2.67133612397359D-15,-6.54121403165269D-16/
      DATA DBK4(1), DBK4(2), DBK4(3), DBK4(4), DBK4(5), DBK4(6),
     1     DBK4(7), DBK4(8), DBK4(9), DBK4(10),DBK4(11),DBK4(12),
     2    DBK4(13),DBK4(14)/ 4.93072999188036D-01, 4.38335419803815D-03,
     3-8.37413882246205D-05, 3.20268810484632D-06,-1.75661979548270D-07,
     4 1.22269906524508D-08,-1.01381314366052D-09, 9.63639784237475D-11,
     5-1.02344993379648D-11, 1.19264576554355D-12,-1.50443899103287D-13,
     6 2.03299052379349D-14,-2.91890652008292D-15, 4.42322081975475D-16/
      DATA DBJP(1), DBJP(2), DBJP(3), DBJP(4), DBJP(5), DBJP(6),
     1     DBJP(7), DBJP(8), DBJP(9), DBJP(10),DBJP(11),DBJP(12),
     2     DBJP(13),DBJP(14),DBJP(15),DBJP(16),DBJP(17),DBJP(18),
     3     DBJP(19)        / 1.13140872390745D-01,-2.08301511416328D-01,
     4 1.69396341953138D-02, 2.90895212478621D-02,-3.41467131311549D-03,
     5-1.46455339197417D-03, 1.63313272898517D-04, 3.91145328922162D-05,
     6-3.96757190808119D-06,-6.51846913772395D-07, 5.98707495269280D-08,
     7 7.44108654536549D-09,-6.21241056522632D-10,-6.18768017313526D-11,
     8 4.72323484752324D-12, 3.91652459802532D-13,-2.74985937845226D-14,
     9-1.95036497762750D-15, 1.26669643809444D-16/
      DATA DBJN(1), DBJN(2), DBJN(3), DBJN(4), DBJN(5), DBJN(6),
     1     DBJN(7), DBJN(8), DBJN(9), DBJN(10),DBJN(11),DBJN(12),
     2     DBJN(13),DBJN(14),DBJN(15),DBJN(16),DBJN(17),DBJN(18),
     3     DBJN(19)        /-1.88091260068850D-02,-1.47798180826140D-01,
     4 5.46075900433171D-01, 1.52146932663116D-01,-9.58260412266886D-02,
     5-1.63102731696130D-02, 5.75364806680105D-03, 7.12145408252655D-04,
     6-1.75452116846724D-04,-1.71063171685128D-05, 3.24435580631680D-06,
     7 2.61190663932884D-07,-4.03026865912779D-08,-2.76435165853895D-09,
     8 3.59687929062312D-10, 2.14953308456051D-11,-2.41849311903901D-12,
     9-1.28068004920751D-13, 1.26939834401773D-14/
      DATA DAA(1),  DAA(2),  DAA(3),  DAA(4),  DAA(5),  DAA(6),
     1     DAA(7),  DAA(8),  DAA(9),  DAA(10), DAA(11), DAA(12),
     2     DAA(13), DAA(14)/ 2.77571356944231D-01,-4.44212833419920D-03,
     3 8.42328522190089D-05, 2.58040318418710D-06,-3.42389720217621D-07,
     4 6.24286894709776D-09, 2.36377836844577D-09,-3.16991042656673D-10,
     5 4.40995691658191D-12, 5.18674221093575D-12,-9.64874015137022D-13,
     6 4.90190576608710D-14, 1.77253430678112D-14,-5.55950610442662D-15/
      DATA DBB(1),  DBB(2),  DBB(3),  DBB(4),  DBB(5),  DBB(6),
     1     DBB(7),  DBB(8),  DBB(9),  DBB(10), DBB(11), DBB(12),
     2     DBB(13), DBB(14)/ 4.91627321104601D-01, 3.11164930427489D-03,
     3 8.23140762854081D-05,-4.61769776172142D-06,-6.13158880534626D-08,
     4 2.87295804656520D-08,-1.81959715372117D-09,-1.44752826642035D-10,
     5 4.53724043420422D-11,-3.99655065847223D-12,-3.24089119830323D-13,
     6 1.62098952568741D-13,-2.40765247974057D-14, 1.69384811284491D-16/
C***FIRST EXECUTABLE STATEMENT  DYAIRY
      AX = ABS(X)
      RX = SQRT(AX)
      C = CON1*AX*RX
      IF (X.LT.0.0D0) GO TO 120
      IF (C.GT.8.0D0) GO TO 60
      IF (X.GT.2.5D0) GO TO 30
      T = (X+X-2.5D0)*0.4D0
      TT = T + T
      J = N1
      F1 = BK1(J)
      F2 = 0.0D0
      DO 10 I=1,M1
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + BK1(J)
        F2 = TEMP1
   10 CONTINUE
      BI = T*F1 - F2 + BK1(1)
      J = N1D
      F1 = DBK1(J)
      F2 = 0.0D0
      DO 20 I=1,M1D
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + DBK1(J)
        F2 = TEMP1
   20 CONTINUE
      DBI = T*F1 - F2 + DBK1(1)
      RETURN
   30 CONTINUE
      RTRX = SQRT(RX)
      T = (X+X-CON2)*CON3
      TT = T + T
      J = N1
      F1 = BK2(J)
      F2 = 0.0D0
      DO 40 I=1,M1
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + BK2(J)
        F2 = TEMP1
   40 CONTINUE
      BI = (T*F1-F2+BK2(1))/RTRX
      EX = EXP(C)
      BI = BI*EX
      J = N2D
      F1 = DBK2(J)
      F2 = 0.0D0
      DO 50 I=1,M2D
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + DBK2(J)
        F2 = TEMP1
   50 CONTINUE
      DBI = (T*F1-F2+DBK2(1))*RTRX
      DBI = DBI*EX
      RETURN
C
   60 CONTINUE
      RTRX = SQRT(RX)
      T = 16.0D0/C - 1.0D0
      TT = T + T
      J = N1
      F1 = BK3(J)
      F2 = 0.0D0
      DO 70 I=1,M1
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + BK3(J)
        F2 = TEMP1
   70 CONTINUE
      S1 = T*F1 - F2 + BK3(1)
      J = N2D
      F1 = DBK3(J)
      F2 = 0.0D0
      DO 80 I=1,M2D
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + DBK3(J)
        F2 = TEMP1
   80 CONTINUE
      D1 = T*F1 - F2 + DBK3(1)
      TC = C + C
      EX = EXP(C)
      IF (TC.GT.35.0D0) GO TO 110
      T = 10.0D0/C - 1.0D0
      TT = T + T
      J = N3
      F1 = BK4(J)
      F2 = 0.0D0
      DO 90 I=1,M3
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + BK4(J)
        F2 = TEMP1
   90 CONTINUE
      S2 = T*F1 - F2 + BK4(1)
      BI = (S1+EXP(-TC)*S2)/RTRX
      BI = BI*EX
      J = N4D
      F1 = DBK4(J)
      F2 = 0.0D0
      DO 100 I=1,M4D
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + DBK4(J)
        F2 = TEMP1
  100 CONTINUE
      D2 = T*F1 - F2 + DBK4(1)
      DBI = RTRX*(D1+EXP(-TC)*D2)
      DBI = DBI*EX
      RETURN
  110 BI = EX*S1/RTRX
      DBI = EX*RTRX*D1
      RETURN
C
  120 CONTINUE
      IF (C.GT.5.0D0) GO TO 150
      T = 0.4D0*C - 1.0D0
      TT = T + T
      J = N2
      F1 = BJP(J)
      E1 = BJN(J)
      F2 = 0.0D0
      E2 = 0.0D0
      DO 130 I=1,M2
        J = J - 1
        TEMP1 = F1
        TEMP2 = E1
        F1 = TT*F1 - F2 + BJP(J)
        E1 = TT*E1 - E2 + BJN(J)
        F2 = TEMP1
        E2 = TEMP2
  130 CONTINUE
      BI = (T*E1-E2+BJN(1)) - AX*(T*F1-F2+BJP(1))
      J = N3D
      F1 = DBJP(J)
      E1 = DBJN(J)
      F2 = 0.0D0
      E2 = 0.0D0
      DO 140 I=1,M3D
        J = J - 1
        TEMP1 = F1
        TEMP2 = E1
        F1 = TT*F1 - F2 + DBJP(J)
        E1 = TT*E1 - E2 + DBJN(J)
        F2 = TEMP1
        E2 = TEMP2
  140 CONTINUE
      DBI = X*X*(T*F1-F2+DBJP(1)) + (T*E1-E2+DBJN(1))
      RETURN
C
  150 CONTINUE
      RTRX = SQRT(RX)
      T = 10.0D0/C - 1.0D0
      TT = T + T
      J = N3
      F1 = AA(J)
      E1 = BB(J)
      F2 = 0.0D0
      E2 = 0.0D0
      DO 160 I=1,M3
        J = J - 1
        TEMP1 = F1
        TEMP2 = E1
        F1 = TT*F1 - F2 + AA(J)
        E1 = TT*E1 - E2 + BB(J)
        F2 = TEMP1
        E2 = TEMP2
  160 CONTINUE
      TEMP1 = T*F1 - F2 + AA(1)
      TEMP2 = T*E1 - E2 + BB(1)
      CV = C - FPI12
      BI = (TEMP1*COS(CV)+TEMP2*SIN(CV))/RTRX
      J = N4D
      F1 = DAA(J)
      E1 = DBB(J)
      F2 = 0.0D0
      E2 = 0.0D0
      DO 170 I=1,M4D
        J = J - 1
        TEMP1 = F1
        TEMP2 = E1
        F1 = TT*F1 - F2 + DAA(J)
        E1 = TT*E1 - E2 + DBB(J)
        F2 = TEMP1
        E2 = TEMP2
  170 CONTINUE
      TEMP1 = T*F1 - F2 + DAA(1)
      TEMP2 = T*E1 - E2 + DBB(1)
      CV = C - SPI12
      DBI = (TEMP1*COS(CV)-TEMP2*SIN(CV))*RTRX
      RETURN
      END