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 |
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 |