GNU Octave  4.2.1
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
i e
Definition: data.cc:2724
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
double complex cmplx
Definition: Faddeeva.cc:217
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
Definition: Quad-opts.cc:233
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
subroutine cunik(ZR, FNU, IKFLG, IPMTR, TOL, INIT, PHI, ZETA1, ZETA2, SUM, CWRK)
Definition: cunik.f:1
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
b
Definition: cellfun.cc:398