1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
|
*DECK SYB31C
SUBROUTINE SYB31C (PPLUS,TAUP,XOPJ,XOPI)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Evaluation of the $C_{ij}$ function in 1D cylindrical and 2D geometry.
*
*Copyright:
* Copyright (C) 2002 Ecole Polytechnique de Montreal
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version
*
*Author(s): A. Hebert
*
*Parameters: input
* TAUP side to side optical path.
* XOPJ optical path in volume j (or volume i).
* XOPI optical path in volume i (or volume j).
*
*Parameters: output
* PPLUS value of the probability.
*
*-----------------------------------------------------------------------
*
*----
* SUBROUTINE ARGUMENTS
*----
REAL PPLUS,TAUP,XOPJ,XOPI
*----
* LOCAL VARIABLES
*----
PARAMETER (MKI3=600)
REAL T(2,2),B(2,2)
COMMON /BICKL3/BI3(0:MKI3),BI31(0:MKI3),BI32(0:MKI3),PAS3,XLIM3,L3
*----
* ASYMPTOTIC VALUE
*----
IF(TAUP.GE.XLIM3) THEN
PPLUS = 0.
RETURN
ENDIF
*----
* SYMMETRIC FORMULA IN (I,J). WE SET POPI <= POPJ
*----
IF(XOPJ.LT.XOPI) THEN
POPI = XOPJ
POPJ = XOPI
ELSE
POPI = XOPI
POPJ = XOPJ
ENDIF
*----
* GENERAL CASE
*
* POPI < POPJ < 1.
* VOID (I,J) UP TO THE END OF PROGRAM
*----
T(2,2) = TAUP + POPJ + POPI
T(1,2) = TAUP + POPJ
T(2,1) = TAUP + POPI
T(1,1) = TAUP
*
DO 15 I=1,2
DO 10 J=1,2
B(I,J) = TABKI(3,T(I,J))
10 CONTINUE
15 CONTINUE
*----
* GENERAL DIFFERENCE FORMULA. THIS FORMULAS SHOULD NOT BE APPLIED TO
* VOIDED (I,J) VOLUMES
*----
IF(POPI.GE.0.004) THEN
* LARGE POPI AND POPJ (MOST GENERAL CASE)
PNLUS = B(2,2) + B(1,1) - (B(1,2) + B(2,1))
PPLUS = PNLUS
ELSE IF(POPJ.GT.0.002) THEN
* SMALL POPI, LARGE POPJ. USE DERIVATIVE DIFFERENCES.
IF(TAUP.GE.XLIM3) THEN
APLUS=0.
ELSE IF(TAUP+POPI.GE.XLIM3) THEN
APLUS=TABKI(3,TAUP)
ELSE IF(POPI.LE.0.002) THEN
APLUS=(TABKI(2,TAUP)+TABKI(2,TAUP+POPI))*POPI*0.5
ELSE IF(POPI.LT.0.004) THEN
PQLUS=(TABKI(2,TAUP)+TABKI(2,TAUP+POPI))*POPI*0.5
PRLUS=TABKI(3,TAUP)-TABKI(3,TAUP+POPI)
FACT=500.0*(POPI-0.002)
APLUS=PRLUS*FACT+PQLUS*(1.0-FACT)
ELSE
APLUS=TABKI(3,TAUP)-TABKI(3,TAUP+POPI)
ENDIF
IF(TAUP+POPJ.GE.XLIM3) THEN
BPLUS=0.
ELSE IF(TAUP+POPI+POPJ.GE.XLIM3) THEN
BPLUS=TABKI(3,TAUP+POPJ)
ELSE IF(POPI.LE.0.002) THEN
BPLUS=(TABKI(2,TAUP+POPJ)+TABKI(2,TAUP+POPI+POPJ))*POPI*0.5
ELSE IF(POPI.LT.0.004) THEN
PQLUS=(TABKI(2,TAUP+POPJ)+TABKI(2,TAUP+POPI+POPJ))*POPI*0.5
PRLUS=TABKI(3,TAUP+POPJ)-TABKI(3,TAUP+POPI+POPJ)
FACT=500.0*(POPI-0.002)
BPLUS=PRLUS*FACT+PQLUS*(1.0-FACT)
ELSE
BPLUS=TABKI(3,TAUP+POPJ)-TABKI(3,TAUP+POPI+POPJ)
ENDIF
PPLUS = APLUS - BPLUS
ELSE IF(T(2,2).GE.XLIM3) THEN
* SIMILAR TO A SECOND DERIVATIVE. ASYMPTOTIC TAUP.
PNLUS = B(2,2) + B(1,1) - (B(1,2) + B(2,1))
PPLUS = PNLUS
ELSE
* SIMILAR TO A SECOND DERIVATIVE.
IF(T(2,2).EQ.0.) THEN
PPLUS = 0.
ELSE IF(TAUP.EQ.0.) THEN
PPLUS = 1.57079632679489 + TABKI(1,TAUP+POPI+POPJ)
PPLUS = PPLUS * 0.5 * POPI*POPJ
ELSE IF(TAUP.LE.0.002) THEN
PPLUS = TABKI(1,TAUP) + TABKI(1,TAUP+POPI+POPJ)
PPLUS = PPLUS * 0.5 * POPI*POPJ
ELSE IF(TAUP.LT.0.004) THEN
PQLUS = TABKI(1,TAUP) + TABKI(1,TAUP+POPI+POPJ)
PRLUS = TABKI(1,TAUP)
TAUP1 = TAUP + POPI + POPJ
PRLUS = PRLUS + TABKI(1,TAUP1)
FACT=500.0*(TAUP-0.002)
PPLUS = PRLUS * FACT + PQLUS * (1.0-FACT)
PPLUS = PPLUS * 0.5 * POPI*POPJ
ELSE
* VOIDED VOLUMES (I,J) SEPARATED WITH NON-VOIDED VOLUMES.
PRLUS = TABKI(1,TAUP)
TAUP1 = TAUP + POPI + POPJ
PRLUS = PRLUS + TABKI(1,TAUP1)
PPLUS = PRLUS * 0.5 * POPI*POPJ
ENDIF
ENDIF
RETURN
END
|
