GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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
Generated on Mon Dec 30 2013 03:04:45 for GNU Octave by   doxygen 1.8.1.2