1c 2c CMLIB - public domain 3c 4 subroutine dqag(f,a,b,epsabs,epsrel,key,result,abserr,neval,ier, 5 * limit,lenw,last,iwork,work,phi,lambda1,zk0,Pup,Tup,rurd,xflow, 6 * kup) 7c***begin prologue dqag 8c***date written 800101 (yymmdd) 9c***revision date 830518 (yymmdd) 10c***category no. h2a1a1 11c***keywords automatic integrator, general-purpose, 12c integrand examinator, globally adaptive, 13c gauss-kronrod 14c***author piessens,robert,appl. math. & progr. div - k.u.leuven 15c de doncker,elise,appl. math. & progr. div. - k.u.leuven 16c***purpose the routine calculates an approximation result to a given 17c definite integral i = integral of f over (a,b), 18c hopefully satisfying following claim for accuracy 19c abs(i-result)le.max(epsabs,epsrel*abs(i)). 20c***description 21c 22c computation of a definite integral 23c standard fortran subroutine 24c double precision version 25c 26c f - double precision 27c function subprogam defining the integrand 28c function f(x). the actual name for f needs to be 29c declared e x t e r n a l in the driver program. 30c 31c a - double precision 32c lower limit of integration 33c 34c b - double precision 35c upper limit of integration 36c 37c epsabs - double precision 38c absolute accoracy requested 39c epsrel - double precision 40c relative accuracy requested 41c if epsabs.le.0 42c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), 43c the routine will end with ier = 6. 44c 45c key - integer 46c key for choice of local integration rule 47c a gauss-kronrod pair is used with 48c 7 - 15 points if key.lt.2, 49c 10 - 21 points if key = 2, 50c 15 - 31 points if key = 3, 51c 20 - 41 points if key = 4, 52c 25 - 51 points if key = 5, 53c 30 - 61 points if key.gt.5. 54c 55c on return 56c result - double precision 57c approximation to the integral 58c 59c abserr - double precision 60c estimate of the modulus of the absolute error, 61c which should equal or exceed abs(i-result) 62c 63c neval - integer 64c number of integrand evaluations 65c 66c ier - integer 67c ier = 0 normal and reliable termination of the 68c routine. it is assumed that the requested 69c accuracy has been achieved. 70c ier.gt.0 abnormal termination of the routine 71c the estimates for result and error are 72c less reliable. it is assumed that the 73c requested accuracy has not been achieved. 74c error messages 75c ier = 1 maximum number of subdivisions allowed 76c has been achieved. one can allow more 77c subdivisions by increasing the value of 78c limit (and taking the according dimension 79c adjustments into account). however, if 80c this yield no improvement it is advised 81c to analyze the integrand in order to 82c determine the integration difficulaties. 83c if the position of a local difficulty can 84c be determined (i.e.singularity, 85c discontinuity within the interval) one 86c will probably gain from splitting up the 87c interval at this point and calling the 88c integrator on the subranges. if possible, 89c an appropriate special-purpose integrator 90c should be used which is designed for 91c handling the type of difficulty involved. 92c = 2 the occurrence of roundoff error is 93c detected, which prevents the requested 94c tolerance from being achieved. 95c = 3 extremely bad integrand behaviour occurs 96c at some points of the integration 97c interval. 98c = 6 the input is invalid, because 99c (epsabs.le.0 and 100c epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) 101c or limit.lt.1 or lenw.lt.limit*4. 102c result, abserr, neval, last are set 103c to zero. 104c except when lenw is invalid, iwork(1), 105c work(limit*2+1) and work(limit*3+1) are 106c set to zero, work(1) is set to a and 107c work(limit+1) to b. 108c 109c dimensioning parameters 110c limit - integer 111c dimensioning parameter for iwork 112c limit determines the maximum number of subintervals 113c in the partition of the given integration interval 114c (a,b), limit.ge.1. 115c if limit.lt.1, the routine will end with ier = 6. 116c 117c lenw - integer 118c dimensioning parameter for work 119c lenw must be at least limit*4. 120c if lenw.lt.limit*4, the routine will end with 121c ier = 6. 122c 123c last - integer 124c on return, last equals the number of subintervals 125c produced in the subdiviosion process, which 126c determines the number of significant elements 127c actually in the work arrays. 128c 129c work arrays 130c iwork - integer 131c vector of dimension at least limit, the first k 132c elements of which contain pointers to the error 133c estimates over the subintervals, such that 134c work(limit*3+iwork(1)),... , work(limit*3+iwork(k)) 135c form a decreasing sequence with k = last if 136c last.le.(limit/2+2), and k = limit+1-last otherwise 137c 138c work - double precision 139c vector of dimension at least lenw 140c on return 141c work(1), ..., work(last) contain the left end 142c points of the subintervals in the partition of 143c (a,b), 144c work(limit+1), ..., work(limit+last) contain the 145c right end points, 146c work(limit*2+1), ..., work(limit*2+last) contain 147c the integral approximations over the subintervals, 148c work(limit*3+1), ..., work(limit*3+last) contain 149c the error estimates. 150c 151c***references (none) 152c***routines called dqage,xerror 153c***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 157c 158 dimension iwork(limit),work(lenw) 159c 160 external f 161c 162c check validity of lenw. 163c 164 d1mach(1)=1E21 165 d1mach(2)=0d0 166 d1mach(3)=0d0 167 d1mach(4)=1E-21 168c 169c***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 176c 177c prepare call for dqage. 178c 179 l1 = limit+1 180 l2 = limit+l1 181 l3 = limit+l2 182c 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) 186c 187c call error handler if necessary. 188c 189 lvl = 0 19010 if(ier.eq.6) lvl = 1 191! if(ier.ne.0) call xerror(26habnormal return from dqag ,26,ier,lvl) 192 return 193 end 194 subroutine dqage(f,a,b,epsabs,epsrel,key,limit,result,abserr, 195 * neval,ier,alist,blist,rlist,elist,iord,last,phi,lambda1,zk0, 196 * Pup,Tup,rurd,xflow,kup) 197c***begin prologue dqage 198c***date written 800101 (yymmdd) 199c***revision date 830518 (yymmdd) 200c***category no. h2a1a1 201c***keywords automatic integrator, general-purpose, 202c integrand examinator, globally adaptive, 203c gauss-kronrod 204c***author piessens,robert,appl. math. & progr. div. - k.u.leuven 205c de doncker,elise,appl. math. & progr. div. - k.u.leuven 206c***purpose the routine calculates an approximation result to a given 207c definite integral i = integral of f over (a,b), 208c hopefully satisfying following claim for accuracy 209c abs(i-reslt).le.max(epsabs,epsrel*abs(i)). 210c***description 211c 212c computation of a definite integral 213c standard fortran subroutine 214c double precision version 215c 216c parameters 217c on entry 218c f - double precision 219c function subprogram defining the integrand 220c function f(x). the actual name for f needs to be 221c declared e x t e r n a l in the driver program. 222c 223c a - double precision 224c lower limit of integration 225c 226c b - double precision 227c upper limit of integration 228c 229c epsabs - double precision 230c absolute accuracy requested 231c epsrel - double precision 232c relative accuracy requested 233c if epsabs.le.0 234c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), 235c the routine will end with ier = 6. 236c 237c key - integer 238c key for choice of local integration rule 239c a gauss-kronrod pair is used with 240c 7 - 15 points if key.lt.2, 241c 10 - 21 points if key = 2, 242c 15 - 31 points if key = 3, 243c 20 - 41 points if key = 4, 244c 25 - 51 points if key = 5, 245c 30 - 61 points if key.gt.5. 246c 247c limit - integer 248c gives an upperbound on the number of subintervals 249c in the partition of (a,b), limit.ge.1. 250c 251c on return 252c result - double precision 253c approximation to the integral 254c 255c abserr - double precision 256c estimate of the modulus of the absolute error, 257c which should equal or exceed abs(i-result) 258c 259c neval - integer 260c number of integrand evaluations 261c 262c ier - integer 263c ier = 0 normal and reliable termination of the 264c routine. it is assumed that the requested 265c accuracy has been achieved. 266c ier.gt.0 abnormal termination of the routine 267c the estimates for result and error are 268c less reliable. it is assumed that the 269c requested accuracy has not been achieved. 270c error messages 271c ier = 1 maximum number of subdivisions allowed 272c has been achieved. one can allow more 273c subdivisions by increasing the value 274c of limit. 275c however, if this yields no improvement it 276c is rather advised to analyze the integrand 277c in order to determine the integration 278c difficulties. if the position of a local 279c difficulty can be determined(e.g. 280c singularity, discontinuity within the 281c interval) one will probably gain from 282c splitting up the interval at this point 283c and calling the integrator on the 284c subranges. if possible, an appropriate 285c special-purpose integrator should be used 286c which is designed for handling the type of 287c difficulty involved. 288c = 2 the occurrence of roundoff error is 289c detected, which prevents the requested 290c tolerance from being achieved. 291c = 3 extremely bad integrand behaviour occurs 292c at some points of the integration 293c interval. 294c = 6 the input is invalid, because 295c (epsabs.le.0 and 296c epsrel.lt.max(50*rel.mach.acc.,0.5d-28), 297c result, abserr, neval, last, rlist(1) , 298c elist(1) and iord(1) are set to zero. 299c alist(1) and blist(1) are set to a and b 300c respectively. 301c 302c alist - double precision 303c vector of dimension at least limit, the first 304c last elements of which are the left 305c end points of the subintervals in the partition 306c of the given integration range (a,b) 307c 308c blist - double precision 309c vector of dimension at least limit, the first 310c last elements of which are the right 311c end points of the subintervals in the partition 312c of the given integration range (a,b) 313c 314c rlist - double precision 315c vector of dimension at least limit, the first 316c last elements of which are the 317c integral approximations on the subintervals 318c 319c elist - double precision 320c vector of dimension at least limit, the first 321c last elements of which are the moduli of the 322c absolute error estimates on the subintervals 323c 324c iord - integer 325c vector of dimension at least limit, the first k 326c elements of which are pointers to the 327c error estimates over the subintervals, 328c such that elist(iord(1)), ..., 329c elist(iord(k)) form a decreasing sequence, 330c with k = last if last.le.(limit/2+2), and 331c k = limit+1-last otherwise 332c 333c last - integer 334c number of subintervals actually produced in the 335c subdivision process 336c 337c***references (none) 338c***routines called d1mach,dqk15,dqk21,dqk31, 339c dqk41,dqk51,dqk61,dqpsrt 340c***end prologue dqage 341c 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 348c 349 dimension alist(limit),blist(limit),elist(limit),iord(limit), 350 * rlist(limit) 351c 352 external f 353 354 355 d1mach(1)=1E21 356 d1mach(2)=0d0 357 d1mach(3)=0d0 358 d1mach(4)=1E-21 359c 360c list of major variables 361c ----------------------- 362c 363c alist - list of left end points of all subintervals 364c considered up to now 365c blist - list of right end points of all subintervals 366c considered up to now 367c rlist(i) - approximation to the integral over 368c (alist(i),blist(i)) 369c elist(i) - error estimate applying to rlist(i) 370c maxerr - pointer to the interval with largest 371c error estimate 372c errmax - elist(maxerr) 373c area - sum of the integrals over the subintervals 374c errsum - sum of the errors over the subintervals 375c errbnd - requested accuracy max(epsabs,epsrel* 376c abs(result)) 377c *****1 - variable for the left subinterval 378c *****2 - variable for the right subinterval 379c last - index for subdivision 380c 381c 382c machine dependent constants 383c --------------------------- 384c 385c epmach is the largest relative spacing. 386c uflow is the smallest positive magnitude. 387c 388c***first executable statement dqage 389 epmach = d1mach(4) 390 uflow = d1mach(1) 391c 392c test on validity of parameters 393c ------------------------------ 394c 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 408c 409c first approximation to the integral 410c ----------------------------------- 411c 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 432c 433c test on accuracy. 434c 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 440c 441c initialization 442c -------------- 443c 444c 445 errmax = abserr 446 maxerr = 1 447 area = result 448 errsum = abserr 449 nrmax = 1 450 iroff1 = 0 451 iroff2 = 0 452c 453c main do-loop 454c ------------ 455c 456 do 30 last = 2,limit 457c 458c bisect the subinterval with the largest error estimate. 459c 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 ) 488c 489c improve previous approximations to integral 490c and error and test for accuracy. 491c 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 505c 506c test for roundoff error and eventually set error flag. 507c 508 if(iroff1.ge.6.or.iroff2.ge.20) ier = 2 509c 510c set error flag in the case that the number of subintervals 511c equals limit. 512c 513 if(last.eq.limit) ier = 1 514c 515c set error flag in the case of bad integrand behaviour 516c at a point of the integration range. 517c 518 if(max(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03* 519 * epmach)*(dabs(a2)+0.1d+04*uflow)) ier = 3 520c 521c append the newly-created intervals to the list. 522c 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 537c 538c call subroutine dqpsrt to maintain the descending ordering 539c in the list of error estimates and select the subinterval 540c with the largest error estimate (to be bisected next). 541c 542 20 call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax,phi, 543 * lambda1,zk0,Pup,Tup,rurd,xflow,kup) 544c ***jump out of do-loop 545 if(ier.ne.0.or.errsum.le.errbnd) go to 40 546 30 continue 547c 548c compute final result. 549c --------------------- 550c 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 559 end 560 subroutine dqk15(f,a,b,result,abserr,resabs,resasc,phi,lambda1, 561 * zk0,Pup,Tup,rurd,xflow,kup) 562c***begin prologue dqk15 563c***date written 800101 (yymmdd) 564c***revision date 830518 (yymmdd) 565c***category no. h2a1a2 566c***keywords 15-point gauss-kronrod rules 567c***author piessens,robert,appl. math. & progr. div. - k.u.leuven 568c de doncker,elise,appl. math. & progr. div - k.u.leuven 569c***purpose to compute i = integral of f over (a,b), with error 570c estimate 571c j = integral of abs(f) over (a,b) 572c***description 573c 574c integration rules 575c standard fortran subroutine 576c double precision version 577c 578c parameters 579c on entry 580c f - double precision 581c function subprogram defining the integrand 582c function f(x). the actual name for f needs to be 583c declared e x t e r n a l in the calling program. 584c 585c a - double precision 586c lower limit of integration 587c 588c b - double precision 589c upper limit of integration 590c 591c on return 592c result - double precision 593c approximation to the integral i 594c result is computed by applying the 15-point 595c kronrod rule (resk) obtained by optimal addition 596c of abscissae to the7-point gauss rule(resg). 597c 598c abserr - double precision 599c estimate of the modulus of the absolute error, 600c which should not exceed abs(i-result) 601c 602c resabs - double precision 603c approximation to the integral j 604c 605c resasc - double precision 606c approximation to the integral of abs(f-i/(b-a)) 607c over (a,b) 608c 609c***references (none) 610c***routines called d1mach 611c***end prologue dqk15 612c 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 619c 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 628c 629c the abscissae and weights are given for the interval (-1,1). 630c because of symmetry only the positive abscissae and their 631c corresponding weights are given. 632c 633c xgk - abscissae of the 15-point kronrod rule 634c xgk(2), xgk(4), ... abscissae of the 7-point 635c gauss rule 636c xgk(1), xgk(3), ... abscissae which are optimally 637c added to the 7-point gauss rule 638c 639c wgk - weights of the 15-point kronrod rule 640c 641c wg - weights of the 7-point gauss rule 642c 643c 644c gauss quadrature weights and kronron quadrature abscissae and weights 645c as evaluated with 80 decimal digit arithmetic by l. w. fullerton, 646c bell labs, nov. 1981. 647c 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 / 652c 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 / 661c 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 / 670c 671c 672c list of major variables 673c ----------------------- 674c 675c centr - mid point of the interval 676c hlgth - half-length of the interval 677c absc - abscissa 678c fval* - function value 679c resg - result of the 7-point gauss formula 680c resk - result of the 15-point kronrod formula 681c reskh - approximation to the mean value of f over (a,b), 682c i.e. to i/(b-a) 683c 684c machine dependent constants 685c --------------------------- 686c 687c epmach is the largest relative spacing. 688c uflow is the smallest positive magnitude. 689c 690c***first executable statement dqk15 691 epmach = d1mach(4) 692 uflow = d1mach(1) 693c 694 centr = 0.5d+00*(a+b) 695 hlgth = 0.5d+00*(b-a) 696 dhlgth = dabs(hlgth) 697c 698c compute the 15-point kronrod approximation to 699c the integral, and estimate the absolute error. 700c 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 742 end 743 subroutine dqk21(f,a,b,result,abserr,resabs,resasc,phi,lambda1, 744 & zk0,Pup,Tup,rurd,xflow,kup) 745c***begin prologue dqk21 746c***date written 800101 (yymmdd) 747c***revision date 830518 (yymmdd) 748c***category no. h2a1a2 749c***keywords 21-point gauss-kronrod rules 750c***author piessens,robert,appl. math. & progr. div. - k.u.leuven 751c de doncker,elise,appl. math. & progr. div. - k.u.leuven 752c***purpose to compute i = integral of f over (a,b), with error 753c estimate 754c j = integral of abs(f) over (a,b) 755c***description 756c 757c integration rules 758c standard fortran subroutine 759c double precision version 760c 761c parameters 762c on entry 763c f - double precision 764c function subprogram defining the integrand 765c function f(x). the actual name for f needs to be 766c declared e x t e r n a l in the driver program. 767c 768c a - double precision 769c lower limit of integration 770c 771c b - double precision 772c upper limit of integration 773c 774c on return 775c result - double precision 776c approximation to the integral i 777c result is computed by applying the 21-point 778c kronrod rule (resk) obtained by optimal addition 779c of abscissae to the 10-point gauss rule (resg). 780c 781c abserr - double precision 782c estimate of the modulus of the absolute error, 783c which should not exceed abs(i-result) 784c 785c resabs - double precision 786c approximation to the integral j 787c 788c resasc - double precision 789c approximation to the integral of abs(f-i/(b-a)) 790c over (a,b) 791c 792c***references (none) 793c***routines called d1mach 794c***end prologue dqk21 795c 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 802c 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 809c 810c the abscissae and weights are given for the interval (-1,1). 811c because of symmetry only the positive abscissae and their 812c corresponding weights are given. 813c 814c xgk - abscissae of the 21-point kronrod rule 815c xgk(2), xgk(4), ... abscissae of the 10-point 816c gauss rule 817c xgk(1), xgk(3), ... abscissae which are optimally 818c added to the 10-point gauss rule 819c 820c wgk - weights of the 21-point kronrod rule 821c 822c wg - weights of the 10-point gauss rule 823c 824c 825c gauss quadrature weights and kronron quadrature abscissae and weights 826c as evaluated with 80 decimal digit arithmetic by l. w. fullerton, 827c bell labs, nov. 1981. 828c 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 / 834c 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 / 846c 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 / 858c 859c 860c list of major variables 861c ----------------------- 862c 863c centr - mid point of the interval 864c hlgth - half-length of the interval 865c absc - abscissa 866c fval* - function value 867c resg - result of the 10-point gauss formula 868c resk - result of the 21-point kronrod formula 869c reskh - approximation to the mean value of f over (a,b), 870c i.e. to i/(b-a) 871c 872c 873c machine dependent constants 874c --------------------------- 875c 876c epmach is the largest relative spacing. 877c uflow is the smallest positive magnitude. 878c 879c***first executable statement dqk21 880 epmach = d1mach(4) 881 uflow = d1mach(1) 882c 883 centr = 0.5d+00*(a+b) 884 hlgth = 0.5d+00*(b-a) 885 dhlgth = dabs(hlgth) 886c 887c compute the 21-point kronrod approximation to 888c the integral, and estimate the absolute error. 889c 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 931 end 932 subroutine dqk31(f,a,b,result,abserr,resabs,resasc,phi,lambda1, 933 * zk0,Pup,Tup,rurd,xflow,kup) 934c***begin prologue dqk31 935c***date written 800101 (yymmdd) 936c***revision date 830518 (yymmdd) 937c***category no. h2a1a2 938c***keywords 31-point gauss-kronrod rules 939c***author piessens,robert,appl. math. & progr. div. - k.u.leuven 940c de doncker,elise,appl. math. & progr. div. - k.u.leuven 941c***purpose to compute i = integral of f over (a,b) with error 942c estimate 943c j = integral of abs(f) over (a,b) 944c***description 945c 946c integration rules 947c standard fortran subroutine 948c double precision version 949c 950c parameters 951c on entry 952c f - double precision 953c function subprogram defining the integrand 954c function f(x). the actual name for f needs to be 955c declared e x t e r n a l in the calling program. 956c 957c a - double precision 958c lower limit of integration 959c 960c b - double precision 961c upper limit of integration 962c 963c on return 964c result - double precision 965c approximation to the integral i 966c result is computed by applying the 31-point 967c gauss-kronrod rule (resk), obtained by optimal 968c addition of abscissae to the 15-point gauss 969c rule (resg). 970c 971c abserr - double precison 972c estimate of the modulus of the modulus, 973c which should not exceed abs(i-result) 974c 975c resabs - double precision 976c approximation to the integral j 977c 978c resasc - double precision 979c approximation to the integral of abs(f-i/(b-a)) 980c over (a,b) 981c 982c***references (none) 983c***routines called d1mach 984c***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 991c 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 998c 999c the abscissae and weights are given for the interval (-1,1). 1000c because of symmetry only the positive abscissae and their 1001c corresponding weights are given. 1002c 1003c xgk - abscissae of the 31-point kronrod rule 1004c xgk(2), xgk(4), ... abscissae of the 15-point 1005c gauss rule 1006c xgk(1), xgk(3), ... abscissae which are optimally 1007c added to the 15-point gauss rule 1008c 1009c wgk - weights of the 31-point kronrod rule 1010c 1011c wg - weights of the 15-point gauss rule 1012c 1013c 1014c gauss quadrature weights and kronron quadrature abscissae and weights 1015c as evaluated with 80 decimal digit arithmetic by l. w. fullerton, 1016c bell labs, nov. 1981. 1017c 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 / 1026c 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 / 1043c 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 / 1060c 1061c 1062c list of major variables 1063c ----------------------- 1064c centr - mid point of the interval 1065c hlgth - half-length of the interval 1066c absc - abscissa 1067c fval* - function value 1068c resg - result of the 15-point gauss formula 1069c resk - result of the 31-point kronrod formula 1070c reskh - approximation to the mean value of f over (a,b), 1071c i.e. to i/(b-a) 1072c 1073c machine dependent constants 1074c --------------------------- 1075c epmach is the largest relative spacing. 1076c uflow is the smallest positive magnitude. 1077c***first executable statement dqk31 1078 epmach = d1mach(4) 1079 uflow = d1mach(1) 1080c 1081 centr = 0.5d+00*(a+b) 1082 hlgth = 0.5d+00*(b-a) 1083 dhlgth = dabs(hlgth) 1084c 1085c compute the 31-point kronrod approximation to 1086c the integral, and estimate the absolute error. 1087c 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 1129 end 1130 subroutine dqk41(f,a,b,result,abserr,resabs,resasc,phi,lambda1, 1131 * zk0,Pup,Tup,rurd,xflow,kup) 1132c***begin prologue dqk41 1133c***date written 800101 (yymmdd) 1134c***revision date 830518 (yymmdd) 1135c***category no. h2a1a2 1136c***keywords 41-point gauss-kronrod rules 1137c***author piessens,robert,appl. math. & progr. div. - k.u.leuven 1138c de doncker,elise,appl. math. & progr. div. - k.u.leuven 1139c***purpose to compute i = integral of f over (a,b), with error 1140c estimate 1141c j = integral of abs(f) over (a,b) 1142c***description 1143c 1144c integration rules 1145c standard fortran subroutine 1146c double precision version 1147c 1148c parameters 1149c on entry 1150c f - double precision 1151c function subprogram defining the integrand 1152c function f(x). the actual name for f needs to be 1153c declared e x t e r n a l in the calling program. 1154c 1155c a - double precision 1156c lower limit of integration 1157c 1158c b - double precision 1159c upper limit of integration 1160c 1161c on return 1162c result - double precision 1163c approximation to the integral i 1164c result is computed by applying the 41-point 1165c gauss-kronrod rule (resk) obtained by optimal 1166c addition of abscissae to the 20-point gauss 1167c rule (resg). 1168c 1169c abserr - double precision 1170c estimate of the modulus of the absolute error, 1171c which should not exceed abs(i-result) 1172c 1173c resabs - double precision 1174c approximation to the integral j 1175c 1176c resasc - double precision 1177c approximation to the integal of abs(f-i/(b-a)) 1178c over (a,b) 1179c 1180c***references (none) 1181c***routines called d1mach 1182c***end prologue dqk41 1183c 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 1190c 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 1196c 1197c the abscissae and weights are given for the interval (-1,1). 1198c because of symmetry only the positive abscissae and their 1199c corresponding weights are given. 1200c 1201c xgk - abscissae of the 41-point gauss-kronrod rule 1202c xgk(2), xgk(4), ... abscissae of the 20-point 1203c gauss rule 1204c xgk(1), xgk(3), ... abscissae which are optimally 1205c added to the 20-point gauss rule 1206c 1207c wgk - weights of the 41-point gauss-kronrod rule 1208c 1209c wg - weights of the 20-point gauss rule 1210c 1211c 1212c gauss quadrature weights and kronron quadrature abscissae and weights 1213c as evaluated with 80 decimal digit arithmetic by l. w. fullerton, 1214c bell labs, nov. 1981. 1215c 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 / 1226c 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 / 1248c 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 / 1270c 1271c 1272c list of major variables 1273c ----------------------- 1274c 1275c centr - mid point of the interval 1276c hlgth - half-length of the interval 1277c absc - abscissa 1278c fval* - function value 1279c resg - result of the 20-point gauss formula 1280c resk - result of the 41-point kronrod formula 1281c reskh - approximation to mean value of f over (a,b), i.e. 1282c to i/(b-a) 1283c 1284c machine dependent constants 1285c --------------------------- 1286c 1287c epmach is the largest relative spacing. 1288c uflow is the smallest positive magnitude. 1289c 1290c***first executable statement dqk41 1291 epmach = d1mach(4) 1292 uflow = d1mach(1) 1293c 1294 centr = 0.5d+00*(a+b) 1295 hlgth = 0.5d+00*(b-a) 1296 dhlgth = dabs(hlgth) 1297c 1298c compute the 41-point gauss-kronrod approximation to 1299c the integral, and estimate the absolute error. 1300c 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 1342 end 1343 subroutine dqk51(f,a,b,result,abserr,resabs,resasc,phi,lambda1, 1344 * zk0,Pup,Tup,rurd,xflow,kup) 1345c***begin prologue dqk51 1346c***date written 800101 (yymmdd) 1347c***revision date 830518 (yymmdd) 1348c***category no. h2a1a2 1349c***keywords 51-point gauss-kronrod rules 1350c***author piessens,robert,appl. math. & progr. div. - k.u.leuven 1351c de doncker,elise,appl. math & progr. div. - k.u.leuven 1352c***purpose to compute i = integral of f over (a,b) with error 1353c estimate 1354c j = integral of abs(f) over (a,b) 1355c***description 1356c 1357c integration rules 1358c standard fortran subroutine 1359c double precision version 1360c 1361c parameters 1362c on entry 1363c f - double precision 1364c function subroutine defining the integrand 1365c function f(x). the actual name for f needs to be 1366c declared e x t e r n a l in the calling program. 1367c 1368c a - double precision 1369c lower limit of integration 1370c 1371c b - double precision 1372c upper limit of integration 1373c 1374c on return 1375c result - double precision 1376c approximation to the integral i 1377c result is computed by applying the 51-point 1378c kronrod rule (resk) obtained by optimal addition 1379c of abscissae to the 25-point gauss rule (resg). 1380c 1381c abserr - double precision 1382c estimate of the modulus of the absolute error, 1383c which should not exceed abs(i-result) 1384c 1385c resabs - double precision 1386c approximation to the integral j 1387c 1388c resasc - double precision 1389c approximation to the integral of abs(f-i/(b-a)) 1390c over (a,b) 1391c 1392c***references (none) 1393c***routines called d1mach 1394c***end prologue dqk51 1395c 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 1402c 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 1408c 1409c the abscissae and weights are given for the interval (-1,1). 1410c because of symmetry only the positive abscissae and their 1411c corresponding weights are given. 1412c 1413c xgk - abscissae of the 51-point kronrod rule 1414c xgk(2), xgk(4), ... abscissae of the 25-point 1415c gauss rule 1416c xgk(1), xgk(3), ... abscissae which are optimally 1417c added to the 25-point gauss rule 1418c 1419c wgk - weights of the 51-point kronrod rule 1420c 1421c wg - weights of the 25-point gauss rule 1422c 1423c 1424c gauss quadrature weights and kronron quadrature abscissae and weights 1425c as evaluated with 80 decimal digit arithmetic by l. w. fullerton, 1426c bell labs, nov. 1981. 1427c 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 / 1441c 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 / 1468c 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 / 1494c note: wgk (26) was calculated from the values of wgk(1..25) 1495 data wgk ( 26) / 0.0615808180 6783293507 8759824240 066 d0 / 1496c 1497c 1498c list of major variables 1499c ----------------------- 1500c 1501c centr - mid point of the interval 1502c hlgth - half-length of the interval 1503c absc - abscissa 1504c fval* - function value 1505c resg - result of the 25-point gauss formula 1506c resk - result of the 51-point kronrod formula 1507c reskh - approximation to the mean value of f over (a,b), 1508c i.e. to i/(b-a) 1509c 1510c machine dependent constants 1511c --------------------------- 1512c 1513c epmach is the largest relative spacing. 1514c uflow is the smallest positive magnitude. 1515c 1516c***first executable statement dqk51 1517 epmach = d1mach(4) 1518 uflow = d1mach(1) 1519c 1520 centr = 0.5d+00*(a+b) 1521 hlgth = 0.5d+00*(b-a) 1522 dhlgth = dabs(hlgth) 1523c 1524c compute the 51-point kronrod approximation to 1525c the integral, and estimate the absolute error. 1526c 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 1568 end 1569 subroutine dqk61(f,a,b,result,abserr,resabs,resasc,phi,lambda1, 1570 * zk0,Pup,Tup,rurd,xflow,kup) 1571c***begin prologue dqk61 1572c***date written 800101 (yymmdd) 1573c***revision date 830518 (yymmdd) 1574c***category no. h2a1a2 1575c***keywords 61-point gauss-kronrod rules 1576c***author piessens,robert,appl. math. & progr. div. - k.u.leuven 1577c de doncker,elise,appl. math. & progr. div. - k.u.leuven 1578c***purpose to compute i = integral of f over (a,b) with error 1579c estimate 1580c j = integral of dabs(f) over (a,b) 1581c***description 1582c 1583c integration rule 1584c standard fortran subroutine 1585c double precision version 1586c 1587c 1588c parameters 1589c on entry 1590c f - double precision 1591c function subprogram defining the integrand 1592c function f(x). the actual name for f needs to be 1593c declared e x t e r n a l in the calling program. 1594c 1595c a - double precision 1596c lower limit of integration 1597c 1598c b - double precision 1599c upper limit of integration 1600c 1601c on return 1602c result - double precision 1603c approximation to the integral i 1604c result is computed by applying the 61-point 1605c kronrod rule (resk) obtained by optimal addition of 1606c abscissae to the 30-point gauss rule (resg). 1607c 1608c abserr - double precision 1609c estimate of the modulus of the absolute error, 1610c which should equal or exceed dabs(i-result) 1611c 1612c resabs - double precision 1613c approximation to the integral j 1614c 1615c resasc - double precision 1616c approximation to the integral of dabs(f-i/(b-a)) 1617c 1618c 1619c***references (none) 1620c***routines called d1mach 1621c***end prologue dqk61 1622c 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 1629c 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 1635c 1636c the abscissae and weights are given for the 1637c interval (-1,1). because of symmetry only the positive 1638c abscissae and their corresponding weights are given. 1639c 1640c xgk - abscissae of the 61-point kronrod rule 1641c xgk(2), xgk(4) ... abscissae of the 30-point 1642c gauss rule 1643c xgk(1), xgk(3) ... optimally added abscissae 1644c to the 30-point gauss rule 1645c 1646c wgk - weights of the 61-point kronrod rule 1647c 1648c wg - weigths of the 30-point gauss rule 1649c 1650c 1651c gauss quadrature weights and kronron quadrature abscissae and weights 1652c as evaluated with 80 decimal digit arithmetic by l. w. fullerton, 1653c bell labs, nov. 1981. 1654c 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 / 1670c 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 / 1702c 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 / 1734c 1735c list of major variables 1736c ----------------------- 1737c 1738c centr - mid point of the interval 1739c hlgth - half-length of the interval 1740c dabsc - abscissa 1741c fval* - function value 1742c resg - result of the 30-point gauss rule 1743c resk - result of the 61-point kronrod rule 1744c reskh - approximation to the mean value of f 1745c over (a,b), i.e. to i/(b-a) 1746c 1747c machine dependent constants 1748c --------------------------- 1749c 1750c epmach is the largest relative spacing. 1751c uflow is the smallest positive magnitude. 1752c 1753 epmach = d1mach(4) 1754 uflow = d1mach(1) 1755c 1756 centr = 0.5d+00*(b+a) 1757 hlgth = 0.5d+00*(b-a) 1758 dhlgth = dabs(hlgth) 1759c 1760c compute the 61-point kronrod approximation to the 1761c integral, and estimate the absolute error. 1762c 1763c***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 1805 end 1806 subroutine dqpsrt(limit,last,maxerr,ermax,elist,iord,nrmax, 1807 * phi,lambda1,zk0,Pup,Tup,rurd,xflow,kup) 1808c***begin prologue dqpsrt 1809c***refer to dqage,dqagie,dqagpe,dqawse 1810c***routines called (none) 1811c***revision date 810101 (yymmdd) 1812c***keywords sequential sorting 1813c***author piessens,robert,appl. math. & progr. div. - k.u.leuven 1814c de doncker,elise,appl. math. & progr. div. - k.u.leuven 1815c***purpose this routine maintains the descending ordering in the 1816c list of the local error estimated resulting from the 1817c interval subdivision process. at each call two error 1818c estimates are inserted using the sequential search 1819c method, top-down for the largest error estimate and 1820c bottom-up for the smallest error estimate. 1821c***description 1822c 1823c ordering routine 1824c standard fortran subroutine 1825c double precision version 1826c 1827c parameters (meaning at output) 1828c limit - integer 1829c maximum number of error estimates the list 1830c can contain 1831c 1832c last - integer 1833c number of error estimates currently in the list 1834c 1835c maxerr - integer 1836c maxerr points to the nrmax-th largest error 1837c estimate currently in the list 1838c 1839c ermax - double precision 1840c nrmax-th largest error estimate 1841c ermax = elist(maxerr) 1842c 1843c elist - double precision 1844c vector of dimension last containing 1845c the error estimates 1846c 1847c iord - integer 1848c vector of dimension last, the first k elements 1849c of which contain pointers to the error 1850c estimates, such that 1851c elist(iord(1)),..., elist(iord(k)) 1852c form a decreasing sequence, with 1853c k = last if last.le.(limit/2+2), and 1854c k = limit+1-last otherwise 1855c 1856c nrmax - integer 1857c maxerr = iord(nrmax) 1858c 1859c***end prologue dqpsrt 1860c 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) 1866c 1867c check whether the list contains more than 1868c two error estimates. 1869c 1870c***first executable statement dqpsrt 1871 if(last.gt.2) go to 10 1872 iord(1) = 1 1873 iord(2) = 2 1874 go to 90 1875c 1876c this part of the routine is only executed if, due to a 1877c difficult integrand, subdivision increased the error 1878c estimate. in the normal case the insert procedure should 1879c start after the nrmax-th largest error estimate. 1880c 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) 1886c ***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 1891c 1892c compute the number of elements in the list to be maintained 1893c in descending order. this number depends on the number of 1894c subdivisions still allowed. 1895c 1896 30 jupbn = last 1897 if(last.gt.(limit/2+2)) jupbn = limit+3-last 1898 errmin = elist(last) 1899c 1900c insert errmax by traversing the list top-down, 1901c starting comparison from the element elist(iord(nrmax+1)). 1902c 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) 1908c ***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 1915c 1916c insert errmin by traversing the list bottom-up. 1917c 1918 60 iord(i-1) = maxerr 1919 k = jbnd 1920 do 70 j=i,jbnd 1921 isucc = iord(k) 1922c ***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 1930c 1931c set maxerr and ermax. 1932c 1933 90 maxerr = iord(nrmax) 1934 ermax = elist(maxerr) 1935 return 1936 end 1937