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
|
*DECK MOCSCAL
SUBROUTINE MOCSCAL(N,NREG,NSOUT,M,NOM,NZON,H,SIGANG,DSIG,EXPT,
1 EXP2,NMU,ZMU)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Calculate coefficients of a track for the cyclic characteristics
* integration: Linear-Discontinuous-Characteristics scheme with
* tabulated exponential and
* 'source term isolation' option turned off.
*
*Copyright:
* Copyright (C) 2015 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
* N number of elements in the current track.
* NREG number of volumes.
* NSOUT number of surfaces.
* M number of material mixtures.
* NOM vector containing the region number of the different segments
* of this track.
* NZON index-number of the mixture type assigned to each volume.
* H vector containing the lenght of the different segments of this
* track.
* SIGANG macroscopic total cross sections and albedos.
* NMU order of the polar quadrature set.
* ZMU inverse of polar quadrature cosines.
*
*Parameters: output
* DSIG macroscopic total cross sections.
* EXPT track coefficient.
* EXP2 quadratic expansion of (1-exp(-a*L))/L with small argument.
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER N,NREG,NSOUT,M,NOM(N),NZON(-NSOUT:NREG),NMU
REAL SIGANG(-6:M),ZMU(NMU)
DOUBLE PRECISION H(N),DSIG(N),EXPT(NMU,N),EXP2(5,NMU,N)
*----
* LOCAL VARIABLES
*----
INTEGER I,NOMI,NUMOLD,NZI,IMU
DOUBLE PRECISION TAUDMIN,SQ3,HID,TAUD,TEMP,C1,C2,H2,H3
* tabulated exponential common block
REAL E0, E1, PAS1, DX1, XLIM1
INTEGER MEX1, LAU
PARAMETER ( MEX1=7936, TAUDMIN=2.0D-2 )
COMMON /EXP1/ E0(0:MEX1),E1(0:MEX1),PAS1,DX1,XLIM1
*
SQ3=SQRT(3.0D0)
NUMOLD=NOM(1)
DO I=1,N
NOMI=NOM(I)
NZI=NZON(NOMI)
DSIG(I)=SIGANG(NZI)
IF(NZI.LT.0) THEN
DO IMU=1,NMU
EXP2(2,IMU,I)=0.D0
EXP2(3,IMU,I)=0.D0
EXP2(5,IMU,I)=0.D0
ENDDO
IF(NUMOLD.NE.NOMI) THEN
DO IMU=1,NMU
EXP2(1,IMU,I)=SIGANG(NZI)
EXP2(4,IMU,I)=EXP2(1,IMU,I)
EXPT(IMU,I)=EXP2(1,IMU,I)
ENDDO
ELSE
DO IMU=1,NMU
EXP2(1,IMU,I)=1.D0
EXP2(4,IMU,I)=1.D0
EXPT(IMU,I)=EXP2(1,IMU,I)
ENDDO
ENDIF
ELSE
DO IMU=1,NMU
HID=DBLE(H(I)*ZMU(IMU))
TAUD=SIGANG(NZI)*HID
IF(TAUD.LE.TAUDMIN) THEN
* Use Taylor series expansions
H2=HID*HID
H3=H2*HID
EXPT(IMU,I)=TAUD*(0.5D0*TAUD-1.0D0)+1.0D0
EXP2(1,IMU,I)=HID*(TAUD*(TAUD/6.0D0-0.5D0)+1.0D0)
EXP2(2,IMU,I)=H2*(TAUD*(TAUD-4.0D0)+12.0D0)/24.0D0
EXP2(3,IMU,I)=-SQ3*H3*(TAUD*(TAUD-2.0D0)+4.0D0)/24.0D0
EXP2(4,IMU,I)=-SQ3*TAUD*(TAUD*(TAUD-2.0D0)+4.0D0)
1 /24.0D0
EXP2(5,IMU,I)=H3*(TAUD*TAUD-TAUD+4.0D0)/40.0D0
ELSE IF(TAUD.GE.XLIM1) THEN
* Out of the table range
EXPT(IMU,I)=0.D0
TEMP=1.D0/TAUD
EXP2(1,IMU,I)=TEMP*HID
EXP2(2,IMU,I)=HID*(1.D0-TEMP)/DSIG(I)
EXP2(3,IMU,I)=-SQ3*HID*(2.0D0-(TAUD+2.0D0)*TEMP)
1 /DSIG(I)**2
EXP2(4,IMU,I)=-SQ3*(2.0D0-(TAUD+2.0D0)*TEMP)/TAUD
C1=TAUD*(TAUD-6.0D0)-12.0D0
C2=TAUD*(3.0D0*TAUD+12.0D0)+12.0D0
EXP2(5,IMU,I)=(C1+C2*TEMP)/(TAUD*DSIG(I)**3)
ELSE
* Use tabulated exponential
LAU=INT(REAL(TAUD)*PAS1)
TEMP=DBLE(E0(LAU)+E1(LAU)*TAUD)
EXPT(IMU,I)=1.0D0-TEMP*TAUD
EXP2(1,IMU,I)=TEMP*HID
EXP2(2,IMU,I)=HID*(1.D0-TEMP)/DSIG(I)
EXP2(3,IMU,I)=-SQ3*HID*(2.0D0-(TAUD+2.0D0)*TEMP)
1 /DSIG(I)**2
EXP2(4,IMU,I)=-SQ3*(2.0D0-(TAUD+2.0D0)*TEMP)/TAUD
C1=TAUD*(TAUD-6.0D0)-12.0D0
C2=TAUD*(3.0D0*TAUD+12.0D0)+12.0D0
EXP2(5,IMU,I)=(C1+C2*TEMP)/(TAUD*DSIG(I)**3)
ENDIF
ENDDO
ENDIF
NUMOLD=NOMI
ENDDO
*
RETURN
END
|
