subroutine GAUSS(npts, kode, a, b, xs, wts) c *** see r h landau, program lpott, 1981 c c subroutine to return a scalled set of gaussian-legendre points c *npts* is the number of points and must be .gt. 1 c *kode* is a three digit (decimal) number -- htu c *u* = 0 means *npts* must be in the set *numbrs*, else syserr is c called. *npts* may be either a constant or variable. c = 1 means if *npts* isn*t in the set *numbrs*, use the larges c number .lt. *npts* and change *npts* accordingly. c *npts* must be a variable. c = 2 means if *npts* isn*t in the set *numbrs*, use the c smallest member .gt. *npts* and change *npts* accordingal c however if *npts* .gt. numbrs(istop) then use numbrs(ist c *npts* must be a variable. c = 5 means use equal spaced points. for odd *npts* they c will be simpson points and weights. for even *npts* they c will have equal weights. c note that in this case *npts* may be any value .ge. 1. c note that they are still scalled according to *t*. c *t* determines how the internal points (on c (-1, +1)) are scalled. see the computed go to c after stmnt no 109 for all the details. c *h* = 0 means scale as per t above. c = 1 means scale as per *t* above but also remove square c root at lower bound. this removes integrable c 1/sqrt(x-lower) discontinuities and gives better c values for sqrt(x-lower) f(x) type integrands. c *a* and *b* are used in the scalling c *xs* will be the set of points, must be real dimensioned .gt. n c *ws* will be the set of weights must be real dimensioned .gt. n c implicit double precision (a-h,o-z) integer t, u, base dimension xs(npts), wts(npts) dimension hh(2) dimension gx(118), wt(118) dimension x1(59), w1(59), x2(59), w2(59) equivalence (x1(1),gx(1)),(w1(1),wt(1)),(x2(1),gx(60)) 1 ,(w2(1),wt(60)) integer numbrs(13) data istop /13/ data numbrs / 2, 4, 6, 8, 10, 12, 14, 16, 20, 24, 32, 40, 48 / data x1 / 0.5773502691896261d0, 0.8611363115940530d0, & 0.3399810435848560d0, 0.9324695142031520d0, 0.6612093864662651d0, & 0.2386191860831970d0, 0.9602898564975360d0, 0.7966664774136270d0, & 0.5255324099163290d0, 0.1834346424956500d0, 0.9739065285171720d0, & 0.8650633666889850d0, 0.6794095682990241d0, 0.4333953941292470d0, & 0.1488743389816310d0, 0.9815606342467190d0, 0.9041172563704750d0, & 0.7699026741943050d0, 0.5873179542866171d0, 0.3678314989981800d0, & 0.1252334085114690d0, 0.9862838086968120d0, 0.9284348836635741d0, & 0.8272013150697650d0, 0.6872929048116851d0, 0.5152486363581540d0, & 0.3191123689278900d0, 0.1080549487073440d0, 0.9894009349916501d0, & 0.9445750230732330d0, 0.8656312023878320d0, 0.7554044083550030d0, & 0.6178762444026441d0, 0.4580167776572270d0, 0.2816035507792590d0, & 0.0950125098376370d0, 0.9931285991850949d0, 0.9639719272779138d0, & 0.9122344282513259d0, 0.8391169718222189d0, 0.7463319064601507d0, & 0.6360536807265151d0, 0.5108670019508270d0, 0.3737060887154195d0, & 0.2277858511416450d0, 0.0765265211334973d0, 0.9951872199970213d0, & 0.9747285559713094d0, 0.9382745520027326d0, 0.8864155270044010d0, & 0.8200019859739029d0, 0.7401241915785542d0, 0.6480936519369756d0, & 0.5454214713888396d0, 0.4337935076260451d0, 0.3150426796961633d0, & 0.1911188674736163d0, 0.0640568928626056d0, 0.9972638618494816d0/ data x2 / 0.9856115115452683d0, 0.9647622555875064d0, & 0.9349060759377397d0, 0.8963211557660522d0, 0.8493676137325699d0, & 0.7944837959679424d0, 0.7321821187402896d0, 0.6630442669302153d0, & 0.5877157572407623d0, 0.5068999089322293d0, 0.4213512761306353d0, & 0.3318686022821276d0, 0.2392873622521370d0, 0.1444719615827964d0, & 0.0483076656877383d0, 0.9982377097105592d0, 0.9907262386994570d0, & 0.9772599499837742d0, 0.9579168192137917d0, 0.9328128082786765d0, & 0.9020988069688742d0, 0.8659595032122595d0, 0.8246122308333115d0, & 0.7783056514265194d0, 0.7273182551899270d0, 0.6719566846141796d0, & 0.6125538896679803d0, 0.5494671250951283d0, 0.4830758016861787d0, & 0.4137792043716050d0, 0.3419940908257584d0, 0.2681521850072536d0, & 0.1926975807013710d0, 0.1160840706752552d0, 0.0387724175060508d0, & 0.9987710072524261d0, 0.9935301722663507d0, 0.9841245837228269d0, & 0.9705915925462472d0, 0.9529877031604309d0, 0.9313866907065542d0, & 0.9058791367155696d0, 0.8765720202742478d0, 0.8435882616243934d0, & 0.8070662040294426d0, 0.7671590325157403d0, 0.7240341309238146d0, & 0.6778723796326640d0, 0.6288673967765137d0, 0.5772247260839728d0, & 0.5231609747222331d0, 0.4669029047509584d0, 0.4086864819907167d0, & 0.3487558862921607d0, 0.2873624873554555d0, 0.2247637903946890d0, & 0.1612223560688917d0, 0.0970046992094627d0, 0.0323801709628694d0/ data w1 / 1.0000000000000000d0, 0.3478548451374540d0, & 0.6521451548625461d0, 0.1713244923791700d0, 0.3607615730481390d0, & 0.4679139345726910d0, 0.1012285362903760d0, 0.2223810344533740d0, & 0.3137066458778870d0, 0.3626837833783620d0, 0.0666713443086880d0, & 0.1494513491505810d0, 0.2190863625159820d0, 0.2692667193099960d0, & 0.2955242247147530d0, 0.0471753363865120d0, 0.1069393259953180d0, & 0.1600783285433460d0, 0.2031674267230660d0, 0.2334925365383550d0, & 0.2491470458134030d0, 0.0351194603317520d0, 0.0801580871597600d0, & 0.1215185706879030d0, 0.1572031671581940d0, 0.1855383974779380d0, & 0.2051984637212960d0, 0.2152638534631580d0, 0.0271524594117540d0, & 0.0622535239386480d0, 0.0951585116824930d0, 0.1246289712555340d0, & 0.1495959888165770d0, 0.1691565193950030d0, 0.1826034150449240d0, & 0.1894506104550690d0, 0.0176140071391521d0, 0.0406014298003869d0, & 0.0626720483341091d0, 0.0832767415767047d0, 0.1019301198172404d0, & 0.1181945319615184d0, 0.1316886384491766d0, 0.1420961093183820d0, & 0.1491729864726037d0, 0.1527533871307258d0, 0.0123412297999872d0, & 0.0285313886289337d0, 0.0442774388174198d0, 0.0592985849154368d0, & 0.0733464814110803d0, 0.0861901615319533d0, 0.0976186521041139d0, & 0.1074442701159656d0, 0.1155056680537256d0, 0.1216704729278033d0, & 0.1258374563468282d0, 0.1279381953467521d0, 0.0070186100094701d0/ data w2 / 0.0162743947309057d0, 0.0253920653092620d0, & 0.0342738629130214d0, 0.0428358980222267d0, 0.0509980592623762d0, & 0.0586840934785355d0, 0.0658222227763618d0, 0.0723457941088485d0, & 0.0781938957870703d0, 0.0833119242269467d0, 0.0876520930044038d0, & 0.0911738786957639d0, 0.0938443990808046d0, 0.0956387200792748d0, & 0.0965400885147278d0, 0.0045212770985332d0, 0.0104982845311528d0, & 0.0164210583819079d0, 0.0222458491941669d0, 0.0279370069800234d0, & 0.0334601952825478d0, 0.0387821679744720d0, 0.0438709081856733d0, & 0.0486958076350722d0, 0.0532278469839368d0, 0.0574397690993916d0, & 0.0613062424929289d0, 0.0648040134566010d0, 0.0679120458152339d0, & 0.0706116473912868d0, 0.0728865823958040d0, 0.0747231690579683d0, & 0.0761103619006262d0, 0.0770398181642480d0, 0.0775059479784248d0, & 0.0031533460523058d0, 0.0073275539012763d0, 0.0114772345792345d0, & 0.0155793157229438d0, 0.0196161604573555d0, 0.0235707608393244d0, & 0.0274265097083569d0, 0.0311672278327981d0, 0.0347772225647704d0, & 0.0382413510658307d0, 0.0415450829434647d0, 0.0446745608566943d0, & 0.0476166584924905d0, 0.0503590355538545d0, 0.0528901894851937d0, & 0.0551995036999842d0, 0.0572772921004032d0, 0.0591148396983956d0, & 0.0607044391658939d0, 0.0620394231598927d0, 0.0631141922862540d0, & 0.0639242385846482d0, 0.0644661644359501d0, 0.0647376968126839d0/ c c this *1 - epsilon* is for simpson points so that when we scale c to infinity we won*t divide by zero. c data one / 1.00000000001d0 / c c FIRST EXECUTABLE >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> t = kode/10 u = kode-10*t h = t/10 t = t-10*h if ((h .lt. 0) .or. (h .gt. 1)) go to 90 if (npts .lt. 1) go to 95 if (u .eq. 5) go to 400 if ((u .lt. 0) .or. (u .gt. 2)) go to 90 base = 0 do 30 i = 1, istop if (numbrs(i) .eq. npts) go to 100 if (numbrs(i) .gt. npts) go to 50 30 base = base + (numbrs(i)+1)/2 c npts .gr. last if (u .eq. 0) go to 95 npts = numbrs(istop) base = base - (npts+1)/2 go to 100 c 50 if ((npts .le. 1) .or. (u .eq. 0)) go to 95 go to (55, 60), u 55 npts = numbrs(i-1) base = base - (npts+1)/2 go to 100 60 npts = numbrs(i) go to 100 c 90 write(6,91) kode 91 format (32h0***gauss*** invalid ""kode"" -- , i20) stop 95 write(6,96) npts, kode 96 format(28h0gauss,invalid npts,ok kode= ,2i20) stop c 100 do 300 n = 1, npts if (u .eq. 5) go to 108 if (n .gt. (npts+1)/2) go to 105 x = -gx(n+base) w = wt(n+base) go to 109 105 nn = npts-n+1 x = gx(nn+base) w = wt(nn+base) go to 109 108 x = xs(n) w = wts(n) c are we removeing sqrae root at lower. 109 if (h .eq. 1) w = (x+1) * w if (h .eq. 1) x = .5 * (x*(x+2)-1) go to(110, 120, 130, 140, 150,160 ), t 160 go to 90 c c 1, scale to (a, b) uniformly 110 xs(n) = (b+a)/2 + (b-a)/2*x wts(n) = (b-a)/2 * w go to 300 c c 2, scale to (0, infinity) 120 xs(n) = a*(1+x)/(1-x) wts(n) = 2*a*w/(1-x)**2 go to 300 c c 3, scale to (-infinity, infinity) 130 xs(n) = a*x/(1-x**2) wts(n) = a*(1+x**2)/(1-x**2)**2 go to 300 c c 4 scale to (b, infinity) so that for b = 0 we have case 2 140 xs(n) = (a+2*b + a*x)/(1-x) wts(n) = 2*(b+a) * w/(1-x)**2 go to 300 c c 5, scale to (0,b) so that for b = infinity we have case 2 c (almost) 150 xs(n) = a*b * (1+x)/(b+a - (b-a)*x) wts(n) = 2*a*b**2 * w/(b+a - (b-a)*x)**2 go to 300 300 continue return c c uniform spaced points for any n. c 400 if ((npts/2)*2 .eq. npts) go to 450 c c odd number = sumpsons if (npts .gt. 1) go to 410 xs(1) = 0 wts(1) = 2 go to 100 410 h = 2. / (npts-1) hb3 = h/3 xs(1) = -1 xs(npts) = 1 if (t .eq. 3) xs(1) = -one if ((t .ge. 2) .and. (t .le. 4)) xs(npts) = one wts(1) = hb3 wts(npts) = hb3 hh(1) = 4*hb3 hh(2) = 2*hb3 nb2 = npts/2 ih = 1 do 430 i = 1, nb2 xs(i+1) = -1 + i*h xs(npts-i) = 1-i*h wts(i+1) = hh(ih) wts(npts-i) = hh(ih) 430 ih = 3-ih go to 100 c c even number = equal weight 450 h = 2./npts do 480 i = 1, npts xs(i) = -1-h/2 + i*h 480 wts(i) = h go to 100 c end