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

Go to the documentation of this file.
00001       SUBROUTINE ZUNIK(ZRR, ZRI, FNU, IKFLG, IPMTR, TOL, INIT, PHIR,
00002      * PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
00003 C***BEGIN PROLOGUE  ZUNIK
00004 C***REFER TO  ZBESI,ZBESK
00005 C
00006 C        ZUNIK 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  ZDIV,XZLOG,XZSQRT,D1MACH
00022 C***END PROLOGUE  ZUNIK
00023 C     COMPLEX CFN,CON,CONE,CRFN,CWRK,CZERO,PHI,S,SR,SUM,T,T2,ZETA1,
00024 C    *ZETA2,ZN,ZR
00025       DOUBLE PRECISION AC, C, CON, CONEI, CONER, CRFNI, CRFNR, CWRKI,
00026      * CWRKR, FNU, PHII, PHIR, RFN, SI, SR, SRI, SRR, STI, STR, SUMI,
00027      * SUMR, TEST, TI, TOL, TR, T2I, T2R, ZEROI, ZEROR, ZETA1I, ZETA1R,
00028      * ZETA2I, ZETA2R, ZNI, ZNR, ZRI, ZRR, D1MACH
00029       INTEGER I, IDUM, IKFLG, INIT, IPMTR, J, K, L
00030       DIMENSION C(120), CWRKR(16), CWRKI(16), CON(2)
00031       DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 /
00032       DATA CON(1), CON(2)  /
00033      1 3.98942280401432678D-01,  1.25331413731550025D+00 /
00034       DATA C(1), C(2), C(3), C(4), C(5), C(6), C(7), C(8), C(9), C(10),
00035      1     C(11), C(12), C(13), C(14), C(15), C(16), C(17), C(18),
00036      2     C(19), C(20), C(21), C(22), C(23), C(24)/
00037      3     1.00000000000000000D+00,    -2.08333333333333333D-01,
00038      4     1.25000000000000000D-01,     3.34201388888888889D-01,
00039      5    -4.01041666666666667D-01,     7.03125000000000000D-02,
00040      6    -1.02581259645061728D+00,     1.84646267361111111D+00,
00041      7    -8.91210937500000000D-01,     7.32421875000000000D-02,
00042      8     4.66958442342624743D+00,    -1.12070026162229938D+01,
00043      9     8.78912353515625000D+00,    -2.36408691406250000D+00,
00044      A     1.12152099609375000D-01,    -2.82120725582002449D+01,
00045      B     8.46362176746007346D+01,    -9.18182415432400174D+01,
00046      C     4.25349987453884549D+01,    -7.36879435947963170D+00,
00047      D     2.27108001708984375D-01,     2.12570130039217123D+02,
00048      E    -7.65252468141181642D+02,     1.05999045252799988D+03/
00049       DATA C(25), C(26), C(27), C(28), C(29), C(30), C(31), C(32),
00050      1     C(33), C(34), C(35), C(36), C(37), C(38), C(39), C(40),
00051      2     C(41), C(42), C(43), C(44), C(45), C(46), C(47), C(48)/
00052      3    -6.99579627376132541D+02,     2.18190511744211590D+02,
00053      4    -2.64914304869515555D+01,     5.72501420974731445D-01,
00054      5    -1.91945766231840700D+03,     8.06172218173730938D+03,
00055      6    -1.35865500064341374D+04,     1.16553933368645332D+04,
00056      7    -5.30564697861340311D+03,     1.20090291321635246D+03,
00057      8    -1.08090919788394656D+02,     1.72772750258445740D+00,
00058      9     2.02042913309661486D+04,    -9.69805983886375135D+04,
00059      A     1.92547001232531532D+05,    -2.03400177280415534D+05,
00060      B     1.22200464983017460D+05,    -4.11926549688975513D+04,
00061      C     7.10951430248936372D+03,    -4.93915304773088012D+02,
00062      D     6.07404200127348304D+00,    -2.42919187900551333D+05,
00063      E     1.31176361466297720D+06,    -2.99801591853810675D+06/
00064       DATA C(49), C(50), C(51), C(52), C(53), C(54), C(55), C(56),
00065      1     C(57), C(58), C(59), C(60), C(61), C(62), C(63), C(64),
00066      2     C(65), C(66), C(67), C(68), C(69), C(70), C(71), C(72)/
00067      3     3.76327129765640400D+06,    -2.81356322658653411D+06,
00068      4     1.26836527332162478D+06,    -3.31645172484563578D+05,
00069      5     4.52187689813627263D+04,    -2.49983048181120962D+03,
00070      6     2.43805296995560639D+01,     3.28446985307203782D+06,
00071      7    -1.97068191184322269D+07,     5.09526024926646422D+07,
00072      8    -7.41051482115326577D+07,     6.63445122747290267D+07,
00073      9    -3.75671766607633513D+07,     1.32887671664218183D+07,
00074      A    -2.78561812808645469D+06,     3.08186404612662398D+05,
00075      B    -1.38860897537170405D+04,     1.10017140269246738D+02,
00076      C    -4.93292536645099620D+07,     3.25573074185765749D+08,
00077      D    -9.39462359681578403D+08,     1.55359689957058006D+09,
00078      E    -1.62108055210833708D+09,     1.10684281682301447D+09/
00079       DATA C(73), C(74), C(75), C(76), C(77), C(78), C(79), C(80),
00080      1     C(81), C(82), C(83), C(84), C(85), C(86), C(87), C(88),
00081      2     C(89), C(90), C(91), C(92), C(93), C(94), C(95), C(96)/
00082      3    -4.95889784275030309D+08,     1.42062907797533095D+08,
00083      4    -2.44740627257387285D+07,     2.24376817792244943D+06,
00084      5    -8.40054336030240853D+04,     5.51335896122020586D+02,
00085      6     8.14789096118312115D+08,    -5.86648149205184723D+09,
00086      7     1.86882075092958249D+10,    -3.46320433881587779D+10,
00087      8     4.12801855797539740D+10,    -3.30265997498007231D+10,
00088      9     1.79542137311556001D+10,    -6.56329379261928433D+09,
00089      A     1.55927986487925751D+09,    -2.25105661889415278D+08,
00090      B     1.73951075539781645D+07,    -5.49842327572288687D+05,
00091      C     3.03809051092238427D+03,    -1.46792612476956167D+10,
00092      D     1.14498237732025810D+11,    -3.99096175224466498D+11,
00093      E     8.19218669548577329D+11,    -1.09837515608122331D+12/
00094       DATA C(97), C(98), C(99), C(100), C(101), C(102), C(103), C(104),
00095      1     C(105), C(106), C(107), C(108), C(109), C(110), C(111),
00096      2     C(112), C(113), C(114), C(115), C(116), C(117), C(118)/
00097      3     1.00815810686538209D+12,    -6.45364869245376503D+11,
00098      4     2.87900649906150589D+11,    -8.78670721780232657D+10,
00099      5     1.76347306068349694D+10,    -2.16716498322379509D+09,
00100      6     1.43157876718888981D+08,    -3.87183344257261262D+06,
00101      7     1.82577554742931747D+04,     2.86464035717679043D+11,
00102      8    -2.40629790002850396D+12,     9.10934118523989896D+12,
00103      9    -2.05168994109344374D+13,     3.05651255199353206D+13,
00104      A    -3.16670885847851584D+13,     2.33483640445818409D+13,
00105      B    -1.23204913055982872D+13,     4.61272578084913197D+12,
00106      C    -1.19655288019618160D+12,     2.05914503232410016D+11,
00107      D    -2.18229277575292237D+10,     1.24700929351271032D+09/
00108       DATA C(119), C(120)/
00109      1    -2.91883881222208134D+07,     1.18838426256783253D+05/
00110 C
00111       IF (INIT.NE.0) GO TO 40
00112 C-----------------------------------------------------------------------
00113 C     INITIALIZE ALL VARIABLES
00114 C-----------------------------------------------------------------------
00115       RFN = 1.0D0/FNU
00116 C-----------------------------------------------------------------------
00117 C     OVERFLOW TEST (ZR/FNU TOO SMALL)
00118 C-----------------------------------------------------------------------
00119       TEST = D1MACH(1)*1.0D+3
00120       AC = FNU*TEST
00121       IF (DABS(ZRR).GT.AC .OR. DABS(ZRI).GT.AC) GO TO 15
00122       ZETA1R = 2.0D0*DABS(DLOG(TEST))+FNU
00123       ZETA1I = 0.0D0
00124       ZETA2R = FNU
00125       ZETA2I = 0.0D0
00126       PHIR = 1.0D0
00127       PHII = 0.0D0
00128       RETURN
00129    15 CONTINUE
00130       TR = ZRR*RFN
00131       TI = ZRI*RFN
00132       SR = CONER + (TR*TR-TI*TI)
00133       SI = CONEI + (TR*TI+TI*TR)
00134       CALL XZSQRT(SR, SI, SRR, SRI)
00135       STR = CONER + SRR
00136       STI = CONEI + SRI
00137       CALL ZDIV(STR, STI, TR, TI, ZNR, ZNI)
00138       CALL XZLOG(ZNR, ZNI, STR, STI, IDUM)
00139       ZETA1R = FNU*STR
00140       ZETA1I = FNU*STI
00141       ZETA2R = FNU*SRR
00142       ZETA2I = FNU*SRI
00143       CALL ZDIV(CONER, CONEI, SRR, SRI, TR, TI)
00144       SRR = TR*RFN
00145       SRI = TI*RFN
00146       CALL XZSQRT(SRR, SRI, CWRKR(16), CWRKI(16))
00147       PHIR = CWRKR(16)*CON(IKFLG)
00148       PHII = CWRKI(16)*CON(IKFLG)
00149       IF (IPMTR.NE.0) RETURN
00150       CALL ZDIV(CONER, CONEI, SR, SI, T2R, T2I)
00151       CWRKR(1) = CONER
00152       CWRKI(1) = CONEI
00153       CRFNR = CONER
00154       CRFNI = CONEI
00155       AC = 1.0D0
00156       L = 1
00157       DO 20 K=2,15
00158         SR = ZEROR
00159         SI = ZEROI
00160         DO 10 J=1,K
00161           L = L + 1
00162           STR = SR*T2R - SI*T2I + C(L)
00163           SI = SR*T2I + SI*T2R
00164           SR = STR
00165    10   CONTINUE
00166         STR = CRFNR*SRR - CRFNI*SRI
00167         CRFNI = CRFNR*SRI + CRFNI*SRR
00168         CRFNR = STR
00169         CWRKR(K) = CRFNR*SR - CRFNI*SI
00170         CWRKI(K) = CRFNR*SI + CRFNI*SR
00171         AC = AC*RFN
00172         TEST = DABS(CWRKR(K)) + DABS(CWRKI(K))
00173         IF (AC.LT.TOL .AND. TEST.LT.TOL) GO TO 30
00174    20 CONTINUE
00175       K = 15
00176    30 CONTINUE
00177       INIT = K
00178    40 CONTINUE
00179       IF (IKFLG.EQ.2) GO TO 60
00180 C-----------------------------------------------------------------------
00181 C     COMPUTE SUM FOR THE I FUNCTION
00182 C-----------------------------------------------------------------------
00183       SR = ZEROR
00184       SI = ZEROI
00185       DO 50 I=1,INIT
00186         SR = SR + CWRKR(I)
00187         SI = SI + CWRKI(I)
00188    50 CONTINUE
00189       SUMR = SR
00190       SUMI = SI
00191       PHIR = CWRKR(16)*CON(1)
00192       PHII = CWRKI(16)*CON(1)
00193       RETURN
00194    60 CONTINUE
00195 C-----------------------------------------------------------------------
00196 C     COMPUTE SUM FOR THE K FUNCTION
00197 C-----------------------------------------------------------------------
00198       SR = ZEROR
00199       SI = ZEROI
00200       TR = CONER
00201       DO 70 I=1,INIT
00202         SR = SR + TR*CWRKR(I)
00203         SI = SI + TR*CWRKI(I)
00204         TR = -TR
00205    70 CONTINUE
00206       SUMR = SR
00207       SUMI = SI
00208       PHIR = CWRKR(16)*CON(2)
00209       PHII = CWRKI(16)*CON(2)
00210       RETURN
00211       END