MathLib.f

00001 *//////////////////////////////////////////////////////////////////////////////////////////////////
00002 *//                                                                                              //
00003 *//                                                                                              //
00004 *//                     Pseudo-CLASS  Mathlib                                                    //
00005 *//                                                                                              //
00006 *//      Purpose:  library of math utilies                                                       //
00007 *//                                                                                              //
00008 *//      SUBROUTINE Mathlib_GausJad(fun,aa,bb,eeps,result)     : Gauss integration               //
00009 *//      DOUBLE PRECISION FUNCTION Mathlib_Gauss(f,a,b,eeps)   : Gauss integration               //
00010 *//      DOUBLE PRECISION FUNCTION Mathlib_dilogy(x)           : Dilog function Li_2             //
00011 *//      DOUBLE PRECISION FUNCTION Mathlib_dpgamm(z)           : Euler Gamma function            //
00012 *//                                                                                              //
00013 *//////////////////////////////////////////////////////////////////////////////////////////////////
00014 
00015 
00016       SUBROUTINE Mathlib_GausJad(fun,aa,bb,eeps,result)
00017 *//////////////////////////////////////////////////////////////////////////////
00018 *//                                                                          //
00019 *//  Gauss-type integration by S. Jadach, Oct. 1990, June 1997               //
00020 *//  this is non-adaptive (!!!!) unoptimized (!!!) integration subprogram.   //
00021 *//                                                                          //
00022 *//  Eeps>0 treated as ABSOLUTE requested error                              //
00023 *//  Eeps<0 treated as RELATIVE requested error                              //
00024 *//                                                                          //
00025 *//////////////////////////////////////////////////////////////////////////////
00026       IMPLICIT NONE
00027 *
00028       DOUBLE PRECISION  fun,aa,bb,eeps,result
00029       EXTERNAL fun
00030 *
00031       DOUBLE PRECISION  a,b,xplus,sum16,range,sum8,erabs,erela,fminu,xminu
00032       DOUBLE PRECISION  fplus,xmidle,calk8,eps,x1,x2,delta,calk16
00033       INTEGER           iter,ndivi,itermx,k,i
00034       DOUBLE PRECISION  wg(12),xx(12)
00035       DATA wg
00036      $/0.101228536290376d0, 0.222381034453374d0, 0.313706645877887d0,
00037      $ 0.362683783378362d0, 0.027152459411754d0, 0.062253523938648d0,
00038      $ 0.095158511682493d0, 0.124628971255534d0, 0.149595988816577d0,
00039      $ 0.169156519395003d0, 0.182603415044924d0, 0.189450610455069d0/
00040       DATA xx
00041      $/0.960289856497536d0, 0.796666477413627d0, 0.525532409916329d0,
00042      $ 0.183434642495650d0, 0.989400934991650d0, 0.944575023073233d0,
00043      $ 0.865631202387832d0, 0.755404408355003d0, 0.617876244402644d0,
00044      $ 0.458016777657227d0, 0.281603550779259d0, 0.095012509837637d0/
00045       DATA itermx / 15/
00046 *-------------------------------------------------------------------------------
00047       a  = aa
00048       b  = bb
00049       eps= ABS(eeps)
00050       ndivi=1
00051 *** iteration over subdivisions terminated by precision requirement
00052       DO iter=1,itermx
00053          calk8  =0d0
00054          calk16 =0d0
00055 ***   sum over delta subintegrals
00056          DO k = 1,ndivi
00057             delta = (b-a)/ndivi
00058             x1    =  a + (k-1)*delta
00059             x2    =  x1+ delta
00060             xmidle= 0.5d0*(x2+x1)
00061             range = 0.5d0*(x2-x1)
00062             sum8 =0d0
00063             sum16=0d0
00064 ***   8- and 12-point   gauss integration over single delta subinterval
00065             DO i=1,12
00066                xplus= xmidle+range*xx(i)
00067                xminu= xmidle-range*xx(i)
00068                fplus=fun(xplus)
00069                fminu=fun(xminu)
00070                IF(i .LE. 4) THEN
00071                   sum8 =sum8  +(fplus+fminu)*wg(i)/2d0
00072                ELSE
00073                   sum16=sum16 +(fplus+fminu)*wg(i)/2d0
00074                ENDIF
00075             ENDDO
00076             calk8 = calk8 + sum8 *(x2-x1)
00077             calk16= calk16+ sum16*(x2-x1)
00078          ENDDO
00079          erabs = ABS(calk16-calk8)
00080          erela = 0d0
00081          IF(calk16 .NE. 0d0) erela= erabs/ABS(calk16)
00082 ***   WRITE(*,*) 'gausjad: calk8,calk16=',iter,calk8,calk16,erela
00083 ***   precision check to terminate integration
00084          IF(eeps .GT. 0d0) THEN
00085             IF(erabs .LT.  eps) GOTO 800
00086          ELSE
00087             IF(erela .LT.  eps) GOTO 800
00088          ENDIF
00089          ndivi=ndivi*2
00090       ENDDO
00091       WRITE(*,*) '++++ Mathlib_GausJad: required precision to high!'
00092       WRITE(*,*) '++++ Mathlib_GausJad: eeps,erela,erabs=',eeps,erela,erabs
00093   800 CONTINUE
00094       result = calk16
00095       END
00096 
00097 
00098       DOUBLE PRECISION FUNCTION Mathlib_Gauss(f,a,b,eeps)
00099 *//////////////////////////////////////////////////////////////////////////////
00100 *//                                                                          //
00101 *// this is iterative integration procedure                                  //
00102 *// originates  probably from cern library                                   //
00103 *// it subdivides inegration range until required PRECISION is reached       //
00104 *// PRECISION is a difference from 8 and 16 point gauss itegr. result        //
00105 *// eeps positive treated as absolute PRECISION                              //
00106 *// eeps negative treated as relative PRECISION                              //
00107 *//                                                                          //
00108 *//////////////////////////////////////////////////////////////////////////////
00109       IMPLICIT NONE
00110       DOUBLE PRECISION  f,a,b,eeps
00111 *
00112       DOUBLE PRECISION  c1,c2,bb,s8,s16,y,aa,const,delta,eps,u
00113       INTEGER           i
00114 *
00115       DOUBLE PRECISION  w(12),x(12)
00116       EXTERNAL f
00117       DATA const /1.0d-19/
00118       DATA w
00119      1/0.10122 85362 90376, 0.22238 10344 53374, 0.31370 66458 77887,
00120      2 0.36268 37833 78362, 0.02715 24594 11754, 0.06225 35239 38648,
00121      3 0.09515 85116 82493, 0.12462 89712 55534, 0.14959 59888 16577,
00122      4 0.16915 65193 95003, 0.18260 34150 44924, 0.18945 06104 55069/
00123       DATA x
00124      1/0.96028 98564 97536, 0.79666 64774 13627, 0.52553 24099 16329,
00125      2 0.18343 46424 95650, 0.98940 09349 91650, 0.94457 50230 73233,
00126      3 0.86563 12023 87832, 0.75540 44083 55003, 0.61787 62444 02644,
00127      4 0.45801 67776 57227, 0.28160 35507 79259, 0.09501 25098 37637/
00128 *-----------------------------------------------------------------------------
00129       eps=abs(eeps)
00130       delta=const*abs(a-b)
00131       Mathlib_Gauss=0d0
00132       aa=a
00133     5 y=b-aa
00134       IF(abs(y)  .LE.  delta) RETURN
00135     2 bb=aa+y
00136       c1=0.5d0*(aa+bb)
00137       c2=c1-aa
00138       s8=0d0
00139       s16=0d0
00140       DO 1 i=1,4
00141       u=x(i)*c2
00142     1 s8=s8+w(i)*(f(c1+u)+f(c1-u))
00143       DO 3 i=5,12
00144       u=x(i)*c2
00145     3 s16=s16+w(i)*(f(c1+u)+f(c1-u))
00146       s8=s8*c2
00147       s16=s16*c2
00148       IF(eeps .LT. 0d0) THEN
00149         IF(abs(s16-s8)  .GT.  eps*abs(s16)) GOTO 4
00150       ELSE
00151         IF(abs(s16-s8)  .GT.  eps) GOTO 4
00152       ENDIF
00153       Mathlib_Gauss=Mathlib_Gauss+s16
00154       aa=bb
00155       GOTO 5
00156     4 y=0.5d0*y
00157       IF(abs(y)  .GT.  delta) GOTO 2
00158       WRITE(*,7)
00159       Mathlib_Gauss=0d0
00160       RETURN
00161     7 FORMAT(1x,36hgaus  ... too high accuracy required)
00162       END
00163 
00164 
00165       DOUBLE PRECISION FUNCTION Mathlib_dilogy(x)
00166 *//////////////////////////////////////////////////////////////////////////////
00167 *//                                                                          //
00168 *// dilogarithm FUNCTION: dilog(x)=int( -ln(1-z)/z ) , 0 < z < x .           //
00169 *// this is the cernlib version.                                             //
00170 *//                                                                          //
00171 *//////////////////////////////////////////////////////////////////////////////
00172       IMPLICIT NONE
00173       DOUBLE PRECISION x
00174 * locals
00175       DOUBLE PRECISION a,b,s,t,y,z
00176 *------------------------------------------------------------------------------
00177       z=-1.644934066848226d0
00178       IF(x  .LT. -1.d0) go to 1
00179       IF(x  .LE.  0.5d0) go to 2
00180       IF(x  .EQ.  1.d0) go to 3
00181       IF(x  .LE.  2.d0) go to 4
00182       z=3.289868133696453d0
00183     1 t=1.d0/x
00184       s=-0.5d0
00185       z=z-0.5d0*dlog(dabs(x))**2
00186       go to 5
00187     2 t=x
00188       s=0.5d0
00189       z=0.d0
00190       go to 5
00191     3 Mathlib_Dilogy=1.644934066848226d0
00192       RETURN
00193     4 t=1.d0-x
00194       s=-0.5d0
00195       z=1.644934066848226d0-dlog(x)*dlog(dabs(t))
00196     5 y=2.666666666666667d0*t+0.666666666666667d0
00197       b=      0.000000000000001d0
00198       a=y*b  +0.000000000000004d0
00199       b=y*a-b+0.000000000000011d0
00200       a=y*b-a+0.000000000000037d0
00201       b=y*a-b+0.000000000000121d0
00202       a=y*b-a+0.000000000000398d0
00203       b=y*a-b+0.000000000001312d0
00204       a=y*b-a+0.000000000004342d0
00205       b=y*a-b+0.000000000014437d0
00206       a=y*b-a+0.000000000048274d0
00207       b=y*a-b+0.000000000162421d0
00208       a=y*b-a+0.000000000550291d0
00209       b=y*a-b+0.000000001879117d0
00210       a=y*b-a+0.000000006474338d0
00211       b=y*a-b+0.000000022536705d0
00212       a=y*b-a+0.000000079387055d0
00213       b=y*a-b+0.000000283575385d0
00214       a=y*b-a+0.000001029904264d0
00215       b=y*a-b+0.000003816329463d0
00216       a=y*b-a+0.000014496300557d0
00217       b=y*a-b+0.000056817822718d0
00218       a=y*b-a+0.000232002196094d0
00219       b=y*a-b+0.001001627496164d0
00220       a=y*b-a+0.004686361959447d0
00221       b=y*a-b+0.024879322924228d0
00222       a=y*b-a+0.166073032927855d0
00223       a=y*a-b+1.935064300869969d0
00224       Mathlib_Dilogy = s*t*(a-b)+z
00225       END
00226 
00227 
00228       DOUBLE PRECISION FUNCTION Mathlib_dpgamm(z)
00229 *//////////////////////////////////////////////////////////////////////////////
00230 *//                                                                          //
00231 *//            Double precision gamma function                               //
00232 *//                                                                          //
00233 *//////////////////////////////////////////////////////////////////////////////
00234       DOUBLE PRECISION z,z1,x,x1,x2,d1,d2,s1,s2,s3,pi,c(20),const
00235       SAVE c,pi,const
00236       DATA c( 1) / 8.3333333333333333333333333332d-02/
00237       DATA c( 2) /-2.7777777777777777777777777777d-03/
00238       DATA c( 3) / 7.9365079365079365079365079364d-04/
00239       DATA c( 4) /-5.9523809523809523809523809523d-04/
00240       DATA c( 5) / 8.4175084175084175084175084175d-04/
00241       DATA c( 6) /-1.9175269175269175269175269175d-03/
00242       DATA c( 7) / 6.4102564102564102564102564102d-03/
00243       DATA c( 8) /-2.9550653594771241830065359477d-02/
00244       DATA c( 9) / 1.7964437236883057316493849001d-01/
00245       DATA c(10) /-1.3924322169059011164274322169d+00/
00246       DATA c(11) / 1.3402864044168391994478951001d+01/
00247       DATA c(12) /-1.5684828462600201730636513245d+02/
00248       DATA c(13) / 2.1931033333333333333333333333d+03/
00249       DATA c(14) /-3.6108771253724989357173265219d+04/
00250       DATA c(15) / 6.9147226885131306710839525077d+05/
00251       DATA c(16) /-1.5238221539407416192283364959d+07/
00252       DATA c(17) / 3.8290075139141414141414141414d+08/
00253       DATA c(18) /-1.0882266035784391089015149165d+10/
00254       DATA c(19) / 3.4732028376500225225225225224d+11/
00255       DATA c(20) /-1.2369602142269274454251710349d+13/
00256       DATA pi    / 3.1415926535897932384626433832d+00/
00257       DATA const / 9.1893853320467274178032973641d-01/
00258       IF(z .GT. 5.75d 1)                                     GOTO  6666
00259       nn = z
00260       IF (z  -  dble(float(nn)))                 3,1,3
00261     1 IF (z     .LE.     0.d 0)                    GOTO 6667
00262       Mathlib_dpgamm = 1.d 0
00263       IF (z     .LE.     2.d 0)                    RETURN
00264       z1 = z
00265     2 z1 = z1  -  1.d 0
00266       Mathlib_dpgamm = Mathlib_dpgamm * z1
00267       IF (z1  -  2.d 0)                          61,61,2
00268     3 IF (dabs(z)     .LT.     1.d-29)             GOTO 60
00269       IF (z     .LT.     0.d 0)                    GOTO 4
00270       x  = z
00271       kk = 1
00272       GOTO 10
00273     4 x  = 1.d 0  -  z
00274       kk = 2
00275    10 x1 = x
00276       IF (x     .GT.     19.d 0)                   GOTO 13
00277       d1 = x
00278    11 x1 = x1  +  1.d 0
00279       IF (x1     .GE.     19.d 0)                  GOTO 12
00280       d1 = d1 * x1
00281       GOTO 11
00282    12 s3 = -dlog(d1)
00283       GOTO 14
00284    13 s3 = 0.d 0
00285    14 d1 = x1 * x1
00286       s1 = (x1  -  5.d-1) * dlog(x1)  -  x1  +  const
00287       DO 20                  k=1,20
00288       s2 = s1  +  c(k)/x1
00289       IF (dabs(s2  -  s1)     .LT.     1.d-28)     GOTO 21
00290       x1 = x1 * d1
00291    20 s1 = s2
00292    21 s3 = s3  +  s2
00293       GOTO (50,22),    kk
00294    22 d2 = dabs(z  -  nn)
00295       d1 = d2 * pi
00296       IF (d1     .LT.     1.d-15)                  GOTO 31
00297    30 x2 =  dlog(pi/dsin(d1))  -  s3
00298       GOTO 40
00299    31 x2 = -dlog(d2)
00300    40 mm = dabs(z)
00301       IF(x2       .GT.       1.74d2)                  GOTO 6666
00302       Mathlib_dpgamm = dexp(x2)
00303       IF (mm    .ne.    (mm/2) * 2)              RETURN
00304       Mathlib_dpgamm = -Mathlib_dpgamm
00305       RETURN
00306    50 IF(s3       .GT.       1.74d2)                  GOTO 6666
00307       Mathlib_dpgamm = dexp(s3)
00308       RETURN
00309  6666 print *, 2000
00310       RETURN
00311  6667 print *, 2001
00312       RETURN
00313    60 Mathlib_dpgamm = 0.d 0
00314       IF(dabs(z)    .LT.    1.d-77)   RETURN
00315       Mathlib_dpgamm = 1.d 0/z
00316    61 RETURN
00317  2000 FORMAT (/////, 2x, 32hdpgamm ..... argument too large., /////)
00318  2001 FORMAT (/////, 2x, 32hdpgamm ..... argument is a pole., /////)
00319       END
00320 
00321 
00322 
00323       SUBROUTINE Mathlib_Gaus8(fun,aa,bb,result)
00324 *//////////////////////////////////////////////////////////////////////////////
00325 *//   8-point Gauss                                                          //
00326 *//////////////////////////////////////////////////////////////////////////////
00327       IMPLICIT NONE
00328       DOUBLE PRECISION  fun,aa,bb,result
00329       EXTERNAL fun
00330       DOUBLE PRECISION  a,b,sum8,xmidle,range,xplus,xminu
00331       INTEGER           k,i
00332       DOUBLE PRECISION  wg(4),xx(4)
00333       DATA wg /0.101228536290376d0, 0.222381034453374d0, 0.313706645877887d0, 0.362683783378362d0/
00334       DATA xx /0.960289856497536d0, 0.796666477413627d0, 0.525532409916329d0, 0.183434642495650d0/
00335 *-------------------------------------------------------------------------------
00336       a  = aa
00337       b  = bb
00338       xmidle= 0.5d0*(a+b)
00339       range = 0.5d0*(b-a)
00340       sum8 =0d0
00341       DO i=1,4
00342          xplus= xmidle+range*xx(i)
00343          xminu= xmidle-range*xx(i)
00344          sum8 =sum8  +(fun(xplus)+fun(xminu))*wg(i)/2d0
00345       ENDDO
00346       result = sum8*(b-a)
00347       END
00348 
00349       SUBROUTINE Mathlib_Gaus16(fun,aa,bb,result)
00350 *//////////////////////////////////////////////////////////////////////////////
00351 *//   12-point Gauss                                                          //
00352 *//////////////////////////////////////////////////////////////////////////////
00353       IMPLICIT NONE
00354       DOUBLE PRECISION  fun,aa,bb,result
00355       EXTERNAL fun
00356       DOUBLE PRECISION  a,b,sum16,xmidle,range,xplus,xminu
00357       INTEGER           k,i
00358       DOUBLE PRECISION  wg(8),xx(8)
00359       DATA wg              /0.027152459411754d0, 0.062253523938648d0,
00360      $ 0.095158511682493d0, 0.124628971255534d0, 0.149595988816577d0,
00361      $ 0.169156519395003d0, 0.182603415044924d0, 0.189450610455069d0/
00362       DATA xx              /0.989400934991650d0, 0.944575023073233d0,
00363      $ 0.865631202387832d0, 0.755404408355003d0, 0.617876244402644d0,
00364      $ 0.458016777657227d0, 0.281603550779259d0, 0.095012509837637d0/
00365 *-------------------------------------------------------------------------------
00366       a  = aa
00367       b  = bb
00368       xmidle= 0.5d0*(a+b)
00369       range = 0.5d0*(b-a)
00370       sum16 =0d0
00371       DO i=1,8
00372          xplus= xmidle+range*xx(i)
00373          xminu= xmidle-range*xx(i)
00374          sum16 =sum16  +(fun(xplus)+fun(xminu))*wg(i)/2d0
00375       ENDDO
00376       result = sum16*(b-a)
00377       END
00378 
00379 
00380 
00381 
00382 *//////////////////////////////////////////////////////////////////////////////
00383 *//                                                                          //
00384 *//                End  Pseudo-CLASS  Mathlib                                //
00385 *//////////////////////////////////////////////////////////////////////////////