PDF/spdf/whit1q.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.2  1996/10/30 08:31:59  cernlib
   6 * Version 7.04
   7 *
   8 * Revision 1.1.1.1  1996/04/12 15:30:12  plothow
   9 * Version 7.01
  10 *
  11 *
  12 #include "pdf/pilot.h"
  13 c
  14 cccccccccccccccccccccccccccccc
  15 c   QPM calculation
  16       subroutine WHIT1Q(x,r,cv,cs)
  17 ccc INPUTS : x,r=mc^2/Q^2
  18 ccc OUTPUTS:   cv,cs (valence- and sea- charm quark dist)
  19 ccc                  cv <-- cv / ( alpha / 2PI)
  20 ccc                  cs <-- cs / ( alpha_s/2PI)
  21 c
  22       implicit none
  23 c arg
  24       DOUBLE PRECISION x,r
  25       DOUBLE PRECISION cv,cs
  26 c CONST
  27       DOUBLE PRECISION ec,mc
  28       parameter(ec=2.0d0/3.0d0,mc=1.5d0)
  29 c  N=15 Gauss int. weights and points
  30       integer GN,i
  31       parameter(GN=15)
  32       DOUBLE PRECISION XG(GN), XW(GN)
  33       DATA (XG(i),i=1,GN)/6.003741d-03,
  34      $ 3.136330d-02, 7.589671d-02, 1.377911d-01, 2.145139d-01,
  35      $ 3.029243d-01, 3.994030d-01, 5.000000d-01, 6.005970d-01,
  36      $ 6.970757d-01, 7.854861d-01, 8.622089d-01, 9.241033d-01,
  37      $ 9.686367d-01, 9.939963d-01/
  38       DATA (XW(i),i=1,GN)/1.537662d-02,
  39      $ 3.518302d-02, 5.357961d-02, 6.978534d-02, 8.313460d-02,
  40      $ 9.308050d-02, 9.921574d-02, 1.012891d-01, 9.921574d-02,
  41      $ 9.308050d-02, 8.313460d-02, 6.978534d-02, 5.357961d-02,
  42      $ 3.518302d-02, 1.537662d-02/
  43 c local
  44       DOUBLE PRECISION sum, y,z,beta,w,WHIT1G,L
  45       parameter(L=4.0d0)
  46       DOUBLE PRECISION x1,rx,z1,rz
  47 c
  48 c begin
  49       x1=1.0d0-x
  50       rx=4.0d0*r*x
  51 c
  52 c direct
  53          beta=dsqrt(1.0d0-rx/x1)
  54          w=x*( beta*(-1.0d0+8.0d0*x*x1-rx*x1)
  55      $        +(x**2+x1**2+rx*(1.0d0-3.0d0*x)-0.5d0*rx**2)
  56      $         *  log( (1.0d0+beta)/(1.0d0-beta) ))
  57          cv = 3.0d0*ec**2 * w / x
  58 c
  59 c resolved
  60       sum=0.0d0
  61       do 10 i=1,GN
  62          y= x+rx + (x1-rx)*XG(i)**L
  63          z=x/y
  64          z1=1.0d0-z
  65          rz=4.0d0*r*z
  66          beta=dsqrt(1.0d0-rz/z1)
  67          w=z*( beta*(-1.0d0+8.0d0*z*z1-rz*z1)
  68      $        +(z**2+z1**2+rz*(1.0d0-3.0d0*z)-0.5d0*rz**2)
  69      $         *  log( (1.0d0+beta)/(1.0d0-beta) ))
  70 c
  71          sum= sum + w * WHIT1G(y,mc**2/r)* L*XG(i)**(L-1.0d0)*XW(i)
  72 c
  73    10 continue
  74 c
  75          cs = 0.5d0/x * (x1-rx) * sum
  76 c
  77       return
  78       end