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
|
*DECK SYB4VO
SUBROUTINE SYB4VO(NSECT,NR,DX,DY,RAD,VOLINT)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute the volumes of a Cartesian sectorized cell.
*
*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
* NSECT number of sectors (multiple of 4).
* NR one plus the number of tubes.
* DX X-oriented side.
* DY Y-oriented side.
* RAD radius of the tubes.
*
*Parameters: output
* VOLINT volumes.
*
*-----------------------------------------------------------------------
*
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NSECT,NR
REAL DX,DY,RAD(NR-1),VOLINT(NSECT,NR)
*
IF(MOD(4+NSECT,4).NE.0) CALL XABORT('SYB4VO: INVALID NSECT.')
DDX = DX / 2.
DDY = DY / 2.
RADDD = SQRT(DDX*DDX + DDY*DDY)
TETA0 = 3.14159265358979 * 2. / REAL(NSECT)
*----
* ITER=1 IS FOR (0, PI/4) AND ITER=2 IS FOR (PI/2, PI/4)
*----
DO 15 IR=1,NR
DO 10 IS=1,NSECT
VOLINT(IS,IR)=0.0
10 CONTINUE
15 CONTINUE
DO 60 ITER=1,2
IF(ITER.EQ.1) THEN
ISA = 1
ISB = 0
DDXY=DDX
ELSE
ISA = - 1
ISB = NSECT/4 + 1
DDXY=DDY
ENDIF
*
TETAC = ACOS(DDXY / RADDD)
TETA2 = 0.
DO 20 IR = 1, NR-1
IF (RAD(IR) .GT. DDXY) THEN
IRNEXT = IR
GOTO 30
ENDIF
20 CONTINUE
IRNEXT = NR
*----
* IRNEXT : NEXT RADIUS INTERCEPTING A SIDE
* DDY2 : LAST PROCESSED COORDINATE
* ISNEXT : NEXT SECTOR
*----
30 DDY2 = 0.
ISNEXT = 1
*----
* NEXT SECTOR
*----
40 IS = ISNEXT
IF (IS .NE. 0) THEN
ISV = ISA * IS + ISB
IR0 = IRNEXT
DDY1 = DDY2
TETA1 = TETA2
TETA2 = IS * TETA0
*
* THE ANGLE IS LIMITED BY THE DIAGONAL.
IF (TETA2 .GE. (TETAC - 1.E-6)) THEN
ISNEXT = 0
TETA2 = TETAC
RAD2 = RADDD
ELSE
ISNEXT = IS + 1
RAD2 = DDXY / COS(TETA2)
ENDIF
*
* THE NEXT RADIUS IS INTERCEPTING THE SECTOR.
IF (IR0 .LT. NR) THEN
RADIR = RAD(IR0)
IF (RADIR .LT. (RAD2 * (1. - 1.E-6))) THEN
RAD2 = RADIR
TETA2 = ACOS(DDXY / RAD2)
IRNEXT = IR0 + 1
ISNEXT = IS
ELSE IF (RADIR .LE. (RAD2 * (1. + 1.E-6))) THEN
* THE NEXT RADIUS IS EQUAL TO THE SECTOR.
IRNEXT = IR0 + 1
ENDIF
ENDIF
*
* DDY2 IS THE NEXT COORDINATE AND DT IS HALF THE ANGLE
* INCREMENT FOR THE SECTOR.
DDY2 = RAD2 * SIN(TETA2)
DT = (TETA2 - TETA1) * 0.5
*
* COMPLETE TUBES.
R1 = 0.
DO 50 IR = 1, IR0 - 1
R0 = R1
R1 = RAD(IR)
DR = (R1 - R0) * (R1 + R0)
VOLINT(ISV, IR) = VOLINT(ISV, IR) + DT * DR
50 CONTINUE
*
* LAST SIDE-INTERCEPTED TUBE.
VOLMAX = DDXY * 0.5 * (DDY2 - DDY1)
VOLMAX = VOLMAX - DT * R1 * R1
VOLINT(ISV, IR0) = VOLINT(ISV, IR0) + VOLMAX
GOTO 40
ENDIF
60 CONTINUE
*
DO 90 IR=1,NR
DO 70 IS=NSECT/4+1,NSECT/2
VOLINT(IS,IR)=VOLINT(NSECT/2-IS+1,IR)
70 CONTINUE
DO 80 IS=NSECT/2+1,NSECT
VOLINT(IS,IR)=VOLINT(NSECT-IS+1,IR)
80 CONTINUE
90 CONTINUE
RETURN
END
|
