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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
|
*DECK MCGPTN
SUBROUTINE MCGPTN(IMPX,NREG,NSOU,NANGL,NMU,VOLSUR,VOLNUM,SURNUM,
1 DENSTY,WZMU)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute errors and normalization factors for 3D prismatic extended
* tracking (adapted from NXTTLS.f).
*
*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): R. Le Tellier
*
*Parameters: input
* IMPX print flag.
* NREG number of regions.
* NSOU number of external surfaces.
* NANGL number of plan tracking angles.
* NMU number of polar angles.
* SURNUM numerical surfaces.
* DENSTY plan tracking track density per angle.
* WZMU polar quadrature weights.
*
*Parameters: input/output
* VOLSUR analytical volume and surfaces
* / analytical volumes and numerical surfaces.
* VOLNUM numerical volumes per angle / normalization factors per angle.
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER IMPX,NREG,NSOU,NANGL,NMU
REAL VOLSUR(NREG+NSOU),WZMU(NMU)
DOUBLE PRECISION DENSTY(NANGL),VOLNUM(NREG,NANGL,NMU,2),
> SURNUM(NSOU)
*----
* LOCAL VARIABLES
*----
INTEGER IOUT
CHARACTER NAMSBR*6,CSGMU(2)*1
PARAMETER (IOUT=6,NAMSBR='MCGPTN')
INTEGER IR,IS,IANGL,IMU,ISGMU,NBVERG,NBVERA,NBSERR,NBV0,ITDIR
REAL RTEMP
DOUBLE PRECISION DCERR,DSVERG,DMVERG,DAVERG,DSVERA,DMVERA,DAVERA,
1 DSSERR,DMSERR,DASERR
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: VGLNUM
DATA CSGMU / '+','-' /
*----
* SCRATCH STORAGE ALLOCATION
*----
ALLOCATE(VGLNUM(NREG))
VGLNUM(:NREG)=0.0D0
*
ITDIR=0
DO 10 IANGL=1,NANGL
IF(DENSTY(IANGL).EQ.0.0D0) GOTO 10
ITDIR=ITDIR+1
DO 20 IMU=1,NMU
DO 30 ISGMU=1,2
NBVERA=0
DSVERA=0.D0
DMVERA=0.D0
DAVERA=0.D0
NBV0=0
DO IR=1,NREG
IF(VOLNUM(IR,ITDIR,IMU,ISGMU).EQ.0.D0) THEN
IF(IMPX.GE.10) WRITE(IOUT,300) NAMSBR,IR,ITDIR,IMU,ISGMU
VOLNUM(IR,ITDIR,IMU,ISGMU)=1.D0
NBV0=NBV0+1
ELSE
VGLNUM(IR)=VGLNUM(IR)
1 +2.0*WZMU(IMU)*VOLNUM(IR,ITDIR,IMU,ISGMU)
VOLNUM(IR,ITDIR,IMU,ISGMU)=VOLNUM(IR,ITDIR,IMU,ISGMU)
1 *DENSTY(IANGL)
VOLNUM(IR,ITDIR,IMU,ISGMU)=DBLE(VOLSUR(IR))
1 /VOLNUM(IR,ITDIR,IMU,ISGMU)
NBVERA=NBVERA+1
ENDIF
DCERR=100.0D0*(1.D0-VOLNUM(IR,ITDIR,IMU,ISGMU))
DMVERA=MAX(DMVERA,ABS(DCERR))
DSVERA=DSVERA+DCERR*DCERR
DAVERA=DAVERA+DCERR
ENDDO
DSVERA=SQRT(DSVERA/DBLE(NBVERA))
DAVERA=DAVERA/DBLE(NBVERA)
IF(NBV0.GT.0) THEN
WRITE(IOUT,500) NAMSBR,NBV0,ITDIR,IMU,CSGMU(ISGMU)
ELSE
IF(IMPX.GE.2) WRITE(IOUT,100) DSVERA,DMVERA,DAVERA,ITDIR,
1 IMU,CSGMU(ISGMU)
ENDIF
30 CONTINUE
20 CONTINUE
10 CONTINUE
*
NBVERG=0
DSVERG=0.D0
DMVERG=0.D0
DAVERG=0.D0
NBV0=0
DO IR=1,NREG
IF(VGLNUM(IR).EQ.0.D0) THEN
WRITE(IOUT,300) NAMSBR,IR
VGLNUM(IR)=1.D0
NBV0=NBV0+1
ELSE
VGLNUM(IR)=DBLE(VOLSUR(IR))/VGLNUM(IR)
NBVERG=NBVERG+1
ENDIF
DCERR=100.0D0*(1.D0-VGLNUM(IR))
DMVERG=MAX(DMVERG,ABS(DCERR))
DSVERG=DSVERG+DCERR*DCERR
DAVERG=DAVERG+DCERR
ENDDO
IF(NBV0.GT.0) THEN
WRITE(IOUT,500) NAMSBR,NBV0
ENDIF
DSVERG=SQRT(DSVERG/DBLE(NBVERG))
DAVERG=DAVERG/DBLE(NBVERG)
IF(IMPX.GE.1) WRITE(IOUT,150) DSVERG,DMVERG,DAVERG
*
NBSERR=0
DSSERR=0.D0
DMSERR=0.D0
DASERR=0.D0
DO IR=1,NSOU
IS=NREG+IR
IF(SURNUM(IR).EQ.0.D0) THEN
WRITE(IOUT,400) NAMSBR,-IR
SURNUM(IR)=1.D0
ELSE
RTEMP=VOLSUR(IS)
VOLSUR(IS)=REAL(SURNUM(IR))
SURNUM(IR)=RTEMP/VOLSUR(IS)
NBSERR=NBSERR+1
ENDIF
DCERR=100.0D0*(1.D0-SURNUM(IR))
DMSERR=MAX(DMSERR,ABS(DCERR))
DSSERR=DSSERR+DCERR*DCERR
DASERR=DASERR+DCERR
ENDDO
DSSERR=SQRT(DSSERR/DBLE(NBSERR))
DASERR=DASERR/DBLE(NBSERR)
IF(IMPX.GE.1) WRITE(IOUT,200) DSSERR,DMSERR,DASERR
*----
* SCRATCH STORAGE DEALLOCATION
*----
DEALLOCATE(VGLNUM)
RETURN
*
100 FORMAT(' Angular RMS, maximum and average errors (%) ',
1 'on region volumes :',3(2X,F10.5),
2 ' for plan tracking angle',I3,' for polar angle',I3,
3 '(',A1,')')
150 FORMAT(' Global RMS, maximum and average errors (%) ',
1 'on region volumes :',3(2X,F10.5))
200 FORMAT(' Global RMS, maximum and average errors (%) ',
1 'on surface areas :',3(2X,F10.5))
300 FORMAT(1X,'***** Warning in ',A6,'*****'/
1 7X,'For region ',I8,
2 1X,'no crossing by angle ',I8,I8,'(',I1,')')
400 FORMAT(1X,'***** Warning in ',A6,'*****'/
1 7X,'For surface ',I8,
2 1X,'no crossing by any angle ')
500 FORMAT(1X,'***** Warning in ',A6,'*****'/
1 7X,I8,' regions not tracked for direction ',I8,I8,
2 '(',A1,')')
END
|
