GNU Octave: libcruft/amos/cunik.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CUNIK(ZR, FNU, IKFLG, IPMTR, TOL, INIT, PHI, ZETA1,
00002      * ZETA2, SUM, CWRK)
00003 C***BEGIN PROLOGUE  CUNIK
00004 C***REFER TO  CBESI,CBESK
00005 C
00006 C        CUNIK COMPUTES PARAMETERS FOR THE UNIFORM ASYMPTOTIC
00007 C        EXPANSIONS OF THE I AND K FUNCTIONS ON IKFLG= 1 OR 2
00008 C        RESPECTIVELY BY
00009 C
00010 C        W(FNU,ZR) = PHI*EXP(ZETA)*SUM
00011 C
00012 C        WHERE       ZETA=-ZETA1 + ZETA2       OR
00013 C                          ZETA1 - ZETA2
00014 C
00015 C        THE FIRST CALL MUST HAVE INIT=0. SUBSEQUENT CALLS WITH THE
00016 C        SAME ZR AND FNU WILL RETURN THE I OR K FUNCTION ON IKFLG=
00017 C        1 OR 2 WITH NO CHANGE IN INIT. CWRK IS A COMPLEX WORK
00018 C        ARRAY. IPMTR=0 COMPUTES ALL PARAMETERS. IPMTR=1 COMPUTES PHI,
00019 C        ZETA1,ZETA2.
00020 C
00021 C***ROUTINES CALLED  (NONE)
00022 C***END PROLOGUE  CUNIK
00023       COMPLEX CFN, CON, CONE, CRFN, CWRK, CZERO, PHI, S, SR, SUM, T,
00024      * T2, ZETA1, ZETA2, ZN, ZR
00025       REAL AC, C, FNU, RFN, TEST, TOL, TSTR, TSTI
00026       INTEGER I, IKFLG, INIT, IPMTR, J, K, L
00027       DIMENSION C(120), CWRK(16), CON(2)
00028       DATA CZERO, CONE / (0.0E0,0.0E0), (1.0E0,0.0E0) /
00029       DATA CON(1), CON(2)  /
00030      1(3.98942280401432678E-01,0.0E0),(1.25331413731550025E+00,0.0E0)/
00031       DATA C(1), C(2), C(3), C(4), C(5), C(6), C(7), C(8), C(9), C(10),
00032      1     C(11), C(12), C(13), C(14), C(15), C(16), C(17), C(18),
00033      2     C(19), C(20), C(21), C(22), C(23), C(24)/
00034      3     1.00000000000000000E+00,    -2.08333333333333333E-01,
00035      4     1.25000000000000000E-01,     3.34201388888888889E-01,
00036      5    -4.01041666666666667E-01,     7.03125000000000000E-02,
00037      6    -1.02581259645061728E+00,     1.84646267361111111E+00,
00038      7    -8.91210937500000000E-01,     7.32421875000000000E-02,
00039      8     4.66958442342624743E+00,    -1.12070026162229938E+01,
00040      9     8.78912353515625000E+00,    -2.36408691406250000E+00,
00041      A     1.12152099609375000E-01,    -2.82120725582002449E+01,
00042      B     8.46362176746007346E+01,    -9.18182415432400174E+01,
00043      C     4.25349987453884549E+01,    -7.36879435947963170E+00,
00044      D     2.27108001708984375E-01,     2.12570130039217123E+02,
00045      E    -7.65252468141181642E+02,     1.05999045252799988E+03/
00046       DATA C(25), C(26), C(27), C(28), C(29), C(30), C(31), C(32),
00047      1     C(33), C(34), C(35), C(36), C(37), C(38), C(39), C(40),
00048      2     C(41), C(42), C(43), C(44), C(45), C(46), C(47), C(48)/
00049      3    -6.99579627376132541E+02,     2.18190511744211590E+02,
00050      4    -2.64914304869515555E+01,     5.72501420974731445E-01,
00051      5    -1.91945766231840700E+03,     8.06172218173730938E+03,
00052      6    -1.35865500064341374E+04,     1.16553933368645332E+04,
00053      7    -5.30564697861340311E+03,     1.20090291321635246E+03,
00054      8    -1.08090919788394656E+02,     1.72772750258445740E+00,
00055      9     2.02042913309661486E+04,    -9.69805983886375135E+04,
00056      A     1.92547001232531532E+05,    -2.03400177280415534E+05,
00057      B     1.22200464983017460E+05,    -4.11926549688975513E+04,
00058      C     7.10951430248936372E+03,    -4.93915304773088012E+02,
00059      D     6.07404200127348304E+00,    -2.42919187900551333E+05,
00060      E     1.31176361466297720E+06,    -2.99801591853810675E+06/
00061       DATA C(49), C(50), C(51), C(52), C(53), C(54), C(55), C(56),
00062      1     C(57), C(58), C(59), C(60), C(61), C(62), C(63), C(64),
00063      2     C(65), C(66), C(67), C(68), C(69), C(70), C(71), C(72)/
00064      3     3.76327129765640400E+06,    -2.81356322658653411E+06,
00065      4     1.26836527332162478E+06,    -3.31645172484563578E+05,
00066      5     4.52187689813627263E+04,    -2.49983048181120962E+03,
00067      6     2.43805296995560639E+01,     3.28446985307203782E+06,
00068      7    -1.97068191184322269E+07,     5.09526024926646422E+07,
00069      8    -7.41051482115326577E+07,     6.63445122747290267E+07,
00070      9    -3.75671766607633513E+07,     1.32887671664218183E+07,
00071      A    -2.78561812808645469E+06,     3.08186404612662398E+05,
00072      B    -1.38860897537170405E+04,     1.10017140269246738E+02,
00073      C    -4.93292536645099620E+07,     3.25573074185765749E+08,
00074      D    -9.39462359681578403E+08,     1.55359689957058006E+09,
00075      E    -1.62108055210833708E+09,     1.10684281682301447E+09/
00076       DATA C(73), C(74), C(75), C(76), C(77), C(78), C(79), C(80),
00077      1     C(81), C(82), C(83), C(84), C(85), C(86), C(87), C(88),
00078      2     C(89), C(90), C(91), C(92), C(93), C(94), C(95), C(96)/
00079      3    -4.95889784275030309E+08,     1.42062907797533095E+08,
00080      4    -2.44740627257387285E+07,     2.24376817792244943E+06,
00081      5    -8.40054336030240853E+04,     5.51335896122020586E+02,
00082      6     8.14789096118312115E+08,    -5.86648149205184723E+09,
00083      7     1.86882075092958249E+10,    -3.46320433881587779E+10,
00084      8     4.12801855797539740E+10,    -3.30265997498007231E+10,
00085      9     1.79542137311556001E+10,    -6.56329379261928433E+09,
00086      A     1.55927986487925751E+09,    -2.25105661889415278E+08,
00087      B     1.73951075539781645E+07,    -5.49842327572288687E+05,
00088      C     3.03809051092238427E+03,    -1.46792612476956167E+10,
00089      D     1.14498237732025810E+11,    -3.99096175224466498E+11,
00090      E     8.19218669548577329E+11,    -1.09837515608122331E+12/
00091       DATA C(97), C(98), C(99), C(100), C(101), C(102), C(103), C(104),
00092      1     C(105), C(106), C(107), C(108), C(109), C(110), C(111),
00093      2     C(112), C(113), C(114), C(115), C(116), C(117), C(118)/
00094      3     1.00815810686538209E+12,    -6.45364869245376503E+11,
00095      4     2.87900649906150589E+11,    -8.78670721780232657E+10,
00096      5     1.76347306068349694E+10,    -2.16716498322379509E+09,
00097      6     1.43157876718888981E+08,    -3.87183344257261262E+06,
00098      7     1.82577554742931747E+04,     2.86464035717679043E+11,
00099      8    -2.40629790002850396E+12,     9.10934118523989896E+12,
00100      9    -2.05168994109344374E+13,     3.05651255199353206E+13,
00101      A    -3.16670885847851584E+13,     2.33483640445818409E+13,
00102      B    -1.23204913055982872E+13,     4.61272578084913197E+12,
00103      C    -1.19655288019618160E+12,     2.05914503232410016E+11,
00104      D    -2.18229277575292237E+10,     1.24700929351271032E+09/
00105       DATA C(119), C(120)/
00106      1    -2.91883881222208134E+07,     1.18838426256783253E+05/
00107 C
00108       IF (INIT.NE.0) GO TO 40
00109 C-----------------------------------------------------------------------
00110 C     INITIALIZE ALL VARIABLES
00111 C-----------------------------------------------------------------------
00112       RFN = 1.0E0/FNU
00113       CRFN = CMPLX(RFN,0.0E0)
00114 C     T = ZR*CRFN
00115 C-----------------------------------------------------------------------
00116 C     OVERFLOW TEST (ZR/FNU TOO SMALL)
00117 C-----------------------------------------------------------------------
00118       TSTR = REAL(ZR)
00119       TSTI = AIMAG(ZR)
00120       TEST = R1MACH(1)*1.0E+3
00121       AC = FNU*TEST
00122       IF (ABS(TSTR).GT.AC .OR. ABS(TSTI).GT.AC) GO TO 15
00123       AC = 2.0E0*ABS(ALOG(TEST))+FNU
00124       ZETA1 = CMPLX(AC,0.0E0)
00125       ZETA2 = CMPLX(FNU,0.0E0)
00126       PHI=CONE
00127       RETURN
00128    15 CONTINUE
00129       T=ZR*CRFN
00130       S = CONE + T*T
00131       SR = CSQRT(S)
00132       CFN = CMPLX(FNU,0.0E0)
00133       ZN = (CONE+SR)/T
00134       ZETA1 = CFN*CLOG(ZN)
00135       ZETA2 = CFN*SR
00136       T = CONE/SR
00137       SR = T*CRFN
00138       CWRK(16) = CSQRT(SR)
00139       PHI = CWRK(16)*CON(IKFLG)
00140       IF (IPMTR.NE.0) RETURN
00141       T2 = CONE/S
00142       CWRK(1) = CONE
00143       CRFN = CONE
00144       AC = 1.0E0
00145       L = 1
00146       DO 20 K=2,15
00147         S = CZERO
00148         DO 10 J=1,K
00149           L = L + 1
00150           S = S*T2 + CMPLX(C(L),0.0E0)
00151    10   CONTINUE
00152         CRFN = CRFN*SR
00153         CWRK(K) = CRFN*S
00154         AC = AC*RFN
00155         TSTR = REAL(CWRK(K))
00156         TSTI = AIMAG(CWRK(K))
00157         TEST = ABS(TSTR) + ABS(TSTI)
00158         IF (AC.LT.TOL .AND. TEST.LT.TOL) GO TO 30
00159    20 CONTINUE
00160       K = 15
00161    30 CONTINUE
00162       INIT = K
00163    40 CONTINUE
00164       IF (IKFLG.EQ.2) GO TO 60
00165 C-----------------------------------------------------------------------
00166 C     COMPUTE SUM FOR THE I FUNCTION
00167 C-----------------------------------------------------------------------
00168       S = CZERO
00169       DO 50 I=1,INIT
00170         S = S + CWRK(I)
00171    50 CONTINUE
00172       SUM = S
00173       PHI = CWRK(16)*CON(1)
00174       RETURN
00175    60 CONTINUE
00176 C-----------------------------------------------------------------------
00177 C     COMPUTE SUM FOR THE K FUNCTION
00178 C-----------------------------------------------------------------------
00179       S = CZERO
00180       T = CONE
00181       DO 70 I=1,INIT
00182         S = S + T*CWRK(I)
00183         T = -T
00184    70 CONTINUE
00185       SUM = S
00186       PHI = CWRK(16)*CON(2)
00187       RETURN
00188       END