EPEMGEN/diffcross.f

   1 C**********************************************************************
   2 C*                                                                    *
   3 C* Exact calculation of the total differential e+ e- -pair production *
   4 C* in Relativistic Heavy Ion Collisions for a point particle in an    *
   5 C* external field approach.                                           *
   6 C*                                                                    *
   7 C* For details see the publication:                                   *
   8 C* "Multiple electromagnetic electron positron pair production in     *
   9 C*  relativistic heavy ion collisions"                                *
  10 C*  Adrian Alscher, Kai Hencken, Dirk Trautmann, and Gerhard Baur     *
  11 C*  Phys. Rev. A55 (1997) 396.                                        *
  12 C*                                                                    *
  13 C* Copyright (c) 1997, 1998, 2002, Adrian Alscher and Kai Hencken     *
  14 C*                                                                    *
  15 C* Permission to use, copy and distribute this software and its       *
  16 C* documentation strictly for non-commercial purposes is hereby       *
  17 C* granted without fee, provided that the above copyright notice      *
  18 C* appears in all copies, that both the copyright notice and this     *
  19 C* permission notice appear in the supporting documentation and that  *
  20 C* the use of this program is acknowledged in scientific              *
  21 C* publications (see reference above). The authors make no claims     *
  22 C* about the suitability of this software for any purpose. It is      *
  23 C* provided "as is" without express or implied warranty. Any change   *
  24 C* of the code should be submitted to the authors                     *
  25 C*                                                                    *
  26 C* To use this program at LHC energies, please make sure that         *
  27 C* "double precision" variables should better be real*16              *
  28 C**********************************************************************
  29
  30 C======================================================================
  31 C
  32 C call this routine first to initialize some parameter needed in the
  33 C function.
  34 C
  35 C gm = Gamma_cm, that is, gm of each ion (~100 for RHIC ~3000 for LHC)
  36 C mass = mass of the produced particle in MeV (~0.511 for e,~100 for mu)
  37 C======================================================================
  38
  39       SUBROUTINE Initdiffcross (gm,mass)
  40       IMPLICIT NONE
  41       DOUBLE PRECISION gm,mass
  42
  43       DOUBLE PRECISION gamma,beta,m,w1xw1,w1xw2,w2xw2,wl,wy
  44       COMMON/PHYSPARAM/gamma,beta,m,w1xw1,w1xw2,w2xw2,wl,wy
  45
  46       DOUBLE PRECISION ARCOSH,x
  47       ARCOSH(x)=LOG(x+SQRT(x**2-1D0))
  48
  49       gamma=gm
  50       beta=sqrt(1D0-1D0/gamma**2)
  51       m=mass
  52
  53       w1xw1 = 1D0/gamma**2
  54       w2xw2 = 1D0/gamma**2
  55       w1xw2 = 2D0-1D0/gamma**2
  56       wl=SQRT (w1xw1)
  57       wy=ARCOSH (gamma)
  58
  59       RETURN
  60       END
  61
  62 C ============================================================================
  63 C
  64 C Diffcross calculates the fivefold differential cross section
  65 C
  66 C dsigma/ dp+t dp-t dy+ dy- ddeltaphi
  67 C
  68 C the trivial integration over the total phi-dependence is not included
  69 C
  70 C ppvt= absolute value of the positron transverse momentum (MeV)
  71 C pmvt=  - " -                electron - " -
  72 C dphi= phi-angle between the electron and positron transverse momentum
  73 C yp = rapidity of the positron
  74 C ym = rapidity of the electron
  75 C      defined to make E = sqrt(pt^2 + m^2) cosh (y)
  76 C      and             Pz= sqrt(pt^2 + m^2) sinh (y)
  77 C
  78 C dsigma = differential cross section in kbarn/MeV^4/ (Zalpha)**4
  79 C
  80 C to get the total cross section, you have to integrate over
  81 C
  82 C Integral dsigma dyp dym ddphi 2 pi ppt dppt pmt dpmt
  83 C
  84 C (see also source for sigma.f)
  85 C
  86 C======================================================================
  87
  88       SUBROUTINE Diffcross(ppvt,yp,pmvt,ym,dphi,dsigma)
  89
  90       IMPLICIT NONE
  91
  92       DOUBLE PRECISION ppvt,pmvt,dphi,dsigma,yp,ym
  93       DOUBLE PRECISION ppvl,pmvl
  94
  95       DOUBLE PRECISION k1(2)
  96       DOUBLE PRECISION kd(2)
  97       DOUBLE PRECISION kx(2)
  98       DOUBLE PRECISION pmt(2)
  99       DOUBLE PRECISION ppt(2)
 100       DOUBLE PRECISION qb
 101       DOUBLE PRECISION mk1(2),mk1d(2),mk1x(2)
 102       DOUBLE PRECISION mkd(2),mkd1(2)
 103       DOUBLE PRECISION mkx(2),mkx1(2)
 104
 105       DOUBLE PRECISION Iz2,Id1,Id2,Id3,Iv1,Iv2
 106
 107       DOUBLE PRECISION PI
 108       PARAMETER (PI=3.141592653589793238462643D0)
 109
 110       DOUBLE PRECISION gamma,beta,m,w1xw1,w1xw2,w2xw2,wl,wy
 111       COMMON/PHYSPARAM/gamma,beta,m,w1xw1,w1xw2,w2xw2,wl,wy
 112
 113       DOUBLE PRECISION m0,m1,md,mx
 114       DOUBLE PRECISION N1,N2,N3,N4,N5
 115       DOUBLE PRECISION N6,N7,N8,N9,N10
 116       DOUBLE PRECISION N11,N12,N13,N14,N15
 117       DOUBLE PRECISION N16,N17,N18,NT
 118
 119       DOUBLE PRECISION pmlxpml,pmlxppl,pmlxql
 120       DOUBLE PRECISION pmtxpmt,pmtxppt,pptxppt
 121       DOUBLE PRECISION pplxppl,pplxql,qlxql
 122       DOUBLE PRECISION w1xpml,w1xppl,w1xql
 123       DOUBLE PRECISION w2xpml,w2xppl,w2xql
 124
 125       DOUBLE PRECISION r1,rd,rx
 126
 127       INTEGER setzero
 128       COMMON /SETZPARAM/ setzero
 129       INTEGER badcount
 130       COMMON/badpar/badcount
 131
 132       setzero=0
 133
 134       ppt(1)=ppvt
 135       ppt(2)=0D0
 136       ppvl=SQRT(ppvt**2+m**2)
 137
 138       pmt(1)=pmvt*cos(dphi)
 139       pmt(2)=pmvt*sin(dphi)
 140       pmvl=SQRT(pmvt**2+m**2)
 141
 142       pmlxpml = pmvl**2
 143       pplxppl = ppvl**2
 144       pmlxppl = pmvl*ppvl*COSH(yp-ym)
 145       w2xpml = wl*pmvl*COSH(ym+wy)
 146       w2xppl = wl*ppvl*COSH(yp+wy)
 147       w1xpml = wl*pmvl*COSH(ym-wy)
 148       w1xppl = wl*ppvl*COSH(yp-wy)
 149       w1xql=0D0
 150       w2xql=w2xpml+w2xppl
 151       qb=gamma*(w2xppl+w2xpml)/SINH(2D0*wy)
 152       qlxql=-qb**2
 153       pmlxql = qb*pmvl*SINH(wy-ym)
 154       pplxql = qb*ppvl*SINH(wy-yp)
 155
 156       pmtxpmt =-pmvt**2
 157       pmtxppt =-pmvt*ppvt*cos(dphi)
 158       pptxppt =-ppvt**2
 159
 160 C======================================================================
 161 C
 162 C Definition of propagatorterms
 163 C
 164       m0 = -qlxql
 165
 166       k1(1) = -ppt(1) -pmt(1)
 167       k1(2) = -pmt(2)
 168       m1=-qlxql-pplxppl-pmlxpml+2D0*(pplxql+pmlxql-pmlxppl)
 169
 170       kd(1)= -pmt(1)
 171       kd(2)= -pmt(2)
 172       md= m**2 - (qlxql+pmlxpml-2D0*pmlxql)
 173
 174       kx(1)= -ppt(1)
 175       kx(2)=0D0
 176       mx= m**2 - (qlxql+pplxppl-2D0*pplxql)
 177
 178       r1 = m1 - m0 + k1(1)*k1(1)+k1(2)*k1(2)
 179       rd = md - m0 + kd(1)*kd(1)+kd(2)*kd(2)
 180       rx = mx - m0 + kx(1)*kx(1)+kx(2)*kx(2)
 181
 182       mk1(1)= -k1(1)
 183       mk1(2)= -k1(2)
 184       mkd(1)= -kd(1)
 185       mkd(2)= -kd(2)
 186       mkx(1)= -kx(1)
 187       mkx(2)= -kx(2)
 188       mk1d(1)= k1(1)-kd(1)
 189       mk1d(2)= k1(2)-kd(2)
 190       mk1x(1)= k1(1)-kx(1)
 191       mk1x(2)= k1(2)-kx(2)
 192       mkd1(1)= kd(1)-k1(1)
 193       mkd1(2)= kd(2)-k1(2)
 194       mkx1(1)= kx(1)-k1(1)
 195       mkx1(2)= kx(2)-k1(2)
 196
 197 C
 198 C calculate all different integrals
 199 C
 200
 201       N1 = Id1(mkd,mk1d,md,m0,m1) * (  - 2*w1xw1*w2xw2 )
 202
 203       N2 = Id1(mkx,mk1x,mx,m0,m1) * (  - 2*w1xw1*w2xw2 )
 204
 205       N3 = Iv1(mk1,mkd1,mkx1,m1,m0,md,mx) * ( 4*w1xw1*w2xw2*
 206      &  pmlxpml + 8*w1xw1*w2xw2*pmlxppl - 8*w1xw1*w2xw2*pmlxql + 4*
 207      &  w1xw1*w2xw2*pmtxpmt + 8*w1xw1*w2xw2*pmtxppt + 4*w1xw1*w2xw2*
 208      &  pplxppl - 8*w1xw1*w2xw2*pplxql + 4*w1xw1*w2xw2*pptxppt + 4*
 209      &  w1xw1*w2xw2*rd + 4*w1xw1*w2xw2*rx - 16*w1xw1*w2xpml*w2xppl +
 210      &  8*w1xw1*w2xpml*w2xql + 8*w1xw1*w2xppl*w2xql - 8*w1xw2**2*
 211      &  pmlxpml - 16*w1xw2**2*pmlxppl + 16*w1xw2**2*pmlxql - 8*
 212      &  w1xw2**2*pmtxpmt - 16*w1xw2**2*pmtxppt - 8*w1xw2**2*pplxppl
 213      &  + 16*w1xw2**2*pplxql - 8*w1xw2**2*pptxppt - 8*w1xw2**2*rd -
 214      &  8*w1xw2**2*rx + 16*w1xw2*w1xpml*w2xpml + 16*w1xw2*w1xpml*
 215      &  w2xppl - 16*w1xw2*w1xpml*w2xql + 16*w1xw2*w1xppl*w2xpml + 16*
 216      &  w1xw2*w1xppl*w2xppl - 16*w1xw2*w1xppl*w2xql - 16*w1xw2*w1xql*
 217      &  w2xpml - 16*w1xw2*w1xql*w2xppl - 8*w1xpml**2*w2xw2 + 8*w1xpml
 218      &  *w1xql*w2xw2 - 8*w1xppl**2*w2xw2 + 8*w1xppl*w1xql*w2xw2 )
 219
 220       N4 = Id1(mk1,mkd1,m1,m0,md) * (  - 2*w1xw1*w2xw2 + 8*
 221      &  w1xw2**2 )
 222
 223       N5 = Id2(mk1d,mkd,md,m1,m0) * ( 4*w1xw1*w2xw2*pmlxppl + 4*
 224      &  w1xw1*w2xw2*pmtxppt - 4*w1xw1*w2xw2*pplxql + 4*w1xw1*w2xw2*
 225      &  m**2 + 2*w1xw1*w2xw2*r1 - 2*w1xw1*w2xw2*rd - 8*w1xw1*w2xpml*
 226      &  w2xppl + 8*w1xw1*w2xppl*w2xql )
 227
 228       N6 = Id1(mk1,mkx1,m1,m0,mx) * (  - 2*w1xw1*w2xw2 + 8*
 229      &  w1xw2**2 )
 230
 231       N7 = Id2(mk1x,mkx,mx,m1,m0) * ( 4*w1xw1*w2xw2*pmlxppl - 4*
 232      &  w1xw1*w2xw2*pmlxql + 4*w1xw1*w2xw2*pmtxppt + 4*w1xw1*w2xw2*
 233      &  m**2 + 2*w1xw1*w2xw2*r1 - 2*w1xw1*w2xw2*rx - 8*w1xw1*w2xpml*
 234      &  w2xppl + 8*w1xw1*w2xpml*w2xql )
 235
 236       N8 = Id1(k1,kd,m0,m1,md) * ( 2*w1xw1*w2xw2 )
 237
 238       N9 = Id2(kd,k1,m0,md,m1) * (  - 4*w1xw1*w2xw2*pmlxpml
 239      &  + 4*w1xw1*w2xw2*pmlxql - 4*w1xw1*w2xw2*pmtxpmt + 4*w1xw1*
 240      &  w2xw2*m**2 - 2*w1xw1*w2xw2*rd + 8*w1xpml**2*w2xw2 - 8*w1xpml*
 241      &  w1xql*w2xw2 )
 242
 243       N10 = Id1(k1,kx,m0,m1,mx) * ( 2*w1xw1*w2xw2 )
 244
 245       N11 = Id2(kx,k1,m0,mx,m1) * (  - 4*w1xw1*w2xw2*pplxppl
 246      &  + 4*w1xw1*w2xw2*pplxql - 4*w1xw1*w2xw2*pptxppt + 4*w1xw1*
 247      &  w2xw2*m**2 - 2*w1xw1*w2xw2*rx + 8*w1xppl**2*w2xw2 - 8*w1xppl*
 248      &  w1xql*w2xw2 )
 249
 250       N12 = Iv2(k1,kd,kx,m0,m1,md,mx) * ( 8*w1xw1*w2xw2*
 251      &  pmlxpml*pplxppl - 8*w1xw1*w2xw2*pmlxpml*pplxql + 8*w1xw1*
 252      &  w2xw2*pmlxpml*pptxppt - 8*w1xw1*w2xw2*pmlxpml*m**2 + 4*w1xw1*
 253      &  w2xw2*pmlxpml*rx - 8*w1xw1*w2xw2*pmlxppl*qlxql - 8*w1xw1*
 254      &  w2xw2*pmlxppl*m0 - 8*w1xw1*w2xw2*pmlxql*pplxppl + 16*w1xw1
 255      &  *w2xw2*pmlxql*pplxql - 8*w1xw1*w2xw2*pmlxql*pptxppt + 8*w1xw1
 256      &  *w2xw2*pmlxql*m**2 - 8*w1xw1*w2xw2*pmlxql*rx + 8*w1xw1*w2xw2*
 257      &  pmtxpmt*pplxppl - 8*w1xw1*w2xw2*pmtxpmt*pplxql + 8*w1xw1*
 258      &  w2xw2*pmtxpmt*pptxppt - 8*w1xw1*w2xw2*pmtxpmt*m**2 + 4*w1xw1*
 259      &  w2xw2*pmtxpmt*rx - 8*w1xw1*w2xw2*pmtxppt*qlxql - 8*w1xw1*
 260      &  w2xw2*pmtxppt*m0 - 8*w1xw1*w2xw2*pplxppl*m**2 + 4*w1xw1*
 261      &  w2xw2*pplxppl*rd + 8*w1xw1*w2xw2*pplxql*m**2 - 8*w1xw1*w2xw2*
 262      &  pplxql*rd - 8*w1xw1*w2xw2*pptxppt*m**2 + 4*w1xw1*w2xw2*
 263      &  pptxppt*rd - 8*w1xw1*w2xw2*qlxql*m**2 + 8*w1xw1*w2xw2*m**4 -
 264      &  8*w1xw1*w2xw2*m**2*m0 - 4*w1xw1*w2xw2*m**2*rd - 4*w1xw1*
 265      &  w2xw2*m**2*rx + 4*w1xw1*w2xw2*rd*rx + 16*w1xw1*w2xpml*w2xppl*
 266      &  qlxql +  16*w1xw1*w2xpml*
 267      &  w2xppl*m0 - 16*w1xw1*w2xpml*w2xql*pplxql + 8*w1xw1*w2xpml*
 268      &  w2xql*rx - 16*w1xw1*w2xppl*w2xql*pmlxql + 8*w1xw1*w2xppl*
 269      &  w2xql*rd + 16*w1xw1*w2xql**2*pmlxppl + 16*w1xw1*w2xql**2*
 270      &  pmtxppt + 16*w1xw1*w2xql**2*m**2 - 16*w1xw2**2*pmlxpml*
 271      &  pplxppl + 16*w1xw2**2*pmlxpml*pplxql - 16*w1xw2**2*pmlxpml*
 272      &  pptxppt + 16*w1xw2**2*pmlxpml*m**2 - 8*w1xw2**2*pmlxpml*rx +
 273      &  16*w1xw2**2*pmlxppl*qlxql + 16*w1xw2**2*pmlxppl*m0 + 16*
 274      &  w1xw2**2*pmlxql*pplxppl - 32*w1xw2**2*pmlxql*pplxql + 16*
 275      &  w1xw2**2*pmlxql*pptxppt - 16*w1xw2**2*pmlxql*m**2 + 16*
 276      &  w1xw2**2*pmlxql*rx - 16*w1xw2**2*pmtxpmt*pplxppl + 16*
 277      &  w1xw2**2*pmtxpmt*pplxql - 16*w1xw2**2*pmtxpmt*pptxppt + 16*
 278      &  w1xw2**2*pmtxpmt*m**2 - 8*w1xw2**2*pmtxpmt*rx + 16*w1xw2**2*
 279      &  pmtxppt*qlxql + 16*w1xw2**2*pmtxppt*m0 + 16*w1xw2**2*
 280      &  pplxppl*m**2 - 8*w1xw2**2*pplxppl*rd - 16*w1xw2**2*pplxql*
 281      &  m**2 + 16*w1xw2**2*pplxql*rd + 16*w1xw2**2*pptxppt*m**2 - 8*
 282      &  w1xw2**2*pptxppt*rd + 16*w1xw2**2*qlxql
 283      &  *m**2 - 16*w1xw2**2*m**4 + 16*w1xw2**2*m**2*m0 + 8*
 284      &  w1xw2**2*m**2*rd + 8*w1xw2**2*m**2*rx - 8*w1xw2**2*rd*rx + 32
 285      &  *w1xw2*w1xpml*w2xpml*pplxppl - 32*w1xw2*w1xpml*w2xpml*pplxql
 286      &  + 32*w1xw2*w1xpml*w2xpml*pptxppt - 32*w1xw2*w1xpml*w2xpml*
 287      &  m**2 + 16*w1xw2*w1xpml*w2xpml*rx - 16*w1xw2*w1xpml*w2xppl*
 288      &  qlxql - 16*w1xw2*w1xpml*w2xppl*m0 - 16*w1xw2*w1xpml*w2xql*
 289      &  pplxppl + 32*w1xw2*w1xpml*w2xql*pplxql - 16*w1xw2*w1xpml*
 290      &  w2xql*pptxppt + 16*w1xw2*w1xpml*w2xql*m**2 - 16*w1xw2*w1xpml*
 291      &  w2xql*rx - 16*w1xw2*w1xppl*w2xpml*qlxql - 16*w1xw2*w1xppl*
 292      &  w2xpml*m0 + 32*w1xw2*w1xppl*w2xppl*pmlxpml - 32*w1xw2*
 293      &  w1xppl*w2xppl*pmlxql + 32*w1xw2*w1xppl*w2xppl*pmtxpmt - 32*
 294      &  w1xw2*w1xppl*w2xppl*m**2 + 16*w1xw2*w1xppl*w2xppl*rd - 16*
 295      &  w1xw2*w1xppl*w2xql*pmlxpml + 32*w1xw2*w1xppl*w2xql*pmlxql -
 296      &  16*w1xw2*w1xppl*w2xql*pmtxpmt + 16*w1xw2*w1xppl*w2xql*m**2 -
 297      &  16*w1xw2*w1xppl*w2xql*rd - 16*w1xw2*w1xql*w2xpml*pplxppl + 32
 298      &  *w1xw2*w1xql*w2xpml*pplxql - 16*w1xw2*w1xql
 299      &  *w2xpml*pptxppt + 16*w1xw2*w1xql*w2xpml*m**2 - 16*w1xw2*w1xql
 300      &  *w2xpml*rx - 16*w1xw2*w1xql*w2xppl*pmlxpml + 32*w1xw2*w1xql*
 301      &  w2xppl*pmlxql - 16*w1xw2*w1xql*w2xppl*pmtxpmt + 16*w1xw2*
 302      &  w1xql*w2xppl*m**2 - 16*w1xw2*w1xql*w2xppl*rd - 32*w1xw2*w1xql
 303      &  *w2xql*pmlxppl - 32*w1xw2*w1xql*w2xql*pmtxppt - 32*w1xw2*
 304      &  w1xql*w2xql*m**2 - 16*w1xpml**2*w2xw2*pplxppl + 16*w1xpml**2*
 305      &  w2xw2*pplxql - 16*w1xpml**2*w2xw2*pptxppt + 16*w1xpml**2*
 306      &  w2xw2*m**2 - 8*w1xpml**2*w2xw2*rx + 32*w1xpml*w1xppl*w2xw2*
 307      &  pmlxppl - 16*w1xpml*w1xppl*w2xw2*pmlxql + 32*w1xpml*w1xppl*
 308      &  w2xw2*pmtxppt - 16*w1xpml*w1xppl*w2xw2*pplxql + 16*w1xpml*
 309      &  w1xppl*w2xw2*qlxql + 32*w1xpml*w1xppl*w2xw2*m**2 + 16*w1xpml*
 310      &  w1xppl*w2xw2*m0 + 8*w1xpml*w1xppl*w2xw2*rd + 8*w1xpml*
 311      &  w1xppl*w2xw2*rx - 64*w1xpml*w1xppl*w2xpml*w2xppl + 32*w1xpml*
 312      &  w1xppl*w2xpml*w2xql + 32*w1xpml*w1xppl*w2xppl*w2xql - 32*
 313      &  w1xpml*w1xppl*w2xql**2 - 16*w1xpml*w1xql*w2xw2*pmlxppl - 16*
 314      &  w1xpml*w1xql*w2xw2*pmtxppt + 16*w1xpml*w1xql*
 315      &  w2xw2*pplxppl - 16*w1xpml*w1xql*w2xw2*pplxql + 16*w1xpml*
 316      &  w1xql*w2xw2*pptxppt - 32*w1xpml*w1xql*w2xw2*m**2 + 8*w1xpml*
 317      &  w1xql*w2xw2*rx + 32*w1xpml*w1xql*w2xpml*w2xppl - 16*w1xppl**2
 318      &  *w2xw2*pmlxpml + 16*w1xppl**2*w2xw2*pmlxql - 16*w1xppl**2*
 319      &  w2xw2*pmtxpmt + 16*w1xppl**2*w2xw2*m**2 - 8*w1xppl**2*w2xw2*
 320      &  rd + 16*w1xppl*w1xql*w2xw2*pmlxpml - 16*w1xppl*w1xql*w2xw2*
 321      &  pmlxppl - 16*w1xppl*w1xql*w2xw2*pmlxql + 16*w1xppl*w1xql*
 322      &  w2xw2*pmtxpmt - 16*w1xppl*w1xql*w2xw2*pmtxppt - 32*w1xppl*
 323      &  w1xql*w2xw2*m**2 + 8*w1xppl*w1xql*w2xw2*rd + 32*w1xppl*w1xql*
 324      &  w2xpml*w2xppl + 16*w1xql**2*w2xw2*pmlxppl + 16*w1xql**2*w2xw2
 325      &  *pmtxppt + 16*w1xql**2*w2xw2*m**2 - 32*w1xql**2*w2xpml*w2xppl
 326      &  )
 327
 328       N13 = Id2(k1,kd,m0,m1,md) * ( 4*w1xw1*w2xw2*pmlxql + 4*
 329      &  w1xw1*w2xw2*pplxql - 2*w1xw1*w2xw2*r1 - 8*w1xw1*w2xpml*w2xql
 330      &  - 8*w1xw1*w2xppl*w2xql + 8*w1xw2**2*pmlxpml - 16*w1xw2**2*
 331      &  pmlxql + 8*w1xw2**2*pmtxpmt - 8*w1xw2**2*m**2 + 8*w1xw2**2*rd
 332      &  - 16*w1xw2*w1xpml*w2xpml + 16*w1xw2*w1xpml*w2xppl + 16*w1xw2
 333      &  *w1xpml*w2xql + 16*w1xw2*w1xql*w2xpml - 16*w1xpml*w1xppl*
 334      &  w2xw2 )
 335
 336       N14 = Id3(k1,kd,m0,m1,md) * ( 4*w1xw1*w2xw2*pmlxpml*
 337      &  pmlxppl + 4*w1xw1*w2xw2*pmlxpml*pmtxppt - 8*w1xw1*w2xw2*
 338      &  pmlxpml*pplxql + 4*w1xw1*w2xw2*pmlxpml*m**2 + 4*w1xw1*w2xw2*
 339      &  pmlxpml*r1 - 4*w1xw1*w2xw2*pmlxpml*rd + 4*w1xw1*w2xw2*pmlxppl
 340      &  *pmtxpmt - 4*w1xw1*w2xw2*pmlxppl*qlxql - 4*w1xw1*w2xw2*
 341      &  pmlxppl*m**2 - 4*w1xw1*w2xw2*pmlxppl*m0 + 8*w1xw1*w2xw2*
 342      &  pmlxql*pplxql - 4*w1xw1*w2xw2*pmlxql*r1 + 4*w1xw1*w2xw2*
 343      &  pmlxql*rd + 4*w1xw1*w2xw2*pmtxpmt*pmtxppt - 8*w1xw1*w2xw2*
 344      &  pmtxpmt*pplxql + 4*w1xw1*w2xw2*pmtxpmt*m**2 + 4*w1xw1*w2xw2*
 345      &  pmtxpmt*r1 - 4*w1xw1*w2xw2*pmtxpmt*rd - 4*w1xw1*w2xw2*pmtxppt
 346      &  *qlxql - 4*w1xw1*w2xw2*pmtxppt*m**2 - 4*w1xw1*w2xw2*pmtxppt*
 347      &  m0 + 8*w1xw1*w2xw2*pplxql*m**2 - 4*w1xw1*w2xw2*pplxql*rd
 348      &  - 4*w1xw1*w2xw2*qlxql*m**2 - 4*w1xw1*w2xw2*m**4 - 4*w1xw1*
 349      &  w2xw2*m**2*m0 - 4*w1xw1*w2xw2*m**2*r1 + 4*w1xw1*w2xw2*m**2
 350      &  *rd + 2*w1xw1*w2xw2*r1*rd - 2*w1xw1*w2xw2*rd**2 - 8*w1xw1*
 351      &  w2xpml*w2xppl*pmlxpml - 8*w1xw1*w2xpml*w2xppl*pmtxpmt + 8*
 352      &  w1xw1*w2xpml*w2xppl*qlxql + 8*w1xw1*w2xpml*w2xppl*m**2
 353      &  + 8*w1xw1*w2xpml*w2xppl*m0 + 16*w1xw1*w2xppl*w2xql*
 354      &  pmlxpml - 16*w1xw1*w2xppl*w2xql*pmlxql + 16*w1xw1*w2xppl*
 355      &  w2xql*pmtxpmt - 16*w1xw1*w2xppl*w2xql*m**2 + 8*w1xw1*w2xppl*
 356      &  w2xql*rd - 16*w1xw2*w1xpml*w2xppl*pmlxpml + 32*w1xw2*w1xpml*
 357      &  w2xppl*pmlxql - 16*w1xw2*w1xpml*w2xppl*pmtxpmt - 16*w1xw2*
 358      &  w1xpml*w2xppl*qlxql + 16*w1xw2*w1xpml*w2xppl*m**2 - 16*w1xw2*
 359      &  w1xpml*w2xppl*m0 - 16*w1xw2*w1xpml*w2xppl*rd - 16*
 360      &  w1xpml**2*w2xw2*pmlxppl - 16*w1xpml**2*w2xw2*pmtxppt + 16*
 361      &  w1xpml**2*w2xw2*pplxql - 16*w1xpml**2*w2xw2*m**2 - 8*
 362      &  w1xpml**2*w2xw2*r1 + 8*w1xpml**2*w2xw2*rd + 32*w1xpml**2*
 363      &  w2xpml*w2xppl - 32*w1xpml**2*w2xppl*w2xql + 8*w1xpml*w1xppl*
 364      &  w2xw2*pmlxpml - 16*w1xpml*w1xppl*w2xw2*pmlxql + 8*w1xpml*
 365      &  w1xppl*w2xw2*pmtxpmt + 8*w1xpml*w1xppl*w2xw2*qlxql - 8*w1xpml
 366      &  *w1xppl*w2xw2*m**2 + 8*w1xpml*w1xppl*w2xw2*m0 + 8*w1xpml*
 367      &  w1xppl*w2xw2*rd + 16*w1xpml*w1xql*w2xw2*pmlxppl + 16*w1xpml*
 368      &  w1xql*w2xw2*pmtxppt - 16*w1xpml*w1xql*w2xw2*
 369      &  pplxql + 16*w1xpml*w1xql*w2xw2*m**2 + 8*w1xpml*w1xql*w2xw2*r1
 370      &  - 8*w1xpml*w1xql*w2xw2*rd - 32*w1xpml*w1xql*w2xpml*w2xppl +
 371      &  32*w1xpml*w1xql*w2xppl*w2xql )
 372
 373       N15 = Id2(k1,kx,m0,m1,mx) * ( 4*w1xw1*w2xw2*pmlxql + 4*
 374      &  w1xw1*w2xw2*pplxql - 2*w1xw1*w2xw2*r1 - 8*w1xw1*w2xpml*w2xql
 375      &  - 8*w1xw1*w2xppl*w2xql + 8*w1xw2**2*pplxppl - 16*w1xw2**2*
 376      &  pplxql + 8*w1xw2**2*pptxppt - 8*w1xw2**2*m**2 + 8*w1xw2**2*rx
 377      &  + 16*w1xw2*w1xppl*w2xpml - 16*w1xw2*w1xppl*w2xppl + 16*w1xw2
 378      &  *w1xppl*w2xql + 16*w1xw2*w1xql*w2xppl - 16*w1xpml*w1xppl*
 379      &  w2xw2 )
 380
 381       N16 = Id3(k1,kx,m0,m1,mx) * ( 4*w1xw1*w2xw2*pmlxppl*
 382      &  pplxppl + 4*w1xw1*w2xw2*pmlxppl*pptxppt - 4*w1xw1*w2xw2*
 383      &  pmlxppl*qlxql - 4*w1xw1*w2xw2*pmlxppl*m**2 - 4*w1xw1*w2xw2*
 384      &  pmlxppl*m0 - 8*w1xw1*w2xw2*pmlxql*pplxppl + 8*w1xw1*w2xw2*
 385      &  pmlxql*pplxql - 8*w1xw1*w2xw2*pmlxql*pptxppt + 8*w1xw1*w2xw2*
 386      &  pmlxql*m**2 - 4*w1xw1*w2xw2*pmlxql*rx + 4*w1xw1*w2xw2*pmtxppt
 387      &  *pplxppl + 4*w1xw1*w2xw2*pmtxppt*pptxppt - 4*w1xw1*w2xw2*
 388      &  pmtxppt*qlxql - 4*w1xw1*w2xw2*pmtxppt*m**2 - 4*w1xw1*w2xw2*
 389      &  pmtxppt*m0 + 4*w1xw1*w2xw2*pplxppl*m**2 + 4*w1xw1*w2xw2*
 390      &  pplxppl*r1 - 4*w1xw1*w2xw2*pplxppl*rx - 4*w1xw1*w2xw2*pplxql*
 391      &  r1 + 4*w1xw1*w2xw2*pplxql*rx + 4*w1xw1*w2xw2*pptxppt*m**2 + 4
 392      &  *w1xw1*w2xw2*pptxppt*r1 - 4*w1xw1*w2xw2*pptxppt*rx - 4*w1xw1*
 393      &  w2xw2*qlxql*m**2 - 4*w1xw1*w2xw2*m**4 - 4*w1xw1*w2xw2*m**2*
 394      &  m0 - 4*w1xw1*w2xw2*m**2*r1 + 4*w1xw1*w2xw2*m**2*rx + 2*
 395      &  w1xw1*w2xw2*r1*rx - 2*w1xw1*w2xw2*rx**2 - 8*w1xw1*w2xpml*
 396      &  w2xppl*pplxppl - 8*w1xw1*w2xpml*w2xppl*pptxppt + 8*w1xw1*
 397      &  w2xpml*w2xppl*qlxql + 8*w1xw1*w2xpml*w2xppl*m**2
 398      &  + 8*w1xw1*w2xpml*w2xppl*m0 + 16*w1xw1*w2xpml*w2xql*
 399      &  pplxppl - 16*w1xw1*w2xpml*w2xql*pplxql + 16*w1xw1*w2xpml*
 400      &  w2xql*pptxppt - 16*w1xw1*w2xpml*w2xql*m**2 + 8*w1xw1*w2xpml*
 401      &  w2xql*rx - 16*w1xw2*w1xppl*w2xpml*pplxppl + 32*w1xw2*w1xppl*
 402      &  w2xpml*pplxql - 16*w1xw2*w1xppl*w2xpml*pptxppt - 16*w1xw2*
 403      &  w1xppl*w2xpml*qlxql + 16*w1xw2*w1xppl*w2xpml*m**2 - 16*w1xw2*
 404      &  w1xppl*w2xpml*m0 - 16*w1xw2*w1xppl*w2xpml*rx + 8*w1xpml*
 405      &  w1xppl*w2xw2*pplxppl - 16*w1xpml*w1xppl*w2xw2*pplxql + 8*
 406      &  w1xpml*w1xppl*w2xw2*pptxppt + 8*w1xpml*w1xppl*w2xw2*qlxql - 8
 407      &  *w1xpml*w1xppl*w2xw2*m**2 + 8*w1xpml*w1xppl*w2xw2*m0 + 8*
 408      &  w1xpml*w1xppl*w2xw2*rx - 16*w1xppl**2*w2xw2*pmlxppl + 16*
 409      &  w1xppl**2*w2xw2*pmlxql - 16*w1xppl**2*w2xw2*pmtxppt - 16*
 410      &  w1xppl**2*w2xw2*m**2 - 8*w1xppl**2*w2xw2*r1 + 8*w1xppl**2*
 411      &  w2xw2*rx + 32*w1xppl**2*w2xpml*w2xppl - 32*w1xppl**2*w2xpml*
 412      &  w2xql + 16*w1xppl*w1xql*w2xw2*pmlxppl - 16*w1xppl*w1xql*w2xw2
 413      &  *pmlxql + 16*w1xppl*w1xql*w2xw2*
 414      &  pmtxppt + 16*w1xppl*w1xql*w2xw2*m**2 + 8*w1xppl*w1xql*w2xw2*
 415      &  r1 - 8*w1xppl*w1xql*w2xw2*rx - 32*w1xppl*w1xql*w2xpml*w2xppl
 416      &  + 32*w1xppl*w1xql*w2xpml*w2xql )
 417
 418       N17 = Iz2(k1,m0,m1) * (  - 8*w1xw2**2 )
 419
 420       N18 = Id1(mkd1,mkx1,m1,md,mx) * ( 4*w1xw1*w2xw2 - 8*w1xw2**2 )
 421
 422 C
 423 C dsigma is summ of all terms
 424 C
 425       NT=N1+N2+N3+N4+N5+N6+N7+N8+N9+N10+N11+N12+N13+N14+
 426      &  N15+N16+N17+N18
 427
 428 C
 429 C correction from w/u
 430       NT=NT*4D0/beta**2
 431 C 1/(2pi)**6 from d3p, (2*pi)**2 from F.T.
 432       NT=NT/(2*pi)**6*(2*pi)**2
 433 C from 1/2E+ 1/2E-
 434       NT=NT/4D0
 435 C transform from MeV^-2 to kbarn
 436       NT=NT*(1.9733D0)**2/10D0
 437
 438       dsigma=NT
 439       IF((setzero.EQ.1).OR.(dsigma.LT.0)) THEN
 440       dsigma=0D0
 441       badcount=badcount+1
 442       ENDIF
 443       END
 444
 445 C========================================================================
 446 C All the differential integral forms are calculated here
 447 C
 448       DOUBLE PRECISION FUNCTION Iz0(x,u,v)
 449       IMPLICIT NONE
 450
 451       DOUBLE PRECISION x(2)
 452       DOUBLE PRECISION u,v
 453       DOUBLE PRECISION s,arg
 454
 455       INTEGER setzero
 456       COMMON /SETZPARAM/ setzero
 457
 458       DOUBLE PRECISION tepxx
 459
 460       DOUBLE PRECISION PI
 461       PARAMETER (PI=3.141592653589793238462643D0)
 462
 463       tepxx=x(1)*x(1)+x(2)*x(2)
 464       s=SQRT((tepxx+u+v)**2 - 4*u*v)
 465       arg=(tepxx+u+v+s)**2/(4*u*v)
 466       IF(arg.lt.0D0) THEN
 467       setzero=1
 468       Iz0=0D0
 469       ELSE
 470       Iz0=pi * LOG((tepxx+u+v+s)**2/(4*u*v)) / s
 471       ENDIF
 472       END
 473
 474 C     ---------------------------------------------------------------
 475
 476       DOUBLE PRECISION FUNCTION Iz1(x,u,v)
 477       IMPLICIT NONE
 478       DOUBLE PRECISION Iz0
 479       EXTERNAL         Iz0
 480
 481       DOUBLE PRECISION x(2)
 482       DOUBLE PRECISION u,v
 483       DOUBLE PRECISION a,b,s
 484
 485       DOUBLE PRECISION tepxx
 486       DOUBLE PRECISION PI
 487       PARAMETER (PI=3.141592653589793238462643D0)
 488
 489       tepxx=x(1)*x(1)+x(2)*x(2)
 490
 491       s=SQRT((tepxx+u+v)**2 - 4*u*v)
 492       A= (2*(tepxx+u-v+s)) / (s**2*(tepxx+u+v+s))
 493      &  - 1/(s*u)
 494       B= -Iz0(x,u,v)/pi * (tepxx+u-v)/s**2
 495       Iz1=(-1)*pi*(A + B)
 496
 497       END
 498
 499 C     -----------------------------------------------------------
 500
 501       DOUBLE PRECISION FUNCTION Iz2(x,u,v)
 502       IMPLICIT NONE
 503
 504       DOUBLE PRECISION Iz0,Iz1
 505       EXTERNAL         Iz0,Iz1
 506
 507       DOUBLE PRECISION x(2)
 508       DOUBLE PRECISION u,v
 509       DOUBLE PRECISION s,A,B,C,D,E
 510
 511       DOUBLE PRECISION tepxx
 512
 513       DOUBLE PRECISION PI
 514       PARAMETER (PI=3.141592653589793238462643D0)
 515
 516       tepxx=x(1)*x(1)+x(2)*x(2)
 517
 518       s=SQRT((tepxx+u+v)**2 - 4*u*v)
 519       A= (2*(tepxx-u+v-s))/(s**3*(tepxx+u+v+s))
 520       B= -2*(tepxx+u-v+s)*(
 521      &  2*(tepxx-u+v)*(tepxx+u+v+s) +
 522      &  s*(tepxx-u+v+s)  ) /
 523      &  (s**4*(tepxx+u+v+s)**2)
 524       C= (tepxx-u+v)/(u*s**3)
 525       D= (Iz1(x,v,u)/pi)*(tepxx+u-v)/(s**2)
 526       E= (Iz0(x,u,v)/pi)*(
 527      &  1/s**2 + (tepxx+u-v)*2*(tepxx-u+v) /
 528      &  (s**4) )
 529       Iz2= pi*(A + B + C +D + E)
 530
 531       END
 532
 533 C     -------------------------------------------------------------
 534
 535       DOUBLE PRECISION FUNCTION Id0(x,y,u,v,w)
 536       IMPLICIT NONE
 537       DOUBLE PRECISION Iz0,Iz1
 538       EXTERNAL         Iz0,Iz1
 539       DOUBLE PRECISION x(2),y(2),dyx(2)
 540       DOUBLE PRECISION u,v,w
 541       DOUBLE PRECISION A,AXY,B,C,D,E,RX,RY
 542       DOUBLE PRECISION tepxx,tepxy,tepyy
 543
 544       tepxx=x(1)*x(1)+x(2)*x(2)
 545       tepyy=y(1)*y(1)+y(2)*y(2)
 546       tepxy=x(1)*y(1)+x(2)*y(2)
 547
 548       rx= v-u+tepxx
 549       ry= w-u+tepyy
 550       axy= x(1)*y(2)-x(2)*y(1)
 551
 552       A= rx**2*tepyy - 2*rx*ry*tepxy + ry**2*tepxx
 553       B= 4*u*axy**2 + A
 554       C= 2*axy**2 + (rx+ry)*tepxy - rx*tepyy - ry*tepxx
 555       D= rx*tepyy - ry*tepxy
 556       E= ry*tepxx - rx*tepxy
 557       dyx(1)= y(1)-x(1)
 558       dyx(2)= y(2)-x(2)
 559       Id0=1/B*(C*Iz0(dyx,v,w) + D*Iz0(y,u,w) + E*Iz0(x,u,v))
 560
 561       END
 562
 563 C     -------------------------------------------------------------
 564
 565       DOUBLE PRECISION FUNCTION Id1(x,y,u,v,w)
 566       IMPLICIT NONE
 567       DOUBLE PRECISION Id0,Iz0,Iz1
 568       EXTERNAL         Id0,Iz0,Iz1
 569
 570       DOUBLE PRECISION x(2),y(2),dyx(2)
 571       DOUBLE PRECISION u,v,w
 572
 573       DOUBLE PRECISION A,AXY,B,C,D,E,F,G,H,RX,RY
 574       DOUBLE PRECISION tepxx,tepxy,tepyy
 575
 576       tepxx=x(1)*x(1)+x(2)*x(2)
 577       tepxy=x(1)*y(1)+x(2)*y(2)
 578       tepyy=y(1)*y(1)+y(2)*y(2)
 579
 580       rx= v-u+tepxx
 581       ry= w-u+tepyy
 582       axy= x(1)*y(2)-x(2)*y(1)
 583       A= rx**2*tepyy - 2*rx*ry*tepxy + ry**2*tepxx
 584       B= 4*u*axy**2 + A
 585       C= -4*axy**2 -2*(-rx*tepyy+(rx+ry)*tepxy-ry*tepxx)
 586       D= -2*tepxy + tepxx + tepyy
 587       E= -tepyy + tepxy
 588       F= rx*tepyy - ry*tepxy
 589       G= -tepxx + tepxy
 590       H= ry*tepxx - rx*tepxy
 591
 592       dyx(1)= y(1)-x(1)
 593       dyx(2)= y(2)-x(2)
 594
 595       Id1= -( (C/B)*Id0(x,y,u,v,w) + (1/B)*(
 596      &  D*Iz0(dyx,v,w) +
 597      &  E*Iz0(y,u,w) - F*Iz1(y,u,w) +
 598      &  G*Iz0(x,u,v) - H*Iz1(x,u,v)   ) )
 599
 600       END
 601
 602 C     --------------------------------------------------------
 603
 604       DOUBLE PRECISION FUNCTION Id2(x,y,u,v,w)
 605       IMPLICIT NONE
 606       DOUBLE PRECISION Id1,Id0,Iz0,Iz1,Iz2
 607       EXTERNAL         Id1,Id0,Iz0,Iz1,Iz2
 608
 609       DOUBLE PRECISION x(2),y(2),mx(2),dyx(2)
 610       DOUBLE PRECISION u,v,w
 611
 612       DOUBLE PRECISION A0,A1,A2,A3,A4,AXY
 613       DOUBLE PRECISION B,C1,C2,C3,D1,D2,D3
 614       DOUBLE PRECISION E1,E2,E3,F1,F2,F3
 615       DOUBLE PRECISION G1,G2,G3,RX,RY
 616       DOUBLE PRECISION tepxx,tepxy,tepyy
 617
 618       tepxx=x(1)*x(1)+x(2)*x(2)
 619       tepxy=x(1)*y(1)+x(2)*y(2)
 620       tepyy=y(1)*y(1)+y(2)*y(2)
 621
 622       mx(1)=-x(1)
 623       mx(2)=-x(2)
 624       dyx(1)= y(1)-x(1)
 625       dyx(2)= y(2)-x(2)
 626       rx= v-u+tepxx
 627       ry= w-u+tepyy
 628       axy= x(1)*y(2)-x(2)*y(1)
 629
 630       A0= rx**2*tepyy - 2*rx*ry*tepxy + ry**2*tepxx
 631       B= 4*u*axy**2 + A0
 632       A1= 2*(tepyy-tepxy)*B
 633       A2= (-4*axy**2 +2*(rx*tepyy-(rx+ry)*tepxy+ry*tepxx))*
 634      &  2*2*(rx*tepyy-ry*tepxy)
 635       A3= -4*axy**2 +2*(rx*tepyy-(rx+ry)*tepxy+ry*tepxx)
 636       A4= -2*(rx*tepyy-ry*tepxy)
 637       C1= tepxy-tepyy
 638       D1= 2*axy**2+(rx+ry)*tepxy-rx*tepyy-ry*tepxx
 639       E1= tepyy
 640       F1= -tepxy
 641       G1= ry*tepxx-rx*tepxy
 642       C2= -2*tepxy+tepxx+tepyy
 643       D2= -tepyy + tepxy
 644       E2= rx*tepyy - ry*tepxy
 645       F2= -tepxx + tepxy
 646       G2= ry*tepxx - rx*tepxy
 647       C3= -2*tepxy+tepxx+tepyy
 648       D3= -tepyy
 649       E3= -tepxx+tepxy
 650       F3= tepxy
 651       G3= -(ry*tepxx-rx*tepxy)
 652
 653       Id2= ((A1-A2)/B**2)*Id0(x,y,u,v,w) +
 654      &  (A3/B**2)*(
 655      &  C1*Iz0(dyx,v,w) - D1*Iz1(dyx,v,w) +
 656      &  E1*Iz0(y,u,w) +
 657      &  F1*Iz0(x,u,v) - G1*Iz1(mx,v,u) )  +
 658      &  (A4/B**2)*(
 659      &  C2*Iz0(dyx,v,w) +
 660      &  D2*Iz0(y,u,w) - E2*Iz1(y,u,w) +
 661      &  F2*Iz0(x,u,v) - G2*Iz1(x,u,v)  ) +
 662      &  (1/B)*(
 663      &  - C3*Iz1(dyx,v,w)
 664      &  + D3*Iz1(y,u,w)
 665      &  - E3*Iz1(mx,v,u) + F3*Iz1(x,u,v) - G3*Iz2(mx,v,u) )
 666
 667       END
 668
 669 C     --------------------------------------------------------
 670
 671       DOUBLE PRECISION FUNCTION Id3(x,y,u,v,w)
 672       IMPLICIT NONE
 673       DOUBLE PRECISION Id2,Id1,Id0,Iz0,Iz1,Iz2
 674       EXTERNAL         Id2,Id1,Id0,Iz0,Iz1,Iz2
 675
 676       DOUBLE PRECISION x(2),y(2),dyx(2),mdyx(2),mx(2),my(2)
 677       DOUBLE PRECISION u,v,w
 678
 679       DOUBLE PRECISION A0,A1,A2,A3A,A3B,A3C,A3D
 680       DOUBLE PRECISION A4,A5A,A5B,A5C,A5D,A5E,A5F
 681       DOUBLE PRECISION A6,A8A,A8B,A8C,A9,A10,AXY
 682       DOUBLE PRECISION B,C2,C5,C6,C7,C8,C9,C10,C11
 683       DOUBLE PRECISION D2,D5,D6,D7,D8,D9,D10,D11
 684       DOUBLE PRECISION E2,E5,E6,E7,E8,E9,E10,E11
 685       DOUBLE PRECISION F2,F5,F6,F7,F8,F9,F10
 686       DOUBLE PRECISION G2,G5,G6,G7,G8,G9,G10
 687       DOUBLE PRECISION RX,RY
 688       DOUBLE PRECISION tepxx,tepxy,tepyy
 689
 690       tepxx=x(1)*x(1)+x(2)*x(2)
 691       tepxy=x(1)*y(1)+x(2)*y(2)
 692       tepyy=y(1)*y(1)+y(2)*y(2)
 693
 694       dyx(1)= y(1)-x(1)
 695       dyx(2)= y(2)-x(2)
 696       mdyx(1)=-dyx(1)
 697       mdyx(2)=-dyx(2)
 698       mx(1)=-x(1)
 699       mx(2)=-x(2)
 700       my(1)=-y(1)
 701       my(2)=-y(2)
 702
 703       rx= v-u+tepxx
 704       ry= w-u+tepyy
 705       axy= x(1)*y(2)-x(2)*y(1)
 706
 707       A0= rx**2*tepyy - 2*rx*ry*tepxy + ry**2*tepxx
 708       B= 4*u*axy**2 + A0
 709
 710       A1= 2*(tepyy-tepxy)*(-2)*2*(-rx*tepxy+ry*tepxx)
 711       A2= 2*(tepyy-tepxy)
 712       C2= tepxy-tepxx
 713       D2= 2*axy**2+(rx+ry)*tepxy-rx*tepyy-ry*tepxx
 714       E2= -tepxy
 715       F2= rx*tepyy-ry*tepxy
 716       G2= tepxx
 717
 718       A3a= 2*(-tepxy+tepxx) * 2*(rx*tepyy-ry*tepxy)
 719       A3b= (-4*axy**2 +2*(rx*tepyy-(rx+ry)*tepxy+ry*tepxx))
 720      +  * 2*(-tepxy)
 721       A3c= (-4*axy**2 +2*(rx*tepyy-(rx+ry)*tepxy+ry*tepxx))
 722      &  * 2*(rx*tepyy-ry*tepxy)
 723       A3d= 3*2*(-rx*tepxy+ry*tepxx)
 724
 725       A4= A3c
 726       A5a= 2*(-tepxy+tepxx)
 727       A5b= -4*axy**2 +2*(rx*tepyy-(rx+ry)*tepxy+ry*tepxx)
 728       A5c= 2*2*(-rx*tepxy+ ry*tepxx)
 729       A5d= -2*(rx*tepyy-ry*tepxy)
 730       A5e= -2*(-tepxy)
 731       A5f= -2*(rx*tepyy-ry*tepxy)*(2)*(-rx*tepxy+ry*tepxx)
 732       C5= tepxy-tepyy
 733       D5= 2*axy**2+(rx+ry)*tepxy-rx*tepyy-ry*tepxx
 734       E5= tepyy
 735       F5= -tepxy
 736       G5= ry*tepxx-rx*tepxy
 737
 738       A6= -2*(rx*tepyy- ry*tepxy)
 739       C6= tepxy-tepxx
 740       D6= 2*axy**2+(rx+ry)*tepxy-rx*tepyy-ry*tepxx
 741       E6= -tepxy
 742       F6= rx*tepyy-ry*tepxy
 743       G6= tepxx
 744
 745       C7= tepxy-tepyy
 746       D7= -(tepxy-tepxx)
 747       E7= 2*axy**2+(rx+ry)*tepxy-rx*tepyy-ry*tepxx
 748       F7= tepyy
 749       G7= -tepxx
 750
 751       A8a= -2*(-tepxy)
 752       A8b= -2*(rx*tepyy-ry*tepxy)
 753       A8c=  2*2*(-rx*tepxy+ry*tepxx)
 754       C8= -2*tepxy+tepxx+tepyy
 755       D8= -tepyy + tepxy
 756       E8= rx*tepyy - ry*tepxy
 757       F8= -tepxx + tepxy
 758       G8= ry*tepxx - rx*tepxy
 759
 760       A9= -2*(rx*tepyy-ry*tepxy)
 761       C9= -2*tepxy+tepxx+tepyy
 762       D9= -tepyy+tepxy
 763       E9= tepxy
 764       F9= -(rx*tepyy-ry*tepxy)
 765       G9= -tepxx
 766
 767       A10= -2*(-rx*tepxy+ry*tepxx)
 768       C10= -2*tepxy+tepxx+tepyy
 769       D10= -tepyy
 770       E10= -tepxx+tepxy
 771       F10= tepxy
 772       G10= -(ry*tepxx-rx*tepxy)
 773
 774       C11= -(-2*tepxy+tepyy+tepxx)
 775       D11= -tepyy
 776       E11= tepxx
 777
 778         Id3= -( (A1/B**2)*Id0(x,y,u,v,w) + (A2/B**2)*(
 779      &  C2*Iz0(dyx,v,w) - D2*Iz1(mdyx,w,v)+
 780      &  E2*Iz0(y,u,w) - F2*Iz1(my,w,u) +
 781      &  G2*Iz0(x,u,v) )
 782      &  -( ((A3a+A3b)*B - A3c*A3d)/B**3 * Id0(x,y,u,v,w) +
 783      &  (A4/B**3)*(
 784      &  C2*Iz0(dyx,v,w) - D2*Iz1(mdyx,w,v)+
 785      &  E2*Iz0(y,u,w) - F2*Iz1(my,w,u)+
 786      &  G2*Iz0(x,u,v) ) ) +
 787      &  (A5a*B - A5b*A5c)/B**3*( A5d*Id0(x,y,u,v,w) +
 788      &  C5*Iz0(dyx,v,w) - D5*Iz1(dyx,v,w) +
 789      &  E5*Iz0(y,u,w) +
 790      &  F5*Iz0(x,u,v) - G5*Iz1(mx,v,u) ) +
 791      &  A5b/B**2 * ( (A5e - A5f/B) * Id0(x,y,u,v,w) +
 792      &  (A6/B)*(
 793      &  C6*Iz0(dyx,v,w) - D6*Iz1(mdyx,w,v) +
 794      &  E6*Iz0(y,u,w) - F6*Iz1(my,w,u)+
 795      &  G6*Iz0(x,u,v)  )
 796      &  - C7*Iz1(mdyx,w,v) + D7*Iz1(dyx,v,w)
 797      &  + E7*Iz2(mdyx,w,v)
 798      &  - F7*Iz1(my,w,u)
 799      &  + G7*Iz1(mx,v,u) )+
 800      &  (A8a*B - A8b*A8c)/B**3*(
 801      &  C8*Iz0(dyx,v,w) +
 802      &  D8*Iz0(y,u,w) - E8*Iz1(y,u,w) +
 803      &  F8*Iz0(x,u,v) - G8*Iz1(x,u,v)  ) +
 804      &  A9/B**2*(
 805      &  - C9*Iz1(mdyx,w,v)
 806      &  - D9*Iz1(my,w,u) + E9*Iz1(y,u,w) - F9*Iz2(y,u,w)
 807      &  + G9*Iz1(x,u,v)  ) +
 808      &  A10/B**2*(
 809      &  - C10*Iz1(dyx,v,w)
 810      &  + D10*Iz1(y,u,w)
 811      &  - E10*Iz1(mx,v,u) + F10*Iz1(x,u,v) - G10*Iz2(mx,v,u) )+
 812      &  1/B*(
 813      &  - C11*Iz2(dyx,v,w)
 814      &  - D11*Iz2(y,u,w)
 815      &  + E11*Iz2(x,u,v)  ))
 816
 817       END
 818
 819 C     --------------------------------------------------------
 820
 821       DOUBLE PRECISION FUNCTION Iv0(x,y,z,u,v,w,w2)
 822       IMPLICIT NONE
 823       DOUBLE PRECISION Id0
 824       EXTERNAL         Id0
 825       DOUBLE PRECISION x(2),y(2),z(2),dyx(2),dzx(2)
 826       DOUBLE PRECISION u,v,w,w2
 827       DOUBLE PRECISION rx,ry,rz
 828       DOUBLE PRECISION axy,ayz,azx
 829       DOUBLE PRECISION A
 830
 831       rx= v-u+x(1)*x(1)+x(2)*x(2)
 832       ry= w-u+y(1)*y(1)+y(2)*y(2)
 833       rz= w2-u+z(1)*z(1)+z(2)*z(2)
 834
 835       axy= x(1)*y(2)-x(2)*y(1)
 836       ayz= y(1)*z(2)-y(2)*z(1)
 837       azx= z(1)*x(2)-z(2)*x(1)
 838
 839       dyx(1)= y(1)-x(1)
 840       dyx(2)= y(2)-x(2)
 841       dzx(1)= z(1)-x(1)
 842       dzx(2)= z(2)-x(2)
 843
 844       A= 1/(axy*rz + ayz*rx + azx*ry)
 845
 846       Iv0= A*(-(axy+ayz+azx)*Id0(dyx,dzx,v,w,w2)
 847      &  +axy*Id0(x,y,u,v,w)
 848      &  +ayz*Id0(y,z,u,w,w2)
 849      &  +azx*Id0(x,z,u,v,w2) )
 850
 851       END
 852
 853 C     --------------------------------------------------------
 854
 855       DOUBLE PRECISION FUNCTION Iv1(x,y,z,u,v,w,w2)
 856       IMPLICIT NONE
 857       DOUBLE PRECISION Iv0,Id0,Id1
 858       EXTERNAL         Iv0,Id0,Id1
 859
 860       DOUBLE PRECISION x(2),y(2),z(2),dyx(2),dzx(2)
 861       DOUBLE PRECISION u,v,w,w2
 862       DOUBLE PRECISION A,B,RX,RY,RZ
 863       DOUBLE PRECISION AXY,AYZ,AZX
 864
 865       dyx(1)= y(1)-x(1)
 866       dyx(2)= y(2)-x(2)
 867       dzx(1)= z(1)-x(1)
 868       dzx(2)= z(2)-x(2)
 869
 870       rx= v-u+x(1)*x(1)+x(2)*x(2)
 871       ry= w-u+y(1)*y(1)+y(2)*y(2)
 872       rz= w2-u+z(1)*z(1)+z(2)*z(2)
 873
 874       axy= x(1)*y(2)-x(2)*y(1)
 875       ayz= y(1)*z(2)-y(2)*z(1)
 876       azx= z(1)*x(2)-z(2)*x(1)
 877
 878       A= 1/(axy*rz + ayz*rx + azx*ry)
 879       B= axy+ayz+azx
 880       Iv1=-((B*A**2)*(-(axy+ayz+azx)*Id0(dyx,dzx,v,w,w2)
 881      &  +axy*Id0(x,y,u,v,w)
 882      &  +ayz*Id0(y,z,u,w,w2)
 883      &  +azx*Id0(x,z,u,v,w2) )+
 884      &  A*( -axy*Id1(x,y,u,v,w)
 885      &  -ayz*Id1(y,z,u,w,w2)
 886      &  -azx*Id1(x,z,u,v,w2) ))
 887
 888       END
 889
 890 C     --------------------------------------------------------
 891
 892       DOUBLE PRECISION FUNCTION Iv2(x,y,z,u,v,w,w2)
 893       IMPLICIT NONE
 894       DOUBLE PRECISION Iv1,Iv0,Id0,Id1,Id2
 895       EXTERNAL         Iv1,Iv0,Id0,Id1,Id2
 896       DOUBLE PRECISION x(2),mx(2),y(2),z(2),dyx(2)
 897       DOUBLE PRECISION dzx(2)
 898       DOUBLE PRECISION u,v,w,w2
 899       DOUBLE PRECISION axy,azx,ayz
 900       DOUBLE PRECISION rx,ry,rz,A,B
 901
 902       mx(1)=-x(1)
 903       mx(2)=-x(2)
 904       dyx(1)= y(1)-x(1)
 905       dyx(2)= y(2)-x(2)
 906       dzx(1)= z(1)-x(1)
 907       dzx(2)= z(2)-x(2)
 908
 909       rx= v-u+x(1)*x(1)+x(2)*x(2)
 910       ry= w-u+y(1)*y(1)+y(2)*y(2)
 911       rz= w2-u+z(1)*z(1)+z(2)*z(2)
 912
 913       axy= x(1)*y(2)-x(2)*y(1)
 914       azx= z(1)*x(2)-z(2)*x(1)
 915       ayz= y(1)*z(2)-y(2)*z(1)
 916
 917       A= 1/(axy*rz + ayz*rx + azx*ry)
 918       B= axy+ayz+azx
 919
 920       Iv2=(-B*2D0*ayz*A**3)*(
 921      &  -B*Id0(dyx,dzx,v,w,w2)
 922      &  +axy*Id0(x,y,u,v,w)
 923      &  +ayz*Id0(y,z,u,w,w2)
 924      &  +azx*Id0(x,z,u,v,w2) )+
 925      &  (B*A**2)*(B*Id1(dyx,dzx,v,w,w2)
 926      &  -axy*Id1(mx,dyx,v,u,w)
 927      &  -azx*Id1(mx,dzx,v,u,w2) )+
 928      &  (-ayz*A**2)*(-axy*Id1(x,y,u,v,w)
 929      &  -ayz*Id1(y,z,u,w,w2)
 930      &  -azx*Id1(x,z,u,v,w2) )+
 931      &  A*(axy*Id2(x,y,u,v,w)
 932      &  +azx*Id2(x,z,u,v,w2) )
 933
 934       END