Branch data Line data Source code
1 : : !!****m* ABINIT/m_simtet
2 : : !! NAME
3 : : !! m_simtet
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !!
9 : : !! SOURCE
10 : :
11 : : #if defined HAVE_CONFIG_H
12 : : #include "config.h"
13 : : #endif
14 : :
15 : : #include "abi_common.h"
16 : :
17 : : module m_simtet
18 : :
19 : : use defs_basis
20 : : use m_errors, only: unused_var
21 : :
22 : : implicit none
23 : :
24 : : private
25 : :
26 : : public :: sim0onei
27 : : public :: sim0twoi
28 : : !!***
29 : :
30 : : contains
31 : :
32 : : !C file: sim0onei.f
33 : : !C date: 2011-04-01
34 : : !C who: S.Kaprzyk
35 : : !C what: Complex linear form integral over standard terahedron
36 : : !C ------------------------------------------------------------
37 : : !C - * -
38 : : !C - * 1 -
39 : : !C - SIM0=6 * dt1*dt2*dt3 ---------------------------------- -
40 : : !C - * verm(4)+(verm(1)-verm(4))*t1+ ..... -
41 : : !C - * -
42 : : !C - 0<t1+t2+t3<1 -
43 : : !C ------------------------------------------------------------
44 : 0 : SUBROUTINE SIM0ONEI(SIM0, SIM0I, VERM)
45 : :
46 : : !IMPLICIT NONE
47 : : INTEGER iuerr
48 : : PARAMETER (iuerr=6)
49 : : DOUBLE COMPLEX SIM0, SIM0I
50 : : DOUBLE COMPLEX VERM(4)
51 : : !C
52 : : DOUBLE COMPLEX SIM, SIMI
53 : : LOGICAL LONE(4),LEPS(3)
54 : : !c
55 : : DOUBLE COMPLEX VERM2D(3)
56 : : DOUBLE COMPLEX U(3)
57 : : DOUBLE PRECISION AS(3), AL(4)
58 : : INTEGER I, I1, I2, I3, K, N
59 : : DOUBLE PRECISION ZERO, EPS, SMALL
60 : : DATA EPS/1.0D-6/, SMALL/1.0D-5/
61 : : DATA ZERO/0.0D0/
62 : : !C
63 [ # # ]: 0 : DO 100 I = 1, 3
64 [ # # ]: 0 : IF (DIMAG(VERM(4))*DIMAG(VERM(I)).LT.ZERO) THEN
65 [ # # # # ]: 0 : WRITE (iuerr,9010) (DIMAG(VERM(K)),K=1,4)
66 : : 9010 FORMAT (' ***sim0onei: not signed ImgVERM()=',4(D13.6,1X))
67 : : !c STOP ' ***sim0onei: '
68 : : END IF
69 : 0 : 100 CONTINUE
70 [ # # ]: 0 : DO 200 I = 1, 4
71 : 0 : AL(I) = CDABS(VERM(I))
72 : 0 : 200 CONTINUE
73 [ # # ]: 0 : DO 300 N = 1, 4
74 [ # # ]: 0 : IF (AL(N).LT.SMALL) THEN
75 : 0 : I1 = MOD(N+0,4) + 1
76 : 0 : I2 = MOD(N+1,4) + 1
77 : 0 : I3 = MOD(N+2,4) + 1
78 : 0 : VERM2D(1) = VERM(I1)
79 : 0 : VERM2D(2) = VERM(I2)
80 : 0 : VERM2D(3) = VERM(I3)
81 : 0 : CALL S2D0ONEI(SIM,SIMI,VERM2D)
82 : 0 : SIM0 = 3*SIM/2
83 : 0 : SIM0I = SIMI + 1.75D0
84 : 0 : RETURN
85 : : END IF
86 : 0 : 300 CONTINUE
87 : : !c
88 : 0 : N = 4 ! N = 3, 2, 1
89 [ # # ]: 0 : DO 201 I = 1, 3
90 [ # # ]: 0 : IF (CDABS(VERM(I)).GT.CDABS(VERM(N))) N = I
91 : 0 : 201 CONTINUE
92 : : !C
93 : : !C Here are 4 cases a), b), c), d); see, Fig.3.1
94 : : !C LONE(1..4) [w1<w2<w3<ONE]; [w1<w2<ONE<w3];
95 : : !C [w1<ONE<w2<w3]; [ONE<w1<w2<w3]
96 : : !C LEPS(1..3) [|w1-w2|<e,|w3-w2|>e]; [|w1-w2|<e,|w3-w2|<e];
97 : : !C [|w1-w2|>e,|w3-w2|<e]
98 : : !*-ASIS
99 : 0 : CALL SIM0LEPS(VERM, N, U, AS, LONE, EPS, LEPS, iuerr)
100 : : !c U(1) = W(1); U(2) = W(2); U(3) = W(3)
101 [ # # ]: 0 : IF (LONE(1)) THEN
102 : 0 : CALL SIM1ONEI(U,AS,LEPS,SIM,SIMI)
103 : : END IF
104 : : !c U(1) = W(1); U(2) = W(2); U(3) = CONE/W(3)
105 [ # # ]: 0 : IF (LONE(2)) THEN
106 : 0 : CALL SIM2ONEI(U,AS,LEPS,SIM,SIMI)
107 : : END IF
108 : : !c U(1) = W(1); U(2) = CONE/W(2); U(3) = CONE/W(3)
109 [ # # ]: 0 : IF (LONE(3)) THEN
110 : 0 : CALL SIM3ONEI(U,AS,LEPS,SIM,SIMI)
111 : : END IF
112 : : !c U(1) = CONE/W(1); U(2) = CONE/W(2); U(3) = CONE/W(3)
113 [ # # ]: 0 : IF (LONE(4)) THEN
114 : 0 : CALL SIM4ONEI(U,AS,LEPS,SIM,SIMI)
115 : : END IF
116 : 0 : SIM0 = 0.75D0*SIM/VERM(N)
117 : 0 : SIM0I = log(VERM(N)) + 0.25D0*SIMI
118 : 0 : SIM0I = SIM0I + (0.25D0,0.0D0)
119 : 0 : RETURN
120 : : END SUBROUTINE SIM0ONEI
121 : : !C
122 : : !C -----------------------------------------------------------------
123 : : !C file: sim0twoi.f
124 : : !C date: 2011-04-02
125 : : !C who: S.Kaprzyk
126 : : !C what: Complex linear form integrals over standard tetrahedron
127 : : !C ------------------(i=1,2,3)--------------------------------------
128 : : !C * -
129 : : !C * t_i -
130 : : !C VERL(i)=6*dt1dt2dt3 ------------------------------------------ -
131 : : !C * VERM(4)+[VERM(1)-VERM(4)]*t1+... -
132 : : !C * -
133 : : !C 0<t1+t2+t3<1 -
134 : : !C -------------------(i=4)-----------------------------------------
135 : : !C * -
136 : : !C * 1- t1 -t2 - t3 -
137 : : !C VERL(4)=6*dt1dt2dt3 ------------------------------------------ -
138 : : !C * VERM(4)+[VERM(1)-VERM(4)]*t1+... -
139 : : !C * -
140 : : !C 0<t1+t2+t3<1 -
141 : : !C -----------------------------------------------------------------
142 : 55164096 : SUBROUTINE SIM0TWOI(VERL, VERLI, VERM)
143 : :
144 : : !IMPLICIT NONE
145 : : INTEGER iuerr
146 : : PARAMETER (iuerr=6)
147 : : DOUBLE COMPLEX VERL(4), VERLI(4), VERM(4)
148 : : !C
149 : : DOUBLE COMPLEX SIM, SIMI
150 : : DOUBLE COMPLEX VERL2D(3), VERL2DI(3), VERM2D(3)
151 : : LOGICAL LONE(4), LEPS(3)
152 : : !c
153 : : INTEGER I, I1, I2, I3, K, N
154 : : DOUBLE PRECISION AS(3), AL(4)
155 : : DOUBLE COMPLEX U(3)
156 : : DOUBLE COMPLEX CZERO
157 : : DOUBLE PRECISION EPS, SMALL
158 : : DOUBLE PRECISION ZERO
159 : : DATA EPS/1.0D-6/, SMALL/1.0D-5/
160 : : DATA ZERO/0.0D0/
161 : : !C
162 : 55164096 : CZERO = (0.0D0,0.0D0)
163 [ + + ]: 220656384 : DO 50 I = 1, 3
164 [ - + ]: 165492288 : IF (DIMAG(VERM(4))*DIMAG(VERM(I)).LT.ZERO) THEN
165 [ # # # # ]: 0 : WRITE (iuerr,9010) (DIMAG(VERM(K)),K=1,4)
166 : : 9010 FORMAT (' ***sim0twoi: not signed ImgVERM()=',4(D13.6,1X))
167 : : !c STOP ' ***sim0twoi: '
168 : : END IF
169 : 55164096 : 50 CONTINUE
170 [ + + ]: 275820480 : DO 200 I = 1, 4
171 : 220656384 : AL(I) = CDABS(VERM(I))
172 : 55164096 : 200 CONTINUE
173 : : !c
174 [ + + ]: 275820480 : DO 100 N = 1, 4
175 : 220656384 : VERL(N) = CZERO
176 : 220656384 : VERLI(N) = CZERO
177 : 220656384 : I1 = MOD(N+0,4) + 1
178 : 220656384 : I2 = MOD(N+1,4) + 1
179 : 220656384 : I3 = MOD(N+2,4) + 1
180 [ - + ]: 220656384 : IF (AL(I1).LT.SMALL) THEN
181 : 0 : VERM2D(1) = VERM(I2)
182 : 0 : VERM2D(2) = VERM(I3)
183 : 0 : VERM2D(3) = VERM(N)
184 : 0 : CALL S2D0TWOI(VERL2D,VERL2DI,VERM2D)
185 : 0 : VERL(N) = VERL2D(3)
186 : 0 : VERLI(N) = VERL2DI(3)*3/4
187 : 0 : GO TO 100
188 : : END IF
189 [ - + ]: 220656384 : IF (AL(I2).LT.SMALL) THEN
190 : 0 : VERM2D(1) = VERM(I3)
191 : 0 : VERM2D(2) = VERM(N)
192 : 0 : VERM2D(3) = VERM(I1)
193 : 0 : CALL S2D0TWOI(VERL2D,VERL2DI,VERM2D)
194 : 0 : VERL(N) = VERL2D(2)
195 : 0 : VERLI(N) = VERL2DI(2)*3/4
196 : 0 : GO TO 100
197 : : END IF
198 [ - + ]: 220656384 : IF (AL(I3).LT.SMALL) THEN
199 : 0 : VERM2D(1) = VERM(N)
200 : 0 : VERM2D(2) = VERM(I1)
201 : 0 : VERM2D(3) = VERM(I2)
202 : 0 : CALL S2D0TWOI(VERL2D,VERL2DI,VERM2D)
203 : 0 : VERL(N) = VERL2D(1)
204 : 0 : VERLI(N) = VERL2DI(1)*3/4
205 : 0 : GO TO 100
206 : : END IF
207 [ - + ]: 220656384 : IF (AL(N).LT.SMALL) THEN
208 : 0 : VERM2D(1) = VERM(I1)
209 : 0 : VERM2D(2) = VERM(I2)
210 : 0 : VERM2D(3) = VERM(I3)
211 : 0 : CALL S2D0ONEI(SIM,SIMI,VERM2D)
212 : 0 : VERL(N) = SIM/2
213 : 0 : VERLI(N) = SIMI/4 +1.75D0
214 : 0 : GO TO 100
215 : : END IF
216 : : !C Here are 4 cases a), b), c), d); see, Fig.3.1
217 : : !C LONE(1..4) [w1<w2<w3<ONE]; [w1<w2<ONE<w3];
218 : : !C [w1<ONE<w2<w3]; [ONE<w1<w2<w3]
219 : : !C LEPS(1..3) [|w1-w2|<e,|w3-w2|>e]; [|w1-w2|<e,|w3-w2|<e];
220 : : !C [|w1-w2|>e,|w3-w2|<e]
221 : : !*-ASIS
222 : 220656384 : CALL SIM0LEPS(VERM, N, U, AS, LONE, EPS, LEPS, iuerr)
223 : : !c U(1) = W(1); U(2) = W(2); U(3) = W(3)
224 [ + + ]: 220656384 : IF (LONE(1)) THEN
225 : 209614532 : CALL SIM1TWOI(U,AS,LEPS,SIM,SIMI)
226 : : END IF
227 : : !c U(1) = W(1); U(2) = W(2); U(3) = CONE/W(3)
228 [ + + ]: 220656384 : IF (LONE(2)) THEN
229 : 6337479 : CALL SIM2TWOI(U,AS,LEPS,SIM,SIMI)
230 : : END IF
231 : : !c U(1) = W(1); U(2) = CONE/W(2); U(3) = CONE/W(3)
232 [ + + ]: 220656384 : IF (LONE(3)) THEN
233 : 2626200 : CALL SIM3TWOI(U,AS,LEPS,SIM,SIMI)
234 : : END IF
235 : : !c U(1) = CONE/W(1); U(2) = CONE/W(2); U(3) = CONE/W(3)
236 [ + + ]: 220656384 : IF (LONE(4)) THEN
237 : 2078173 : CALL SIM4TWOI(U,AS,LEPS,SIM,SIMI)
238 : : END IF
239 : : !c
240 : 220656384 : VERL(N) = SIM/VERM(N)/8
241 : 220656384 : VERLI(N) = log(VERM(N))/4 + SIMI/32
242 : 220656384 : VERLI(N) = VERLI(N) + (0.25D0,0.0D0)
243 : 55164096 : 100 CONTINUE ! end DO ... N = 1, 4
244 : 55164096 : RETURN
245 : : END SUBROUTINE SIM0TWOI
246 : : !C
247 : : !C ---------------------------------------------------------------------
248 : : !C file: sim0ur0.f
249 : : !C date: 2011-03-16
250 : : !C who: S.Kaprzyk
251 : : !C what: R0(U)=Arth(w)/w=(ln(1+U)-ln(1-U))/(2*U)
252 : : !C what: R0(U)= 1 + U^2/3 + U^4/5 + ...
253 : 0 : DOUBLE COMPLEX FUNCTION SIM0UR0(U)
254 : :
255 : : !IMPLICIT NONE
256 : : DOUBLE COMPLEX U
257 : : !c
258 : : !c DOUBLE PRECISION DLAMCH
259 : : !c EXTERNAL DLAMCH
260 : : !c
261 : : INTEGER K
262 : : DOUBLE COMPLEX AV, BV, CV, UU
263 : : DOUBLE COMPLEX CONE
264 : : DOUBLE PRECISION DEV, TOL
265 : : DATA TOL /1.0D-13/
266 : : !c
267 : : !c TOL = DLAMCH('e')*10
268 : 0 : CONE = (1.0D0,0.0D0)
269 : 0 : UU = U*U
270 [ # # ]: 0 : IF (CDABS(U).GT.0.27D0) THEN
271 : 0 : SIM0UR0 = log((CONE+U)/(CONE-U))/(2*U)
272 : : !c SIM0UR0 = (log(CONE+U)-log(CONE-U))/(2*U)
273 : : ELSE
274 : : CV = UU
275 : : AV = CONE
276 : 0 : DO 50 K = 3, 13, 2
277 : 0 : AV = AV + CV/DBLE(K)
278 [ # # ]: 0 : CV = CV*UU
279 : 0 : 50 CONTINUE
280 : 0 : DO 100 K = 15, 49, 2
281 : 0 : BV = CV/DBLE(K)
282 : 0 : AV = AV + BV
283 : 0 : DEV = CDABS(BV)
284 [ # # ]: 0 : IF (DEV.LE.TOL*CDABS(AV)) GO TO 150
285 [ # # ]: 0 : CV = CV*UU
286 : 0 : 100 CONTINUE
287 : 0 : WRITE (std_out,9010) U, DEV, TOL
288 : : 9010 FORMAT ('***sim0ur0: U=',2(D15.8,1X),' DEV=',D14.7,' TOL=',D14.7)
289 : 0 : STOP '***sim0ur0: accuracy not reached'
290 : : 150 CONTINUE
291 : : SIM0UR0 = AV
292 : : END IF
293 : : !C
294 : : RETURN
295 : : END FUNCTION SIM0UR0
296 : : !C
297 : : !C ----------------------------------------------------------------------
298 : : !C file: sim0ux0.f
299 : : !C date: 2010-08-17
300 : : !C who: S.Kaprzyk
301 : : !C what: Complex version of X0(U)=(R0(U)/U-1-U^2/3-U^4/5)/U^6
302 : : !C R0(U)=Arth(w)/w=Log((1+U)/(1-U))/(2*U)
303 : 661969152 : DOUBLE COMPLEX FUNCTION SIM0UX0(U)
304 : :
305 : : !IMPLICIT NONE
306 : : DOUBLE COMPLEX U
307 : : !c
308 : : !c DOUBLE PRECISION DLAMCH
309 : : !c EXTERNAL DLAMCH
310 : : !c
311 : : INTEGER K
312 : : DOUBLE COMPLEX AV, CV, UU
313 : : DOUBLE COMPLEX CONE
314 : : DOUBLE PRECISION DEV, TOL
315 : : DATA TOL /1.0D-13/
316 : : !c TOL = DLAMCH('e')*10
317 : 661969152 : CONE = (1.0D0,0.0D0)
318 : 661969152 : UU = U*U
319 [ + + ]: 661969152 : IF (CDABS(U).GT.0.27D0) THEN
320 : 52466800 : SIM0UX0 =(log((CONE+U)/(CONE-U))/(2*U) - (CONE+UU*(CONE/3.0D0+UU*CONE/5.0D0)))/(UU*UU*UU)
321 : : ELSE
322 : : CV = UU
323 : : AV = CONE/7.0D0
324 : 1828507056 : DO 50 K = 9, 13, 2
325 : 1828507056 : AV = AV + CV/DBLE(K)
326 [ + + ]: 1828507056 : CV = CV*UU
327 : 609502352 : 50 CONTINUE
328 [ + - ]: 867229108 : DO 100 K = 15, 49, 2
329 : 1476731460 : AV = AV + CV/DBLE(K)
330 : 1476731460 : DEV = CDABS(CV/AV)/DBLE(K)
331 : 1476731460 : CV = CV*UU
332 [ + + ]: 1476731460 : IF (DEV.LE.TOL) GO TO 150
333 : 0 : 100 CONTINUE
334 : 0 : WRITE (std_out,9010) U, DEV, TOL
335 : : 9010 FORMAT ('***sim0ux0: U=',2(D15.8,1X), ' DEV=',D14.7,' TOL=',D14.7)
336 : 0 : STOP '***sim0ux0: accuracy not reached'
337 : : 150 CONTINUE
338 : : SIM0UX0 = AV
339 : : END IF
340 : : !C
341 : : RETURN
342 : : END FUNCTION SIM0UX0
343 : : !C
344 : : !C ----------------------------------------------------------------------
345 : : !C file: sim0udx0.f
346 : : !C date: 2010-08-17
347 : : !C who: S.Kaprzyk
348 : : !C what: Complex version of DX0(u)=dX0/du
349 : 55560960 : SUBROUTINE SIM0UDX0(U, X0, DX0)
350 : :
351 : : !IMPLICIT NONE
352 : : DOUBLE COMPLEX U, X0, DX0(4)
353 : : !c
354 : : INTEGER I
355 : : DOUBLE COMPLEX AV, BV, CV, DV
356 : : DOUBLE COMPLEX CZERO, CONE
357 : : !C
358 : 55560960 : CZERO = (0.0D0,0.0D0)
359 : 55560960 : CONE = (1.0D0,0.0D0)
360 : : !c
361 : 55560960 : AV = U
362 : 55560960 : BV = AV*AV
363 [ + + ]: 55560960 : IF (CDABS(BV).GE.0.37D0) THEN
364 : 1358638 : BV = CONE/(CONE-BV)
365 : 1358638 : CV = BV*BV
366 : 1358638 : DX0(1) = (BV-7.0D0*X0)/AV
367 : 1358638 : DX0(2) = 2.0D0*(CV-4.0D0*DX0(1)/AV)
368 : 1358638 : DX0(3) = 8.0D0*AV*CV*BV + (2.0D0*CV-9.0D0*DX0(2))/AV
369 : 1358638 : DX0(4) = 24.0D0*CV*BV*(CONE+2.0D0*AV*AV*BV) - 10.0D0*DX0(3)/AV
370 : : ELSE
371 : 54202322 : DX0(3) = CZERO
372 : 54202322 : DX0(2) = (2.0D0,0.0D0)/9.0D0
373 : 54202322 : DX0(1) = DX0(2)*AV
374 : 54202322 : DX0(4) = (24.0D0,0.0D0)/11.0D0
375 : 54202322 : CV = AV
376 : 54202322 : BV = AV*AV
377 : 54202322 : DO 50 I = 2, 48, 2
378 : 1300855728 : DV = DBLE(2+I)/DBLE(9+I)*CV
379 : 1300855728 : DX0(4) = DX0(4) + DBLE((4+I)*(3+I)*(2+I)*(1+I))/DBLE(11+I)*CV*AV
380 : 1300855728 : DX0(3) = DX0(3) + DBLE((1+I)*I)*DV
381 : 1300855728 : DX0(2) = DX0(2) + DBLE(1+I)*DV*AV
382 : 1300855728 : DX0(1) = DX0(1) + DV*BV
383 [ + + ]: 1300855728 : CV = CV*BV
384 : 54202322 : 50 CONTINUE
385 : : END IF
386 : : !C
387 : 55560960 : RETURN
388 : : END SUBROUTINE SIM0UDX0
389 : : !C
390 : : !C ----------------------------------------------------------------------
391 : : !C file: sim1onei.f
392 : : !C date: 2010-09-22
393 : : !C who: S.Kaprzyk
394 : : !C what: CASE a): [w1<w2<w3<ONE]
395 : : !C U(1)=W(1); U(2)=W(2); U(3)=W(3)
396 : : !C LEPS(1..3) [|w1-w2|<e,|w3-w2|>e];
397 : : !C [|w1-w2|<e,|w3-w2|<e];
398 : : !C [|w1-w2|>e,|w3-w2|<e]
399 : 0 : SUBROUTINE SIM1ONEI(U, AS, LEPS, SIM1, SIM1I)
400 : :
401 : : !IMPLICIT NONE
402 : : DOUBLE COMPLEX U(3)
403 : : DOUBLE PRECISION AS(3)
404 : : LOGICAL LEPS(3)
405 : : DOUBLE COMPLEX SIM1, SIM1I
406 : : !C
407 : : !DOUBLE COMPLEX SIM0UX0
408 : : !EXTERNAL SIM0UX0
409 : : DOUBLE COMPLEX X0(3), QX0(3), DX0(4)
410 : : DOUBLE COMPLEX R0(3), QR0(3)
411 : : DOUBLE COMPLEX S0(3), QS0(3)
412 : : DOUBLE COMPLEX AV, BV, CV, DU, G(5)
413 : : INTEGER I, K !, N
414 : : DOUBLE COMPLEX CZERO, CONE
415 : : DOUBLE PRECISION ZERO, ONE
416 : : DATA ZERO/0.0D0/, ONE/1.0D0/
417 : : !C
418 : : ABI_UNUSED(as)
419 : :
420 : 0 : CZERO = (0.0D0,0.0D0)
421 : 0 : CONE = (1.0D0,0.0D0)
422 : : !C ------------------------------------------------------------------
423 : : !C R0(U) = ARTH(U)/U
424 : : !C R0(U) = 1+U**2/3+U**6*X0(U)
425 : : !C ------------------------------------------------------------------
426 [ # # ]: 0 : DO 100 I = 1, 3
427 : 0 : AV = U(I)*U(I)
428 : 0 : X0(I) = SIM0UX0(U(I))
429 : 0 : R0(I) = CONE + AV*(CONE/3.0D0+AV*(CONE/5.0D0+AV*X0(I)))
430 : 0 : S0(I) = (CONE-U(I))/(CONE+U(I))*R0(I)
431 : 0 : 100 CONTINUE
432 : : !C -----------------------------------------------------------------
433 : : !C - The DX0(I) contains derivatives of the function X0(u) -
434 : : !C -----------------------------------------------------------------
435 [ # # # # : 0 : IF (LEPS(1).OR.LEPS(2).OR.LEPS(3)) THEN
# # ]
436 : 0 : CALL SIM0UDX0(U(2),X0(2),DX0)
437 : : END IF
438 : 0 : DO 200 I = 1, 3, 2
439 : 0 : DU = U(I) - U(2)
440 [ # # ]: 0 : IF (.NOT.LEPS(I)) THEN
441 : 0 : QX0(I) = (X0(I)-X0(2))/DU
442 : : ELSE
443 : 0 : QX0(I) = DX0(1) + (DX0(2)/2.0D0+(DX0(3)/6.0D0+DX0(4)/24.0D0*DU) *DU)*DU
444 : : END IF
445 : 0 : AV = U(I)
446 : 0 : G(1) = U(I) + U(2)
447 [ # # ]: 0 : DO 150 K = 2, 5
448 : 0 : AV = AV*U(I)
449 : 0 : G(K) = G(K-1)*U(2) + AV
450 : 0 : 150 CONTINUE
451 : 0 : QR0(I) = G(1)/3.0D0 + G(3)/5.0D0 + G(5)*X0(2) + AV*U(I)*QX0(I)
452 [ # # ]: 0 : QS0(I) = (-S0(2)-R0(2)+(CONE-U(I))*QR0(I))/(CONE+U(I))
453 : 0 : 200 CONTINUE ! end DO ... I = 1, 3, 2
454 [ # # ]: 0 : IF (LEPS(2)) THEN
455 : : !*-ASIS
456 : : QX0(2) = DX0(2)/2.0D0 + DX0(3)*(U(1)+U(3)-2.0D0*U(2))/6.0D0+ &
457 : : & DX0(4)*((U(3)-U(2))*(U(3)-U(2))+(U(1)-U(2))*(U(1)-U(2))+ &
458 : 0 : & (U(3)-U(2))*(U(1)-U(2)))/24.0D0
459 : : ELSE
460 : 0 : QX0(2) = (QX0(3)-QX0(1))/(U(3)-U(1))
461 : : END IF
462 : 0 : AV = U(1)
463 : 0 : G(1) = U(1) + U(3)
464 [ # # ]: 0 : DO 300 K = 2, 5
465 : 0 : AV = AV*U(1)
466 : 0 : G(K) = G(K-1)*U(3) + AV
467 : 0 : 300 CONTINUE
468 : : !*-ASIS
469 : : QR0(2) = CONE/3.0D0 + (G(2)+U(2)*(U(2)+G(1)))/5.0D0 + &
470 : : & (G(4)+U(2)*(G(3)+U(2)*(G(2)+U(2)*(G(1)+U(2)))))*X0(2)+ &
471 : 0 : & G(5)*QX0(3) + AV*U(1)*QX0(2)
472 : 0 : QS0(2) = (-QS0(3)-QR0(3)+(CONE-U(1))*QR0(2))/(CONE+U(1))
473 : : !c
474 : 0 : AV = (CONE+U(1))*(CONE+U(2))*(CONE+U(3))
475 : 0 : BV = U(3) + U(1) - (2.0D0,0.0D0)
476 : 0 : CV = (CONE-U(1))*(CONE-U(1))
477 : 0 : SIM1 = AV*(R0(2)+QR0(3)*BV+QR0(2)*CV)
478 : 0 : SIM1I = AV*(S0(2)+QS0(3)*BV+QS0(2)*CV)
479 : 0 : RETURN
480 : : END SUBROUTINE SIM1ONEI
481 : : !C
482 : : !C ----------------------------------------------------------------------
483 : : !C file: sim1twoi.f
484 : : !C date: 2010-09-22
485 : : !C who: S.Kaprzyk
486 : : !C what: CASE a): [w1<w2<w3<ONE]
487 : : !C U(1)=W(1); U(2)=W(2); U(3)=W(3)
488 : : !C LEPS(1..3) [|w1-w2|<e,|w3-w2|>e];
489 : : !C [|w1-w2|<e,|w3-w2|<e];
490 : : !C [|w1-w2|>e,|w3-w2|<e]
491 : 209614532 : SUBROUTINE SIM1TWOI(U, AS, LEPS, SIM, SIMI)
492 : :
493 : : !IMPLICIT NONE
494 : : DOUBLE COMPLEX U(3)
495 : : DOUBLE PRECISION AS(3)
496 : : LOGICAL LEPS(3)
497 : : DOUBLE COMPLEX SIM, SIMI
498 : : !C
499 : : !DOUBLE COMPLEX SIM0UX0
500 : : !EXTERNAL SIM0UX0
501 : : !c
502 : : DOUBLE COMPLEX R4(3), X0(3), DX0(4)
503 : : DOUBLE COMPLEX QR4(3), QX0(3)
504 : : DOUBLE COMPLEX S4(3), QS4(3)
505 : : DOUBLE COMPLEX AV, BV, CV, UV, G(5)
506 : : INTEGER I, K
507 : : DOUBLE COMPLEX CZERO, CONE, CTWO
508 : : DOUBLE PRECISION ZERO, ONE
509 : : DATA ZERO/0.0D0/, ONE/1.0D0/
510 : : !C
511 : : ABI_UNUSED(as)
512 : 209614532 : CZERO = (0.0D0,0.0D0)
513 : 209614532 : CONE = (1.0D0,0.0D0)
514 : 209614532 : CTWO = (2.0D0,0.0D0)
515 : : !C -----------------------------------------------------------------
516 : : !C - -
517 : : !C - R4(W)=1+(W-1)*(ARTH(W)/W-1)/W -
518 : : !C - R4(W)=1-W/3+W**2/3-W**3/5+W**4/5+(W-1)*W**5*X0(W) -
519 : : !C - -
520 : : !C -----------------------------------------------------------------
521 [ + + ]: 838458128 : DO 100 I = 1, 3
522 : 628843596 : AV = U(I)
523 : 628843596 : X0(I) = SIM0UX0(AV)
524 : : !*-ASIS
525 : : R4(I) = CONE + AV*(-CONE/3.0D0 + AV*(CONE/3.0D0 + &
526 : 628843596 : & AV*((-0.2D0,0.0D0)+AV*((0.2D0,0.0D0)+AV*(AV-CONE)*X0(I)))))
527 : 628843596 : S4(I) = (CONE-U(I))/(CONE+U(I))*R4(I)
528 : 209614532 : 100 CONTINUE
529 : : !C -----------------------------------------------------------------
530 : : !C - The DX0(I) Contains Derivatives Of The Function X0(W) -
531 : : !C -----------------------------------------------------------------
532 [ + + + - : 209614532 : IF (LEPS(1).OR.LEPS(2).OR.LEPS(3)) THEN
+ + ]
533 : 53619796 : CALL SIM0UDX0(U(2),X0(2),DX0)
534 : : END IF
535 : : !C
536 : 419229064 : DO 200 I = 1, 3, 2
537 : 419229064 : UV = U(I) - U(2)
538 [ + + ]: 419229064 : IF (LEPS(I)) THEN
539 : 53619796 : QX0(I) = DX0(1) + (DX0(2)/2.0D0+(DX0(3)/6.0D0+DX0(4)/24.0D0*UV) *UV)*UV
540 : : ELSE
541 : 365609268 : QX0(I) = (X0(I)-X0(2))/UV
542 : : END IF
543 : 419229064 : AV = U(I)
544 : 419229064 : G(1) = U(I) + U(2)
545 [ + + ]: 2096145320 : DO 150 K = 2, 5
546 : 1676916256 : AV = AV*U(I)
547 : 1676916256 : G(K) = G(K-1)*U(2) + AV
548 : 419229064 : 150 CONTINUE
549 : : QR4(I) = -CONE/3.0D0 + G(1)/3.0D0 - G(2)/5.0D0 + G(3) &
550 : 419229064 : & /5.0D0 + (G(5)-G(4))*X0(2) + (U(I)-CONE)*AV*QX0(I)
551 [ + + ]: 419229064 : QS4(I) = (-S4(2)-R4(2)+(CONE-U(I))*QR4(I))/(CONE+U(I))
552 : 209614532 : 200 CONTINUE ! DO 750 I = 1, 3, 2
553 [ - + ]: 209614532 : IF (LEPS(2)) THEN
554 : : !*-ASIS
555 : : QX0(2) = DX0(2)/2.0D0 + DX0(3)*(U(1)+U(3)-2.0D0*U(2))/6.0D0 + &
556 : : & DX0(4)*((U(3)-U(2))*(U(3)-U(2))+(U(1)-U(2))*(U(1)-U(2))+ &
557 : 0 : & (U(3)-U(2))*(U(1)-U(2)))/24.0D0
558 : : ELSE
559 : 209614532 : QX0(2) = (QX0(3)-QX0(1))/(U(3)-U(1))
560 : : END IF
561 : 209614532 : AV = U(1)
562 : 209614532 : G(1) = U(1) + U(3)
563 [ + + ]: 1048072660 : DO 300 K = 2, 5
564 : 838458128 : AV = AV*U(1)
565 : 838458128 : G(K) = G(K-1)*U(3) + AV
566 : 209614532 : 300 CONTINUE
567 : : !*-ASIS
568 : : QR4(2) = CONE/3.0D0- (U(2)+G(1))/5.0D0 + &
569 : : & (G(2)+U(2)*(U(2)+G(1)))/5.0D0 + &
570 : : & (G(4)+(U(2)-CONE)*(G(3) + G(2)*U(2) + G(1)*U(2)*U(2) + &
571 : : & U(2)*U(2)*U(2)))*X0(2) + (G(5)-G(4))*QX0(3) + &
572 : 209614532 : & (U(1)-CONE)*AV*QX0(2)
573 : 209614532 : QS4(2) = (-QS4(3)-QR4(3)+(CONE-U(1))*QR4(2))/(CONE+U(1))
574 : : !c
575 : 209614532 : AV = (CONE+U(1))*(CONE+U(2))*(CONE+U(3))
576 : 209614532 : BV = U(3) + U(1) - (2.0D0,0.0D0)
577 : 209614532 : CV = (CONE-U(1))*(CONE-U(1))
578 : 209614532 : SIM = AV*(R4(2)+QR4(3)*BV+QR4(2)*CV)
579 : 209614532 : SIMI = AV*(S4(2)+QS4(3)*BV+QS4(2)*CV)
580 : 209614532 : RETURN
581 : : END SUBROUTINE SIM1TWOI
582 : : !C
583 : : !C ----------------------------------------------------------------------
584 : : !C file: sim2onei.f
585 : : !C date: 2010-08-19
586 : : !C who: S.Kaprzyk
587 : : !C what: CASE b) : [w1<w2<ONE<w3]
588 : : !C U(1)=W(1); U(2)=W(2); U(3)=CONE/W(3)
589 : : !C LEPS(1..3) [|w1-w2|<e,|w3-w2|>e];
590 : : !C [|w1-w2|<e,|w3-w2|<e];
591 : : !C [|w1-w2|>e,|w3-w2|<e]
592 : 0 : SUBROUTINE SIM2ONEI(U,AS,LEPS, SIM2, SIM2I)
593 : :
594 : : !IMPLICIT NONE
595 : : DOUBLE COMPLEX U(3)
596 : : DOUBLE PRECISION AS(3)
597 : : LOGICAL LEPS(3)
598 : : DOUBLE COMPLEX SIM2, SIM2I
599 : : !c
600 : : !DOUBLE COMPLEX SIM0UX0
601 : : !EXTERNAL SIM0UX0
602 : : DOUBLE COMPLEX X0(3), QX0(3), DX0(4)
603 : : DOUBLE COMPLEX R0(3), QR0(3)
604 : : DOUBLE COMPLEX S0(3), QS0(3)
605 : : DOUBLE COMPLEX AV, BV, CV, DV, G(5)
606 : : DOUBLE COMPLEX UD
607 : : INTEGER I, K !, N
608 : : DOUBLE COMPLEX CZERO, CONE, CTWO
609 : : DOUBLE PRECISION ZERO, ONE
610 : : DATA ZERO/0.0D0/, ONE/1.0D0/
611 : : !C
612 : 0 : CZERO = (0.0D0,0.0D0)
613 : 0 : CONE = (1.0D0,0.0D0)
614 : 0 : CTWO = (2.0D0,0.0D0)
615 : : !C ------------------------------------------------------------------
616 : : !C R0(U) = ARTH(U)/U
617 : : !C R0(U) = 1+U**2/3+U**6*X0(U)
618 : : !C ------------------------------------------------------------------
619 [ # # ]: 0 : DO 100 I = 1, 3
620 : 0 : AV = U(I)*U(I)
621 : 0 : X0(I) = SIM0UX0(U(I))
622 : 0 : R0(I) = CONE + AV*(CONE/3.0D0+AV*(CONE/5.0D0+AV*X0(I)))
623 : 0 : S0(I) = (CONE-U(I))/(CONE+U(I))*R0(I)
624 : 0 : 100 CONTINUE
625 : : !C -----------------------------------------------------------------
626 : : !C - THE DX0(I) CONTAINS DERIVATIVES OF THE FUNCTION X0(W) -
627 : : !C -----------------------------------------------------------------
628 [ # # # # : 0 : IF (LEPS(1).OR.LEPS(2).OR.LEPS(3)) THEN
# # ]
629 : 0 : CALL SIM0UDX0(U(2),X0(2),DX0)
630 : : END IF
631 : : !c DO 400 I = 1, 3, 2
632 : 0 : I = 1
633 : 0 : UD = U(I) - U(2)
634 [ # # ]: 0 : IF (.NOT.LEPS(I)) THEN
635 : 0 : QX0(I) = (X0(I)-X0(2))/UD
636 : : ELSE
637 : 0 : QX0(I) = DX0(1) + (DX0(2)/2.0D0+(DX0(3)/6.0D0+DX0(4)/24.0D0*UD) *UD)*UD
638 : : END IF
639 : 0 : AV = U(I)
640 : 0 : G(1) = U(I) + U(2)
641 [ # # ]: 0 : DO 200 K = 2, 5
642 : 0 : AV = AV*U(I)
643 : 0 : G(K) = G(K-1)*U(2) + AV
644 : 0 : 200 CONTINUE
645 : 0 : QR0(I) = G(1)/3.0D0 + G(3)/5.0D0 + G(5)*X0(2) + AV*U(I)*QX0(I)
646 : 0 : QS0(I) = (-S0(2)-R0(2)+(CONE-U(I))*QR0(I))/(CONE+U(I))
647 : : !c 400 CONTINUE ! end DO ... I = 1, 3, 2
648 : 0 : AV = (CONE+U(1))*(CONE+U(2))*(U(3)+CONE) /((U(3)*U(1)-CONE)*(U(3)*U(2)-CONE))
649 : 0 : BV = CTWO - (U(1)+U(2)) - U(3)*(CONE-U(1)*U(2))
650 : 0 : CV = (U(3)*U(2)-CONE)*(CONE-U(1))**2
651 : 0 : DV = (U(3)-CONE)**2
652 : 0 : SIM2 = AV*(BV*R0(2)+CV*QR0(1)+DV*(DCMPLX(ZERO,AS(3))+U(3)*R0(3)))
653 : 0 : DV = DV*(U(3)-CONE)/(U(3)+CONE)
654 : 0 : SIM2I = AV*(BV*S0(2)+CV*QS0(1)+DV*(DCMPLX(ZERO,AS(3))+U(3)*R0(3)))
655 : : !c
656 : 0 : RETURN
657 : : END SUBROUTINE SIM2ONEI
658 : : !C
659 : : !C ----------------------------------------------------------------------
660 : : !C file: sim2twoi.f
661 : : !C date: 2010-08-19
662 : : !C who: S.Kaprzyk
663 : : !C what: CASE b) : [w1<w2<ONE<w3]
664 : : !C U(1)=W(1); U(2)=W(2); U(3)=CONE/W(3)
665 : : !C LEPS(1..3) [|w1-w2|<e,|w3-w2|>e];
666 : : !C [|w1-w2|<e,|w3-w2|<e];
667 : : !C [|w1-w2|>e,|w3-w2|<e]
668 : 6337479 : SUBROUTINE SIM2TWOI(U,AS,LEPS, SIM, SIMI)
669 : :
670 : : !IMPLICIT NONE
671 : : DOUBLE COMPLEX U(3)
672 : : DOUBLE PRECISION AS(3)
673 : : LOGICAL LEPS(3)
674 : : DOUBLE COMPLEX SIM, SIMI
675 : : !C
676 : : !DOUBLE COMPLEX SIM0UX0
677 : : !EXTERNAL SIM0UX0
678 : : !c
679 : : DOUBLE COMPLEX R4(3), QR4(3)
680 : : DOUBLE COMPLEX S4(3), QS4(3)
681 : : DOUBLE COMPLEX X0(3), QX0(3), DX0(4), G(5)
682 : : DOUBLE COMPLEX AV, BV, CV, DV, UV
683 : : INTEGER I, K
684 : : DOUBLE COMPLEX CZERO, CONE, CTWO
685 : : DOUBLE PRECISION ZERO, ONE
686 : : DATA ZERO/0.0D0/, ONE/1.0D0/
687 : : !C
688 : 6337479 : CZERO = (0.0D0,0.0D0)
689 : 6337479 : CONE = (1.0D0,0.0D0)
690 : 6337479 : CTWO = (2.0D0,0.0D0)
691 : : !C -----------------------------------------------------------------
692 : : !C - -
693 : : !C - R4(W)=1+(W-1)*(ARTH(W)/W-1)/W -
694 : : !C - R4(W)=1-W/3+W**2/3-W**3/5+W**4/5+(W-1)*W**5*X0(W) -
695 : : !C - -
696 : : !C -----------------------------------------------------------------
697 [ + + ]: 25349916 : DO 100 I = 1, 3
698 : 19012437 : AV = U(I)
699 : 19012437 : X0(I) = SIM0UX0(AV)
700 : : !*-ASIS
701 : : R4(I) = CONE + AV*(-CONE/3.0D0 + AV*(CONE/3.0D0 + &
702 : 19012437 : & AV*((-0.2D0,0.0D0)+AV*((0.2D0,0.0D0)+AV*(AV-CONE)*X0(I)))))
703 : 19012437 : S4(I) = (CONE-U(I))/(CONE+U(I))*R4(I)
704 : 6337479 : 100 CONTINUE
705 : : !C -----------------------------------------------------------------
706 : : !C - The DX0(I) Contains Derivatives Of The Function X0(W) -
707 : : !C -----------------------------------------------------------------
708 [ + + + - : 6337479 : IF (LEPS(1).OR.LEPS(2).OR.LEPS(3)) THEN
- + ]
709 : 716892 : CALL SIM0UDX0(U(2),X0(2),DX0)
710 : : END IF
711 : : !c DO 200 I = 1, 3, 2
712 : 6337479 : I = 1
713 : 6337479 : UV = U(I) - U(2)
714 [ + + ]: 6337479 : IF (LEPS(I)) THEN
715 : 716892 : QX0(I) = DX0(1) + (DX0(2)/2.0D0+(DX0(3)/6.0D0+DX0(4)/24.0D0*UV) *UV)*UV
716 : : ELSE
717 : 5620587 : QX0(I) = (X0(I)-X0(2))/UV
718 : : END IF
719 : 6337479 : AV = U(I)
720 : 6337479 : G(1) = U(I) + U(2)
721 [ + + ]: 31687395 : DO 200 K = 2, 5
722 : 25349916 : AV = AV*U(I)
723 : 25349916 : G(K) = G(K-1)*U(2) + AV
724 : 6337479 : 200 CONTINUE
725 : : !*-ASIS
726 : : QR4(I) = -CONE/3.0D0 + G(1)/3.0D0 - G(2)/5.0D0 + &
727 : 6337479 : & G(3) /5.0D0 + (G(5)-G(4))*X0(2) + (U(I)-CONE)*AV*QX0(I)
728 : 6337479 : QS4(I) = (-S4(2)-R4(2)+(CONE-U(I))*QR4(I))/(CONE+U(I))
729 : : !c 200 CONTINUE ! DO ... I = 1, 3, 2
730 : : AV = (CONE+U(1))*(CONE+U(2))*(U(3)+CONE) &
731 : 6337479 : & /((U(3)*U(1)-CONE)*(U(3)*U(2)-CONE))
732 : 6337479 : BV = CTWO - (U(1)+U(2)) - U(3)*(CONE-U(1)*U(2))
733 : 6337479 : CV = (U(3)*U(2)-CONE)*(CONE-U(1))**2
734 : 6337479 : DV = (U(3)-CONE)**2
735 : : SIM = AV*(BV*R4(2)+CV*QR4(1) &
736 : 6337479 : & +DV*((CONE-U(3))*DCMPLX(ZERO,AS(3))+CONE+U(3)-U(3)*U(3)*R4(3)))
737 : 6337479 : DV = DV*(U(3)-CONE)/(U(3)+CONE)
738 : : SIMI = AV*(BV*S4(2)+CV*QS4(1) &
739 : 6337479 : & +DV*((CONE-U(3))*DCMPLX(ZERO,AS(3))+CONE+U(3)-U(3)*U(3)*R4(3)))
740 : : !C
741 : 6337479 : RETURN
742 : : END SUBROUTINE SIM2TWOI
743 : : !C
744 : : !C ----------------------------------------------------------------------
745 : : !C file: sim3onei.f
746 : : !C date: 2010-08-19
747 : : !C who: S.Kaprzyk
748 : : !C what: CASE c) : [w1<ONE<w2<w3]
749 : : !C U(1)=W(1); U(2)=CONE/W(2); U(3)=CONE/W(3)
750 : : !C LEPS(1..3) [|w1-w2|<e,|w3-w2|>e];
751 : : !C [|w1-w2|<e,|w3-w2|<e];
752 : : !C [|w1-w2|>e,|w3-w2|<e]
753 : 0 : SUBROUTINE SIM3ONEI(U,AS,LEPS,SIM3,SIM3I)
754 : :
755 : : !IMPLICIT NONE
756 : : DOUBLE COMPLEX U(3)
757 : : DOUBLE PRECISION AS(3)
758 : : LOGICAL LEPS(3)
759 : : DOUBLE COMPLEX SIM3, SIM3I
760 : : !C
761 : : !DOUBLE COMPLEX SIM0UX0
762 : : !EXTERNAL SIM0UX0
763 : : DOUBLE COMPLEX X0(3), QX0(3), DX0(4)
764 : : DOUBLE COMPLEX R0(3), QR0(3)
765 : : DOUBLE COMPLEX S0(3), QS0(3)
766 : : DOUBLE COMPLEX AV, BV, CV, DV, UD, G(5)
767 : : INTEGER I, K !, N
768 : : DOUBLE COMPLEX CZERO, CONE, CTWO, CIMAG
769 : : DOUBLE PRECISION ZERO, ONE
770 : : DATA ZERO/0.0D0/, ONE/1.0D0/
771 : : !C
772 : 0 : CZERO = (0.0D0,0.0D0)
773 : 0 : CONE = (1.0D0,0.0D0)
774 : 0 : CTWO = (2.0D0,0.0D0)
775 : 0 : CIMAG = (0.0D0,1.0D0)
776 : : !C ------------------------------------------------------------------
777 : : !C R0(U) = ARTH(U)/U
778 : : !C R0(U) = 1+U**2/3+U**6*X0(U)
779 : : !C ------------------------------------------------------------------
780 [ # # ]: 0 : DO 100 I = 1, 3
781 : 0 : AV = U(I)*U(I)
782 : 0 : X0(I) = SIM0UX0(U(I))
783 : 0 : R0(I) = CONE + AV*(CONE/3.0D0+AV*(CONE/5.0D0+AV*X0(I)))
784 : : !c S0(I) = (CONE-U(I))/(CONE+U(I))*R0(I)
785 : 0 : 100 CONTINUE
786 : : !C -----------------------------------------------------------------
787 : : !C - THE DX0(I) CONTAINS DERIVATIVES OF THE FUNCTION X0(U) -
788 : : !C -----------------------------------------------------------------
789 [ # # # # : 0 : IF (LEPS(1).OR.LEPS(2).OR.LEPS(3)) THEN
# # ]
790 : 0 : CALL SIM0UDX0(U(2),X0(2),DX0)
791 : : END IF
792 : : !c DO 400 I = 1, 3, 2
793 : 0 : I = 3
794 : 0 : UD = U(I) - U(2)
795 [ # # ]: 0 : IF (.NOT.LEPS(I)) THEN
796 : 0 : QX0(I) = (X0(I)-X0(2))/UD
797 : : ELSE
798 : 0 : QX0(I) = DX0(1) + (DX0(2)/2.0D0+(DX0(3)/6.0D0+DX0(4)/24.0D0*UD) *UD)*UD
799 : : END IF
800 : 0 : AV = U(I)
801 : 0 : G(1) = U(I) + U(2)
802 [ # # ]: 0 : DO 200 K = 2, 5
803 : 0 : AV = AV*U(I)
804 : 0 : G(K) = G(K-1)*U(2) + AV
805 : 0 : 200 CONTINUE
806 : 0 : QR0(I) = G(1)/3.0D0 + G(3)/5.0D0 + G(5)*X0(2) + AV*U(I)*QX0(I)
807 : : !c 400 CONTINUE ! end DO ... I = 1, 3, 2
808 : 0 : QR0(3) = R0(2) + U(3)*QR0(3)
809 [ # # ]: 0 : IF (.NOT.LEPS(3)) THEN
810 : 0 : QR0(3) = QR0(3) + DCMPLX(ZERO,(AS(3)-AS(2)))/(U(3)-U(2))
811 : : END IF
812 : 0 : R0(2) = DCMPLX(ZERO,AS(2)) + U(2)*R0(2)
813 : 0 : AV = (ONE+U(1))*(U(2)+CONE)*(U(3)+CONE) /((U(1)*U(2)-CONE)*(U(1)*U(3)-CONE))
814 : 0 : BV = (CONE-U(1))**2
815 : 0 : CV = CTWO - (U(2)+U(3)) - U(1)*(CONE-U(2)*U(3))
816 : 0 : DV = (U(1)*U(2)-CONE)*(U(3)-CONE)**2
817 : 0 : SIM3 = AV*(BV*R0(1)+CV*R0(2)+DV*QR0(3))
818 : : !c
819 : 0 : S0(1) = (CONE-U(1))/(CONE+U(1))*R0(1)
820 : 0 : S0(2) = (U(2)-CONE)/(U(2)+CONE)*R0(2)
821 : 0 : QS0(3) = (-S0(2)+R0(2)+(U(3)-CONE)*QR0(3))/(U(3)+CONE)
822 : 0 : SIM3I = AV*(BV*S0(1)+CV*S0(2)+DV*QS0(3))
823 : : !c
824 : 0 : RETURN
825 : : END SUBROUTINE SIM3ONEI
826 : : !C
827 : : !C ----------------------------------------------------------------------
828 : : !C file: sim3twoi.f
829 : : !C date: 2010-08-19
830 : : !C who: S.Kaprzyk
831 : : !C what: CASE c) : [w1<ONE<w2<w3]
832 : : !C U(1)=W(1); U(2)=CONE/W(2); U(3)=CONE/W(3)
833 : : !C LEPS(1..3) [|w1-w2|<e,|w3-w2|>e];
834 : : !C [|w1-w2|<e,|w3-w2|<e];
835 : : !C [|w1-w2|>e,|w3-w2|<e]
836 : 2626200 : SUBROUTINE SIM3TWOI(U,AS,LEPS,SIM,SIMI)
837 : :
838 : : !IMPLICIT NONE
839 : : DOUBLE COMPLEX U(3)
840 : : DOUBLE PRECISION AS(3)
841 : : LOGICAL LEPS(3)
842 : : DOUBLE COMPLEX SIM, SIMI
843 : : !C
844 : : !DOUBLE COMPLEX SIM0UX0
845 : : !EXTERNAL SIM0UX0
846 : : !C
847 : : DOUBLE COMPLEX X0(3), QX0(3), DX0(4)
848 : : DOUBLE COMPLEX R4(3), QR4(3)
849 : : DOUBLE COMPLEX S4(3), QS4(3)
850 : : DOUBLE COMPLEX AV, BV, CV, DV, UV, G(5)
851 : : INTEGER I, K
852 : : DOUBLE COMPLEX CZERO, CONE, CIMAG, CTWO
853 : : DOUBLE PRECISION ZERO, ONE
854 : : DATA ZERO/0.0D0/, ONE/1.0D0/
855 : : !C
856 : 2626200 : CZERO = (0.0D0,0.0D0)
857 : 2626200 : CIMAG = (0.0D0,1.0D0)
858 : 2626200 : CONE = (1.0D0,0.0D0)
859 : 2626200 : CTWO = (2.0D0,0.0D0)
860 : : !C -----------------------------------------------------------------
861 : : !C - -
862 : : !C - R4(W)=1+(W-1)*(ARTH(W)/W-1)/W -
863 : : !C - R4(W)=1-W/3+W**2/3-W**3/5+W**4/5+(W-1)*W**5*X0(W) -
864 : : !C - -
865 : : !C -----------------------------------------------------------------
866 [ + + ]: 10504800 : DO 100 I = 1, 3
867 : 7878600 : AV = U(I)
868 : 7878600 : X0(I) = SIM0UX0(AV)
869 : : !*-ASIS
870 : : R4(I) = CONE + AV*(-CONE/3.0D0 + AV*(CONE/3.0D0 + &
871 : 7878600 : & AV*((-0.2D0,0.0D0)+AV*((0.2D0,0.0D0)+AV*(AV-CONE)*X0(I)))))
872 : 2626200 : 100 CONTINUE ! end DO 350 I = 1, 3
873 : : !C -----------------------------------------------------------------
874 : : !C - The DX0(I) Contains Derivatives Of The Function X0(W) -
875 : : !C -----------------------------------------------------------------
876 [ + - + - : 2626200 : IF (LEPS(1).OR.LEPS(2).OR.LEPS(3)) THEN
+ + ]
877 : 513524 : CALL SIM0UDX0(U(2),X0(2),DX0)
878 : : END IF
879 : : !c DO 750 I = 1, 3, 2
880 : 2626200 : I = 3
881 : 2626200 : UV = U(I) - U(2)
882 [ + + ]: 2626200 : IF (.NOT.LEPS(I)) THEN
883 : 2112676 : QX0(I) = (X0(I)-X0(2))/UV
884 : : ELSE
885 : 513524 : QX0(I) = DX0(1) + (DX0(2)/2.0D0+(DX0(3)/6.0D0+DX0(4)/24.0D0*UV) *UV)*UV
886 : : END IF
887 : 2626200 : AV = U(I)
888 : 2626200 : G(1) = U(I) + U(2)
889 [ + + ]: 13131000 : DO 200 K = 2, 5
890 : 10504800 : AV = AV*U(I)
891 : 10504800 : G(K) = G(K-1)*U(2) + AV
892 : 2626200 : 200 CONTINUE
893 : : !*-ASIS
894 : : QR4(I) = -CONE/3.0D0 + G(1)/3.0D0 - G(2)/5.0D0 + G(3)/5.0D0 + &
895 : 2626200 : & (G(5)-G(4))*X0(2) + (U(I)-CONE)*AV*QX0(I)
896 : : !C 750 CONTINUE ! DO .. I = 1, 3, 2
897 : 2626200 : QR4(3) = -DCMPLX(ZERO,AS(2)) + CONE - (U(3)+U(2))*R4(2) - U(3)*U(3)*QR4(3)
898 [ + + ]: 2626200 : IF (.NOT.LEPS(3)) THEN
899 : 2112676 : QR4(3) = QR4(3) + (CONE-U(3))*DCMPLX(ZERO,(AS(3)-AS(2))) /(U(3)-U(2))
900 : : END IF
901 : 2626200 : R4(2) = (CONE-U(2))*DCMPLX(ZERO,AS(2)) + CONE + U(2) - U(2)*U(2)*R4(2)
902 : : !c
903 : 2626200 : AV = (ONE+U(1))*(U(2)+CONE)*(U(3)+CONE)/((U(1)*U(2)-CONE)*(U(1)*U(3)-CONE))
904 : 2626200 : BV = (CONE-U(1))**2
905 : 2626200 : CV = CTWO - (U(2)+U(3)) - U(1)*(CONE-U(2)*U(3))
906 : 2626200 : DV = (U(1)*U(2)-CONE)*(U(3)-CONE)**2
907 : 2626200 : SIM = AV*(BV*R4(1)+CV*R4(2)+DV*QR4(3))
908 : : !c
909 : 2626200 : S4(1) = (CONE-U(1))/(CONE+U(1))*R4(1)
910 : 2626200 : S4(2) = (U(2)-CONE)/(U(2)+CONE)*R4(2)
911 : 2626200 : QS4(3) = (-S4(2)+R4(2)+(U(3)-CONE)*QR4(3))/(U(3)+CONE)
912 : 2626200 : SIMI = AV*(BV*S4(1)+CV*S4(2)+DV*QS4(3))
913 : : !c
914 : 2626200 : RETURN
915 : : END SUBROUTINE SIM3TWOI
916 : : !C
917 : : !C ----------------------------------------------------------------------
918 : : !C file: sim4onei.f
919 : : !C date: 2010-09-22
920 : : !C who: S.Kaprzyk
921 : : !C what: CASE d): [ ONE<w1<w2<w3 ]
922 : : !C U(1)=CONE/W(1); U(2)=CONE/W(2); U(3)=CONE/W(3)
923 : : !C LEPS(1..3) [|w1-w2|<e,|w3-w2|>e];
924 : : !C [|w1-w2|<e,|w3-w2|<e];
925 : : !C [|w1-w2|>e,|w3-w2|<e]
926 : 0 : SUBROUTINE SIM4ONEI(U,AS,LEPS,SIM4,SIM4I)
927 : :
928 : : !IMPLICIT NONE
929 : : DOUBLE COMPLEX U(3)
930 : : DOUBLE PRECISION AS(3)
931 : : LOGICAL LEPS(3)
932 : : DOUBLE COMPLEX SIM4, SIM4I
933 : : !c
934 : : !DOUBLE COMPLEX SIM0UX0
935 : : !EXTERNAL SIM0UX0
936 : : DOUBLE COMPLEX X0(3), QX0(3), DX0(4)
937 : : DOUBLE COMPLEX R0(3), QR0(3)
938 : : DOUBLE COMPLEX S0(3), QS0(3)
939 : : DOUBLE COMPLEX QAS(3)
940 : : DOUBLE COMPLEX AV, BV, CV, UD, G(5)
941 : : INTEGER I, K !, N
942 : : DOUBLE COMPLEX CZERO, CONE, CTWO, CIMAG
943 : : DOUBLE PRECISION ZERO, ONE
944 : : DATA ZERO/0.0D0/, ONE/1.0D0/
945 : : !C
946 : 0 : CZERO = (0.0D0,0.0D0)
947 : 0 : CONE = (1.0D0,0.0D0)
948 : 0 : CTWO = (2.0D0,0.0D0)
949 : 0 : CIMAG = (0.0D0,1.0D0)
950 : : !C ------------------------------------------------------------------
951 : : !C R0(U)=ARTH(U)/U
952 : : !C R0(U)=1+U**2/3+U**6*X0(U)
953 : : !C ------------------------------------------------------------------
954 [ # # ]: 0 : DO 100 I = 1, 3
955 : 0 : AV = U(I)*U(I)
956 : 0 : X0(I) = SIM0UX0(U(I))
957 : 0 : R0(I) = CONE + AV*(CONE/3.0D0+AV*(CONE/5.0D0+AV*X0(I)))
958 : : !c S0(I) = (CONE-U(I))/(CONE+U(I))*R0(I)
959 : 0 : 100 CONTINUE
960 : : !C -----------------------------------------------------------------
961 : : !C - THE DX0(I) CONTAINS DERIVATIVES OF THE FUNCTION X0(U)
962 : : !C -----------------------------------------------------------------
963 [ # # # # : 0 : IF (LEPS(1).OR.LEPS(2).OR.LEPS(3)) THEN
# # ]
964 : 0 : CALL SIM0UDX0(U(2),X0(2),DX0)
965 : : END IF
966 : 0 : DO 200 I = 1, 3, 2
967 : 0 : UD = U(I) - U(2)
968 [ # # ]: 0 : IF (.NOT.LEPS(I)) THEN
969 : 0 : QX0(I) = (X0(I)-X0(2))/UD
970 : 0 : QAS(I) = (AS(I)-AS(2))/UD
971 : : ELSE
972 : 0 : QX0(I) = DX0(1) + (DX0(2)/2.0D0+(DX0(3)/6.0D0+DX0(4)/24.0D0*UD) *UD)*UD
973 : 0 : QAS(I) = CZERO
974 : : END IF
975 : 0 : AV = U(I)
976 : 0 : G(1) = U(I) + U(2)
977 [ # # ]: 0 : DO 150 K = 2, 5
978 : 0 : AV = AV*U(I)
979 : 0 : G(K) = G(K-1)*U(2) + AV
980 : 0 : 150 CONTINUE
981 [ # # ]: 0 : QR0(I) = G(1)/3.0D0 + G(3)/5.0D0 + G(5)*X0(2) + AV*U(I)*QX0(I)
982 : 0 : 200 CONTINUE ! end DO ... I = 1, 3, 2
983 [ # # ]: 0 : IF (LEPS(2)) THEN
984 : 0 : QAS(2) = CZERO
985 : : !*-ASIS
986 : : QX0(2) = DX0(2)/2.0D0 + DX0(3)*(U(1)+U(3)-2.0D0*U(2))/6.0D0 + &
987 : : & DX0(4)*((U(3)-U(2))*(U(3)-U(2))+(U(1)-U(2))*(U(1)-U(2)) + &
988 : 0 : & (U(3)-U(2))*(U(1)-U(2)))/24.0D0
989 : : ELSE
990 : 0 : QAS(2) = (QAS(3)-QAS(1))/(U(3)-U(1))
991 : 0 : QX0(2) = (QX0(3)-QX0(1))/(U(3)-U(1))
992 : : END IF
993 : 0 : AV = U(1)
994 : 0 : G(1) = U(1) + U(3)
995 [ # # ]: 0 : DO 300 K = 2, 5
996 : 0 : AV = AV*U(1)
997 : 0 : G(K) = G(K-1)*U(3) + AV
998 : 0 : 300 CONTINUE
999 : : !*-ASIS
1000 : : QR0(2) = CONE/3.0D0 + (G(2)+U(2)*(U(2)+G(1)))/5.0D0 + &
1001 : : & (G(4)+U(2)*(G(3)+U(2)*(G(2)+U(2)*(G(1)+U(2)))))*X0(2)+ &
1002 : 0 : & G(5)*QX0(3) + AV*U(1)*QX0(2)
1003 : : !c
1004 : 0 : QR0(2) = CIMAG*QAS(2) + QR0(3) + U(1)*QR0(2)
1005 : 0 : QR0(3) = CIMAG*QAS(3) + R0(2) + U(3)*QR0(3)
1006 : 0 : R0(2) = CIMAG*AS(2) + U(2)*R0(2)
1007 : : !c
1008 : 0 : S0(2) = (U(2)-CONE)/(U(2)+CONE)*R0(2)
1009 : 0 : QS0(3) = (-S0(2)+R0(2)+(U(3)-CONE)*QR0(3))/(U(3)+CONE)
1010 : 0 : QS0(2) = (-QS0(3)+QR0(3)+(U(1)-CONE)*QR0(2))/(U(1)+CONE)
1011 : : !c
1012 : 0 : AV = (CONE+U(1))*(CONE+U(2))*(CONE+U(3))
1013 : 0 : BV = U(3) + U(1) - CTWO
1014 : 0 : CV = (CONE-U(1))**2
1015 : 0 : SIM4 = AV*(R0(2)+BV*QR0(3)+CV*QR0(2))
1016 : 0 : SIM4I = AV*(S0(2)+BV*QS0(3)+CV*QS0(2))
1017 : : !c
1018 : 0 : RETURN
1019 : : END SUBROUTINE SIM4ONEI
1020 : : !C
1021 : : !C ----------------------------------------------------------------------
1022 : : !C file: sim4twoi.f
1023 : : !C date: 2010-09-22
1024 : : !C who: S.Kaprzyk
1025 : : !C what: CASE d): [ ONE<w1<w2<w3 ]
1026 : : !C U(1)=CONE/W(1); U(2)=CONE/W(2); U(3)=CONE/W(3)
1027 : : !C LEPS(1..3) [|w1-w2|<e,|w3-w2|>e];
1028 : : !C [|w1-w2|<e,|w3-w2|<e];
1029 : : !C [|w1-w2|>e,|w3-w2|<e]
1030 : 2078173 : SUBROUTINE SIM4TWOI(U,AS,LEPS,SIM,SIMI)
1031 : :
1032 : : !IMPLICIT NONE
1033 : : DOUBLE COMPLEX U(3)
1034 : : DOUBLE PRECISION AS(3)
1035 : : LOGICAL LEPS(3)
1036 : : DOUBLE COMPLEX SIM, SIMI
1037 : : !C
1038 : : !DOUBLE COMPLEX SIM0UX0
1039 : : !EXTERNAL SIM0UX0
1040 : : !c
1041 : : DOUBLE COMPLEX X0(3), QX0(3), DX0(4)
1042 : : DOUBLE COMPLEX R4(3), QR4(3)
1043 : : DOUBLE COMPLEX S4(3), QS4(3)
1044 : : !c DOUBLE COMPLEX RU(3), TU(3)
1045 : : DOUBLE COMPLEX QAS(3), AV, BV, CV, UV, G(5)
1046 : : INTEGER I, K
1047 : : DOUBLE COMPLEX CZERO, CONE, CTWO, CIMAG
1048 : : DOUBLE PRECISION ZERO, ONE
1049 : : DATA ZERO/0.0D0/, ONE/1.0D0/
1050 : : !C
1051 : 2078173 : CZERO = (0.0D0,0.0D0)
1052 : 2078173 : CONE = (1.0D0,0.0D0)
1053 : 2078173 : CTWO = (2.0D0,0.0D0)
1054 : 2078173 : CIMAG = (0.0D0,1.0D0)
1055 : : !C -----------------------------------------------------------------
1056 : : !C - -
1057 : : !C - R4(U)=1+(U-1)*(ARTH(U)/U-1)/U -
1058 : : !C - R4(U)=1-U/3+U**2/3-U**3/5+U**4/5+(U-1)*U**5*X0(U) -
1059 : : !C - -
1060 : : !C -----------------------------------------------------------------
1061 [ + + ]: 8312692 : DO 100 I = 1, 3
1062 : 6234519 : AV = U(I)
1063 : 6234519 : X0(I) = SIM0UX0(AV)
1064 : : !*-ASIS
1065 : : R4(I) = CONE + AV*(-CONE/3.0D0 + AV*(CONE/3.0D0 + &
1066 : 6234519 : & AV*((-0.2D0,0.0D0)+AV*((0.2D0,0.0D0)+AV*(AV-CONE)*X0(I)))))
1067 : 2078173 : 100 CONTINUE ! end DO ... I = 1, 3
1068 : : !C -----------------------------------------------------------------
1069 : : !C - The DX0(I) Contains Derivatives Of The Function X0(U) -
1070 : : !C -----------------------------------------------------------------
1071 [ + + + - : 2078173 : IF (LEPS(1).OR.LEPS(2).OR.LEPS(3)) THEN
+ + ]
1072 : 710748 : CALL SIM0UDX0(U(2),X0(2),DX0)
1073 : : END IF
1074 : : !C
1075 : 4156346 : DO 200 I = 1, 3, 2
1076 : 4156346 : UV = U(I) - U(2)
1077 [ + + ]: 4156346 : IF (LEPS(I)) THEN
1078 : 710748 : QX0(I) = DX0(1) + (DX0(2)/2.0D0+(DX0(3)/6.0D0+DX0(4)/24.0D0*UV) *UV)*UV
1079 : 710748 : QAS(I) = CZERO
1080 : : ELSE
1081 : 3445598 : QX0(I) = (X0(I)-X0(2))/UV
1082 : 3445598 : QAS(I) = (AS(I)-AS(2))/UV
1083 : : END IF
1084 : 4156346 : AV = U(I)
1085 : 4156346 : G(1) = U(I) + U(2)
1086 [ + + ]: 20781730 : DO 150 K = 2, 5
1087 : 16625384 : AV = AV*U(I)
1088 : 16625384 : G(K) = G(K-1)*U(2) + AV
1089 : 4156346 : 150 CONTINUE
1090 : : QR4(I) = -CONE/3.0D0 + G(1)/3.0D0 - G(2)/5.0D0 + G(3) &
1091 [ + + ]: 4156346 : & /5.0D0 + (G(5)-G(4))*X0(2) + (U(I)-CONE)*AV*QX0(I)
1092 : 2078173 : 200 CONTINUE ! end DO ... I = 1, 3, 2
1093 : : !C
1094 [ - + ]: 2078173 : IF (LEPS(2)) THEN
1095 : 0 : QAS(2) = CZERO
1096 : : !*-ASIS
1097 : : QX0(2) = DX0(2)/2.0D0 + DX0(3)*(U(1)+U(3)-2.0D0*U(2))/6.0D0 + &
1098 : : & DX0(4)*((U(3)-U(2))*(U(3)-U(2))+(U(1)-U(2))*(U(1)-U(2)) + &
1099 : 0 : & (U(3)-U(2))*(U(1)-U(2)))/24.0D0
1100 : : ELSE
1101 : 2078173 : QAS(2) = (QAS(3)-QAS(1))/(U(3)-U(1))
1102 : 2078173 : QX0(2) = (QX0(3)-QX0(1))/(U(3)-U(1))
1103 : : END IF
1104 : 2078173 : AV = U(1)
1105 : 2078173 : G(1) = U(1) + U(3)
1106 [ + + ]: 10390865 : DO 300 K = 2, 5
1107 : 8312692 : AV = AV*U(1)
1108 : 8312692 : G(K) = G(K-1)*U(3) + AV
1109 : 2078173 : 300 CONTINUE
1110 : : !*-ASIS
1111 : : QR4(2) = CONE/3.0D0- (U(2)+G(1))/5.0D0 + &
1112 : : & (G(2)+U(2)*(U(2)+G(1)))/5.0D0 + &
1113 : : & (G(4)+(U(2)-CONE)*(G(3) + G(2)*U(2) + G(1)*U(2)*U(2) + &
1114 : : & U(2)*U(2)*U(2)))*X0(2) + (G(5)-G(4))*QX0(3) + &
1115 : 2078173 : & (U(1)-CONE)*AV*QX0(2)
1116 : : !C
1117 : : QR4(2) = CIMAG*(-QAS(3)+(CONE-U(1))*QAS(2)) &
1118 : 2078173 : & - (R4(2)+(U(3)+U(1))*QR4(3)+U(1)*U(1)*QR4(2))
1119 : : QR4(3) = CIMAG*(-AS(2)+(CONE-U(3))*QAS(3)) + CONE - (U(3)+U(2)) &
1120 : 2078173 : & *R4(2) - U(3)*U(3)*QR4(3)
1121 : 2078173 : R4(2) = CIMAG*(CONE-U(2))*AS(2) + CONE + U(2) - U(2)*U(2)*R4(2)
1122 : : !c
1123 : 2078173 : S4(2) = (U(2)-CONE)/(U(2)+CONE)*R4(2)
1124 : 2078173 : QS4(3) = (-S4(2)+R4(2)+(U(3)-CONE)*QR4(3))/(U(3)+CONE)
1125 : 2078173 : QS4(2) = (-QS4(3)+QR4(3)+(U(1)-CONE)*QR4(2))/(U(1)+CONE)
1126 : : !c
1127 : 2078173 : AV = (CONE+U(1))*(CONE+U(2))*(CONE+U(3))
1128 : 2078173 : BV = U(3) + U(1) - CTWO
1129 : 2078173 : CV = (CONE-U(1))**2
1130 : 2078173 : SIM = AV*(R4(2)+BV*QR4(3)+CV*QR4(2))
1131 : 2078173 : SIMI = AV*(S4(2)+BV*QS4(3)+CV*QS4(2))
1132 : : !c
1133 : 2078173 : RETURN
1134 : : END SUBROUTINE SIM4TWOI
1135 : : !C
1136 : : !C ----------------------------------------------------------------------
1137 : : !C file: sim0leps.f
1138 : : !C date: 2011-03-12
1139 : : !C who: S.Kaprzyk
1140 : : !C what: Return a proper case a), b), c) or d): see, Fig.3.1
1141 : : !C INPUT:
1142 : : !C VERM(1..4) - four complex numbers
1143 : : !C N - integer number for W(I)=(VERM(I)-VERM(N))/(VERM(I)+VERM(N))
1144 : : !C EPS - critical distance |w1-w2| etc.
1145 : : !C OUTPUT:
1146 : : !C LONE(1..4) [w1<w2<w3<ONE]; [w1<w2<ONE<w3];
1147 : : !C [w1<ONE<w2<w3]; [ONE<w1<w2<w3]
1148 : : !C LEPS(1..3) [|w1-w2|<e,|w3-w2|>e];
1149 : : !C [|w1-w2|<e,|w3-w2|<e];
1150 : : !C [|w1-w2|>e,|w3-w2|<e]
1151 : 220656384 : SUBROUTINE SIM0LEPS(VERM, N, W, AS, LONE, EPS, LEPS, iuerr)
1152 : :
1153 : : !IMPLICIT NONE
1154 : : DOUBLE COMPLEX VERM(4)
1155 : : INTEGER N
1156 : : DOUBLE COMPLEX W(3)
1157 : : DOUBLE PRECISION AS(3)
1158 : : DOUBLE PRECISION EPS
1159 : : LOGICAL LONE(4), LEPS(3)
1160 : : INTEGER iuerr
1161 : : !C
1162 : : !c DOUBLE PRECISION DLAMCH
1163 : : !c EXTERNAl DLAMCH
1164 : : !C
1165 : : LOGICAL LWONE(3), LV
1166 : : DOUBLE COMPLEX AV
1167 : : DOUBLE PRECISION AL(3), AA, AIW, OZERO
1168 : : INTEGER I, II, K
1169 : : DOUBLE PRECISION ZERO, ONE, PI
1170 : : DATA ZERO/0.0D0/,ONE/1.0D0/,PI/3.141592653589793D0/
1171 : : DATA OZERO/ 1.0D-13/
1172 : : !C
1173 : : !c OZERO = DLAMCH('e')*10
1174 [ - + ]: 220656384 : IF ((N.GT.4).OR.(N.LT.1)) THEN
1175 : 0 : WRITE(iuerr,9003) N
1176 : : 9003 FORMAT(' ***sim0leps: N =',I6,' must be 1,2,3 or 4' )
1177 : 0 : STOP ' ***sim0leps: '
1178 : : ENDIF
1179 : : !c DO 100 I = 1, 4
1180 : : !c IF (I.EQ.N) GO TO 100
1181 : : !c IF (DIMAG(VERM(N))*DIMAG(VERM(I)).LT.ZERO) THEN
1182 : : !c WRITE (iuerr,9010) (DIMAG(VERM(K)),K=1,4)
1183 : : !c 9010 FORMAT (' ***sim0leps: not signed ImgVERM()=',4(D13.6,1X))
1184 : : !c STOP ' ***sim0leps: '
1185 : : !c END IF
1186 : : !c 100 CONTINUE
1187 : : II = 0
1188 [ + + ]: 1103281920 : DO 200 I = 1, 4
1189 [ + + ]: 882625536 : IF (I.EQ.N) GO TO 200
1190 : 661969152 : II = II + 1
1191 [ + + ]: 661969152 : IF (CDABS(VERM(N)-VERM(I)).LT.CDABS(VERM(N)+VERM(I))) THEN
1192 : 644144754 : LWONE(II) = .TRUE.
1193 : 644144754 : W(II) = (VERM(N)-VERM(I))/(VERM(N)+VERM(I))
1194 : 644144754 : AIW = DIMAG(W(II))
1195 [ + + ]: 644144754 : IF (DABS(AIW).GE.OZERO) THEN
1196 : 588583794 : AS(II) = 0.5D0*PI*DSIGN(ONE,AIW)
1197 : : ELSE
1198 : 55560960 : AS(II) = 0.5D0*PI*DSIGN(ONE,DREAL(VERM(I)-VERM(N)))
1199 : : END IF
1200 : : ELSE
1201 : 17824398 : LWONE(II) = .FALSE.
1202 : 17824398 : W(II) = (VERM(N)+VERM(I))/(VERM(N)-VERM(I))
1203 : 17824398 : AIW = DIMAG(W(II))
1204 [ + - ]: 17824398 : IF (DABS(AIW).GE.OZERO) THEN
1205 : 17824398 : AS(II) = -0.5D0*PI*DSIGN(ONE,AIW)
1206 : : ELSE
1207 : 0 : AS(II) = 0.5D0*PI*DSIGN(ONE,DREAL(VERM(I)-VERM(N)))
1208 : : END IF
1209 : : END IF
1210 : 882625536 : AL(II) = CDABS(W(II))
1211 : 220656384 : 200 CONTINUE
1212 [ + + ]: 882625536 : DO 300 I = 1, 3
1213 [ + + ]: 1985907456 : DO 250 K = I, 3
1214 [ + + + + : 1323938304 : IF (LWONE(I).AND.LWONE(K).AND.(AL(I).LT.AL(K))) GO TO 250
+ + ]
1215 [ + + + + : 1039231323 : IF (.NOT.LWONE(I).AND..NOT.LWONE(K).AND.(AL(K).LT.AL(I))) GO TO 250
+ + ]
1216 [ + + + + : 1035147868 : IF (.NOT.LWONE(I).AND.LWONE(K).AND.(ONE.LT.AL(I)*AL(K))) GO TO 250
- + ]
1217 [ + + + + : 1035147868 : IF (LWONE(I).AND..NOT.LWONE(K).AND.(AL(I)*AL(K).LT.ONE)) GO TO 250
+ - ]
1218 : 1025626548 : AA = AL(K)
1219 : 1025626548 : AL(K) = AL(I)
1220 : 1025626548 : AL(I) = AA
1221 : 1025626548 : AA = AS(K)
1222 : 1025626548 : AS(K) = AS(I)
1223 : 1025626548 : AS(I) = AA
1224 : 1025626548 : AV = W(K)
1225 : 1025626548 : W(K) = W(I)
1226 : 1025626548 : W(I) = AV
1227 : 1025626548 : LV = LWONE(K)
1228 : 1025626548 : LWONE(K) = LWONE(I)
1229 : 1323938304 : LWONE(I) = LV
1230 : 661969152 : 250 CONTINUE
1231 : 220656384 : 300 CONTINUE
1232 : : !C
1233 : 220656384 : LONE(1) = .FALSE.
1234 : 220656384 : LONE(2) = .FALSE.
1235 : 220656384 : LONE(3) = .FALSE.
1236 : 220656384 : LONE(4) = .FALSE.
1237 [ + + + + : 220656384 : IF (LWONE(1).AND.LWONE(2).AND.LWONE(3)) THEN
+ + ]
1238 : 209614532 : LONE(1) = .TRUE.
1239 : : END IF
1240 [ + + + + : 220656384 : IF (LWONE(1).AND.LWONE(2).AND.(.NOT.LWONE(3))) THEN
+ + ]
1241 : 6337479 : LONE(2) = .TRUE.
1242 : : END IF
1243 [ + + + + : 220656384 : IF (LWONE(1).AND.(.NOT.LWONE(2)).AND.(.NOT.LWONE(3))) THEN
+ - ]
1244 : 2626200 : LONE(3) = .TRUE.
1245 : : END IF
1246 [ + + + - : 220656384 : IF ((.NOT.LWONE(1)).AND.(.NOT.LWONE(2)).AND.(.NOT.LWONE(3))) THEN
+ - ]
1247 : 2078173 : LONE(4) = .TRUE.
1248 : : END IF
1249 : : !C here are 3-cases, how close are w()
1250 : 220656384 : LEPS(1) = .FALSE.
1251 : 220656384 : LEPS(2) = .FALSE.
1252 : 220656384 : LEPS(3) = .FALSE.
1253 [ + + + + ]: 220656384 : IF ( LONE(1).OR.LONE(4) ) THEN
1254 [ + + ]: 211692705 : IF (ABS(W(1)-W(2)).LE.EPS) THEN
1255 : 23298912 : LEPS(1) = .TRUE.
1256 : : END IF
1257 [ + + ]: 211692705 : IF (ABS(W(3)-W(2)).LE.EPS) THEN
1258 : 31031632 : LEPS(3) = .TRUE.
1259 : : END IF
1260 [ + + - + ]: 211692705 : IF ((ABS(W(1)-W(2)).LE.EPS).AND.(ABS(W(3)-W(2)).LE.EPS)) THEN
1261 : 0 : LEPS(2) = .TRUE.
1262 : : END IF
1263 : : END IF
1264 : : !c
1265 [ + + ]: 220656384 : IF ( LONE(2)) THEN
1266 [ + + ]: 6337479 : IF (ABS(W(1)-W(2)).LE.EPS) THEN
1267 : 716892 : LEPS(1) = .TRUE.
1268 : : END IF
1269 : : END IF
1270 : : !c
1271 [ + + ]: 220656384 : IF ( LONE(3)) THEN
1272 [ + + ]: 2626200 : IF (ABS(W(2)-W(3)).LE.EPS) THEN
1273 : 513524 : LEPS(3) = .TRUE.
1274 : : END IF
1275 : : END IF
1276 : : !C
1277 [ + + + + : 220656384 : IF ((.NOT.LONE(1)).AND.(.NOT.LONE(2)).AND.(.NOT.LONE(3)).AND. (.NOT.LONE(4))) THEN
+ + - + ]
1278 : 0 : WRITE (iuerr,9020) (LONE(I),I=1,4)
1279 : : 9020 FORMAT ('***sim0leps: LONE()=',4(L1,1X),'not in order')
1280 : 0 : STOP '***sim0leps: '
1281 : : END IF
1282 : 220656384 : RETURN
1283 : : END SUBROUTINE SIM0LEPS
1284 : : !C
1285 : : !C file: s2d0onei.f
1286 : : !C date: 2011-04-01
1287 : : !C who: S.Kaprzyk
1288 : : !C what: Complex linear form integral over standard 2-d simplex
1289 : : !C ----------------------------------------------------
1290 : : !C - * -
1291 : : !C - * 1 -
1292 : : !C - SIM0= 2* dt1dt2 ----------------------------------
1293 : : !C - * verm(3)+(verm(1)-verm(3))*t1+.. -
1294 : : !C - * -
1295 : : !C - 0<t1<1 -
1296 : : !C ----------------------------------------------------
1297 : 0 : SUBROUTINE S2D0ONEI(SIM0, SIM0I, VERM)
1298 : :
1299 : : !IMPLICIT NONE
1300 : : INTEGER iuerr
1301 : : PARAMETER (iuerr=6)
1302 : : DOUBLE COMPLEX SIM0, SIM0I
1303 : : DOUBLE COMPLEX VERM (3)
1304 : : !C
1305 : : DOUBLE COMPLEX VERM1D(2)
1306 : : DOUBLE COMPLEX SIM, SIMI
1307 : : LOGICAL LONE(3),LEPS(1)
1308 : : !c
1309 : : DOUBLE COMPLEX U(2)
1310 : : DOUBLE PRECISION AS(2), AL(3)
1311 : : INTEGER I, I1, I2, K, N
1312 : : DOUBLE PRECISION ZERO, EPS, SMALL
1313 : : DATA EPS/1.0D-6/, SMALL/1.0D-5/
1314 : : DATA ZERO/0.0D0/
1315 : : !C
1316 [ # # ]: 0 : DO 100 I = 1, 2
1317 [ # # ]: 0 : IF (DIMAG(VERM(3))*DIMAG(VERM(I)).LT.ZERO) THEN
1318 [ # # # # ]: 0 : WRITE (iuerr,9010) (DIMAG(VERM(K)),K=1,3)
1319 : : 9010 FORMAT (' ***s2d0onei: not signed ImgVERM()=',3(D13.6,1X))
1320 : : !c STOP ' ***s2d0onei: '
1321 : : END IF
1322 : 0 : 100 CONTINUE
1323 : : !C
1324 [ # # ]: 0 : DO 200 I = 1, 3
1325 : 0 : AL(I) = CDABS(VERM(I))
1326 : 0 : 200 CONTINUE
1327 [ # # ]: 0 : DO 300 N = 1, 3
1328 [ # # ]: 0 : IF (AL(N).LT.SMALL) THEN
1329 : 0 : I1 = MOD(N,3) + 1
1330 : 0 : I2 = MOD(N+1,3) + 1
1331 : 0 : VERM1D(1) = VERM(I1)
1332 : 0 : VERM1D(2) = VERM(I2)
1333 : 0 : CALL S1D0ONEI(SIM,SIMI,VERM1D)
1334 : 0 : SIM0 = 2*SIM
1335 : 0 : SIM0I = SIMI - 0.5D0
1336 : 0 : RETURN
1337 : : END IF
1338 : 0 : 300 CONTINUE
1339 : : !C
1340 : 0 : N = 3
1341 [ # # ]: 0 : DO 400 I = 1, 2
1342 [ # # ]: 0 : IF (AL(I).GT.AL(N)) N = I
1343 : 0 : 400 CONTINUE
1344 : : !C Here are 3 cases, how w() are placed on complex-plane
1345 : : !C LONE(1..3) [w1<w2<ONE];[w1<ONE<w2];[ONE<w1<w2]
1346 : : !C LEPS(1) [|w1-w2|<e
1347 : : !*-ASIS
1348 : 0 : CALL S2D0LEPS(VERM, N, U, AS, LONE, EPS, LEPS, iuerr)
1349 : : !c U(1) = W(1); U(2) = W(2);
1350 [ # # ]: 0 : IF (LONE(1)) THEN
1351 : 0 : CALL S2D1ONEI(U,AS,LEPS,SIM,SIMI)
1352 : : END IF
1353 : : !c U(1) = W(1); U(2) = 1/W(2);
1354 [ # # ]: 0 : IF (LONE(2)) THEN
1355 : 0 : CALL S2D2ONEI(U,AS,LEPS,SIM,SIMI)
1356 : : END IF
1357 : : !c U(1) = 1/W(1); U(2) = 1/W(2);
1358 [ # # ]: 0 : IF (LONE(3)) THEN
1359 : 0 : CALL S2D3ONEI(U,AS,LEPS,SIM,SIMI)
1360 : : END IF
1361 : : !C
1362 : 0 : SIM0 = SIM/VERM(N)
1363 : 0 : SIM0I = log(VERM(N))+0.5D0*SIMI - 1.5D0
1364 : 0 : RETURN
1365 : : END SUBROUTINE S2D0ONEI
1366 : : !C
1367 : : !C ----------------------------------------------------------------------
1368 : : !C file: s2d0twoi.f
1369 : : !C date: 2011-04-01
1370 : : !C who: S.Kaprzyk
1371 : : !C what: Complex linear form integrals over standard 2-d simplex
1372 : : !C ------------------(i=1,2)----------------------------------------
1373 : : !C * -
1374 : : !C * t_i -
1375 : : !C VERL(i)=2*dt1dt2 --------------------------------------------- -
1376 : : !C * VERM(3)+[VERM(1)-VERM(3)]*t1+... -
1377 : : !C * -
1378 : : !C 0<t1+t2<1 -
1379 : : !C -------------------(i=3)-----------------------------------------
1380 : : !C * -
1381 : : !C * 1- t1 -t2
1382 : : !C VERL(3)=2*dt1dt2 -----------------------------------------------
1383 : : !C * VERM(3)+[VERM(1)-VERM(3)]*t1+... -
1384 : : !C * -
1385 : : !C 0<t1+t2<1 -
1386 : : !C -----------------------------------------------------------------
1387 : 0 : SUBROUTINE S2D0TWOI(VERL, VERLI, VERM)
1388 : :
1389 : : !IMPLICIT NONE
1390 : : INTEGER iuerr
1391 : : PARAMETER (iuerr=6)
1392 : : DOUBLE COMPLEX VERL(3), VERLI(3), VERM(3)
1393 : : !C
1394 : : DOUBLE COMPLEX SIM, SIMI
1395 : : DOUBLE COMPLEX VERL1D(2), VERL1DI(2), VERM1D(2)
1396 : : LOGICAL LONE(3),LEPS(2)
1397 : : !c
1398 : : INTEGER I, I1, I2, K, N
1399 : : DOUBLE PRECISION AS(2), AL(3)
1400 : : DOUBLE COMPLEX U(2)
1401 : : DOUBLE COMPLEX CZERO
1402 : : DOUBLE PRECISION EPS, SMALL
1403 : : DOUBLE PRECISION ZERO
1404 : : DATA EPS/1.0D-6/, SMALL/1.0D-5/
1405 : : DATA ZERO/0.0D0/
1406 : : !C
1407 : 0 : CZERO = (0.0D0,0.0D0)
1408 [ # # ]: 0 : DO 100 I = 1, 2
1409 [ # # ]: 0 : IF (DIMAG(VERM(3))*DIMAG(VERM(I)).LT.ZERO) THEN
1410 [ # # # # ]: 0 : WRITE (iuerr,9010) (DIMAG(VERM(K)),K=1,3)
1411 : : 9010 FORMAT (' ***s2d0twoi: not signed ImgVERM()=',3(D13.6,1X))
1412 : 0 : STOP ' ***s2d0twoi: '
1413 : : END IF
1414 : 0 : 100 CONTINUE
1415 [ # # ]: 0 : DO 200 I = 1, 3
1416 : 0 : AL(I) = CDABS(VERM(I))
1417 : 0 : 200 CONTINUE
1418 : : !C
1419 [ # # ]: 0 : DO 300 N = 1, 3
1420 : 0 : VERL(N) = CZERO
1421 : 0 : VERLI(N) = CZERO
1422 : 0 : I1 = MOD(N,3) + 1
1423 : 0 : I2 = MOD(N+1,3) + 1
1424 [ # # ]: 0 : IF (AL(I1).LT.SMALL) THEN
1425 : 0 : VERM1D(1) = VERM(I2)
1426 : 0 : VERM1D(2) = VERM(N)
1427 : 0 : CALL S1D0TWOI(VERL1D,VERL1DI,VERM1D)
1428 : 0 : VERL(N) = VERL1D(2)
1429 : 0 : VERLI(N) = VERL1DI(2)*2/3
1430 : 0 : GO TO 300
1431 : : END IF
1432 [ # # ]: 0 : IF (AL(I2).LT.SMALL) THEN
1433 : 0 : VERM1D(1) = VERM(I1)
1434 : 0 : VERM1D(2) = VERM(N)
1435 : 0 : CALL S1D0TWOI(VERL1D,VERL1DI,VERM1D)
1436 : 0 : VERL(N) = VERL1D(2)
1437 : 0 : VERLI(N) = VERL1DI(2)*2/3
1438 : 0 : GO TO 300
1439 : : END IF
1440 [ # # ]: 0 : IF (AL(N).LT.SMALL) THEN
1441 : 0 : VERM1D(1) = VERM(I1)
1442 : 0 : VERM1D(2) = VERM(I2)
1443 : 0 : CALL S1D0ONEI(SIM,SIMI,VERM1D)
1444 : 0 : VERL(N) = SIM
1445 : 0 : VERLI(N) = SIMI/3 - 0.5D0
1446 : 0 : GO TO 300
1447 : : END IF
1448 : : !C Here are 3 cases, how w() are placed on complex-plane
1449 : : !C LONE(1..3) [w1<w2<ONE];[w1<ONE<w2];[ONE<w1<w2]
1450 : : !C LEPS(1) [|w1-w2|<e
1451 : : !*-ASIS
1452 : 0 : CALL S2D0LEPS(VERM, N, U, AS, LONE, EPS, LEPS, iuerr)
1453 : : !c U(1) = W(1); U(2) = W(2)
1454 [ # # ]: 0 : IF (LONE(1)) THEN
1455 : 0 : CALL S2D1TWOI(U,AS,LEPS,SIM,SIMI)
1456 : : END IF
1457 : : !c U(1) = W(1); U(2) = CONE/W(2)
1458 [ # # ]: 0 : IF (LONE(2)) THEN
1459 : 0 : CALL S2D2TWOI(U,AS,LEPS,SIM,SIMI)
1460 : : END IF
1461 : : !c U(1) = CONE/W(1); U(2) = CONE/W(2)
1462 [ # # ]: 0 : IF (LONE(3)) THEN
1463 : 0 : CALL S2D3TWOI(U,AS,LEPS,SIM,SIMI)
1464 : : END IF
1465 : : !c
1466 : 0 : VERL(N) = SIM/VERM(N)/4
1467 : 0 : VERLI(N) = (log(VERM(N))+SIMI/4)/3 - 5.0D0/18.0D0
1468 : 0 : 300 CONTINUE ! end DO ... N = 1, 3
1469 : 0 : RETURN
1470 : : END SUBROUTINE S2D0TWOI
1471 : : !C
1472 : : !C ----------------------------------------------------------------------
1473 : : !C file: s2d1onei.f
1474 : : !C date: 2011-03-23
1475 : : !C who: S.Kaprzyk
1476 : : !C what: CASE a) : [w1<w2<ONE]
1477 : : !C U(1)=W(1); U(2)=W(2)
1478 : : !C LEPS(1) [|w1-w2|<e]
1479 : 0 : SUBROUTINE S2D1ONEI(U,AS,LEPS, SIM1, SIM1I)
1480 : :
1481 : : !IMPLICIT NONE
1482 : : DOUBLE COMPLEX U(2)
1483 : : DOUBLE PRECISION AS(2)
1484 : : LOGICAL LEPS(1)
1485 : : DOUBLE COMPLEX SIM1, SIM1I
1486 : : !C
1487 : : !DOUBLE COMPLEX SIM0UX0
1488 : : !EXTERNAL SIM0UX0
1489 : : DOUBLE COMPLEX X0(2), DX0(4), R0(2), S0(2)
1490 : : DOUBLE COMPLEX QX0(1), QR0(1), QS0(1)
1491 : : DOUBLE COMPLEX AV, BV, DU, G(5)
1492 : : INTEGER I, K
1493 : : DOUBLE COMPLEX CZERO, CONE
1494 : : DOUBLE PRECISION ZERO, ONE
1495 : : DATA ZERO/0.0D0/, ONE/1.0D0/
1496 : :
1497 : : ABI_UNUSED((/as(1)/))
1498 : : !C
1499 : 0 : CZERO = (0.0D0,0.0D0)
1500 : 0 : CONE = (1.0D0,0.0D0)
1501 : : !C ------------------------------------------------------------------
1502 : : !C R0(U) = ARTH(U)/U
1503 : : !C R0(U) = 1+U**2/3+U**6*X0(U)
1504 : : !C ------------------------------------------------------------------
1505 [ # # ]: 0 : DO 100 I = 1, 2
1506 : 0 : AV = U(I)*U(I)
1507 : 0 : X0(I) = SIM0UX0(U(I))
1508 : 0 : R0(I) = CONE + AV*(CONE/3.0D0+AV*(CONE/5.0D0+AV*X0(I)))
1509 : 0 : S0(I) = (CONE-U(I))/(CONE+U(I))*R0(I)
1510 : 0 : 100 CONTINUE
1511 : : !C -----------------------------------------------------------------
1512 : : !C - The DX0(I) contains derivatives of the function X0(u) -
1513 : : !C -----------------------------------------------------------------
1514 : : !c IF (LEPS(1).OR.LEPS(2).OR.LEPS(3)) THEN
1515 [ # # ]: 0 : IF (LEPS(1)) THEN
1516 : 0 : CALL SIM0UDX0(U(2),X0(2),DX0)
1517 : : END IF
1518 : : !c DO 200 I = 1, 3, 2
1519 : 0 : I = 1
1520 : 0 : DU = U(I) - U(2)
1521 [ # # ]: 0 : IF (.NOT.LEPS(I)) THEN
1522 : 0 : QX0(I) = (X0(I)-X0(2))/DU
1523 : : ELSE
1524 : 0 : QX0(I) = DX0(1) + (DX0(2)/2.0D0+(DX0(3)/6.0D0+ DX0(4)/24.0D0*DU)*DU)*DU
1525 : : END IF
1526 : 0 : AV = U(I)
1527 : 0 : G(1) = U(I) + U(2)
1528 [ # # ]: 0 : DO 150 K = 2, 5
1529 : 0 : AV = AV*U(I)
1530 : 0 : G(K) = G(K-1)*U(2) + AV
1531 : 0 : 150 CONTINUE
1532 : 0 : QR0(I) = G(1)/3.0D0 + G(3)/5.0D0 + G(5)*X0(2) + AV*U(I)*QX0(I)
1533 : 0 : QS0(I) = (-S0(2)-R0(2)+(CONE-U(I))*QR0(I))/(CONE+U(I))
1534 : : !C
1535 : 0 : AV = (CONE+U(1))*(CONE+U(2))
1536 : 0 : BV = (U(2)-CONE)
1537 : 0 : SIM1 = AV*(R0(1) + QR0(1)*BV)
1538 : 0 : SIM1I = AV*(S0(1) + QS0(1)*BV)
1539 : 0 : RETURN
1540 : : END SUBROUTINE S2D1ONEI
1541 : : !C
1542 : : !C ----------------------------------------------------------------------
1543 : : !C file: s2d1twoi.f
1544 : : !C date: 2010-09-22
1545 : : !C who: S.Kaprzyk
1546 : : !C what: CASE a) : [w1<w2<ONE]
1547 : : !C U(1)=W(1); U(2)=W(2)
1548 : : !C LEPS(1) [|w1-w2|<e]
1549 : 0 : SUBROUTINE S2D1TWOI(U,AS,LEPS, SIM, SIMI)
1550 : :
1551 : : !IMPLICIT NONE
1552 : : DOUBLE COMPLEX U(2)
1553 : : DOUBLE PRECISION AS(2)
1554 : : LOGICAL LEPS(2)
1555 : : DOUBLE COMPLEX SIM, SIMI
1556 : : !C
1557 : : !DOUBLE COMPLEX SIM0UX0
1558 : : !EXTERNAL SIM0UX0
1559 : : !c
1560 : : DOUBLE COMPLEX X0(2), DX0(4), R3(2), S3(2)
1561 : : DOUBLE COMPLEX QR3(1), QX0(1)
1562 : : DOUBLE COMPLEX QS3(1)
1563 : : DOUBLE COMPLEX AV, BV, UV, G(5)
1564 : : INTEGER I, K
1565 : : DOUBLE COMPLEX CZERO, CONE, CTWO
1566 : : DOUBLE PRECISION ZERO, ONE
1567 : : DATA ZERO/0.0D0/, ONE/1.0D0/
1568 : : !C
1569 : :
1570 : : ABI_UNUSED((/as/))
1571 : 0 : CZERO = (0.0D0,0.0D0)
1572 : 0 : CONE = (1.0D0,0.0D0)
1573 : 0 : CTWO = (2.0D0,0.0D0)
1574 : : !C -----------------------------------------------------------------
1575 : : !C - -
1576 : : !C - R3(W)=1+(W-1)*(ARTH(W)/W-1)/W -
1577 : : !C - R3(W)=1-W/3+W**2/3-W**3/5+W**4/5+(W-1)*W**5*X0(W) -
1578 : : !C - -
1579 : : !C -----------------------------------------------------------------
1580 [ # # ]: 0 : DO 100 I = 1, 2
1581 : 0 : AV = U(I)
1582 : 0 : X0(I) = SIM0UX0(AV)
1583 : : !*-ASIS
1584 : : R3(I) = CONE + AV*(-CONE/3.0D0 + AV*(CONE/3.0D0 + &
1585 : 0 : & AV*((-0.2D0,0.0D0)+AV*((0.2D0,0.0D0)+AV*(AV-CONE)*X0(I)))))
1586 : 0 : S3(I) = (CONE-U(I))/(CONE+U(I))*R3(I)
1587 : 0 : 100 CONTINUE
1588 : : !C -----------------------------------------------------------------
1589 : : !C - The DX0(I) Contains Derivatives Of The Function X0(W) -
1590 : : !C -----------------------------------------------------------------
1591 [ # # ]: 0 : IF (LEPS(1)) THEN
1592 : 0 : CALL SIM0UDX0(U(2),X0(2),DX0)
1593 : : END IF
1594 : : !C
1595 : : !c DO 200 I = 1, 3, 2
1596 : 0 : I = 1
1597 : 0 : UV = U(I) - U(2)
1598 [ # # ]: 0 : IF (LEPS(I)) THEN
1599 : 0 : QX0(I) = DX0(1) + (DX0(2)/2.0D0+(DX0(3)/6.0D0+ DX0(4)/24.0D0*UV)*UV)*UV
1600 : : ELSE
1601 : 0 : QX0(I) = (X0(I)-X0(2))/UV
1602 : : END IF
1603 : 0 : AV = U(I)
1604 : 0 : G(1) = U(I) + U(2)
1605 [ # # ]: 0 : DO 150 K = 2, 5
1606 : 0 : AV = AV*U(I)
1607 : 0 : G(K) = G(K-1)*U(2) + AV
1608 : 0 : 150 CONTINUE
1609 : : QR3(I) = -CONE/3.0D0 + G(1)/3.0D0 - G(2)/5.0D0 + &
1610 : 0 : & G(3)/5.0D0 + (G(5)-G(4))*X0(2) + (U(I)-CONE)*AV*QX0(I)
1611 : 0 : QS3(I) = (-S3(2)-R3(2)+(CONE-U(I))*QR3(I))/(CONE+U(I))
1612 : : !c 200 CONTINUE ! DO 750 I = 1, 3, 2
1613 : : !c
1614 : 0 : AV = (CONE+U(1))*(CONE+U(2))
1615 : 0 : BV = U(2)-CONE
1616 : 0 : SIM = AV*(R3(1) + QR3(1)*BV)
1617 : 0 : SIMI = AV*(S3(1) + QS3(1)*BV)
1618 : 0 : RETURN
1619 : : END SUBROUTINE S2D1TWOI
1620 : : !C
1621 : : !C ----------------------------------------------------------------------
1622 : : !C file: s2d2onei.f
1623 : : !C date: 2011-03-24
1624 : : !C who: S.Kaprzyk
1625 : : !C what: CASE b) : [w1<ONE<w2]
1626 : : !C U(1)=W(1); U(2)= CONE/W(2)
1627 : : !C LEPS(1) [|w1-w2|<e]
1628 : 0 : SUBROUTINE S2D2ONEI(U,AS,LEPS, SIM2, SIM2I)
1629 : :
1630 : : !IMPLICIT NONE
1631 : : DOUBLE COMPLEX U(2)
1632 : : DOUBLE PRECISION AS(2)
1633 : : LOGICAL LEPS(1)
1634 : : DOUBLE COMPLEX SIM2, SIM2I
1635 : : !c
1636 : : !DOUBLE COMPLEX SIM0UX0
1637 : : !EXTERNAL SIM0UX0
1638 : : DOUBLE COMPLEX X0(2), R0(2), S0(2) !, DX0(4)
1639 : : !c DOUBLE COMPLEX QX0(1), QR0(1), QS0(1)
1640 : : DOUBLE COMPLEX AV, BV, CV !, G(5)
1641 : : !DOUBLE COMPLEX UD
1642 : : INTEGER I !, K, N
1643 : : DOUBLE COMPLEX CZERO, CONE, CTWO
1644 : : DOUBLE PRECISION ZERO, ONE
1645 : : DATA ZERO/0.0D0/, ONE/1.0D0/
1646 : : !C
1647 : : ABI_UNUSED(leps)
1648 : 0 : CZERO = (0.0D0,0.0D0)
1649 : 0 : CONE = (1.0D0,0.0D0)
1650 : 0 : CTWO = (2.0D0,0.0D0)
1651 : : !C ------------------------------------------------------------------
1652 : : !C R0(U) = ARTH(U)/U
1653 : : !C R0(U) = 1+U**2/3+U**6*X0(U)
1654 : : !C ------------------------------------------------------------------
1655 [ # # ]: 0 : DO 100 I = 1, 2
1656 : 0 : AV = U(I)*U(I)
1657 : 0 : X0(I) = SIM0UX0(U(I))
1658 : 0 : R0(I) = CONE + AV*(CONE/3.0D0+AV*(CONE/5.0D0+AV*X0(I)))
1659 : 0 : S0(I) = (CONE-U(I))/(CONE+U(I))*R0(I)
1660 : 0 : 100 CONTINUE
1661 : : !c
1662 : 0 : AV = (CONE+U(1))*(U(2)+CONE)/(CONE-U(2)*U(1))
1663 : 0 : BV = -(U(1)-CONE)
1664 : 0 : CV = (CONE-U(2))
1665 : 0 : SIM2 = AV*(BV*R0(1)+CV*(DCMPLX(ZERO,AS(2))+U(2)*R0(2)))
1666 : 0 : CV = CV*(U(2)-CONE)/(U(2)+CONE)
1667 : 0 : SIM2I= AV*(BV*S0(1)+CV*(DCMPLX(ZERO,AS(2))+U(2)*R0(2)))
1668 : : !c
1669 : 0 : RETURN
1670 : : END SUBROUTINE S2D2ONEI
1671 : : !C
1672 : : !C ----------------------------------------------------------------------
1673 : : !C file: s2d2twoi.f
1674 : : !C date: 2011-03-26
1675 : : !C who: S.Kaprzyk
1676 : : !C what: CASE b) : [w1<ONE<w2]
1677 : : !C U(1)=W(1); U(2)= CONE/W(2)
1678 : : !C LEPS(1) [|w1-w2|<e]
1679 : 0 : SUBROUTINE S2D2TWOI(U,AS,LEPS, SIM, SIMI)
1680 : :
1681 : : !IMPLICIT NONE
1682 : : DOUBLE COMPLEX U(2)
1683 : : DOUBLE PRECISION AS(2)
1684 : : LOGICAL LEPS(1)
1685 : : DOUBLE COMPLEX SIM, SIMI
1686 : : !C
1687 : : !DOUBLE COMPLEX SIM0UX0
1688 : : !EXTERNAL SIM0UX0
1689 : : !c
1690 : : DOUBLE COMPLEX X0(2), R3(2), S3(2)
1691 : : DOUBLE COMPLEX AV, BV
1692 : : INTEGER I
1693 : : DOUBLE COMPLEX CZERO, CONE
1694 : : DOUBLE PRECISION ZERO, ONE
1695 : : DATA ZERO/0.0D0/, ONE/1.0D0/
1696 : : !C
1697 : :
1698 : : ABI_UNUSED((/leps/))
1699 : 0 : CZERO = (0.0D0,0.0D0)
1700 : 0 : CONE = (1.0D0,0.0D0)
1701 : : !C -----------------------------------------------------------------
1702 : : !C - -
1703 : : !C - R3(W)=1+(W-1)*(ARTH(W)/W-1)/W -
1704 : : !C - R3(W)=1-W/3+W**2/3-W**3/5+W**4/5+(W-1)*W**5*X0(W) -
1705 : : !C - -
1706 : : !C -----------------------------------------------------------------
1707 [ # # ]: 0 : DO 100 I = 1, 2
1708 : 0 : AV = U(I)
1709 : 0 : X0(I) = SIM0UX0(AV)
1710 : : !*-ASIS
1711 : : R3(I) = CONE + AV*(-CONE/3.0D0 + AV*(CONE/3.0D0 + &
1712 : 0 : & AV*((-0.2D0,0.0D0)+AV*((0.2D0,0.0D0)+AV*(AV-CONE)*X0(I)))))
1713 : 0 : 100 CONTINUE
1714 : 0 : S3(1) = (CONE-U(1))/(CONE+U(1))*R3(1)
1715 : 0 : AV = (CONE+U(1))*(U(2)+CONE)/(U(1)*U(2)-CONE)
1716 : 0 : BV = (CONE-U(2))*DCMPLX(ZERO,AS(2))+CONE+U(2)-U(2)*U(2)*R3(2)
1717 : 0 : SIM = AV*((U(1)-CONE)*R3(1) - (CONE-U(2))*BV)
1718 : 0 : BV = (U(2)-CONE)/(U(2)+CONE)*BV
1719 : 0 : SIMI= AV*((U(1)-CONE)*S3(1) - (CONE-U(2))*BV)
1720 : : !C
1721 : 0 : RETURN
1722 : : END SUBROUTINE S2D2TWOI
1723 : : !C
1724 : : !C ----------------------------------------------------------------------
1725 : : !C file: s2d3onei.f
1726 : : !C date: 2011-04-01
1727 : : !C who: S.Kaprzyk
1728 : : !C what: CASE c) : [ ONE<w1<w2]
1729 : : !C U(1)=CONE/W(1); U(2)=CONE/W(2)
1730 : : !C LEPS(1) [|w1-w2|<e]
1731 : 0 : SUBROUTINE S2D3ONEI(U,AS,LEPS,SIM3,SIM3I)
1732 : :
1733 : : !IMPLICIT NONE
1734 : : DOUBLE COMPLEX U(2)
1735 : : DOUBLE PRECISION AS(2)
1736 : : LOGICAL LEPS(1)
1737 : : DOUBLE COMPLEX SIM3, SIM3I
1738 : : !c
1739 : : !DOUBLE COMPLEX SIM0UX0
1740 : : !EXTERNAL SIM0UX0
1741 : : !c
1742 : : DOUBLE COMPLEX X0(2), QX0(1), DX0(4)
1743 : : DOUBLE COMPLEX R0(2), UR0(2), QR0(1), QUR0(1)
1744 : : DOUBLE COMPLEX US0(2), QUS0(1), QAS(1)
1745 : : DOUBLE COMPLEX AV, UD, G(5)
1746 : : INTEGER I, K
1747 : : DOUBLE COMPLEX CZERO, CONE, CIMAG
1748 : : DOUBLE PRECISION ZERO, ONE
1749 : : DATA ZERO/0.0D0/, ONE/1.0D0/
1750 : : !C
1751 : 0 : CZERO = (0.0D0,0.0D0)
1752 : 0 : CONE = (1.0D0,0.0D0)
1753 : 0 : CIMAG = (0.0D0,1.0D0)
1754 : : !C ------------------------------------------------------------------
1755 : : !C R0(U)=ARTH(U)/U
1756 : : !C R0(U)=1+U**2/3+U**6*X0(U)
1757 : : !C ------------------------------------------------------------------
1758 [ # # ]: 0 : DO 100 I = 1, 2
1759 : 0 : AV = U(I)*U(I)
1760 : 0 : X0(I) = SIM0UX0(U(I))
1761 : 0 : R0(I) = CONE + AV*(CONE/3.0D0+AV*(CONE/5.0D0+AV*X0(I)))
1762 : 0 : UR0(I)= DCMPLX(ZERO,AS(I)) + U(I)*R0(I)
1763 : 0 : 100 CONTINUE
1764 : : !C -----------------------------------------------------------------
1765 : : !C - THE DX0(I) CONTAINS DERIVATIVES OF THE FUNCTION X0(U)
1766 : : !C -----------------------------------------------------------------
1767 [ # # ]: 0 : IF (LEPS(1)) THEN
1768 : 0 : CALL SIM0UDX0(U(2),X0(2),DX0)
1769 : : END IF
1770 : : !c
1771 : 0 : I = 1
1772 : 0 : UD = U(I) - U(2)
1773 [ # # ]: 0 : IF (.NOT.LEPS(I)) THEN
1774 : 0 : QX0(I) = (X0(I)-X0(2))/UD
1775 : 0 : QAS(I) = (AS(I)-AS(2))/UD
1776 : : ELSE
1777 : 0 : QX0(I) = DX0(1) + (DX0(2)/2.0D0+(DX0(3)/6.0D0 + DX0(4)/24.0D0*UD)*UD)*UD
1778 : 0 : QAS(I) = CZERO
1779 : : END IF
1780 : 0 : AV = U(I)
1781 : 0 : G(1) = U(I) + U(2)
1782 [ # # ]: 0 : DO 150 K = 2, 5
1783 : 0 : AV = AV*U(I)
1784 : 0 : G(K) = G(K-1)*U(2) + AV
1785 : 0 : 150 CONTINUE
1786 : 0 : QR0(I) = G(1)/3.0D0 + G(3)/5.0D0 + G(5)*X0(2) + AV*U(I)*QX0(I)
1787 : : !c
1788 : 0 : QUR0(1) = CIMAG*QAS(1) + R0(2) +U(1)*QR0(1)
1789 : 0 : US0(2) = (U(2)-CONE)/(U(2)+CONE)*UR0(2)
1790 : 0 : QUS0(1) = ( -US0(2)+UR0(2)+(U(1)-CONE)*QUR0(1))/(U(1)+CONE)
1791 : : !c
1792 : 0 : AV = (U(1)+CONE)*(U(2)+CONE)
1793 : 0 : SIM3 = AV*(UR0(2) + (U(1)-CONE)*QUR0(1))
1794 : 0 : SIM3I = AV*(US0(2) + (U(1)-CONE)*QUS0(1))
1795 : 0 : RETURN
1796 : : END SUBROUTINE S2D3ONEI
1797 : : !C
1798 : : !C ----------------------------------------------------------------------
1799 : : !C file: s2d3twoi.f
1800 : : !C date: 2011-03-28
1801 : : !C who: S.Kaprzyk
1802 : : !C what: CASE c) : [ ONE<w1<w2]
1803 : : !C U(1)=CONE/W(1); U(2)=CONE/W(2)
1804 : : !C LEPS(1) [|w1-w2|<e]
1805 : 0 : SUBROUTINE S2D3TWOI(U,AS,LEPS,SIM,SIMI)
1806 : :
1807 : : !IMPLICIT NONE
1808 : : DOUBLE COMPLEX U(2)
1809 : : DOUBLE PRECISION AS(2)
1810 : : LOGICAL LEPS(1)
1811 : : DOUBLE COMPLEX SIM, SIMI
1812 : : !C
1813 : : !DOUBLE COMPLEX SIM0UX0
1814 : : !EXTERNAL SIM0UX0
1815 : : !c
1816 : : DOUBLE COMPLEX X0(2), QX0(1), DX0(4)
1817 : : DOUBLE COMPLEX R3(2), UR3(2), QR3(1), QUR3(1)
1818 : : DOUBLE COMPLEX US3(2), QUS3(1)
1819 : : DOUBLE COMPLEX QAS(1), AV, UV, G(5)
1820 : : INTEGER I, K
1821 : : DOUBLE COMPLEX CZERO, CONE, CTWO, CIMAG
1822 : : DOUBLE PRECISION ZERO, ONE
1823 : : DATA ZERO/0.0D0/, ONE/1.0D0/
1824 : : !C
1825 : 0 : CZERO = (0.0D0,0.0D0)
1826 : 0 : CONE = (1.0D0,0.0D0)
1827 : 0 : CTWO = (2.0D0,0.0D0)
1828 : 0 : CIMAG = (0.0D0,1.0D0)
1829 : : !C -----------------------------------------------------------------
1830 : : !C - -
1831 : : !C - R3(U)=1+(U-1)*(ARTH(U)/U-1)/U -
1832 : : !C - R3(U)=1-U/3+U**2/3-U**3/5+U**4/5+(U-1)*U**5*X0(U) -
1833 : : !C - -
1834 : : !C -----------------------------------------------------------------
1835 [ # # ]: 0 : DO 100 I = 1, 2
1836 : 0 : AV = U(I)
1837 : 0 : X0(I) = SIM0UX0(AV)
1838 : : !*-ASIS
1839 : : R3(I) = CONE + AV*(-CONE/3.0D0 + AV*(CONE/3.0D0 + &
1840 : 0 : & AV*((-0.2D0,0.0D0)+AV*((0.2D0,0.0D0)+AV*(AV-CONE)*X0(I)))))
1841 : 0 : UR3(I)=(CONE-U(I))*DCMPLX(ZERO,AS(I))+CONE+U(I)-U(I)*U(I)*R3(I)
1842 : 0 : 100 CONTINUE ! end DO ... I = 1, 2
1843 : : !C -----------------------------------------------------------------
1844 : : !C - The DX0(I) Contains Derivatives Of The Function X0(U) -
1845 : : !C -----------------------------------------------------------------
1846 [ # # ]: 0 : IF (LEPS(1)) THEN
1847 : 0 : CALL SIM0UDX0(U(2),X0(2),DX0)
1848 : : END IF
1849 : : !C
1850 : 0 : I = 1
1851 : 0 : UV = U(I) - U(2)
1852 [ # # ]: 0 : IF (LEPS(I)) THEN
1853 : 0 : QX0(I) = DX0(1) + (DX0(2)/2.0D0+(DX0(3)/6.0D0+DX0(4)/24.0D0*UV)*UV)*UV
1854 : 0 : QAS(I) = CZERO
1855 : : ELSE
1856 : 0 : QX0(I) = (X0(I)-X0(2))/UV
1857 : 0 : QAS(I) = (AS(I)-AS(2))/UV
1858 : : END IF
1859 : 0 : AV = U(I)
1860 : 0 : G(1) = U(I) + U(2)
1861 [ # # ]: 0 : DO 150 K = 2, 5
1862 : 0 : AV = AV*U(I)
1863 : 0 : G(K) = G(K-1)*U(2) + AV
1864 : 0 : 150 CONTINUE
1865 : : QR3(I) = -CONE/3.0D0 + G(1)/3.0D0 - G(2)/5.0D0 + &
1866 : 0 : & G(3)/5.0D0 + (G(5)-G(4))*X0(2) + (U(I)-CONE)*AV*QX0(I)
1867 : : QUR3(1) = -DCMPLX(ZERO,AS(2))+CIMAG*(CONE-U(1))*QAS(1)+ CONE - &
1868 : 0 : & (U(1)+U(2))*R3(2)-U(1)*U(1)*QR3(1)
1869 : : !C
1870 : 0 : US3(2) = (U(2)-CONE)/(U(2)+CONE)*UR3(2)
1871 : 0 : QUS3(1)=( -US3(2)+UR3(2)+(U(1)-CONE)*QUR3(1))/(U(1)+CONE)
1872 : : !c
1873 : 0 : AV = (U(1)+CONE)*(U(2)+CONE)
1874 : 0 : SIM = AV*(UR3(2) + (U(1)-CONE)*QUR3(1))
1875 : 0 : SIMI= AV*(US3(2) + (U(1)-CONE)*QUS3(1))
1876 : : !c
1877 : 0 : RETURN
1878 : : END SUBROUTINE S2D3TWOI
1879 : : !C
1880 : : !C ----------------------------------------------------------------------
1881 : : !C file: s2d0leps.f
1882 : : !C date: 2011-04-02
1883 : : !C who: S.Kaprzyk
1884 : : !C what: Return a proper CASE:
1885 : : !C a) [w1<w2<ONE]; b) [w1<ONE<w2]; c) [ONE<w1<w2]
1886 : : !C INPUT:
1887 : : !C VERM(3) - three complex numbers
1888 : : !C N - integer number for W(I)=(VERM(I)-VERM(N))/(VERM(I)+VERM(N))
1889 : : !C EPS - critical distance |w1-w2| etc.
1890 : : !C OUTPUT:
1891 : : !C LONE(3) [w1<w2<ONE]; [w1<ONE<w2]; [ONE<w1<w2]
1892 : : !C LEPS(1) [|w1-w2|<eps]
1893 : 0 : SUBROUTINE S2D0LEPS(VERM, N, W, AS, LONE, EPS, LEPS, iuerr)
1894 : :
1895 : : !IMPLICIT NONE
1896 : : DOUBLE COMPLEX VERM(3)
1897 : : INTEGER N
1898 : : DOUBLE COMPLEX W(2)
1899 : : DOUBLE PRECISION AS(2)
1900 : : DOUBLE PRECISION EPS
1901 : : LOGICAL LONE(3), LEPS(1)
1902 : : INTEGER iuerr
1903 : : !C
1904 : : !c DOUBLE PRECISION DLAMCH
1905 : : !c EXTERNAl DLAMCH
1906 : : !C
1907 : : LOGICAL LWONE(2), LV
1908 : : DOUBLE COMPLEX AV
1909 : : DOUBLE PRECISION AL(2), AA, AIW, OZERO
1910 : : INTEGER I, II, K
1911 : : DOUBLE COMPLEX CZERO, CONE
1912 : : DOUBLE PRECISION ZERO, ONE, PI
1913 : : DATA ZERO/0.0D0/,ONE/1.0D0/,PI/3.141592653589793D0/
1914 : : DATA OZERO/ 1.0D-13/
1915 : : !c OZERO = DLAMCH('e')*10
1916 : 0 : CZERO = (0.0D0,0.0D0)
1917 : 0 : CONE = (1.0D0,0.0D0)
1918 [ # # ]: 0 : IF ((N.GT.3).OR.(N.LT.1)) THEN
1919 : 0 : WRITE (iuerr,9010) N
1920 : : 9010 FORMAT (' ***s2d0leps: N =',I6,' must be 1, or 2')
1921 : 0 : STOP ' ***s2d0leps: '
1922 : : END IF
1923 : : II = 0
1924 [ # # ]: 0 : DO 200 I = 1, 3
1925 [ # # ]: 0 : IF (I.EQ.N) GO TO 200
1926 : 0 : II = II + 1
1927 [ # # ]: 0 : IF (CDABS(VERM(N)-VERM(I)).LT.CDABS(VERM(N)+VERM(I))) THEN
1928 : 0 : LWONE(II) = .TRUE.
1929 : 0 : W(II) = (VERM(N)-VERM(I))/(VERM(N)+VERM(I))
1930 : 0 : AIW = DIMAG(W(II))
1931 [ # # ]: 0 : IF (DABS(AIW).GE.OZERO) THEN
1932 : 0 : AS(II) = 0.5D0*PI*DSIGN(ONE,AIW)
1933 : : ELSE
1934 : 0 : AS(II) = 0.5D0*PI*DSIGN(ONE,DREAL(VERM(I)-VERM(N)))
1935 : : END IF
1936 : : ELSE
1937 : 0 : LWONE(II) = .FALSE.
1938 : 0 : W(II) = (VERM(N)+VERM(I))/(VERM(N)-VERM(I))
1939 : 0 : AIW = DIMAG(W(II))
1940 [ # # ]: 0 : IF (DABS(AIW).GE.OZERO) THEN
1941 : 0 : AS(II) = -0.5D0*PI*DSIGN(ONE,AIW)
1942 : : ELSE
1943 : 0 : AS(II) = 0.5D0*PI*DSIGN(ONE,DREAL(VERM(I)-VERM(N)))
1944 : : END IF
1945 : : END IF
1946 : 0 : AL(II) = CDABS(W(II))
1947 : 0 : 200 CONTINUE
1948 [ # # ]: 0 : DO 300 I = 1, 2
1949 [ # # ]: 0 : DO 250 K = I, 2
1950 [ # # # # : 0 : IF (LWONE(I).AND.LWONE(K).AND.(AL(I).LT.AL(K))) GO TO 250
# # ]
1951 [ # # # # : 0 : IF (.NOT.LWONE(I).AND..NOT.LWONE(K).AND.(AL(K).LT.AL(I))) GO TO 250
# # ]
1952 [ # # # # : 0 : IF (.NOT.LWONE(I).AND.LWONE(K).AND.(ONE.LT.AL(I)*AL(K))) GO TO 250
# # ]
1953 [ # # # # : 0 : IF (LWONE(I).AND..NOT.LWONE(K).AND.(AL(I)*AL(K).LT.ONE)) GO TO 250
# # ]
1954 : 0 : AA = AL(K)
1955 : 0 : AL(K) = AL(I)
1956 : 0 : AL(I) = AA
1957 : 0 : AA = AS(K)
1958 : 0 : AS(K) = AS(I)
1959 : 0 : AS(I) = AA
1960 : 0 : AV = W(K)
1961 : 0 : W(K) = W(I)
1962 : 0 : W(I) = AV
1963 : 0 : LV = LWONE(K)
1964 : 0 : LWONE(K) = LWONE(I)
1965 : 0 : LWONE(I) = LV
1966 : 0 : 250 CONTINUE
1967 : 0 : 300 CONTINUE
1968 : : !C
1969 : 0 : LONE(1) = .FALSE.
1970 : 0 : LONE(2) = .FALSE.
1971 : 0 : LONE(3) = .FALSE.
1972 [ # # # # ]: 0 : IF (LWONE(1).AND.LWONE(2)) THEN
1973 : 0 : LONE(1) = .TRUE.
1974 : : END IF
1975 [ # # # # ]: 0 : IF (LWONE(1).AND..NOT.LWONE(2)) THEN
1976 : 0 : LONE(2) = .TRUE.
1977 : : END IF
1978 [ # # # # ]: 0 : IF (.NOT.LWONE(1).AND..NOT.LWONE(2)) THEN
1979 : 0 : LONE(3) = .TRUE.
1980 : : END IF
1981 : : !C here are cases, how close are w()
1982 : 0 : LEPS(1) = .FALSE.
1983 [ # # # # : 0 : IF ( (LONE(1).OR.LONE(3)).AND. (CDABS(W(1)-W(2)).LT.EPS) ) LEPS(1)=.TRUE.
# # ]
1984 : : !c
1985 [ # # # # : 0 : IF ((.NOT.LONE(1)).AND.(.NOT.LONE(2)).AND.(.NOT.LONE(3))) THEN
# # ]
1986 : 0 : WRITE (iuerr,9030) (LONE(I),I=1,3)
1987 : : 9030 FORMAT ('***s2d0leps: LONE()=',3(L1,1X),'not in order')
1988 : 0 : STOP '***s2d0leps: '
1989 : : END IF
1990 : : !C
1991 : 0 : RETURN
1992 : : END SUBROUTINE S2D0LEPS
1993 : : !C s2d0leps()
1994 : : !C file: s1d0onei.f
1995 : : !C date: 2011-03-20
1996 : : !C who: S.Kaprzyk
1997 : : !C what: Complex linear form integral over standard 1-d simplex
1998 : : !C ----------------------------------------------------
1999 : : !C - * -
2000 : : !C - * 1 -
2001 : : !C - SIM0= * dt1 ---------------------------------- -
2002 : : !C - * verm(2)+(verm(1)-verm(2))*t1 -
2003 : : !C - * -
2004 : : !C - 0<t1<1 -
2005 : : !C ----------------------------------------------------
2006 : 0 : SUBROUTINE S1D0ONEI(SIM0, SIM0I, VERM)
2007 : :
2008 : : !IMPLICIT NONE
2009 : : INTEGER iuerr
2010 : : PARAMETER (iuerr=6)
2011 : : DOUBLE COMPLEX SIM0, SIM0I
2012 : : DOUBLE COMPLEX VERM(*)
2013 : : !C
2014 : : DOUBLE COMPLEX SIM, SIMI
2015 : : LOGICAL LONE(2),LEPS(1)
2016 : : !c
2017 : : DOUBLE COMPLEX U(1)
2018 : : DOUBLE PRECISION AS(1)
2019 : : INTEGER I, K, N
2020 : : DOUBLE COMPLEX CZERO, CONE
2021 : : DOUBLE PRECISION ZERO, EPS
2022 : : DATA EPS/1.0D-6/
2023 : : DATA ZERO/0.0D0/
2024 : : !C
2025 : 0 : CZERO = (0.0D0,0.0D0)
2026 : 0 : CONE = (1.0D0,0.0D0)
2027 : 0 : N = 2
2028 [ # # ]: 0 : DO 100 I = 1, 1
2029 [ # # ]: 0 : IF (DIMAG(VERM(N))*DIMAG(VERM(I)).LT.ZERO) THEN
2030 [ # # # # ]: 0 : WRITE (iuerr,9010) (DIMAG(VERM(K)),K=1,2)
2031 : : 9010 FORMAT (' ***s1d0onei: not signed ImgVERM()=',4(D13.6,1X))
2032 : : !c STOP ' ***s1d0onei: '
2033 : : END IF
2034 : 0 : 100 CONTINUE
2035 [ # # ]: 0 : DO 101 I = 1, 1
2036 [ # # ]: 0 : IF (CDABS(VERM(I)).GT.CDABS(VERM(N))) N = I
2037 : 0 : 101 CONTINUE
2038 : : !C Here are 2 cases, how w() are placed on complex-plane
2039 : : !C LONE(1..2) [w1<ONE]; [ONE<w1]
2040 : : !C LEPS(1) null
2041 : : !*-ASIS
2042 : 0 : CALL S1D0LEPS(VERM, N, U, AS, LONE, EPS, LEPS, iuerr)
2043 : : !c U(1) = W(1);
2044 [ # # ]: 0 : IF (LONE(1)) THEN
2045 : 0 : CALL S1D1ONEI(U,AS,LEPS,SIM,SIMI)
2046 : : END IF
2047 : : !c U(1) = CONE/W(1)
2048 [ # # ]: 0 : IF (LONE(2)) THEN
2049 : 0 : CALL S1D2ONEI(U,AS,LEPS,SIM,SIMI)
2050 : : END IF
2051 : : !C
2052 : 0 : SIM0 = SIM/VERM(N)
2053 : 0 : SIM0I= log(VERM(N)) + SIMI - CONE
2054 : 0 : RETURN
2055 : : END SUBROUTINE S1D0ONEI
2056 : : !C
2057 : : !C ----------------------------------------------------------------------
2058 : : !C file: s1d0twoi.f
2059 : : !C date: 2011-03-24
2060 : : !C who: S.Kaprzyk
2061 : : !C what: Complex linear form integrals over standard 1d-simplex
2062 : : !C ------------------(i=1)--------------------------------------
2063 : : !C * -
2064 : : !C * t_i -
2065 : : !C VERL(i)= *dt1 ------------------------------------------ -
2066 : : !C * VERM(2)+[VERM(1)-VERM(2)]*t1 -
2067 : : !C * -
2068 : : !C 0<t1<1 -
2069 : : !C -------------------(i=2)-------------------------------------
2070 : : !C * -
2071 : : !C * 1- t1 -
2072 : : !C VERL(2)= *dt1 ------------------------------------------ -
2073 : : !C * VERM(2)+[VERM(1)-VERM(2)]*t1 -
2074 : : !C * -
2075 : : !C 0<t1<1 -
2076 : : !C -------------------------------------------------------------
2077 : 0 : SUBROUTINE S1D0TWOI(VERL, VERLI, VERM)
2078 : :
2079 : : !IMPLICIT NONE
2080 : : INTEGER iuerr
2081 : : PARAMETER (iuerr=6)
2082 : : DOUBLE COMPLEX VERL(*), VERLI(*), VERM(*)
2083 : : !C
2084 : : DOUBLE COMPLEX SIM, SIMI
2085 : : LOGICAL LONE(2),LEPS(1)
2086 : : !c
2087 : : INTEGER I, K, N
2088 : : DOUBLE PRECISION AS(1)
2089 : : DOUBLE COMPLEX U(1)
2090 : : DOUBLE COMPLEX CZERO
2091 : : DOUBLE PRECISION EPS
2092 : : DOUBLE PRECISION ZERO
2093 : : DATA EPS/1.0D-6/
2094 : : DATA ZERO/0.0D0/
2095 : : !C
2096 : 0 : CZERO = (0.0D0,0.0D0)
2097 [ # # ]: 0 : DO 50 I = 1, 1
2098 [ # # ]: 0 : IF (DIMAG(VERM(2))*DIMAG(VERM(I)).LT.ZERO) THEN
2099 [ # # # # ]: 0 : WRITE (iuerr,9010) (DIMAG(VERM(K)),K=1,2)
2100 : : 9010 FORMAT (' ***s1d0twoi: not signed ImgVERM()=',2(D13.6,1X))
2101 : 0 : STOP ' ***s1dt0woi: '
2102 : : END IF
2103 : 0 : 50 CONTINUE
2104 [ # # ]: 0 : DO 100 N = 1, 2
2105 : 0 : VERL(N) = CZERO
2106 : : !C Here are 2 cases, how w() are placed on complex-plane
2107 : : !C LONE(1..2) [w1<ONE]; [ONE<w1]
2108 : : !C LEPS(1) null
2109 : : !*-ASIS
2110 : 0 : CALL S1D0LEPS(VERM, N, U, AS, LONE, EPS, LEPS, iuerr)
2111 : : !c U(1) = W(1);
2112 [ # # ]: 0 : IF (LONE(1)) THEN
2113 : 0 : CALL S1D1TWOI(U,AS,LEPS,SIM,SIMI)
2114 : : END IF
2115 : : !c U(1) = CONE/W(1);
2116 [ # # ]: 0 : IF (LONE(2)) THEN
2117 : 0 : CALL S1D2TWOI(U,AS,LEPS,SIM,SIMI)
2118 : : END IF
2119 : : !c
2120 : 0 : VERL(N) = 0.50D0*SIM/VERM(N)
2121 : 0 : VERLI(N) = 0.50D0*log(VERM(N)) + 0.25D0*SIMI -0.25D0
2122 : 0 : 100 CONTINUE ! end DO .. N = 1, 2
2123 : 0 : RETURN
2124 : : END SUBROUTINE S1D0TWOI
2125 : : !C
2126 : : !C ----------------------------------------------------------------------
2127 : : !C file: s1d1onei.f
2128 : : !C date: 2010-09-22
2129 : : !C who: S.Kaprzyk
2130 : : !C what: CASE a) : [w1<ONE]
2131 : : !C U(1)=W(1)
2132 : : !C LEPS() null
2133 : 0 : SUBROUTINE S1D1ONEI(U,AS,LEPS, SIM1, SIM1I)
2134 : :
2135 : : !IMPLICIT NONE
2136 : : DOUBLE COMPLEX U(1)
2137 : : DOUBLE PRECISION AS(1)
2138 : : LOGICAL LEPS(1)
2139 : : DOUBLE COMPLEX SIM1, SIM1I
2140 : : !C
2141 : : !DOUBLE COMPLEX SIM0UR0
2142 : : !EXTERNAL SIM0UR0
2143 : : DOUBLE COMPLEX R0(1)
2144 : : DOUBLE COMPLEX CZERO, CONE
2145 : : DOUBLE PRECISION ZERO, ONE
2146 : : DATA ZERO/0.0D0/, ONE/1.0D0/
2147 : : !C
2148 : : ABI_UNUSED(as)
2149 : : ABI_UNUSED(leps)
2150 : 0 : CZERO = (0.0D0,0.0D0)
2151 : 0 : CONE = (1.0D0,0.0D0)
2152 : : !C ------------------------------------------------------------------
2153 : : !C R0(U) = ARTH(U)/U
2154 : : !C R0(U) = 1+U**2/3+U**6*X0(U)
2155 : : !C ------------------------------------------------------------------
2156 : 0 : R0(1) = SIM0UR0(U(1))
2157 : 0 : SIM1 = (CONE+U(1))*R0(1)
2158 : 0 : SIM1I = (CONE-U(1))*R0(1)
2159 : 0 : RETURN
2160 : : END SUBROUTINE S1D1ONEI
2161 : : !C
2162 : : !C ----------------------------------------------------------------------
2163 : : !C file: s1d1twoi.f
2164 : : !C date: 2011-03-25
2165 : : !C who: S.Kaprzyk
2166 : : !C what: CASE a) : [w1<ONE]
2167 : : !C U(1)=W(1)
2168 : : !C LEPS() null
2169 : 0 : SUBROUTINE S1D1TWOI(U,AS,LEPS, SIM, SIMI)
2170 : :
2171 : : !IMPLICIT NONE
2172 : : DOUBLE COMPLEX U(1)
2173 : : DOUBLE PRECISION AS(1)
2174 : : LOGICAL LEPS(1)
2175 : : DOUBLE COMPLEX SIM, SIMI
2176 : : !C
2177 : : !DOUBLE COMPLEX SIM0UX0
2178 : : !EXTERNAL SIM0UX0
2179 : : !c
2180 : : DOUBLE COMPLEX X0(1), R2(1), S2(1)
2181 : : DOUBLE COMPLEX AV
2182 : : DOUBLE COMPLEX CZERO, CONE
2183 : : DOUBLE PRECISION ZERO, ONE
2184 : : DATA ZERO/0.0D0/, ONE/1.0D0/
2185 : : !C
2186 : : ABI_UNUSED(as)
2187 : : ABI_UNUSED(leps)
2188 : 0 : CZERO = (0.0D0,0.0D0)
2189 : 0 : CONE = (1.0D0,0.0D0)
2190 : : !C -----------------------------------------------------------------
2191 : : !C - -
2192 : : !C - R2(W)=1+(W-1)*(ARTH(W)/W-1)/W=[1+(W-1)*R0(W)]/W -
2193 : : !C - R2(W)=1-W/3+W**2/3-W**3/5+W**4/5+(W-1)*W**5*X0(W) -
2194 : : !C - -
2195 : : !C -----------------------------------------------------------------
2196 : 0 : AV = U(1)
2197 : 0 : X0(1) = SIM0UX0(AV)
2198 : : !*-ASIS
2199 : : R2(1) = CONE + AV*(-CONE/3.0D0 + AV*(CONE/3.0D0 + &
2200 : 0 : & AV*((-0.2D0,0.0D0)+AV*((0.2D0,0.0D0)+AV*(AV-CONE)*X0(1)))))
2201 : 0 : S2(1) = (CONE-U(1))/(CONE+U(1))*R2(1)
2202 : 0 : SIM = (CONE+U(1))*R2(1)
2203 : 0 : SIMI = (CONE+U(1))*S2(1)
2204 : 0 : RETURN
2205 : : END SUBROUTINE S1D1TWOI
2206 : : !C
2207 : : !C ----------------------------------------------------------------------
2208 : : !C file: s1d2onei.f
2209 : : !C date: 2011-03-20
2210 : : !C who: S.Kaprzyk
2211 : : !C what: CASE b) : [ ONE<w1]
2212 : : !C U(1) = CONE/W(1)
2213 : : !C LEPS(1) null
2214 : 0 : SUBROUTINE S1D2ONEI(U,AS,LEPS,SIM2,SIM2I)
2215 : :
2216 : : !IMPLICIT NONE
2217 : : DOUBLE COMPLEX U(1)
2218 : : DOUBLE PRECISION AS(1)
2219 : : LOGICAL LEPS(1)
2220 : : DOUBLE COMPLEX SIM2, SIM2I
2221 : : !c
2222 : : !DOUBLE COMPLEX SIM0UR0
2223 : : !EXTERNAL SIM0UR0
2224 : : DOUBLE COMPLEX R0(1)
2225 : : DOUBLE COMPLEX CZERO, CONE, CIMAG
2226 : : DOUBLE PRECISION ZERO, ONE
2227 : : DATA ZERO/0.0D0/, ONE/1.0D0/
2228 : : !C
2229 : : ABI_UNUSED(leps)
2230 : 0 : CZERO = (0.0D0,0.0D0)
2231 : 0 : CONE = (1.0D0,0.0D0)
2232 : 0 : CIMAG = (0.0D0,1.0D0)
2233 : : !C ------------------------------------------------------------------
2234 : : !C R0(U)=ARTH(U)/U
2235 : : !C R0(U)=1+U**2/3+U**6*X0(U)
2236 : : !C ------------------------------------------------------------------
2237 : 0 : R0(1) = SIM0UR0(U(1))
2238 : 0 : SIM2 = (U(1)+CONE)*(CIMAG*AS(1)+U(1)*R0(1))
2239 : 0 : SIM2I= (U(1)-CONE)*(CIMAG*AS(1)+U(1)*R0(1))
2240 : : !c
2241 : 0 : RETURN
2242 : : END SUBROUTINE S1D2ONEI
2243 : : !C
2244 : : !C ----------------------------------------------------------------------
2245 : : !C s1d2onei()
2246 : : !C file: s1d2twoi.f
2247 : : !C date: 2011-03-25
2248 : : !C who: S.Kaprzyk
2249 : : !C what: CASE b) : [ ONE<w1]
2250 : : !C U(1) = CONE/W(1)
2251 : : !C LEPS(1) null
2252 : 0 : SUBROUTINE S1D2TWOI(U,AS,LEPS, SIM, SIMI)
2253 : :
2254 : : !IMPLICIT NONE
2255 : : DOUBLE COMPLEX U(1)
2256 : : DOUBLE PRECISION AS(1)
2257 : : LOGICAL LEPS(1)
2258 : : DOUBLE COMPLEX SIM, SIMI
2259 : : !C
2260 : : !DOUBLE COMPLEX SIM0UX0
2261 : : !EXTERNAL SIM0UX0
2262 : : !c
2263 : : DOUBLE COMPLEX X0(2), R2(2), S2(2)
2264 : : DOUBLE COMPLEX AV
2265 : : !INTEGER I
2266 : : DOUBLE COMPLEX CZERO, CONE
2267 : : DOUBLE PRECISION ZERO, ONE
2268 : : DATA ZERO/0.0D0/, ONE/1.0D0/
2269 : : !C
2270 : : ABI_UNUSED(leps)
2271 : 0 : CZERO = (0.0D0,0.0D0)
2272 : 0 : CONE = (1.0D0,0.0D0)
2273 : : !C -----------------------------------------------------------------
2274 : : !C - -
2275 : : !C - R2(W)=1+(W-1)*(ARTH(W)/W-1)/W -
2276 : : !C - R2(W)=1-W/3+W**2/3-W**3/5+W**4/5+(W-1)*W**5*X0(W) -
2277 : : !C - -
2278 : : !C -----------------------------------------------------------------
2279 : 0 : AV = U(1)
2280 : 0 : X0(1) = SIM0UX0(AV)
2281 : : !*-ASIS
2282 : : R2(1) = CONE + AV*(-CONE/3.0D0 + AV*(CONE/3.0D0 + &
2283 : 0 : & AV*((-0.2D0,0.0D0)+AV*((0.2D0,0.0D0)+AV*(AV-CONE)*X0(1)))))
2284 : : !C
2285 : 0 : R2(1)= (CONE-U(1))*DCMPLX(ZERO,AS(1))+CONE+U(1)-U(1)*U(1)*R2(1)
2286 : 0 : S2(1) = (U(1)-CONE)/(U(1)+CONE)*R2(1)
2287 : 0 : SIM = (U(1)+CONE)*R2(1)
2288 : 0 : SIMI = (U(1)+CONE)*S2(1)
2289 : : !C
2290 : 0 : RETURN
2291 : : END SUBROUTINE S1D2TWOI
2292 : : !C
2293 : : !C ----------------------------------------------------------------------
2294 : : !C file: s1d0leps.f
2295 : : !C date: 2011-04-02
2296 : : !C who: S.Kaprzyk
2297 : : !C what: Return a proper case:
2298 : : !C a) [w1<ONE]; b) [ONE<w1]
2299 : : !C INPUT:
2300 : : !C VERM(1,2) - two complex numbers
2301 : : !C N - integer number for W(I)=(VERM(I)-VERM(N))/(VERM(I)+VERM(N))
2302 : : !C EPS - critical distance
2303 : : !C OUTPUT:
2304 : : !C LONE(2) [w1<ONE]; [ONE<w1]
2305 : : !C LEPS(1) null
2306 : 0 : SUBROUTINE S1D0LEPS(VERM, N, W, AS, LONE, EPS, LEPS, iuerr)
2307 : :
2308 : : !IMPLICIT NONE
2309 : : DOUBLE COMPLEX VERM(2)
2310 : : INTEGER N
2311 : : DOUBLE COMPLEX W(1)
2312 : : DOUBLE PRECISION AS(1)
2313 : : DOUBLE PRECISION EPS
2314 : : LOGICAL LONE(2), LEPS(1)
2315 : : INTEGER iuerr
2316 : : !C
2317 : : !c DOUBLE PRECISION DLAMCH
2318 : : !c EXTERNAl DLAMCH
2319 : : !C
2320 : : LOGICAL LWONE(1)
2321 : : DOUBLE PRECISION AL(1), AIW, OZERO
2322 : : INTEGER I, II !, K
2323 : : DOUBLE COMPLEX CZERO, CONE
2324 : : DOUBLE PRECISION ZERO, ONE, PI
2325 : : DATA ZERO/0.0D0/,ONE/1.0D0/,PI/3.141592653589793D0/
2326 : : DATA OZERO/ 1.0D-13/
2327 : : ABI_UNUSED(eps)
2328 : : !c OZERO = DLAMCH('e')*10
2329 : 0 : CZERO = (0.0D0,0.0D0)
2330 : 0 : CONE = (1.0D0,0.0D0)
2331 [ # # ]: 0 : IF ((N.GT.2).OR.(N.LT.1)) THEN
2332 : 0 : WRITE (iuerr,9010) N
2333 : : 9010 FORMAT (' ***s1d0leps: N =',I6,' must be 1, or 2')
2334 : 0 : STOP ' ***s1d0leps: '
2335 : : END IF
2336 : : !C
2337 : : II = 0
2338 [ # # ]: 0 : DO 200 I = 1, 2
2339 [ # # ]: 0 : IF (I.EQ.N) GO TO 200
2340 : 0 : II = II + 1
2341 [ # # ]: 0 : IF (CDABS(VERM(N)-VERM(I)).LT.CDABS(VERM(N)+VERM(I))) THEN
2342 : 0 : LWONE(II) = .TRUE.
2343 : 0 : W(II) = (VERM(N)-VERM(I))/(VERM(N)+VERM(I))
2344 : 0 : AIW = DIMAG(W(II))
2345 [ # # ]: 0 : IF (DABS(AIW).GE.OZERO) THEN
2346 : 0 : AS(II) = 0.5D0*PI*DSIGN(ONE,AIW)
2347 : : ELSE
2348 : 0 : AS(II) = 0.5D0*PI*DSIGN(ONE,DREAL(VERM(I)-VERM(N)))
2349 : : END IF
2350 : : ELSE
2351 : 0 : LWONE(II) = .FALSE.
2352 : 0 : W(II) = (VERM(N)+VERM(I))/(VERM(N)-VERM(I))
2353 : 0 : AIW = DIMAG(W(II))
2354 [ # # ]: 0 : IF (DABS(AIW).GE.OZERO) THEN
2355 : 0 : AS(II) = -0.5D0*PI*DSIGN(ONE,AIW)
2356 : : ELSE
2357 : 0 : AS(II) = 0.5D0*PI*DSIGN(ONE,DREAL(VERM(I)-VERM(N)))
2358 : : END IF
2359 : : END IF
2360 : 0 : AL(II) = CDABS(W(II))
2361 : 0 : 200 CONTINUE
2362 : : !C
2363 : 0 : LONE(1) = .FALSE.
2364 : 0 : LONE(2) = .FALSE.
2365 [ # # ]: 0 : IF (LWONE(1)) THEN
2366 : 0 : LONE(1) = .TRUE.
2367 : : END IF
2368 [ # # ]: 0 : IF (.NOT.LWONE(1)) THEN
2369 : 0 : LONE(2) = .TRUE.
2370 : : END IF
2371 : : !C here are cases, how close are w()
2372 : 0 : LEPS(1) = .FALSE.
2373 : : !c
2374 [ # # # # ]: 0 : IF ((.NOT.LONE(1)).AND.(.NOT.LONE(2))) THEN
2375 : 0 : WRITE (iuerr,9030) (LONE(I),I=1,2)
2376 : : 9030 FORMAT ('***s1d0leps: LONE()=',2(L1,1X),'not in order')
2377 : 0 : STOP '***s1d0leps: '
2378 : : END IF
2379 : : !c
2380 : 0 : RETURN
2381 : : END SUBROUTINE S1D0LEPS
2382 : : !C
2383 : :
2384 : : end module m_simtet
2385 : : !!***
|
