git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Event generator for e+e- pair production (Yu.Kharlov)
[u/mrichter/AliRoot.git] / EPEMGEN / diffcross.f
CommitLineData
c9c6bd3b 1C**********************************************************************
2C* *
3C* Exact calculation of the total differential e+ e- -pair production *
4C* in Relativistic Heavy Ion Collisions for a point particle in an *
5C* external field approach. *
6C* *
7C* For details see the publication: *
8C* "Multiple electromagnetic electron positron pair production in *
9C* relativistic heavy ion collisions" *
10C* Adrian Alscher, Kai Hencken, Dirk Trautmann, and Gerhard Baur *
11C* Phys. Rev. A55 (1997) 396. *
12C* *
13C* Copyright (c) 1997, 1998, 2002, Adrian Alscher and Kai Hencken *
14C* *
15C* Permission to use, copy and distribute this software and its *
16C* documentation strictly for non-commercial purposes is hereby *
17C* granted without fee, provided that the above copyright notice *
18C* appears in all copies, that both the copyright notice and this *
19C* permission notice appear in the supporting documentation and that *
20C* the use of this program is acknowledged in scientific *
21C* publications (see reference above). The authors make no claims *
22C* about the suitability of this software for any purpose. It is *
23C* provided "as is" without express or implied warranty. Any change *
24C* of the code should be submitted to the authors *
25C* *
26C* To use this program at LHC energies, please make sure that *
27C* "double precision" variables should better be real*16 *
28C**********************************************************************
29
30C======================================================================
31C
32C call this routine first to initialize some parameter needed in the
33C function.
34C
35C gm = Gamma_cm, that is, gm of each ion (~100 for RHIC ~3000 for LHC)
36C mass = mass of the produced particle in MeV (~0.511 for e,~100 for mu)
37C======================================================================
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
62C ============================================================================
63C
64C Diffcross calculates the fivefold differential cross section
65C
66C dsigma/ dp+t dp-t dy+ dy- ddeltaphi
67C
68C the trivial integration over the total phi-dependence is not included
69C
70C ppvt= absolute value of the positron transverse momentum (MeV)
71C pmvt= - " - electron - " -
72C dphi= phi-angle between the electron and positron transverse momentum
73C yp = rapidity of the positron
74C ym = rapidity of the electron
75C defined to make E = sqrt(pt^2 + m^2) cosh (y)
76C and Pz= sqrt(pt^2 + m^2) sinh (y)
77C
78C dsigma = differential cross section in kbarn/MeV^4/ (Zalpha)**4
79C
80C to get the total cross section, you have to integrate over
81C
82C Integral dsigma dyp dym ddphi 2 pi ppt dppt pmt dpmt
83C
84C (see also source for sigma.f)
85C
86C======================================================================
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
160C======================================================================
161C
162C Definition of propagatorterms
163C
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
197C
198C calculate all different integrals
199C
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
422C
423C dsigma is summ of all terms
424C
425 NT=N1+N2+N3+N4+N5+N6+N7+N8+N9+N10+N11+N12+N13+N14+
426 & N15+N16+N17+N18
427
428C
429C correction from w/u
430 NT=NT*4D0/beta**2
431C 1/(2pi)**6 from d3p, (2*pi)**2 from F.T.
432 NT=NT/(2*pi)**6*(2*pi)**2
433C from 1/2E+ 1/2E-
434 NT=NT/4D0
435C 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
445C========================================================================
446C All the differential integral forms are calculated here
447C
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
474C ---------------------------------------------------------------
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
499C -----------------------------------------------------------
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
533C -------------------------------------------------------------
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
563C -------------------------------------------------------------
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
602C --------------------------------------------------------
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
669C --------------------------------------------------------
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
819C --------------------------------------------------------
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
853C --------------------------------------------------------
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
890C --------------------------------------------------------
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