CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
dqag.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine dqag (f, a, b, epsabs, epsrel, key, result, abserr, neval, ier, limit, lenw, last, iwork, work, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
 
subroutine dqage (f, a, b, epsabs, epsrel, key, limit, result, abserr, neval, ier, alist, blist, rlist, elist, iord, last, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
 
subroutine dqk15 (f, a, b, result, abserr, resabs, resasc, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
 
subroutine dqk21 (f, a, b, result, abserr, resabs, resasc, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
 
subroutine dqk31 (f, a, b, result, abserr, resabs, resasc, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
 
subroutine dqk41 (f, a, b, result, abserr, resabs, resasc, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
 
subroutine dqk51 (f, a, b, result, abserr, resabs, resasc, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
 
subroutine dqk61 (f, a, b, result, abserr, resabs, resasc, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
 
subroutine dqpsrt (limit, last, maxerr, ermax, elist, iord, nrmax, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
 

Function/Subroutine Documentation

◆ dqag()

subroutine dqag ( external  f,
real*8  a,
real*8  b,
real*8  epsabs,
real*8  epsrel,
integer  key,
real*8  result,
real*8  abserr,
integer  neval,
integer  ier,
integer  limit,
integer  lenw,
integer  last,
integer, dimension(limit)  iwork,
real*8, dimension(lenw)  work,
real*8  phi,
real*8  lambda1,
real*8  zk0,
real*8  Pup,
real*8  Tup,
real*8  rurd,
real*8  xflow,
real*8  kup 
)
7 c***begin prologue dqag
8 c***date written 800101 (yymmdd)
9 c***revision date 830518 (yymmdd)
10 c***category no. h2a1a1
11 c***keywords automatic integrator, general-purpose,
12 c integrand examinator, globally adaptive,
13 c gauss-kronrod
14 c***author piessens,robert,appl. math. & progr. div - k.u.leuven
15 c de doncker,elise,appl. math. & progr. div. - k.u.leuven
16 c***purpose the routine calculates an approximation result to a given
17 c definite integral i = integral of f over (a,b),
18 c hopefully satisfying following claim for accuracy
19 c abs(i-result)le.max(epsabs,epsrel*abs(i)).
20 c***description
21 c
22 c computation of a definite integral
23 c standard fortran subroutine
24 c double precision version
25 c
26 c f - double precision
27 c function subprogam defining the integrand
28 c function f(x). the actual name for f needs to be
29 c declared e x t e r n a l in the driver program.
30 c
31 c a - double precision
32 c lower limit of integration
33 c
34 c b - double precision
35 c upper limit of integration
36 c
37 c epsabs - double precision
38 c absolute accoracy requested
39 c epsrel - double precision
40 c relative accuracy requested
41 c if epsabs.le.0
42 c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
43 c the routine will end with ier = 6.
44 c
45 c key - integer
46 c key for choice of local integration rule
47 c a gauss-kronrod pair is used with
48 c 7 - 15 points if key.lt.2,
49 c 10 - 21 points if key = 2,
50 c 15 - 31 points if key = 3,
51 c 20 - 41 points if key = 4,
52 c 25 - 51 points if key = 5,
53 c 30 - 61 points if key.gt.5.
54 c
55 c on return
56 c result - double precision
57 c approximation to the integral
58 c
59 c abserr - double precision
60 c estimate of the modulus of the absolute error,
61 c which should equal or exceed abs(i-result)
62 c
63 c neval - integer
64 c number of integrand evaluations
65 c
66 c ier - integer
67 c ier = 0 normal and reliable termination of the
68 c routine. it is assumed that the requested
69 c accuracy has been achieved.
70 c ier.gt.0 abnormal termination of the routine
71 c the estimates for result and error are
72 c less reliable. it is assumed that the
73 c requested accuracy has not been achieved.
74 c error messages
75 c ier = 1 maximum number of subdivisions allowed
76 c has been achieved. one can allow more
77 c subdivisions by increasing the value of
78 c limit (and taking the according dimension
79 c adjustments into account). however, if
80 c this yield no improvement it is advised
81 c to analyze the integrand in order to
82 c determine the integration difficulaties.
83 c if the position of a local difficulty can
84 c be determined (i.e.singularity,
85 c discontinuity within the interval) one
86 c will probably gain from splitting up the
87 c interval at this point and calling the
88 c integrator on the subranges. if possible,
89 c an appropriate special-purpose integrator
90 c should be used which is designed for
91 c handling the type of difficulty involved.
92 c = 2 the occurrence of roundoff error is
93 c detected, which prevents the requested
94 c tolerance from being achieved.
95 c = 3 extremely bad integrand behaviour occurs
96 c at some points of the integration
97 c interval.
98 c = 6 the input is invalid, because
99 c (epsabs.le.0 and
100 c epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
101 c or limit.lt.1 or lenw.lt.limit*4.
102 c result, abserr, neval, last are set
103 c to zero.
104 c except when lenw is invalid, iwork(1),
105 c work(limit*2+1) and work(limit*3+1) are
106 c set to zero, work(1) is set to a and
107 c work(limit+1) to b.
108 c
109 c dimensioning parameters
110 c limit - integer
111 c dimensioning parameter for iwork
112 c limit determines the maximum number of subintervals
113 c in the partition of the given integration interval
114 c (a,b), limit.ge.1.
115 c if limit.lt.1, the routine will end with ier = 6.
116 c
117 c lenw - integer
118 c dimensioning parameter for work
119 c lenw must be at least limit*4.
120 c if lenw.lt.limit*4, the routine will end with
121 c ier = 6.
122 c
123 c last - integer
124 c on return, last equals the number of subintervals
125 c produced in the subdiviosion process, which
126 c determines the number of significant elements
127 c actually in the work arrays.
128 c
129 c work arrays
130 c iwork - integer
131 c vector of dimension at least limit, the first k
132 c elements of which contain pointers to the error
133 c estimates over the subintervals, such that
134 c work(limit*3+iwork(1)),... , work(limit*3+iwork(k))
135 c form a decreasing sequence with k = last if
136 c last.le.(limit/2+2), and k = limit+1-last otherwise
137 c
138 c work - double precision
139 c vector of dimension at least lenw
140 c on return
141 c work(1), ..., work(last) contain the left end
142 c points of the subintervals in the partition of
143 c (a,b),
144 c work(limit+1), ..., work(limit+last) contain the
145 c right end points,
146 c work(limit*2+1), ..., work(limit*2+last) contain
147 c the integral approximations over the subintervals,
148 c work(limit*3+1), ..., work(limit*3+last) contain
149 c the error estimates.
150 c
151 c***references (none)
152 c***routines called dqage,xerror
153 c***end prologue dqag
154  real*8 a,abserr,b,epsabs,epsrel,f,result,work,d1mach(4),
155  * phi,lambda1,zk0,pup,tup,rurd,xflow,kup
156  integer ier,iwork,key,last,lenw,limit,lvl,l1,l2,l3,neval
157 c
158  dimension iwork(limit),work(lenw)
159 c
160  external f
161 c
162 c check validity of lenw.
163 c
164  d1mach(1)=1e21
165  d1mach(2)=0d0
166  d1mach(3)=0d0
167  d1mach(4)=1e-21
168 c
169 c***first executable statement dqag
170  ier = 6
171  neval = 0
172  last = 0
173  result = 0.0d+00
174  abserr = 0.0d+00
175  if(limit.lt.1.or.lenw.lt.limit*4) go to 10
176 c
177 c prepare call for dqage.
178 c
179  l1 = limit+1
180  l2 = limit+l1
181  l3 = limit+l2
182 c
183  call dqage(f,a,b,epsabs,epsrel,key,limit,result,abserr,neval,
184  * ier,work(1),work(l1),work(l2),work(l3),iwork,last,phi,lambda1,
185  * zk0,pup,tup,rurd,xflow,kup)
186 c
187 c call error handler if necessary.
188 c
189  lvl = 0
190 10 if(ier.eq.6) lvl = 1
191 ! if(ier.ne.0) call xerror(26habnormal return from dqag ,26,ier,lvl)
192  return
d1mach
double precision function d1mach(I)
Definition: ddeabm.f:2012
dqage
subroutine dqage(f, a, b, epsabs, epsrel, key, limit, result, abserr, neval, ier, alist, blist, rlist, elist, iord, last, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
Definition: dqag.f:197
e21
static double * e21
Definition: radflowload.c:42

◆ dqage()

subroutine dqage ( external  f,
double precision  a,
double precision  b,
double precision  epsabs,
double precision  epsrel,
integer  key,
integer  limit,
double precision  result,
double precision  abserr,
integer  neval,
integer  ier,
double precision, dimension(limit)  alist,
double precision, dimension(limit)  blist,
double precision, dimension(limit)  rlist,
double precision, dimension(limit)  elist,
integer, dimension(limit)  iord,
integer  last,
double precision  phi,
double precision  lambda1,
double precision  zk0,
double precision  Pup,
double precision  Tup,
double precision  rurd,
double precision  xflow,
double precision  kup 
)
197 c***begin prologue dqage
198 c***date written 800101 (yymmdd)
199 c***revision date 830518 (yymmdd)
200 c***category no. h2a1a1
201 c***keywords automatic integrator, general-purpose,
202 c integrand examinator, globally adaptive,
203 c gauss-kronrod
204 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
205 c de doncker,elise,appl. math. & progr. div. - k.u.leuven
206 c***purpose the routine calculates an approximation result to a given
207 c definite integral i = integral of f over (a,b),
208 c hopefully satisfying following claim for accuracy
209 c abs(i-reslt).le.max(epsabs,epsrel*abs(i)).
210 c***description
211 c
212 c computation of a definite integral
213 c standard fortran subroutine
214 c double precision version
215 c
216 c parameters
217 c on entry
218 c f - double precision
219 c function subprogram defining the integrand
220 c function f(x). the actual name for f needs to be
221 c declared e x t e r n a l in the driver program.
222 c
223 c a - double precision
224 c lower limit of integration
225 c
226 c b - double precision
227 c upper limit of integration
228 c
229 c epsabs - double precision
230 c absolute accuracy requested
231 c epsrel - double precision
232 c relative accuracy requested
233 c if epsabs.le.0
234 c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
235 c the routine will end with ier = 6.
236 c
237 c key - integer
238 c key for choice of local integration rule
239 c a gauss-kronrod pair is used with
240 c 7 - 15 points if key.lt.2,
241 c 10 - 21 points if key = 2,
242 c 15 - 31 points if key = 3,
243 c 20 - 41 points if key = 4,
244 c 25 - 51 points if key = 5,
245 c 30 - 61 points if key.gt.5.
246 c
247 c limit - integer
248 c gives an upperbound on the number of subintervals
249 c in the partition of (a,b), limit.ge.1.
250 c
251 c on return
252 c result - double precision
253 c approximation to the integral
254 c
255 c abserr - double precision
256 c estimate of the modulus of the absolute error,
257 c which should equal or exceed abs(i-result)
258 c
259 c neval - integer
260 c number of integrand evaluations
261 c
262 c ier - integer
263 c ier = 0 normal and reliable termination of the
264 c routine. it is assumed that the requested
265 c accuracy has been achieved.
266 c ier.gt.0 abnormal termination of the routine
267 c the estimates for result and error are
268 c less reliable. it is assumed that the
269 c requested accuracy has not been achieved.
270 c error messages
271 c ier = 1 maximum number of subdivisions allowed
272 c has been achieved. one can allow more
273 c subdivisions by increasing the value
274 c of limit.
275 c however, if this yields no improvement it
276 c is rather advised to analyze the integrand
277 c in order to determine the integration
278 c difficulties. if the position of a local
279 c difficulty can be determined(e.g.
280 c singularity, discontinuity within the
281 c interval) one will probably gain from
282 c splitting up the interval at this point
283 c and calling the integrator on the
284 c subranges. if possible, an appropriate
285 c special-purpose integrator should be used
286 c which is designed for handling the type of
287 c difficulty involved.
288 c = 2 the occurrence of roundoff error is
289 c detected, which prevents the requested
290 c tolerance from being achieved.
291 c = 3 extremely bad integrand behaviour occurs
292 c at some points of the integration
293 c interval.
294 c = 6 the input is invalid, because
295 c (epsabs.le.0 and
296 c epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
297 c result, abserr, neval, last, rlist(1) ,
298 c elist(1) and iord(1) are set to zero.
299 c alist(1) and blist(1) are set to a and b
300 c respectively.
301 c
302 c alist - double precision
303 c vector of dimension at least limit, the first
304 c last elements of which are the left
305 c end points of the subintervals in the partition
306 c of the given integration range (a,b)
307 c
308 c blist - double precision
309 c vector of dimension at least limit, the first
310 c last elements of which are the right
311 c end points of the subintervals in the partition
312 c of the given integration range (a,b)
313 c
314 c rlist - double precision
315 c vector of dimension at least limit, the first
316 c last elements of which are the
317 c integral approximations on the subintervals
318 c
319 c elist - double precision
320 c vector of dimension at least limit, the first
321 c last elements of which are the moduli of the
322 c absolute error estimates on the subintervals
323 c
324 c iord - integer
325 c vector of dimension at least limit, the first k
326 c elements of which are pointers to the
327 c error estimates over the subintervals,
328 c such that elist(iord(1)), ...,
329 c elist(iord(k)) form a decreasing sequence,
330 c with k = last if last.le.(limit/2+2), and
331 c k = limit+1-last otherwise
332 c
333 c last - integer
334 c number of subintervals actually produced in the
335 c subdivision process
336 c
337 c***references (none)
338 c***routines called d1mach,dqk15,dqk21,dqk31,
339 c dqk41,dqk51,dqk61,dqpsrt
340 c***end prologue dqage
341 c
342  double precision a,abserr,alist,area,area1,area12,area2,a1,a2,b,
343  * blist,b1,b2,dabs,defabs,defab1,defab2,d1mach(4),elist,
344  * epmach,epsabs,epsrel,errbnd,errmax,error1,error2,erro12,errsum,f,
345  * resabs,result,rlist,uflow,phi,lambda1,zk0,pup,tup,rurd,xflow,kup
346  integer ier,iord,iroff1,iroff2,k,key,keyf,last,limit,maxerr,neval,
347  * nrmax
348 c
349  dimension alist(limit),blist(limit),elist(limit),iord(limit),
350  * rlist(limit)
351 c
352  external f
353 
354 
355  d1mach(1)=1e21
356  d1mach(2)=0d0
357  d1mach(3)=0d0
358  d1mach(4)=1e-21
359 c
360 c list of major variables
361 c -----------------------
362 c
363 c alist - list of left end points of all subintervals
364 c considered up to now
365 c blist - list of right end points of all subintervals
366 c considered up to now
367 c rlist(i) - approximation to the integral over
368 c (alist(i),blist(i))
369 c elist(i) - error estimate applying to rlist(i)
370 c maxerr - pointer to the interval with largest
371 c error estimate
372 c errmax - elist(maxerr)
373 c area - sum of the integrals over the subintervals
374 c errsum - sum of the errors over the subintervals
375 c errbnd - requested accuracy max(epsabs,epsrel*
376 c abs(result))
377 c *****1 - variable for the left subinterval
378 c *****2 - variable for the right subinterval
379 c last - index for subdivision
380 c
381 c
382 c machine dependent constants
383 c ---------------------------
384 c
385 c epmach is the largest relative spacing.
386 c uflow is the smallest positive magnitude.
387 c
388 c***first executable statement dqage
389  epmach = d1mach(4)
390  uflow = d1mach(1)
391 c
392 c test on validity of parameters
393 c ------------------------------
394 c
395  ier = 0
396  neval = 0
397  last = 0
398  result = 0.0d+00
399  abserr = 0.0d+00
400  alist(1) = a
401  blist(1) = b
402  rlist(1) = 0.0d+00
403  elist(1) = 0.0d+00
404  iord(1) = 0
405  if(epsabs.le.0.0d+00.and.
406  * epsrel.lt.max(0.5d+02*epmach,0.5d-28)) ier = 6
407  if(ier.eq.6) go to 999
408 c
409 c first approximation to the integral
410 c -----------------------------------
411 c
412  keyf = key
413  if(key.le.0) keyf = 1
414  if(key.ge.7) keyf = 6
415  neval = 0
416  if(keyf.eq.1) call dqk15(f,a,b,result,abserr,defabs,resabs,
417  & phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
418  if(keyf.eq.2) call dqk21(f,a,b,result,abserr,defabs,resabs,
419  & phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
420  if(keyf.eq.3) call dqk31(f,a,b,result,abserr,defabs,resabs,
421  & phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
422  if(keyf.eq.4) call dqk41(f,a,b,result,abserr,defabs,resabs,
423  & phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
424  if(keyf.eq.5) call dqk51(f,a,b,result,abserr,defabs,resabs,
425  & phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
426  if(keyf.eq.6) call dqk61(f,a,b,result,abserr,defabs,resabs,
427  & phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
428  last = 1
429  rlist(1) = result
430  elist(1) = abserr
431  iord(1) = 1
432 c
433 c test on accuracy.
434 c
435  errbnd = max(epsabs,epsrel*dabs(result))
436  if(abserr.le.0.5d+02*epmach*defabs.and.abserr.gt.errbnd) ier = 2
437  if(limit.eq.1) ier = 1
438  if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs)
439  * .or.abserr.eq.0.0d+00) go to 60
440 c
441 c initialization
442 c --------------
443 c
444 c
445  errmax = abserr
446  maxerr = 1
447  area = result
448  errsum = abserr
449  nrmax = 1
450  iroff1 = 0
451  iroff2 = 0
452 c
453 c main do-loop
454 c ------------
455 c
456  do 30 last = 2,limit
457 c
458 c bisect the subinterval with the largest error estimate.
459 c
460  a1 = alist(maxerr)
461  b1 = 0.5d+00*(alist(maxerr)+blist(maxerr))
462  a2 = b1
463  b2 = blist(maxerr)
464  if(keyf.eq.1) call dqk15(f,a1,b1,area1,error1,resabs,defab1,
465  * phi,lambda1,zk0,pup,tup,rurd,xflow,kup )
466  if(keyf.eq.2) call dqk21(f,a1,b1,area1,error1,resabs,defab1,
467  * phi,lambda1,zk0,pup,tup,rurd,xflow,kup )
468  if(keyf.eq.3) call dqk31(f,a1,b1,area1,error1,resabs,defab1,
469  * phi,lambda1,zk0,pup,tup,rurd,xflow,kup )
470  if(keyf.eq.4) call dqk41(f,a1,b1,area1,error1,resabs,defab1,
471  * phi,lambda1,zk0,pup,tup,rurd,xflow,kup )
472  if(keyf.eq.5) call dqk51(f,a1,b1,area1,error1,resabs,defab1,
473  * phi,lambda1,zk0,pup,tup,rurd,xflow,kup )
474  if(keyf.eq.6) call dqk61(f,a1,b1,area1,error1,resabs,defab1,
475  * phi,lambda1,zk0,pup,tup,rurd,xflow,kup )
476  if(keyf.eq.1) call dqk15(f,a2,b2,area2,error2,resabs,defab2,
477  * phi,lambda1,zk0,pup,tup,rurd,xflow,kup )
478  if(keyf.eq.2) call dqk21(f,a2,b2,area2,error2,resabs,defab2,
479  * phi,lambda1,zk0,pup,tup,rurd,xflow,kup )
480  if(keyf.eq.3) call dqk31(f,a2,b2,area2,error2,resabs,defab2,
481  * phi,lambda1,zk0,pup,tup,rurd,xflow,kup )
482  if(keyf.eq.4) call dqk41(f,a2,b2,area2,error2,resabs,defab2,
483  * phi,lambda1,zk0,pup,tup,rurd,xflow,kup )
484  if(keyf.eq.5) call dqk51(f,a2,b2,area2,error2,resabs,defab2,
485  * phi,lambda1,zk0,pup,tup,rurd,xflow,kup )
486  if(keyf.eq.6) call dqk61(f,a2,b2,area2,error2,resabs,defab2,
487  * phi,lambda1,zk0,pup,tup,rurd,xflow,kup )
488 c
489 c improve previous approximations to integral
490 c and error and test for accuracy.
491 c
492  neval = neval+1
493  area12 = area1+area2
494  erro12 = error1+error2
495  errsum = errsum+erro12-errmax
496  area = area+area12-rlist(maxerr)
497  if(defab1.eq.error1.or.defab2.eq.error2) go to 5
498  if(dabs(rlist(maxerr)-area12).le.0.1d-04*dabs(area12)
499  * .and.erro12.ge.0.99d+00*errmax) iroff1 = iroff1+1
500  if(last.gt.10.and.erro12.gt.errmax) iroff2 = iroff2+1
501  5 rlist(maxerr) = area1
502  rlist(last) = area2
503  errbnd = max(epsabs,epsrel*dabs(area))
504  if(errsum.le.errbnd) go to 8
505 c
506 c test for roundoff error and eventually set error flag.
507 c
508  if(iroff1.ge.6.or.iroff2.ge.20) ier = 2
509 c
510 c set error flag in the case that the number of subintervals
511 c equals limit.
512 c
513  if(last.eq.limit) ier = 1
514 c
515 c set error flag in the case of bad integrand behaviour
516 c at a point of the integration range.
517 c
518  if(max(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03*
519  * epmach)*(dabs(a2)+0.1d+04*uflow)) ier = 3
520 c
521 c append the newly-created intervals to the list.
522 c
523  8 if(error2.gt.error1) go to 10
524  alist(last) = a2
525  blist(maxerr) = b1
526  blist(last) = b2
527  elist(maxerr) = error1
528  elist(last) = error2
529  go to 20
530  10 alist(maxerr) = a2
531  alist(last) = a1
532  blist(last) = b1
533  rlist(maxerr) = area2
534  rlist(last) = area1
535  elist(maxerr) = error2
536  elist(last) = error1
537 c
538 c call subroutine dqpsrt to maintain the descending ordering
539 c in the list of error estimates and select the subinterval
540 c with the largest error estimate (to be bisected next).
541 c
542  20 call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax,phi,
543  * lambda1,zk0,pup,tup,rurd,xflow,kup)
544 c ***jump out of do-loop
545  if(ier.ne.0.or.errsum.le.errbnd) go to 40
546  30 continue
547 c
548 c compute final result.
549 c ---------------------
550 c
551  40 result = 0.0d+00
552  do 50 k=1,last
553  result = result+rlist(k)
554  50 continue
555  abserr = errsum
556  60 if(keyf.ne.1) neval = (10*keyf+1)*(2*neval+1)
557  if(keyf.eq.1) neval = 30*neval+15
558  999 return
max
#define max(a, b)
Definition: cascade.c:32
dqk31
subroutine dqk31(f, a, b, result, abserr, resabs, resasc, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
Definition: dqag.f:934
d1mach
double precision function d1mach(I)
Definition: ddeabm.f:2012
dqk51
subroutine dqk51(f, a, b, result, abserr, resabs, resasc, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
Definition: dqag.f:1345
dqk61
subroutine dqk61(f, a, b, result, abserr, resabs, resasc, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
Definition: dqag.f:1571
dqk41
subroutine dqk41(f, a, b, result, abserr, resabs, resasc, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
Definition: dqag.f:1132
e21
static double * e21
Definition: radflowload.c:42
dqpsrt
subroutine dqpsrt(limit, last, maxerr, ermax, elist, iord, nrmax, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
Definition: dqag.f:1808
b1
static double * b1
Definition: mafillkmain.c:30
dqk21
subroutine dqk21(f, a, b, result, abserr, resabs, resasc, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
Definition: dqag.f:745
area1
static double * area1
Definition: mafillkmain.c:30
dqk15
subroutine dqk15(f, a, b, result, abserr, resabs, resasc, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
Definition: dqag.f:562

◆ dqk15()

subroutine dqk15 ( external  f,
double precision  a,
double precision  b,
double precision  result,
double precision  abserr,
double precision  resabs,
double precision  resasc,
double precision  phi,
double precision  lambda1,
double precision  zk0,
double precision  Pup,
double precision  Tup,
double precision  rurd,
double precision  xflow,
double precision  kup 
)
562 c***begin prologue dqk15
563 c***date written 800101 (yymmdd)
564 c***revision date 830518 (yymmdd)
565 c***category no. h2a1a2
566 c***keywords 15-point gauss-kronrod rules
567 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
568 c de doncker,elise,appl. math. & progr. div - k.u.leuven
569 c***purpose to compute i = integral of f over (a,b), with error
570 c estimate
571 c j = integral of abs(f) over (a,b)
572 c***description
573 c
574 c integration rules
575 c standard fortran subroutine
576 c double precision version
577 c
578 c parameters
579 c on entry
580 c f - double precision
581 c function subprogram defining the integrand
582 c function f(x). the actual name for f needs to be
583 c declared e x t e r n a l in the calling program.
584 c
585 c a - double precision
586 c lower limit of integration
587 c
588 c b - double precision
589 c upper limit of integration
590 c
591 c on return
592 c result - double precision
593 c approximation to the integral i
594 c result is computed by applying the 15-point
595 c kronrod rule (resk) obtained by optimal addition
596 c of abscissae to the7-point gauss rule(resg).
597 c
598 c abserr - double precision
599 c estimate of the modulus of the absolute error,
600 c which should not exceed abs(i-result)
601 c
602 c resabs - double precision
603 c approximation to the integral j
604 c
605 c resasc - double precision
606 c approximation to the integral of abs(f-i/(b-a))
607 c over (a,b)
608 c
609 c***references (none)
610 c***routines called d1mach
611 c***end prologue dqk15
612 c
613  double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
614  * d1mach(4),epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,
615  * resasc,resg,resk,reskh,result,uflow,wg,wgk,xgk,phi,lambda1,
616  * zk0,pup,tup,rurd,xflow,kup
617  integer j,jtw,jtwm1
618  external f
619 c
620  dimension fv1(7),fv2(7),wg(4),wgk(8),xgk(8)
621 
622  d1mach(1)=1e21
623  d1mach(2)=0d0
624  d1mach(3)=0d0
625  d1mach(4)=1e-21
626 
627 
628 c
629 c the abscissae and weights are given for the interval (-1,1).
630 c because of symmetry only the positive abscissae and their
631 c corresponding weights are given.
632 c
633 c xgk - abscissae of the 15-point kronrod rule
634 c xgk(2), xgk(4), ... abscissae of the 7-point
635 c gauss rule
636 c xgk(1), xgk(3), ... abscissae which are optimally
637 c added to the 7-point gauss rule
638 c
639 c wgk - weights of the 15-point kronrod rule
640 c
641 c wg - weights of the 7-point gauss rule
642 c
643 c
644 c gauss quadrature weights and kronron quadrature abscissae and weights
645 c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
646 c bell labs, nov. 1981.
647 c
648  data wg( 1) / 0.1294849661 6886969327 0611432679 082 d0 /
649  data wg( 2) / 0.2797053914 8927666790 1467771423 780 d0 /
650  data wg( 3) / 0.3818300505 0511894495 0369775488 975 d0 /
651  data wg( 4) / 0.4179591836 7346938775 5102040816 327 d0 /
652 c
653  data xgk( 1) / 0.9914553711 2081263920 6854697526 329 d0 /
654  data xgk( 2) / 0.9491079123 4275852452 6189684047 851 d0 /
655  data xgk( 3) / 0.8648644233 5976907278 9712788640 926 d0 /
656  data xgk( 4) / 0.7415311855 9939443986 3864773280 788 d0 /
657  data xgk( 5) / 0.5860872354 6769113029 4144838258 730 d0 /
658  data xgk( 6) / 0.4058451513 7739716690 6606412076 961 d0 /
659  data xgk( 7) / 0.2077849550 0789846760 0689403773 245 d0 /
660  data xgk( 8) / 0.0000000000 0000000000 0000000000 000 d0 /
661 c
662  data wgk( 1) / 0.0229353220 1052922496 3732008058 970 d0 /
663  data wgk( 2) / 0.0630920926 2997855329 0700663189 204 d0 /
664  data wgk( 3) / 0.1047900103 2225018383 9876322541 518 d0 /
665  data wgk( 4) / 0.1406532597 1552591874 5189590510 238 d0 /
666  data wgk( 5) / 0.1690047266 3926790282 6583426598 550 d0 /
667  data wgk( 6) / 0.1903505780 6478540991 3256402421 014 d0 /
668  data wgk( 7) / 0.2044329400 7529889241 4161999234 649 d0 /
669  data wgk( 8) / 0.2094821410 8472782801 2999174891 714 d0 /
670 c
671 c
672 c list of major variables
673 c -----------------------
674 c
675 c centr - mid point of the interval
676 c hlgth - half-length of the interval
677 c absc - abscissa
678 c fval* - function value
679 c resg - result of the 7-point gauss formula
680 c resk - result of the 15-point kronrod formula
681 c reskh - approximation to the mean value of f over (a,b),
682 c i.e. to i/(b-a)
683 c
684 c machine dependent constants
685 c ---------------------------
686 c
687 c epmach is the largest relative spacing.
688 c uflow is the smallest positive magnitude.
689 c
690 c***first executable statement dqk15
691  epmach = d1mach(4)
692  uflow = d1mach(1)
693 c
694  centr = 0.5d+00*(a+b)
695  hlgth = 0.5d+00*(b-a)
696  dhlgth = dabs(hlgth)
697 c
698 c compute the 15-point kronrod approximation to
699 c the integral, and estimate the absolute error.
700 c
701  fc = f(centr,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
702  resg = fc*wg(4)
703  resk = fc*wgk(8)
704  resabs = dabs(resk)
705  do 10 j=1,3
706  jtw = j*2
707  absc = hlgth*xgk(jtw)
708  fval1 = f(centr-absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
709  fval2 = f(centr+absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
710  fv1(jtw) = fval1
711  fv2(jtw) = fval2
712  fsum = fval1+fval2
713  resg = resg+wg(j)*fsum
714  resk = resk+wgk(jtw)*fsum
715  resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
716  10 continue
717  do 15 j = 1,4
718  jtwm1 = j*2-1
719  absc = hlgth*xgk(jtwm1)
720  fval1 = f(centr-absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
721  fval2 = f(centr+absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
722  fv1(jtwm1) = fval1
723  fv2(jtwm1) = fval2
724  fsum = fval1+fval2
725  resk = resk+wgk(jtwm1)*fsum
726  resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
727  15 continue
728  reskh = resk*0.5d+00
729  resasc = wgk(8)*dabs(fc-reskh)
730  do 20 j=1,7
731  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
732  20 continue
733  result = resk*hlgth
734  resabs = resabs*dhlgth
735  resasc = resasc*dhlgth
736  abserr = dabs((resk-resg)*hlgth)
737  if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
738  * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
739  if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
740  * ((epmach*0.5d+02)*resabs,abserr)
741  return
d1mach
double precision function d1mach(I)
Definition: ddeabm.f:2012
e21
static double * e21
Definition: radflowload.c:42

◆ dqk21()

subroutine dqk21 ( external  f,
double precision  a,
double precision  b,
double precision  result,
double precision  abserr,
double precision  resabs,
double precision  resasc,
double precision  phi,
double precision  lambda1,
double precision  zk0,
double precision  Pup,
double precision  Tup,
double precision  rurd,
double precision  xflow,
double precision  kup 
)
745 c***begin prologue dqk21
746 c***date written 800101 (yymmdd)
747 c***revision date 830518 (yymmdd)
748 c***category no. h2a1a2
749 c***keywords 21-point gauss-kronrod rules
750 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
751 c de doncker,elise,appl. math. & progr. div. - k.u.leuven
752 c***purpose to compute i = integral of f over (a,b), with error
753 c estimate
754 c j = integral of abs(f) over (a,b)
755 c***description
756 c
757 c integration rules
758 c standard fortran subroutine
759 c double precision version
760 c
761 c parameters
762 c on entry
763 c f - double precision
764 c function subprogram defining the integrand
765 c function f(x). the actual name for f needs to be
766 c declared e x t e r n a l in the driver program.
767 c
768 c a - double precision
769 c lower limit of integration
770 c
771 c b - double precision
772 c upper limit of integration
773 c
774 c on return
775 c result - double precision
776 c approximation to the integral i
777 c result is computed by applying the 21-point
778 c kronrod rule (resk) obtained by optimal addition
779 c of abscissae to the 10-point gauss rule (resg).
780 c
781 c abserr - double precision
782 c estimate of the modulus of the absolute error,
783 c which should not exceed abs(i-result)
784 c
785 c resabs - double precision
786 c approximation to the integral j
787 c
788 c resasc - double precision
789 c approximation to the integral of abs(f-i/(b-a))
790 c over (a,b)
791 c
792 c***references (none)
793 c***routines called d1mach
794 c***end prologue dqk21
795 c
796  double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
797  * d1mach(4),epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,
798  * resasc,resg,resk,reskh,result,uflow,wg,wgk,xgk,phi,lambda1,
799  * zk0,pup,tup,rurd,xflow,kup
800  integer j,jtw,jtwm1
801  external f
802 c
803  dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11)
804 
805  d1mach(1)=1e21
806  d1mach(2)=0d0
807  d1mach(3)=0d0
808  d1mach(4)=1e-21
809 c
810 c the abscissae and weights are given for the interval (-1,1).
811 c because of symmetry only the positive abscissae and their
812 c corresponding weights are given.
813 c
814 c xgk - abscissae of the 21-point kronrod rule
815 c xgk(2), xgk(4), ... abscissae of the 10-point
816 c gauss rule
817 c xgk(1), xgk(3), ... abscissae which are optimally
818 c added to the 10-point gauss rule
819 c
820 c wgk - weights of the 21-point kronrod rule
821 c
822 c wg - weights of the 10-point gauss rule
823 c
824 c
825 c gauss quadrature weights and kronron quadrature abscissae and weights
826 c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
827 c bell labs, nov. 1981.
828 c
829  data wg( 1) / 0.0666713443 0868813759 3568809893 332 d0 /
830  data wg( 2) / 0.1494513491 5058059314 5776339657 697 d0 /
831  data wg( 3) / 0.2190863625 1598204399 5534934228 163 d0 /
832  data wg( 4) / 0.2692667193 0999635509 1226921569 469 d0 /
833  data wg( 5) / 0.2955242247 1475287017 3892994651 338 d0 /
834 c
835  data xgk( 1) / 0.9956571630 2580808073 5527280689 003 d0 /
836  data xgk( 2) / 0.9739065285 1717172007 7964012084 452 d0 /
837  data xgk( 3) / 0.9301574913 5570822600 1207180059 508 d0 /
838  data xgk( 4) / 0.8650633666 8898451073 2096688423 493 d0 /
839  data xgk( 5) / 0.7808177265 8641689706 3717578345 042 d0 /
840  data xgk( 6) / 0.6794095682 9902440623 4327365114 874 d0 /
841  data xgk( 7) / 0.5627571346 6860468333 9000099272 694 d0 /
842  data xgk( 8) / 0.4333953941 2924719079 9265943165 784 d0 /
843  data xgk( 9) / 0.2943928627 0146019813 1126603103 866 d0 /
844  data xgk( 10) / 0.1488743389 8163121088 4826001129 720 d0 /
845  data xgk( 11) / 0.0000000000 0000000000 0000000000 000 d0 /
846 c
847  data wgk( 1) / 0.0116946388 6737187427 8064396062 192 d0 /
848  data wgk( 2) / 0.0325581623 0796472747 8818972459 390 d0 /
849  data wgk( 3) / 0.0547558965 7435199603 1381300244 580 d0 /
850  data wgk( 4) / 0.0750396748 1091995276 7043140916 190 d0 /
851  data wgk( 5) / 0.0931254545 8369760553 5065465083 366 d0 /
852  data wgk( 6) / 0.1093871588 0229764189 9210590325 805 d0 /
853  data wgk( 7) / 0.1234919762 6206585107 7958109831 074 d0 /
854  data wgk( 8) / 0.1347092173 1147332592 8054001771 707 d0 /
855  data wgk( 9) / 0.1427759385 7706008079 7094273138 717 d0 /
856  data wgk( 10) / 0.1477391049 0133849137 4841515972 068 d0 /
857  data wgk( 11) / 0.1494455540 0291690566 4936468389 821 d0 /
858 c
859 c
860 c list of major variables
861 c -----------------------
862 c
863 c centr - mid point of the interval
864 c hlgth - half-length of the interval
865 c absc - abscissa
866 c fval* - function value
867 c resg - result of the 10-point gauss formula
868 c resk - result of the 21-point kronrod formula
869 c reskh - approximation to the mean value of f over (a,b),
870 c i.e. to i/(b-a)
871 c
872 c
873 c machine dependent constants
874 c ---------------------------
875 c
876 c epmach is the largest relative spacing.
877 c uflow is the smallest positive magnitude.
878 c
879 c***first executable statement dqk21
880  epmach = d1mach(4)
881  uflow = d1mach(1)
882 c
883  centr = 0.5d+00*(a+b)
884  hlgth = 0.5d+00*(b-a)
885  dhlgth = dabs(hlgth)
886 c
887 c compute the 21-point kronrod approximation to
888 c the integral, and estimate the absolute error.
889 c
890  resg = 0.0d+00
891  fc = f(centr,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
892  resk = wgk(11)*fc
893  resabs = dabs(resk)
894  do 10 j=1,5
895  jtw = 2*j
896  absc = hlgth*xgk(jtw)
897  fval1 = f(centr-absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
898  fval2 = f(centr+absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
899  fv1(jtw) = fval1
900  fv2(jtw) = fval2
901  fsum = fval1+fval2
902  resg = resg+wg(j)*fsum
903  resk = resk+wgk(jtw)*fsum
904  resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
905  10 continue
906  do 15 j = 1,5
907  jtwm1 = 2*j-1
908  absc = hlgth*xgk(jtwm1)
909  fval1 = f(centr-absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
910  fval2 = f(centr+absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
911  fv1(jtwm1) = fval1
912  fv2(jtwm1) = fval2
913  fsum = fval1+fval2
914  resk = resk+wgk(jtwm1)*fsum
915  resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
916  15 continue
917  reskh = resk*0.5d+00
918  resasc = wgk(11)*dabs(fc-reskh)
919  do 20 j=1,10
920  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
921  20 continue
922  result = resk*hlgth
923  resabs = resabs*dhlgth
924  resasc = resasc*dhlgth
925  abserr = dabs((resk-resg)*hlgth)
926  if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
927  * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
928  if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
929  * ((epmach*0.5d+02)*resabs,abserr)
930  return
d1mach
double precision function d1mach(I)
Definition: ddeabm.f:2012
e21
static double * e21
Definition: radflowload.c:42

◆ dqk31()

subroutine dqk31 ( external  f,
double precision  a,
double precision  b,
double precision  result,
double precision  abserr,
double precision  resabs,
double precision  resasc,
double precision  phi,
double precision  lambda1,
double precision  zk0,
double precision  Pup,
double precision  Tup,
double precision  rurd,
double precision  xflow,
double precision  kup 
)
934 c***begin prologue dqk31
935 c***date written 800101 (yymmdd)
936 c***revision date 830518 (yymmdd)
937 c***category no. h2a1a2
938 c***keywords 31-point gauss-kronrod rules
939 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
940 c de doncker,elise,appl. math. & progr. div. - k.u.leuven
941 c***purpose to compute i = integral of f over (a,b) with error
942 c estimate
943 c j = integral of abs(f) over (a,b)
944 c***description
945 c
946 c integration rules
947 c standard fortran subroutine
948 c double precision version
949 c
950 c parameters
951 c on entry
952 c f - double precision
953 c function subprogram defining the integrand
954 c function f(x). the actual name for f needs to be
955 c declared e x t e r n a l in the calling program.
956 c
957 c a - double precision
958 c lower limit of integration
959 c
960 c b - double precision
961 c upper limit of integration
962 c
963 c on return
964 c result - double precision
965 c approximation to the integral i
966 c result is computed by applying the 31-point
967 c gauss-kronrod rule (resk), obtained by optimal
968 c addition of abscissae to the 15-point gauss
969 c rule (resg).
970 c
971 c abserr - double precison
972 c estimate of the modulus of the modulus,
973 c which should not exceed abs(i-result)
974 c
975 c resabs - double precision
976 c approximation to the integral j
977 c
978 c resasc - double precision
979 c approximation to the integral of abs(f-i/(b-a))
980 c over (a,b)
981 c
982 c***references (none)
983 c***routines called d1mach
984 c***end prologue dqk31
985  double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
986  * d1mach(4),epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,
987  * resasc,resg,resk,reskh,result,uflow,wg,wgk,xgk,phi,lambda1,
988  * zk0,pup,tup,rurd,xflow,kup
989  integer j,jtw,jtwm1
990  external f
991 c
992  dimension fv1(15),fv2(15),xgk(16),wgk(16),wg(8)
993 
994  d1mach(1)=1e21
995  d1mach(2)=0d0
996  d1mach(3)=0d0
997  d1mach(4)=1e-21
998 c
999 c the abscissae and weights are given for the interval (-1,1).
1000 c because of symmetry only the positive abscissae and their
1001 c corresponding weights are given.
1002 c
1003 c xgk - abscissae of the 31-point kronrod rule
1004 c xgk(2), xgk(4), ... abscissae of the 15-point
1005 c gauss rule
1006 c xgk(1), xgk(3), ... abscissae which are optimally
1007 c added to the 15-point gauss rule
1008 c
1009 c wgk - weights of the 31-point kronrod rule
1010 c
1011 c wg - weights of the 15-point gauss rule
1012 c
1013 c
1014 c gauss quadrature weights and kronron quadrature abscissae and weights
1015 c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
1016 c bell labs, nov. 1981.
1017 c
1018  data wg( 1) / 0.0307532419 9611726835 4628393577 204 d0 /
1019  data wg( 2) / 0.0703660474 8810812470 9267416450 667 d0 /
1020  data wg( 3) / 0.1071592204 6717193501 1869546685 869 d0 /
1021  data wg( 4) / 0.1395706779 2615431444 7804794511 028 d0 /
1022  data wg( 5) / 0.1662692058 1699393355 3200860481 209 d0 /
1023  data wg( 6) / 0.1861610000 1556221102 6800561866 423 d0 /
1024  data wg( 7) / 0.1984314853 2711157645 6118326443 839 d0 /
1025  data wg( 8) / 0.2025782419 2556127288 0620199967 519 d0 /
1026 c
1027  data xgk( 1) / 0.9980022986 9339706028 5172840152 271 d0 /
1028  data xgk( 2) / 0.9879925180 2048542848 9565718586 613 d0 /
1029  data xgk( 3) / 0.9677390756 7913913425 7347978784 337 d0 /
1030  data xgk( 4) / 0.9372733924 0070590430 7758947710 209 d0 /
1031  data xgk( 5) / 0.8972645323 4408190088 2509656454 496 d0 /
1032  data xgk( 6) / 0.8482065834 1042721620 0648320774 217 d0 /
1033  data xgk( 7) / 0.7904185014 4246593296 7649294817 947 d0 /
1034  data xgk( 8) / 0.7244177313 6017004741 6186054613 938 d0 /
1035  data xgk( 9) / 0.6509967412 9741697053 3735895313 275 d0 /
1036  data xgk( 10) / 0.5709721726 0853884753 7226737253 911 d0 /
1037  data xgk( 11) / 0.4850818636 4023968069 3655740232 351 d0 /
1038  data xgk( 12) / 0.3941513470 7756336989 7207370981 045 d0 /
1039  data xgk( 13) / 0.2991800071 5316881216 6780024266 389 d0 /
1040  data xgk( 14) / 0.2011940939 9743452230 0628303394 596 d0 /
1041  data xgk( 15) / 0.1011420669 1871749902 7074231447 392 d0 /
1042  data xgk( 16) / 0.0000000000 0000000000 0000000000 000 d0 /
1043 c
1044  data wgk( 1) / 0.0053774798 7292334898 7792051430 128 d0 /
1045  data wgk( 2) / 0.0150079473 2931612253 8374763075 807 d0 /
1046  data wgk( 3) / 0.0254608473 2671532018 6874001019 653 d0 /
1047  data wgk( 4) / 0.0353463607 9137584622 2037948478 360 d0 /
1048  data wgk( 5) / 0.0445897513 2476487660 8227299373 280 d0 /
1049  data wgk( 6) / 0.0534815246 9092808726 5343147239 430 d0 /
1050  data wgk( 7) / 0.0620095678 0067064028 5139230960 803 d0 /
1051  data wgk( 8) / 0.0698541213 1872825870 9520077099 147 d0 /
1052  data wgk( 9) / 0.0768496807 5772037889 4432777482 659 d0 /
1053  data wgk( 10) / 0.0830805028 2313302103 8289247286 104 d0 /
1054  data wgk( 11) / 0.0885644430 5621177064 7275443693 774 d0 /
1055  data wgk( 12) / 0.0931265981 7082532122 5486872747 346 d0 /
1056  data wgk( 13) / 0.0966427269 8362367850 5179907627 589 d0 /
1057  data wgk( 14) / 0.0991735987 2179195933 2393173484 603 d0 /
1058  data wgk( 15) / 0.1007698455 2387559504 4946662617 570 d0 /
1059  data wgk( 16) / 0.1013300070 1479154901 7374792767 493 d0 /
1060 c
1061 c
1062 c list of major variables
1063 c -----------------------
1064 c centr - mid point of the interval
1065 c hlgth - half-length of the interval
1066 c absc - abscissa
1067 c fval* - function value
1068 c resg - result of the 15-point gauss formula
1069 c resk - result of the 31-point kronrod formula
1070 c reskh - approximation to the mean value of f over (a,b),
1071 c i.e. to i/(b-a)
1072 c
1073 c machine dependent constants
1074 c ---------------------------
1075 c epmach is the largest relative spacing.
1076 c uflow is the smallest positive magnitude.
1077 c***first executable statement dqk31
1078  epmach = d1mach(4)
1079  uflow = d1mach(1)
1080 c
1081  centr = 0.5d+00*(a+b)
1082  hlgth = 0.5d+00*(b-a)
1083  dhlgth = dabs(hlgth)
1084 c
1085 c compute the 31-point kronrod approximation to
1086 c the integral, and estimate the absolute error.
1087 c
1088  fc = f(centr,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1089  resg = wg(8)*fc
1090  resk = wgk(16)*fc
1091  resabs = dabs(resk)
1092  do 10 j=1,7
1093  jtw = j*2
1094  absc = hlgth*xgk(jtw)
1095  fval1 = f(centr-absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1096  fval2 = f(centr+absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1097  fv1(jtw) = fval1
1098  fv2(jtw) = fval2
1099  fsum = fval1+fval2
1100  resg = resg+wg(j)*fsum
1101  resk = resk+wgk(jtw)*fsum
1102  resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
1103  10 continue
1104  do 15 j = 1,8
1105  jtwm1 = j*2-1
1106  absc = hlgth*xgk(jtwm1)
1107  fval1 = f(centr-absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1108  fval2 = f(centr+absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1109  fv1(jtwm1) = fval1
1110  fv2(jtwm1) = fval2
1111  fsum = fval1+fval2
1112  resk = resk+wgk(jtwm1)*fsum
1113  resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
1114  15 continue
1115  reskh = resk*0.5d+00
1116  resasc = wgk(16)*dabs(fc-reskh)
1117  do 20 j=1,15
1118  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
1119  20 continue
1120  result = resk*hlgth
1121  resabs = resabs*dhlgth
1122  resasc = resasc*dhlgth
1123  abserr = dabs((resk-resg)*hlgth)
1124  if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
1125  * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
1126  if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
1127  * ((epmach*0.5d+02)*resabs,abserr)
1128  return
d1mach
double precision function d1mach(I)
Definition: ddeabm.f:2012
e21
static double * e21
Definition: radflowload.c:42

◆ dqk41()

subroutine dqk41 ( external  f,
double precision  a,
double precision  b,
double precision  result,
double precision  abserr,
double precision  resabs,
double precision  resasc,
double precision  phi,
double precision  lambda1,
double precision  zk0,
double precision  Pup,
double precision  Tup,
double precision  rurd,
double precision  xflow,
double precision  kup 
)
1132 c***begin prologue dqk41
1133 c***date written 800101 (yymmdd)
1134 c***revision date 830518 (yymmdd)
1135 c***category no. h2a1a2
1136 c***keywords 41-point gauss-kronrod rules
1137 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
1138 c de doncker,elise,appl. math. & progr. div. - k.u.leuven
1139 c***purpose to compute i = integral of f over (a,b), with error
1140 c estimate
1141 c j = integral of abs(f) over (a,b)
1142 c***description
1143 c
1144 c integration rules
1145 c standard fortran subroutine
1146 c double precision version
1147 c
1148 c parameters
1149 c on entry
1150 c f - double precision
1151 c function subprogram defining the integrand
1152 c function f(x). the actual name for f needs to be
1153 c declared e x t e r n a l in the calling program.
1154 c
1155 c a - double precision
1156 c lower limit of integration
1157 c
1158 c b - double precision
1159 c upper limit of integration
1160 c
1161 c on return
1162 c result - double precision
1163 c approximation to the integral i
1164 c result is computed by applying the 41-point
1165 c gauss-kronrod rule (resk) obtained by optimal
1166 c addition of abscissae to the 20-point gauss
1167 c rule (resg).
1168 c
1169 c abserr - double precision
1170 c estimate of the modulus of the absolute error,
1171 c which should not exceed abs(i-result)
1172 c
1173 c resabs - double precision
1174 c approximation to the integral j
1175 c
1176 c resasc - double precision
1177 c approximation to the integal of abs(f-i/(b-a))
1178 c over (a,b)
1179 c
1180 c***references (none)
1181 c***routines called d1mach
1182 c***end prologue dqk41
1183 c
1184  double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
1185  * d1mach(4),epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,
1186  * resasc,resg,resk,reskh,result,uflow,wg,wgk,xgk,phi,lambda1,
1187  * zk0,pup,tup,rurd,xflow,kup
1188  integer j,jtw,jtwm1
1189  external f
1190 c
1191  dimension fv1(20),fv2(20),xgk(21),wgk(21),wg(10)
1192  d1mach(1)=1e21
1193  d1mach(2)=0d0
1194  d1mach(3)=0d0
1195  d1mach(4)=1e-21
1196 c
1197 c the abscissae and weights are given for the interval (-1,1).
1198 c because of symmetry only the positive abscissae and their
1199 c corresponding weights are given.
1200 c
1201 c xgk - abscissae of the 41-point gauss-kronrod rule
1202 c xgk(2), xgk(4), ... abscissae of the 20-point
1203 c gauss rule
1204 c xgk(1), xgk(3), ... abscissae which are optimally
1205 c added to the 20-point gauss rule
1206 c
1207 c wgk - weights of the 41-point gauss-kronrod rule
1208 c
1209 c wg - weights of the 20-point gauss rule
1210 c
1211 c
1212 c gauss quadrature weights and kronron quadrature abscissae and weights
1213 c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
1214 c bell labs, nov. 1981.
1215 c
1216  data wg( 1) / 0.0176140071 3915211831 1861962351 853 d0 /
1217  data wg( 2) / 0.0406014298 0038694133 1039952274 932 d0 /
1218  data wg( 3) / 0.0626720483 3410906356 9506535187 042 d0 /
1219  data wg( 4) / 0.0832767415 7670474872 4758143222 046 d0 /
1220  data wg( 5) / 0.1019301198 1724043503 6750135480 350 d0 /
1221  data wg( 6) / 0.1181945319 6151841731 2377377711 382 d0 /
1222  data wg( 7) / 0.1316886384 4917662689 8494499748 163 d0 /
1223  data wg( 8) / 0.1420961093 1838205132 9298325067 165 d0 /
1224  data wg( 9) / 0.1491729864 7260374678 7828737001 969 d0 /
1225  data wg( 10) / 0.1527533871 3072585069 8084331955 098 d0 /
1226 c
1227  data xgk( 1) / 0.9988590315 8827766383 8315576545 863 d0 /
1228  data xgk( 2) / 0.9931285991 8509492478 6122388471 320 d0 /
1229  data xgk( 3) / 0.9815078774 5025025919 3342994720 217 d0 /
1230  data xgk( 4) / 0.9639719272 7791379126 7666131197 277 d0 /
1231  data xgk( 5) / 0.9408226338 3175475351 9982722212 443 d0 /
1232  data xgk( 6) / 0.9122344282 5132590586 7752441203 298 d0 /
1233  data xgk( 7) / 0.8782768112 5228197607 7442995113 078 d0 /
1234  data xgk( 8) / 0.8391169718 2221882339 4529061701 521 d0 /
1235  data xgk( 9) / 0.7950414288 3755119835 0638833272 788 d0 /
1236  data xgk( 10) / 0.7463319064 6015079261 4305070355 642 d0 /
1237  data xgk( 11) / 0.6932376563 3475138480 5490711845 932 d0 /
1238  data xgk( 12) / 0.6360536807 2651502545 2836696226 286 d0 /
1239  data xgk( 13) / 0.5751404468 1971031534 2946036586 425 d0 /
1240  data xgk( 14) / 0.5108670019 5082709800 4364050955 251 d0 /
1241  data xgk( 15) / 0.4435931752 3872510319 9992213492 640 d0 /
1242  data xgk( 16) / 0.3737060887 1541956067 2548177024 927 d0 /
1243  data xgk( 17) / 0.3016278681 1491300432 0555356858 592 d0 /
1244  data xgk( 18) / 0.2277858511 4164507808 0496195368 575 d0 /
1245  data xgk( 19) / 0.1526054652 4092267550 5220241022 678 d0 /
1246  data xgk( 20) / 0.0765265211 3349733375 4640409398 838 d0 /
1247  data xgk( 21) / 0.0000000000 0000000000 0000000000 000 d0 /
1248 c
1249  data wgk( 1) / 0.0030735837 1852053150 1218293246 031 d0 /
1250  data wgk( 2) / 0.0086002698 5564294219 8661787950 102 d0 /
1251  data wgk( 3) / 0.0146261692 5697125298 3787960308 868 d0 /
1252  data wgk( 4) / 0.0203883734 6126652359 8010231432 755 d0 /
1253  data wgk( 5) / 0.0258821336 0495115883 4505067096 153 d0 /
1254  data wgk( 6) / 0.0312873067 7703279895 8543119323 801 d0 /
1255  data wgk( 7) / 0.0366001697 5820079803 0557240707 211 d0 /
1256  data wgk( 8) / 0.0416688733 2797368626 3788305936 895 d0 /
1257  data wgk( 9) / 0.0464348218 6749767472 0231880926 108 d0 /
1258  data wgk( 10) / 0.0509445739 2372869193 2707670050 345 d0 /
1259  data wgk( 11) / 0.0551951053 4828599474 4832372419 777 d0 /
1260  data wgk( 12) / 0.0591114008 8063957237 4967220648 594 d0 /
1261  data wgk( 13) / 0.0626532375 5478116802 5870122174 255 d0 /
1262  data wgk( 14) / 0.0658345971 3361842211 1563556969 398 d0 /
1263  data wgk( 15) / 0.0686486729 2852161934 5623411885 368 d0 /
1264  data wgk( 16) / 0.0710544235 5344406830 5790361723 210 d0 /
1265  data wgk( 17) / 0.0730306903 3278666749 5189417658 913 d0 /
1266  data wgk( 18) / 0.0745828754 0049918898 6581418362 488 d0 /
1267  data wgk( 19) / 0.0757044976 8455667465 9542775376 617 d0 /
1268  data wgk( 20) / 0.0763778676 7208073670 5502835038 061 d0 /
1269  data wgk( 21) / 0.0766007119 1799965644 5049901530 102 d0 /
1270 c
1271 c
1272 c list of major variables
1273 c -----------------------
1274 c
1275 c centr - mid point of the interval
1276 c hlgth - half-length of the interval
1277 c absc - abscissa
1278 c fval* - function value
1279 c resg - result of the 20-point gauss formula
1280 c resk - result of the 41-point kronrod formula
1281 c reskh - approximation to mean value of f over (a,b), i.e.
1282 c to i/(b-a)
1283 c
1284 c machine dependent constants
1285 c ---------------------------
1286 c
1287 c epmach is the largest relative spacing.
1288 c uflow is the smallest positive magnitude.
1289 c
1290 c***first executable statement dqk41
1291  epmach = d1mach(4)
1292  uflow = d1mach(1)
1293 c
1294  centr = 0.5d+00*(a+b)
1295  hlgth = 0.5d+00*(b-a)
1296  dhlgth = dabs(hlgth)
1297 c
1298 c compute the 41-point gauss-kronrod approximation to
1299 c the integral, and estimate the absolute error.
1300 c
1301  resg = 0.0d+00
1302  fc = f(centr,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1303  resk = wgk(21)*fc
1304  resabs = dabs(resk)
1305  do 10 j=1,10
1306  jtw = j*2
1307  absc = hlgth*xgk(jtw)
1308  fval1 = f(centr-absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1309  fval2 = f(centr+absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1310  fv1(jtw) = fval1
1311  fv2(jtw) = fval2
1312  fsum = fval1+fval2
1313  resg = resg+wg(j)*fsum
1314  resk = resk+wgk(jtw)*fsum
1315  resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
1316  10 continue
1317  do 15 j = 1,10
1318  jtwm1 = j*2-1
1319  absc = hlgth*xgk(jtwm1)
1320  fval1 = f(centr-absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1321  fval2 = f(centr+absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1322  fv1(jtwm1) = fval1
1323  fv2(jtwm1) = fval2
1324  fsum = fval1+fval2
1325  resk = resk+wgk(jtwm1)*fsum
1326  resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
1327  15 continue
1328  reskh = resk*0.5d+00
1329  resasc = wgk(21)*dabs(fc-reskh)
1330  do 20 j=1,20
1331  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
1332  20 continue
1333  result = resk*hlgth
1334  resabs = resabs*dhlgth
1335  resasc = resasc*dhlgth
1336  abserr = dabs((resk-resg)*hlgth)
1337  if(resasc.ne.0.0d+00.and.abserr.ne.0.d+00)
1338  * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
1339  if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
1340  * ((epmach*0.5d+02)*resabs,abserr)
1341  return
d1mach
double precision function d1mach(I)
Definition: ddeabm.f:2012
e21
static double * e21
Definition: radflowload.c:42

◆ dqk51()

subroutine dqk51 ( external  f,
double precision  a,
double precision  b,
double precision  result,
double precision  abserr,
double precision  resabs,
double precision  resasc,
double precision  phi,
double precision  lambda1,
double precision  zk0,
double precision  Pup,
double precision  Tup,
double precision  rurd,
double precision  xflow,
double precision  kup 
)
1345 c***begin prologue dqk51
1346 c***date written 800101 (yymmdd)
1347 c***revision date 830518 (yymmdd)
1348 c***category no. h2a1a2
1349 c***keywords 51-point gauss-kronrod rules
1350 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
1351 c de doncker,elise,appl. math & progr. div. - k.u.leuven
1352 c***purpose to compute i = integral of f over (a,b) with error
1353 c estimate
1354 c j = integral of abs(f) over (a,b)
1355 c***description
1356 c
1357 c integration rules
1358 c standard fortran subroutine
1359 c double precision version
1360 c
1361 c parameters
1362 c on entry
1363 c f - double precision
1364 c function subroutine defining the integrand
1365 c function f(x). the actual name for f needs to be
1366 c declared e x t e r n a l in the calling program.
1367 c
1368 c a - double precision
1369 c lower limit of integration
1370 c
1371 c b - double precision
1372 c upper limit of integration
1373 c
1374 c on return
1375 c result - double precision
1376 c approximation to the integral i
1377 c result is computed by applying the 51-point
1378 c kronrod rule (resk) obtained by optimal addition
1379 c of abscissae to the 25-point gauss rule (resg).
1380 c
1381 c abserr - double precision
1382 c estimate of the modulus of the absolute error,
1383 c which should not exceed abs(i-result)
1384 c
1385 c resabs - double precision
1386 c approximation to the integral j
1387 c
1388 c resasc - double precision
1389 c approximation to the integral of abs(f-i/(b-a))
1390 c over (a,b)
1391 c
1392 c***references (none)
1393 c***routines called d1mach
1394 c***end prologue dqk51
1395 c
1396  double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
1397  * d1mach(4),epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,
1398  * resasc,resg,resk,reskh,result,uflow,wg,wgk,xgk,phi,lambda1,
1399  * zk0,pup,tup,rurd,xflow,kup
1400  integer j,jtw,jtwm1
1401  external f
1402 c
1403  dimension fv1(25),fv2(25),xgk(26),wgk(26),wg(13)
1404  d1mach(1)=1e21
1405  d1mach(2)=0d0
1406  d1mach(3)=0d0
1407  d1mach(4)=1e-21
1408 c
1409 c the abscissae and weights are given for the interval (-1,1).
1410 c because of symmetry only the positive abscissae and their
1411 c corresponding weights are given.
1412 c
1413 c xgk - abscissae of the 51-point kronrod rule
1414 c xgk(2), xgk(4), ... abscissae of the 25-point
1415 c gauss rule
1416 c xgk(1), xgk(3), ... abscissae which are optimally
1417 c added to the 25-point gauss rule
1418 c
1419 c wgk - weights of the 51-point kronrod rule
1420 c
1421 c wg - weights of the 25-point gauss rule
1422 c
1423 c
1424 c gauss quadrature weights and kronron quadrature abscissae and weights
1425 c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
1426 c bell labs, nov. 1981.
1427 c
1428  data wg( 1) / 0.0113937985 0102628794 7902964113 235 d0 /
1429  data wg( 2) / 0.0263549866 1503213726 1901815295 299 d0 /
1430  data wg( 3) / 0.0409391567 0130631265 5623487711 646 d0 /
1431  data wg( 4) / 0.0549046959 7583519192 5936891540 473 d0 /
1432  data wg( 5) / 0.0680383338 1235691720 7187185656 708 d0 /
1433  data wg( 6) / 0.0801407003 3500101801 3234959669 111 d0 /
1434  data wg( 7) / 0.0910282619 8296364981 1497220702 892 d0 /
1435  data wg( 8) / 0.1005359490 6705064420 2206890392 686 d0 /
1436  data wg( 9) / 0.1085196244 7426365311 6093957050 117 d0 /
1437  data wg( 10) / 0.1148582591 4571164833 9325545869 556 d0 /
1438  data wg( 11) / 0.1194557635 3578477222 8178126512 901 d0 /
1439  data wg( 12) / 0.1222424429 9031004168 8959518945 852 d0 /
1440  data wg( 13) / 0.1231760537 2671545120 3902873079 050 d0 /
1441 c
1442  data xgk( 1) / 0.9992621049 9260983419 3457486540 341 d0 /
1443  data xgk( 2) / 0.9955569697 9049809790 8784946893 902 d0 /
1444  data xgk( 3) / 0.9880357945 3407724763 7331014577 406 d0 /
1445  data xgk( 4) / 0.9766639214 5951751149 8315386479 594 d0 /
1446  data xgk( 5) / 0.9616149864 2584251241 8130033660 167 d0 /
1447  data xgk( 6) / 0.9429745712 2897433941 4011169658 471 d0 /
1448  data xgk( 7) / 0.9207471152 8170156174 6346084546 331 d0 /
1449  data xgk( 8) / 0.8949919978 7827536885 1042006782 805 d0 /
1450  data xgk( 9) / 0.8658470652 9327559544 8996969588 340 d0 /
1451  data xgk( 10) / 0.8334426287 6083400142 1021108693 570 d0 /
1452  data xgk( 11) / 0.7978737979 9850005941 0410904994 307 d0 /
1453  data xgk( 12) / 0.7592592630 3735763057 7282865204 361 d0 /
1454  data xgk( 13) / 0.7177664068 1308438818 6654079773 298 d0 /
1455  data xgk( 14) / 0.6735663684 7346836448 5120633247 622 d0 /
1456  data xgk( 15) / 0.6268100990 1031741278 8122681624 518 d0 /
1457  data xgk( 16) / 0.5776629302 4122296772 3689841612 654 d0 /
1458  data xgk( 17) / 0.5263252843 3471918259 9623778158 010 d0 /
1459  data xgk( 18) / 0.4730027314 4571496052 2182115009 192 d0 /
1460  data xgk( 19) / 0.4178853821 9303774885 1814394594 572 d0 /
1461  data xgk( 20) / 0.3611723058 0938783773 5821730127 641 d0 /
1462  data xgk( 21) / 0.3030895389 3110783016 7478909980 339 d0 /
1463  data xgk( 22) / 0.2438668837 2098843204 5190362797 452 d0 /
1464  data xgk( 23) / 0.1837189394 2104889201 5969888759 528 d0 /
1465  data xgk( 24) / 0.1228646926 1071039638 7359818808 037 d0 /
1466  data xgk( 25) / 0.0615444830 0568507888 6546392366 797 d0 /
1467  data xgk( 26) / 0.0000000000 0000000000 0000000000 000 d0 /
1468 c
1469  data wgk( 1) / 0.0019873838 9233031592 6507851882 843 d0 /
1470  data wgk( 2) / 0.0055619321 3535671375 8040236901 066 d0 /
1471  data wgk( 3) / 0.0094739733 8617415160 7207710523 655 d0 /
1472  data wgk( 4) / 0.0132362291 9557167481 3656405846 976 d0 /
1473  data wgk( 5) / 0.0168478177 0912829823 1516667536 336 d0 /
1474  data wgk( 6) / 0.0204353711 4588283545 6568292235 939 d0 /
1475  data wgk( 7) / 0.0240099456 0695321622 0092489164 881 d0 /
1476  data wgk( 8) / 0.0274753175 8785173780 2948455517 811 d0 /
1477  data wgk( 9) / 0.0307923001 6738748889 1109020215 229 d0 /
1478  data wgk( 10) / 0.0340021302 7432933783 6748795229 551 d0 /
1479  data wgk( 11) / 0.0371162714 8341554356 0330625367 620 d0 /
1480  data wgk( 12) / 0.0400838255 0403238207 4839284467 076 d0 /
1481  data wgk( 13) / 0.0428728450 2017004947 6895792439 495 d0 /
1482  data wgk( 14) / 0.0455029130 4992178890 9870584752 660 d0 /
1483  data wgk( 15) / 0.0479825371 3883671390 6392255756 915 d0 /
1484  data wgk( 16) / 0.0502776790 8071567196 3325259433 440 d0 /
1485  data wgk( 17) / 0.0523628858 0640747586 4366712137 873 d0 /
1486  data wgk( 18) / 0.0542511298 8854549014 4543370459 876 d0 /
1487  data wgk( 19) / 0.0559508112 2041231730 8240686382 747 d0 /
1488  data wgk( 20) / 0.0574371163 6156783285 3582693939 506 d0 /
1489  data wgk( 21) / 0.0586896800 2239420796 1974175856 788 d0 /
1490  data wgk( 22) / 0.0597203403 2417405997 9099291932 562 d0 /
1491  data wgk( 23) / 0.0605394553 7604586294 5360267517 565 d0 /
1492  data wgk( 24) / 0.0611285097 1705304830 5859030416 293 d0 /
1493  data wgk( 25) / 0.0614711898 7142531666 1544131965 264 d0 /
1494 c note: wgk (26) was calculated from the values of wgk(1..25)
1495  data wgk( 26) / 0.0615808180 6783293507 8759824240 066 d0 /
1496 c
1497 c
1498 c list of major variables
1499 c -----------------------
1500 c
1501 c centr - mid point of the interval
1502 c hlgth - half-length of the interval
1503 c absc - abscissa
1504 c fval* - function value
1505 c resg - result of the 25-point gauss formula
1506 c resk - result of the 51-point kronrod formula
1507 c reskh - approximation to the mean value of f over (a,b),
1508 c i.e. to i/(b-a)
1509 c
1510 c machine dependent constants
1511 c ---------------------------
1512 c
1513 c epmach is the largest relative spacing.
1514 c uflow is the smallest positive magnitude.
1515 c
1516 c***first executable statement dqk51
1517  epmach = d1mach(4)
1518  uflow = d1mach(1)
1519 c
1520  centr = 0.5d+00*(a+b)
1521  hlgth = 0.5d+00*(b-a)
1522  dhlgth = dabs(hlgth)
1523 c
1524 c compute the 51-point kronrod approximation to
1525 c the integral, and estimate the absolute error.
1526 c
1527  fc = f(centr,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1528  resg = wg(13)*fc
1529  resk = wgk(26)*fc
1530  resabs = dabs(resk)
1531  do 10 j=1,12
1532  jtw = j*2
1533  absc = hlgth*xgk(jtw)
1534  fval1 = f(centr-absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1535  fval2 = f(centr+absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1536  fv1(jtw) = fval1
1537  fv2(jtw) = fval2
1538  fsum = fval1+fval2
1539  resg = resg+wg(j)*fsum
1540  resk = resk+wgk(jtw)*fsum
1541  resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
1542  10 continue
1543  do 15 j = 1,13
1544  jtwm1 = j*2-1
1545  absc = hlgth*xgk(jtwm1)
1546  fval1 = f(centr-absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1547  fval2 = f(centr+absc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1548  fv1(jtwm1) = fval1
1549  fv2(jtwm1) = fval2
1550  fsum = fval1+fval2
1551  resk = resk+wgk(jtwm1)*fsum
1552  resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
1553  15 continue
1554  reskh = resk*0.5d+00
1555  resasc = wgk(26)*dabs(fc-reskh)
1556  do 20 j=1,25
1557  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
1558  20 continue
1559  result = resk*hlgth
1560  resabs = resabs*dhlgth
1561  resasc = resasc*dhlgth
1562  abserr = dabs((resk-resg)*hlgth)
1563  if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
1564  * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
1565  if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
1566  * ((epmach*0.5d+02)*resabs,abserr)
1567  return
d1mach
double precision function d1mach(I)
Definition: ddeabm.f:2012
e21
static double * e21
Definition: radflowload.c:42

◆ dqk61()

subroutine dqk61 ( external  f,
double precision  a,
double precision  b,
double precision  result,
double precision  abserr,
double precision  resabs,
double precision  resasc,
double precision  phi,
double precision  lambda1,
double precision  zk0,
double precision  Pup,
double precision  Tup,
double precision  rurd,
double precision  xflow,
double precision  kup 
)
1571 c***begin prologue dqk61
1572 c***date written 800101 (yymmdd)
1573 c***revision date 830518 (yymmdd)
1574 c***category no. h2a1a2
1575 c***keywords 61-point gauss-kronrod rules
1576 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
1577 c de doncker,elise,appl. math. & progr. div. - k.u.leuven
1578 c***purpose to compute i = integral of f over (a,b) with error
1579 c estimate
1580 c j = integral of dabs(f) over (a,b)
1581 c***description
1582 c
1583 c integration rule
1584 c standard fortran subroutine
1585 c double precision version
1586 c
1587 c
1588 c parameters
1589 c on entry
1590 c f - double precision
1591 c function subprogram defining the integrand
1592 c function f(x). the actual name for f needs to be
1593 c declared e x t e r n a l in the calling program.
1594 c
1595 c a - double precision
1596 c lower limit of integration
1597 c
1598 c b - double precision
1599 c upper limit of integration
1600 c
1601 c on return
1602 c result - double precision
1603 c approximation to the integral i
1604 c result is computed by applying the 61-point
1605 c kronrod rule (resk) obtained by optimal addition of
1606 c abscissae to the 30-point gauss rule (resg).
1607 c
1608 c abserr - double precision
1609 c estimate of the modulus of the absolute error,
1610 c which should equal or exceed dabs(i-result)
1611 c
1612 c resabs - double precision
1613 c approximation to the integral j
1614 c
1615 c resasc - double precision
1616 c approximation to the integral of dabs(f-i/(b-a))
1617 c
1618 c
1619 c***references (none)
1620 c***routines called d1mach
1621 c***end prologue dqk61
1622 c
1623  double precision a,dabsc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
1624  * d1mach(4),epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,
1625  * resasc,resg,resk,reskh,result,uflow,wg,wgk,xgk,phi,lambda1,
1626  * zk0,pup,tup,rurd,xflow,kup
1627  integer j,jtw,jtwm1
1628  external f
1629 c
1630  dimension fv1(30),fv2(30),xgk(31),wgk(31),wg(15)
1631  d1mach(1)=1e21
1632  d1mach(2)=0
1633  d1mach(3)=0
1634  d1mach(4)=1e-21
1635 c
1636 c the abscissae and weights are given for the
1637 c interval (-1,1). because of symmetry only the positive
1638 c abscissae and their corresponding weights are given.
1639 c
1640 c xgk - abscissae of the 61-point kronrod rule
1641 c xgk(2), xgk(4) ... abscissae of the 30-point
1642 c gauss rule
1643 c xgk(1), xgk(3) ... optimally added abscissae
1644 c to the 30-point gauss rule
1645 c
1646 c wgk - weights of the 61-point kronrod rule
1647 c
1648 c wg - weigths of the 30-point gauss rule
1649 c
1650 c
1651 c gauss quadrature weights and kronron quadrature abscissae and weights
1652 c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
1653 c bell labs, nov. 1981.
1654 c
1655  data wg( 1) / 0.0079681924 9616660561 5465883474 674 d0 /
1656  data wg( 2) / 0.0184664683 1109095914 2302131912 047 d0 /
1657  data wg( 3) / 0.0287847078 8332336934 9719179611 292 d0 /
1658  data wg( 4) / 0.0387991925 6962704959 6801936446 348 d0 /
1659  data wg( 5) / 0.0484026728 3059405290 2938140422 808 d0 /
1660  data wg( 6) / 0.0574931562 1761906648 1721689402 056 d0 /
1661  data wg( 7) / 0.0659742298 8218049512 8128515115 962 d0 /
1662  data wg( 8) / 0.0737559747 3770520626 8243850022 191 d0 /
1663  data wg( 9) / 0.0807558952 2942021535 4694938460 530 d0 /
1664  data wg( 10) / 0.0868997872 0108297980 2387530715 126 d0 /
1665  data wg( 11) / 0.0921225222 3778612871 7632707087 619 d0 /
1666  data wg( 12) / 0.0963687371 7464425963 9468626351 810 d0 /
1667  data wg( 13) / 0.0995934205 8679526706 2780282103 569 d0 /
1668  data wg( 14) / 0.1017623897 4840550459 6428952168 554 d0 /
1669  data wg( 15) / 0.1028526528 9355884034 1285636705 415 d0 /
1670 c
1671  data xgk( 1) / 0.9994844100 5049063757 1325895705 811 d0 /
1672  data xgk( 2) / 0.9968934840 7464954027 1630050918 695 d0 /
1673  data xgk( 3) / 0.9916309968 7040459485 8628366109 486 d0 /
1674  data xgk( 4) / 0.9836681232 7974720997 0032581605 663 d0 /
1675  data xgk( 5) / 0.9731163225 0112626837 4693868423 707 d0 /
1676  data xgk( 6) / 0.9600218649 6830751221 6871025581 798 d0 /
1677  data xgk( 7) / 0.9443744447 4855997941 5831324037 439 d0 /
1678  data xgk( 8) / 0.9262000474 2927432587 9324277080 474 d0 /
1679  data xgk( 9) / 0.9055733076 9990779854 6522558925 958 d0 /
1680  data xgk( 10) / 0.8825605357 9205268154 3116462530 226 d0 /
1681  data xgk( 11) / 0.8572052335 4606109895 8658510658 944 d0 /
1682  data xgk( 12) / 0.8295657623 8276839744 2898119732 502 d0 /
1683  data xgk( 13) / 0.7997278358 2183908301 3668942322 683 d0 /
1684  data xgk( 14) / 0.7677774321 0482619491 7977340974 503 d0 /
1685  data xgk( 15) / 0.7337900624 5322680472 6171131369 528 d0 /
1686  data xgk( 16) / 0.6978504947 9331579693 2292388026 640 d0 /
1687  data xgk( 17) / 0.6600610641 2662696137 0053668149 271 d0 /
1688  data xgk( 18) / 0.6205261829 8924286114 0477556431 189 d0 /
1689  data xgk( 19) / 0.5793452358 2636169175 6024932172 540 d0 /
1690  data xgk( 20) / 0.5366241481 4201989926 4169793311 073 d0 /
1691  data xgk( 21) / 0.4924804678 6177857499 3693061207 709 d0 /
1692  data xgk( 22) / 0.4470337695 3808917678 0609900322 854 d0 /
1693  data xgk( 23) / 0.4004012548 3039439253 5476211542 661 d0 /
1694  data xgk( 24) / 0.3527047255 3087811347 1037207089 374 d0 /
1695  data xgk( 25) / 0.3040732022 7362507737 2677107199 257 d0 /
1696  data xgk( 26) / 0.2546369261 6788984643 9805129817 805 d0 /
1697  data xgk( 27) / 0.2045251166 8230989143 8957671002 025 d0 /
1698  data xgk( 28) / 0.1538699136 0858354696 3794672743 256 d0 /
1699  data xgk( 29) / 0.1028069379 6673703014 7096751318 001 d0 /
1700  data xgk( 30) / 0.0514718425 5531769583 3025213166 723 d0 /
1701  data xgk( 31) / 0.0000000000 0000000000 0000000000 000 d0 /
1702 c
1703  data wgk( 1) / 0.0013890136 9867700762 4551591226 760 d0 /
1704  data wgk( 2) / 0.0038904611 2709988405 1267201844 516 d0 /
1705  data wgk( 3) / 0.0066307039 1593129217 3319826369 750 d0 /
1706  data wgk( 4) / 0.0092732796 5951776342 8441146892 024 d0 /
1707  data wgk( 5) / 0.0118230152 5349634174 2232898853 251 d0 /
1708  data wgk( 6) / 0.0143697295 0704580481 2451432443 580 d0 /
1709  data wgk( 7) / 0.0169208891 8905327262 7572289420 322 d0 /
1710  data wgk( 8) / 0.0194141411 9394238117 3408951050 128 d0 /
1711  data wgk( 9) / 0.0218280358 2160919229 7167485738 339 d0 /
1712  data wgk( 10) / 0.0241911620 7808060136 5686370725 232 d0 /
1713  data wgk( 11) / 0.0265099548 8233310161 0601709335 075 d0 /
1714  data wgk( 12) / 0.0287540487 6504129284 3978785354 334 d0 /
1715  data wgk( 13) / 0.0309072575 6238776247 2884252943 092 d0 /
1716  data wgk( 14) / 0.0329814470 5748372603 1814191016 854 d0 /
1717  data wgk( 15) / 0.0349793380 2806002413 7499670731 468 d0 /
1718  data wgk( 16) / 0.0368823646 5182122922 3911065617 136 d0 /
1719  data wgk( 17) / 0.0386789456 2472759295 0348651532 281 d0 /
1720  data wgk( 18) / 0.0403745389 5153595911 1995279752 468 d0 /
1721  data wgk( 19) / 0.0419698102 1516424614 7147541285 970 d0 /
1722  data wgk( 20) / 0.0434525397 0135606931 6831728117 073 d0 /
1723  data wgk( 21) / 0.0448148001 3316266319 2355551616 723 d0 /
1724  data wgk( 22) / 0.0460592382 7100698811 6271735559 374 d0 /
1725  data wgk( 23) / 0.0471855465 6929915394 5261478181 099 d0 /
1726  data wgk( 24) / 0.0481858617 5708712914 0779492298 305 d0 /
1727  data wgk( 25) / 0.0490554345 5502977888 7528165367 238 d0 /
1728  data wgk( 26) / 0.0497956834 2707420635 7811569379 942 d0 /
1729  data wgk( 27) / 0.0504059214 0278234684 0893085653 585 d0 /
1730  data wgk( 28) / 0.0508817958 9874960649 2297473049 805 d0 /
1731  data wgk( 29) / 0.0512215478 4925877217 0656282604 944 d0 /
1732  data wgk( 30) / 0.0514261285 3745902593 3862879215 781 d0 /
1733  data wgk( 31) / 0.0514947294 2945156755 8340433647 099 d0 /
1734 c
1735 c list of major variables
1736 c -----------------------
1737 c
1738 c centr - mid point of the interval
1739 c hlgth - half-length of the interval
1740 c dabsc - abscissa
1741 c fval* - function value
1742 c resg - result of the 30-point gauss rule
1743 c resk - result of the 61-point kronrod rule
1744 c reskh - approximation to the mean value of f
1745 c over (a,b), i.e. to i/(b-a)
1746 c
1747 c machine dependent constants
1748 c ---------------------------
1749 c
1750 c epmach is the largest relative spacing.
1751 c uflow is the smallest positive magnitude.
1752 c
1753  epmach = d1mach(4)
1754  uflow = d1mach(1)
1755 c
1756  centr = 0.5d+00*(b+a)
1757  hlgth = 0.5d+00*(b-a)
1758  dhlgth = dabs(hlgth)
1759 c
1760 c compute the 61-point kronrod approximation to the
1761 c integral, and estimate the absolute error.
1762 c
1763 c***first executable statement dqk61
1764  resg = 0.0d+00
1765  fc = f(centr,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1766  resk = wgk(31)*fc
1767  resabs = dabs(resk)
1768  do 10 j=1,15
1769  jtw = j*2
1770  dabsc = hlgth*xgk(jtw)
1771  fval1 = f(centr-dabsc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1772  fval2 = f(centr+dabsc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1773  fv1(jtw) = fval1
1774  fv2(jtw) = fval2
1775  fsum = fval1+fval2
1776  resg = resg+wg(j)*fsum
1777  resk = resk+wgk(jtw)*fsum
1778  resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
1779  10 continue
1780  do 15 j=1,15
1781  jtwm1 = j*2-1
1782  dabsc = hlgth*xgk(jtwm1)
1783  fval1 = f(centr-dabsc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1784  fval2 = f(centr+dabsc,phi,lambda1,zk0,pup,tup,rurd,xflow,kup)
1785  fv1(jtwm1) = fval1
1786  fv2(jtwm1) = fval2
1787  fsum = fval1+fval2
1788  resk = resk+wgk(jtwm1)*fsum
1789  resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
1790  15 continue
1791  reskh = resk*0.5d+00
1792  resasc = wgk(31)*dabs(fc-reskh)
1793  do 20 j=1,30
1794  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
1795  20 continue
1796  result = resk*hlgth
1797  resabs = resabs*dhlgth
1798  resasc = resasc*dhlgth
1799  abserr = dabs((resk-resg)*hlgth)
1800  if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
1801  * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
1802  if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
1803  * ((epmach*0.5d+02)*resabs,abserr)
1804  return
d1mach
double precision function d1mach(I)
Definition: ddeabm.f:2012
e21
static double * e21
Definition: radflowload.c:42

◆ dqpsrt()

subroutine dqpsrt ( integer  limit,
integer  last,
integer  maxerr,
double precision  ermax,
double precision, dimension(last)  elist,
integer, dimension(last)  iord,
integer  nrmax,
double precision  phi,
double precision  lambda1,
double precision  zk0,
double precision  Pup,
double precision  Tup,
double precision  rurd,
double precision  xflow,
double precision  kup 
)
1808 c***begin prologue dqpsrt
1809 c***refer to dqage,dqagie,dqagpe,dqawse
1810 c***routines called (none)
1811 c***revision date 810101 (yymmdd)
1812 c***keywords sequential sorting
1813 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
1814 c de doncker,elise,appl. math. & progr. div. - k.u.leuven
1815 c***purpose this routine maintains the descending ordering in the
1816 c list of the local error estimated resulting from the
1817 c interval subdivision process. at each call two error
1818 c estimates are inserted using the sequential search
1819 c method, top-down for the largest error estimate and
1820 c bottom-up for the smallest error estimate.
1821 c***description
1822 c
1823 c ordering routine
1824 c standard fortran subroutine
1825 c double precision version
1826 c
1827 c parameters (meaning at output)
1828 c limit - integer
1829 c maximum number of error estimates the list
1830 c can contain
1831 c
1832 c last - integer
1833 c number of error estimates currently in the list
1834 c
1835 c maxerr - integer
1836 c maxerr points to the nrmax-th largest error
1837 c estimate currently in the list
1838 c
1839 c ermax - double precision
1840 c nrmax-th largest error estimate
1841 c ermax = elist(maxerr)
1842 c
1843 c elist - double precision
1844 c vector of dimension last containing
1845 c the error estimates
1846 c
1847 c iord - integer
1848 c vector of dimension last, the first k elements
1849 c of which contain pointers to the error
1850 c estimates, such that
1851 c elist(iord(1)),..., elist(iord(k))
1852 c form a decreasing sequence, with
1853 c k = last if last.le.(limit/2+2), and
1854 c k = limit+1-last otherwise
1855 c
1856 c nrmax - integer
1857 c maxerr = iord(nrmax)
1858 c
1859 c***end prologue dqpsrt
1860 c
1861  double precision elist,ermax,errmax,errmin,phi,lambda1,zk0,
1862  * pup,tup,rurd,xflow,kup
1863  integer i,ibeg,ido,iord,isucc,j,jbnd,jupbn,k,last,limit,maxerr,
1864  * nrmax
1865  dimension elist(last),iord(last)
1866 c
1867 c check whether the list contains more than
1868 c two error estimates.
1869 c
1870 c***first executable statement dqpsrt
1871  if(last.gt.2) go to 10
1872  iord(1) = 1
1873  iord(2) = 2
1874  go to 90
1875 c
1876 c this part of the routine is only executed if, due to a
1877 c difficult integrand, subdivision increased the error
1878 c estimate. in the normal case the insert procedure should
1879 c start after the nrmax-th largest error estimate.
1880 c
1881  10 errmax = elist(maxerr)
1882  if(nrmax.eq.1) go to 30
1883  ido = nrmax-1
1884  do 20 i = 1,ido
1885  isucc = iord(nrmax-1)
1886 c ***jump out of do-loop
1887  if(errmax.le.elist(isucc)) go to 30
1888  iord(nrmax) = isucc
1889  nrmax = nrmax-1
1890  20 continue
1891 c
1892 c compute the number of elements in the list to be maintained
1893 c in descending order. this number depends on the number of
1894 c subdivisions still allowed.
1895 c
1896  30 jupbn = last
1897  if(last.gt.(limit/2+2)) jupbn = limit+3-last
1898  errmin = elist(last)
1899 c
1900 c insert errmax by traversing the list top-down,
1901 c starting comparison from the element elist(iord(nrmax+1)).
1902 c
1903  jbnd = jupbn-1
1904  ibeg = nrmax+1
1905  if(ibeg.gt.jbnd) go to 50
1906  do 40 i=ibeg,jbnd
1907  isucc = iord(i)
1908 c ***jump out of do-loop
1909  if(errmax.ge.elist(isucc)) go to 60
1910  iord(i-1) = isucc
1911  40 continue
1912  50 iord(jbnd) = maxerr
1913  iord(jupbn) = last
1914  go to 90
1915 c
1916 c insert errmin by traversing the list bottom-up.
1917 c
1918  60 iord(i-1) = maxerr
1919  k = jbnd
1920  do 70 j=i,jbnd
1921  isucc = iord(k)
1922 c ***jump out of do-loop
1923  if(errmin.lt.elist(isucc)) go to 80
1924  iord(k+1) = isucc
1925  k = k-1
1926  70 continue
1927  iord(i) = last
1928  go to 90
1929  80 iord(k+1) = last
1930 c
1931 c set maxerr and ermax.
1932 c
1933  90 maxerr = iord(nrmax)
1934  ermax = elist(maxerr)
1935  return
Hosted by OpenAircraft.com, (Michigan UAV, LLC)