]>
Commit | Line | Data |
---|---|---|
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 |