blob: d87177a66c27db06100a0acb1f78bf2d2ac07031 (
plain)
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
|
*DECK SYB41C
SUBROUTINE SYB41C(PP,TAU0,TAUI,TAUJ)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Evaluation of the $C_{ij}$ function in 1D spherical geometry.
*
*Copyright:
* Copyright (C) 2005 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
* TAU0 side to side optical path.
* TAUI optical path in volume i (or volume j).
* TAUJ optical path in volume j (or volume i).
*
*Parameters: output
* PP value of the expression.
*
*-----------------------------------------------------------------------
*
*----
* SUBROUTINE ARGUMENTS
*----
REAL PP,TAU0,TAUI,TAUJ
*----
* SYMMETRIC FORMULA IN (I,J). WE SET POPI <= POPJ
*----
IF(TAUJ.LT.TAUI) THEN
POPI = TAUJ
POPJ = TAUI
ELSE
POPI = TAUI
POPJ = TAUJ
ENDIF
IF(POPI.GT.0.004) THEN
PP=EXP(-TAU0)-EXP(-(TAU0+POPI))-EXP(-(TAU0+POPJ))+
1 EXP(-(TAU0+POPI+POPJ))
ELSE IF((POPI.GT.0.002).AND.(POPJ.GT.0.004)) THEN
PP=EXP(-TAU0)-EXP(-(TAU0+POPI))-EXP(-(TAU0+POPJ))+
1 EXP(-(TAU0+POPI+POPJ))
PQ=POPI*(EXP(-TAU0)-EXP(-(TAU0+POPJ)))
FACTI=500.0*(POPI-0.002)
PP=PP*FACTI+PQ*(1.0-FACTI)
ELSE IF((POPI.GT.0.002).AND.(POPJ.GT.0.002)) THEN
PP=EXP(-TAU0)-EXP(-(TAU0+POPI))-EXP(-(TAU0+POPJ))+
1 EXP(-(TAU0+POPI+POPJ))
PQ=POPI*(EXP(-TAU0)-EXP(-(TAU0+POPJ)))
PR=POPJ*(EXP(-TAU0)-EXP(-(TAU0+POPI)))
PS=POPI*POPJ*EXP(-TAU0)
FACTI=500.0*(POPI-0.002)
FACTJ=500.0*(POPJ-0.002)
PP=PP*FACTI*FACTJ+PQ*(1.0-FACTI)*FACTJ+PR*FACTI*(1.0-FACTJ)
1 +PS*(1.0-FACTI)*(1.0-FACTJ)
ELSE IF(POPJ.GT.0.004) THEN
PP=POPI*(EXP(-TAU0)-EXP(-(TAU0+POPJ)))
ELSE IF(POPJ.GT.0.002) THEN
PP=POPI*(EXP(-TAU0)-EXP(-(TAU0+POPJ)))
PS=POPI*POPJ*EXP(-TAU0)
FACTJ=500.0*(POPJ-0.002)
PP=PP*FACTJ+PS*(1.0-FACTJ)
ELSE
PP=POPI*POPJ*EXP(-TAU0)
ENDIF
RETURN
END
|
