Empirical cut to increase robustness
[u/mrichter/AliRoot.git] / TEPEMGEN / diffcross.f
1 C**********************************************************************
2 C*                                                                    *
3 C* Exact calculation of the total differential e+ e- -pair production *
4 C* in Relativistic Heavy Ion Collisions for a point particle in an    *
5 C* external field approach.                                           *
6 C*                                                                    *
7 C* For details see the publication:                                   *
8 C* "Multiple electromagnetic electron positron pair production in     *
9 C*  relativistic heavy ion collisions"                                *
10 C*  Adrian Alscher, Kai Hencken, Dirk Trautmann, and Gerhard Baur     *
11 C*  Phys. Rev. A55 (1997) 396.                                        *
12 C*                                                                    *
13 C* Copyright (c) 1997, 1998, 2002, Adrian Alscher and Kai Hencken     *
14 C*                                                                    *
15 C* Permission to use, copy and distribute this software and its       *
16 C* documentation strictly for non-commercial purposes is hereby       *
17 C* granted without fee, provided that the above copyright notice      *
18 C* appears in all copies, that both the copyright notice and this     *
19 C* permission notice appear in the supporting documentation and that  *
20 C* the use of this program is acknowledged in scientific              *
21 C* publications (see reference above). The authors make no claims     *
22 C* about the suitability of this software for any purpose. It is      *
23 C* provided "as is" without express or implied warranty. Any change   *
24 C* of the code should be submitted to the authors                     *
25 C*                                                                    *
26 C* To use this program at LHC energies, please make sure that         *
27 C* "double precision" variables should better be real*16              *
28 C**********************************************************************
29
30 C======================================================================
31 C
32 C call this routine first to initialize some parameter needed in the
33 C function. 
34 C
35 C gm = Gamma_cm, that is, gm of each ion (~100 for RHIC ~3000 for LHC)
36 C mass = mass of the produced particle in MeV (~0.511 for e,~100 for mu)
37 C======================================================================
38
39       SUBROUTINE Initdiffcross (gm,mass)
40       IMPLICIT NONE
41       DOUBLE PRECISION gm,mass
42
43       DOUBLE PRECISION gamma,beta,m,w1xw1,w1xw2,w2xw2,wl,wy
44       COMMON/PHYSPARAM/gamma,beta,m,w1xw1,w1xw2,w2xw2,wl,wy
45
46       DOUBLE PRECISION ARCOSH,x
47       ARCOSH(x)=LOG(x+SQRT((x-1D0)*(x+1D0)))
48
49       gamma=gm
50       beta=sqrt((1D0-1D0/gamma)*(1D0+1D0/gamma))
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