| c9c6bd3b |
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 |
| 60e55aee |
47 | ARCOSH(x)=LOG(x+SQRT((x-1D0)*(x+1D0))) |
| c9c6bd3b |
48 | |
| 49 | gamma=gm |
| 60e55aee |
50 | beta=sqrt((1D0-1D0/gamma)*(1D0+1D0/gamma)) |
| c9c6bd3b |
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 |
git://git.uio.no / u / mrichter / AliRoot.git / blame