GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
cunik.f
Go to the documentation of this file.
1  SUBROUTINE cunik(ZR, FNU, IKFLG, IPMTR, TOL, INIT, PHI, ZETA1,
2  * ZETA2, SUM, CWRK)
3 C***BEGIN PROLOGUE CUNIK
4 C***REFER TO CBESI,CBESK
5 C
6 C CUNIK COMPUTES PARAMETERS FOR THE UNIFORM ASYMPTOTIC
7 C EXPANSIONS OF THE I AND K FUNCTIONS ON IKFLG= 1 OR 2
8 C RESPECTIVELY BY
9 C
10 C W(FNU,ZR) = PHI*EXP(ZETA)*SUM
11 C
12 C WHERE ZETA=-ZETA1 + ZETA2 OR
13 C ZETA1 - ZETA2
14 C
15 C THE FIRST CALL MUST HAVE INIT=0. SUBSEQUENT CALLS WITH THE
16 C SAME ZR AND FNU WILL RETURN THE I OR K FUNCTION ON IKFLG=
17 C 1 OR 2 WITH NO CHANGE IN INIT. CWRK IS A COMPLEX WORK
18 C ARRAY. IPMTR=0 COMPUTES ALL PARAMETERS. IPMTR=1 COMPUTES PHI,
19 C ZETA1,ZETA2.
20 C
21 C***ROUTINES CALLED (NONE)
22 C***END PROLOGUE CUNIK
23  COMPLEX CFN, CON, CONE, CRFN, CWRK, CZERO, PHI, S, SR, SUM, T,
24  * t2, zeta1, zeta2, zn, zr
25  REAL AC, C, FNU, RFN, TEST, TOL, TSTR, TSTI
26  INTEGER I, IKFLG, INIT, IPMTR, J, K, L
27  dimension c(120), cwrk(16), con(2)
28  DATA czero, cone / (0.0e0,0.0e0), (1.0e0,0.0e0) /
29  DATA con(1), con(2) /
30  1(3.98942280401432678e-01,0.0e0),(1.25331413731550025e+00,0.0e0)/
31  DATA c(1), c(2), c(3), c(4), c(5), c(6), c(7), c(8), c(9), c(10),
32  1 c(11), c(12), c(13), c(14), c(15), c(16), c(17), c(18),
33  2 c(19), c(20), c(21), c(22), c(23), c(24)/
34  3 1.00000000000000000e+00, -2.08333333333333333e-01,
35  4 1.25000000000000000e-01, 3.34201388888888889e-01,
36  5 -4.01041666666666667e-01, 7.03125000000000000e-02,
37  6 -1.02581259645061728e+00, 1.84646267361111111e+00,
38  7 -8.91210937500000000e-01, 7.32421875000000000e-02,
39  8 4.66958442342624743e+00, -1.12070026162229938e+01,
40  9 8.78912353515625000e+00, -2.36408691406250000e+00,
41  a 1.12152099609375000e-01, -2.82120725582002449e+01,
42  b 8.46362176746007346e+01, -9.18182415432400174e+01,
43  c 4.25349987453884549e+01, -7.36879435947963170e+00,
44  d 2.27108001708984375e-01, 2.12570130039217123e+02,
45  e -7.65252468141181642e+02, 1.05999045252799988e+03/
46  DATA c(25), c(26), c(27), c(28), c(29), c(30), c(31), c(32),
47  1 c(33), c(34), c(35), c(36), c(37), c(38), c(39), c(40),
48  2 c(41), c(42), c(43), c(44), c(45), c(46), c(47), c(48)/
49  3 -6.99579627376132541e+02, 2.18190511744211590e+02,
50  4 -2.64914304869515555e+01, 5.72501420974731445e-01,
51  5 -1.91945766231840700e+03, 8.06172218173730938e+03,
52  6 -1.35865500064341374e+04, 1.16553933368645332e+04,
53  7 -5.30564697861340311e+03, 1.20090291321635246e+03,
54  8 -1.08090919788394656e+02, 1.72772750258445740e+00,
55  9 2.02042913309661486e+04, -9.69805983886375135e+04,
56  a 1.92547001232531532e+05, -2.03400177280415534e+05,
57  b 1.22200464983017460e+05, -4.11926549688975513e+04,
58  c 7.10951430248936372e+03, -4.93915304773088012e+02,
59  d 6.07404200127348304e+00, -2.42919187900551333e+05,
60  e 1.31176361466297720e+06, -2.99801591853810675e+06/
61  DATA c(49), c(50), c(51), c(52), c(53), c(54), c(55), c(56),
62  1 c(57), c(58), c(59), c(60), c(61), c(62), c(63), c(64),
63  2 c(65), c(66), c(67), c(68), c(69), c(70), c(71), c(72)/
64  3 3.76327129765640400e+06, -2.81356322658653411e+06,
65  4 1.26836527332162478e+06, -3.31645172484563578e+05,
66  5 4.52187689813627263e+04, -2.49983048181120962e+03,
67  6 2.43805296995560639e+01, 3.28446985307203782e+06,
68  7 -1.97068191184322269e+07, 5.09526024926646422e+07,
69  8 -7.41051482115326577e+07, 6.63445122747290267e+07,
70  9 -3.75671766607633513e+07, 1.32887671664218183e+07,
71  a -2.78561812808645469e+06, 3.08186404612662398e+05,
72  b -1.38860897537170405e+04, 1.10017140269246738e+02,
73  c -4.93292536645099620e+07, 3.25573074185765749e+08,
74  d -9.39462359681578403e+08, 1.55359689957058006e+09,
75  e -1.62108055210833708e+09, 1.10684281682301447e+09/
76  DATA c(73), c(74), c(75), c(76), c(77), c(78), c(79), c(80),
77  1 c(81), c(82), c(83), c(84), c(85), c(86), c(87), c(88),
78  2 c(89), c(90), c(91), c(92), c(93), c(94), c(95), c(96)/
79  3 -4.95889784275030309e+08, 1.42062907797533095e+08,
80  4 -2.44740627257387285e+07, 2.24376817792244943e+06,
81  5 -8.40054336030240853e+04, 5.51335896122020586e+02,
82  6 8.14789096118312115e+08, -5.86648149205184723e+09,
83  7 1.86882075092958249e+10, -3.46320433881587779e+10,
84  8 4.12801855797539740e+10, -3.30265997498007231e+10,
85  9 1.79542137311556001e+10, -6.56329379261928433e+09,
86  a 1.55927986487925751e+09, -2.25105661889415278e+08,
87  b 1.73951075539781645e+07, -5.49842327572288687e+05,
88  c 3.03809051092238427e+03, -1.46792612476956167e+10,
89  d 1.14498237732025810e+11, -3.99096175224466498e+11,
90  e 8.19218669548577329e+11, -1.09837515608122331e+12/
91  DATA c(97), c(98), c(99), c(100), c(101), c(102), c(103), c(104),
92  1 c(105), c(106), c(107), c(108), c(109), c(110), c(111),
93  2 c(112), c(113), c(114), c(115), c(116), c(117), c(118)/
94  3 1.00815810686538209e+12, -6.45364869245376503e+11,
95  4 2.87900649906150589e+11, -8.78670721780232657e+10,
96  5 1.76347306068349694e+10, -2.16716498322379509e+09,
97  6 1.43157876718888981e+08, -3.87183344257261262e+06,
98  7 1.82577554742931747e+04, 2.86464035717679043e+11,
99  8 -2.40629790002850396e+12, 9.10934118523989896e+12,
100  9 -2.05168994109344374e+13, 3.05651255199353206e+13,
101  a -3.16670885847851584e+13, 2.33483640445818409e+13,
102  b -1.23204913055982872e+13, 4.61272578084913197e+12,
103  c -1.19655288019618160e+12, 2.05914503232410016e+11,
104  d -2.18229277575292237e+10, 1.24700929351271032e+09/
105  DATA c(119), c(120)/
106  1 -2.91883881222208134e+07, 1.18838426256783253e+05/
107 C
108  IF (init.NE.0) GO TO 40
109 C-----------------------------------------------------------------------
110 C INITIALIZE ALL VARIABLES
111 C-----------------------------------------------------------------------
112  rfn = 1.0e0/fnu
113  crfn = cmplx(rfn,0.0e0)
114 C T = ZR*CRFN
115 C-----------------------------------------------------------------------
116 C OVERFLOW TEST (ZR/FNU TOO SMALL)
117 C-----------------------------------------------------------------------
118  tstr = real(zr)
119  tsti = aimag(zr)
120  test = r1mach(1)*1.0e+3
121  ac = fnu*test
122  IF (abs(tstr).GT.ac .OR. abs(tsti).GT.ac) GO TO 15
123  ac = 2.0e0*abs(alog(test))+fnu
124  zeta1 = cmplx(ac,0.0e0)
125  zeta2 = cmplx(fnu,0.0e0)
126  phi=cone
127  RETURN
128  15 CONTINUE
129  t=zr*crfn
130  s = cone + t*t
131  sr = csqrt(s)
132  cfn = cmplx(fnu,0.0e0)
133  zn = (cone+sr)/t
134  zeta1 = cfn*clog(zn)
135  zeta2 = cfn*sr
136  t = cone/sr
137  sr = t*crfn
138  cwrk(16) = csqrt(sr)
139  phi = cwrk(16)*con(ikflg)
140  IF (ipmtr.NE.0) RETURN
141  t2 = cone/s
142  cwrk(1) = cone
143  crfn = cone
144  ac = 1.0e0
145  l = 1
146  DO 20 k=2,15
147  s = czero
148  DO 10 j=1,k
149  l = l + 1
150  s = s*t2 + cmplx(c(l),0.0e0)
151  10 CONTINUE
152  crfn = crfn*sr
153  cwrk(k) = crfn*s
154  ac = ac*rfn
155  tstr = real(cwrk(k))
156  tsti = aimag(cwrk(k))
157  test = abs(tstr) + abs(tsti)
158  IF (ac.LT.tol .AND. test.LT.tol) GO TO 30
159  20 CONTINUE
160  k = 15
161  30 CONTINUE
162  init = k
163  40 CONTINUE
164  IF (ikflg.EQ.2) GO TO 60
165 C-----------------------------------------------------------------------
166 C COMPUTE SUM FOR THE I FUNCTION
167 C-----------------------------------------------------------------------
168  s = czero
169  DO 50 i=1,init
170  s = s + cwrk(i)
171  50 CONTINUE
172  sum = s
173  phi = cwrk(16)*con(1)
174  RETURN
175  60 CONTINUE
176 C-----------------------------------------------------------------------
177 C COMPUTE SUM FOR THE K FUNCTION
178 C-----------------------------------------------------------------------
179  s = czero
180  t = cone
181  DO 70 i=1,init
182  s = s + t*cwrk(i)
183  t = -t
184  70 CONTINUE
185  sum = s
186  phi = cwrk(16)*con(2)
187  RETURN
188  END
cmplx
double complex cmplx
Definition: Faddeeva.cc:217
cunik
subroutine cunik(ZR, FNU, IKFLG, IPMTR, TOL, INIT, PHI, ZETA1, ZETA2, SUM, CWRK)
Definition: cunik.f:3
real
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:137
abs
static T abs(T x)
Definition: pr-output.cc:1678
r1mach
real function r1mach(i)
Definition: r1mach.f:23