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 : !!***
|
