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
178
179
180
181
182
183
184
185
186
187
188
|
*DECK NXTIAA
FUNCTION NXTIAA(POSANN ,POSPIN,VOLINT)
*
*----------
*
*Purpose:
* Compute the volume of intersection between
* a 2--D annular region and an annular pin.
*
*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): G. Marleau.
*
*Parameters: input
* POSANN spatial description of the annular region with
* POSANN(0) the radius, POSANN(1) the $X$ position
* of center and POSANN(2) the $Y$ position
* of center.
* POSPIN spatial description of the annular pin region with
* POSPIN(0) the radius, POSPIN(1) the $X$ position
* of center and POSPIN(2) the $Y$ position
* of center.
*
*Parameters: output
* NXTIAA type of intersection between annular region and
* annular pin, where:
* = 0 means that there is no intersection
* between the two regions;
* = 1 means that the annular region
* is all located inside the annular pin;
* = 2 means that the annular pin
* is all located inside the annular region;
* =-1 means that the intersection between
* the annular region and the annular pin is partial.
* VOLINT 2-D volume of intersection (area) between annular region and
* annular pin.
*
*Reference:
* G. Marleau,
* New Geometries Processing in DRAGON: The NXT: Module,
* Report IGE-260, Polytechnique Montreal,
* Montreal, 2005.
*
*----
*
IMPLICIT NONE
*----
* Subroutine arguments
*----
INTEGER NXTIAA
DOUBLE PRECISION POSANN(0:2),POSPIN(0:2)
DOUBLE PRECISION VOLINT
*----
* Local parameters
*----
INTEGER IOUT
CHARACTER NAMSBR*6
PARAMETER (IOUT=6,NAMSBR='NXTIAA')
INTEGER IPRINT
PARAMETER (IPRINT=100)
DOUBLE PRECISION DCUTOF
PARAMETER (DCUTOF=1.0D-8)
DOUBLE PRECISION DZERO
PARAMETER (DZERO=0.0D0)
*----
* Functions
*----
DOUBLE PRECISION XDRCST,PI
*----
* Local variables
*----
INTEGER IFACE
DOUBLE PRECISION PX,PY,PX2,PY2,RA2,RP2,DAP2,DAP,XINT,YINT,
> ANGLEP,ANGLEA,SP,SA,VOLANN,VOLPIN,ACARG
DOUBLE PRECISION DT1,DT2,DT3
*----
* Initialize NXTIAA and VOLINT
*----
IF(IPRINT .GE. 200) THEN
WRITE(IOUT,6000) NAMSBR
WRITE(IOUT,6010) (POSANN(IFACE),IFACE=0,2)
WRITE(IOUT,6011) (POSPIN(IFACE),IFACE=0,2)
ENDIF
PI=XDRCST('Pi',' ')
NXTIAA=0
VOLINT=DZERO
*----
* Find distance from center of annulus to pin center
* $d=\sqrt(y^{2}+y^{2})$
* annular region/pin intersection if $d<r_{p}+r_{a}$
* where $r_{a}$ is annular region radius
* and $r_{p}$ is pin radius
* Annular region in pin if $r_{a}<r_{p}$ $d<r_{p}-r_{a}$
* Pin in annular region if $r_{p}<r_{a}$ $d<r_{a}-r_{p}$
*----
PX=POSANN(1)-POSPIN(1)
PY=POSANN(2)-POSPIN(2)
PX2=PX*PX
PY2=PY*PY
RA2=POSANN(0)*POSANN(0)
RP2=POSPIN(0)*POSPIN(0)
VOLANN=PI*RA2
VOLPIN=PI*RP2
DAP2=PX2+PY2
DAP=SQRT(DAP2)
IF(DAP .LE. ABS(POSANN(0)-POSPIN(0)) ) THEN
IF(POSANN(0) .LE. POSPIN(0)) THEN
VOLINT=VOLANN
NXTIAA=1
ELSE
VOLINT=VOLPIN
NXTIAA=2
ENDIF
ELSE IF(DAP .LT. POSANN(0)+POSPIN(0)) THEN
*----
* Find intersection points $(x_{i},y_{i})$ assuming
* the annular region is rotated in such a way that its center
* is located at $(d,0)$.
* Equation for annular region is :
* $(x-d)^{2}+y^{2}=r^{2}_{a}$
* Equation for pin region is :
* $x^{2}+y^{2}=r^{2}_{p}$
* Intersection point is :
* $x_{i}=(d^{2}+r^{2}_{p}-r^{2}_{s})/(2*d)$
* $y_{i}=\sqrt{r^{2}_{p}-x^{2}_{i}}$
*----
XINT=(DAP2+RP2-RA2)/(2*DAP)
YINT=SQRT(RP2-XINT*XINT)
*----
* Find angular opening for pin and annulus
* Pin surface to the right of $x=x_{i}$ plane
* $S_{p}=r_{p}^{2}\arccos(x_{i}/r_{p})-x_{i}y_{i}$
* Annular surface to the left of $x=x_{i}$ plane
* $S_{a}=r_{a}^{2}\arccos((d-x_{i})/r_{a})-(d-x_{i})y_{i}$
*----
ACARG=XINT/POSPIN(0)
IF(ACARG .GE. 1.0D0) THEN
ANGLEP=ACOS(1.0D0)
ELSE IF(ACARG .LE. -1.0D0) THEN
ANGLEP=ACOS(-1.0D0)
ELSE
ANGLEP=ACOS(ACARG)
ENDIF
ACARG=(DAP-XINT)/POSANN(0)
IF(ACARG .GE. 1.0D0) THEN
ANGLEA=ACOS(1.0D0)
ELSE IF(ACARG .LE. -1.0D0) THEN
ANGLEA=ACOS(-1.0D0)
ELSE
ANGLEA=ACOS(ACARG)
ENDIF
SP=RP2*ANGLEP-XINT*YINT
SA=RA2*ANGLEA-(DAP-XINT)*YINT
VOLINT=SP+SA
DT1=ABS(VOLINT-VOLPIN)
DT2=ABS(VOLINT-VOLANN)
DT3=ABS(VOLINT)
IF(DT1 .LT. DCUTOF) THEN
VOLINT=VOLPIN
NXTIAA=2
ELSE IF(DT2 .LT. DCUTOF) THEN
VOLINT=VOLANN
NXTIAA=1
ELSE IF(DT3 .GT. DCUTOF) THEN
NXTIAA=-1
ELSE
VOLINT=DZERO
ENDIF
ENDIF
IF(IPRINT .GE. 200) THEN
WRITE(IOUT,6012) NAMSBR,NXTIAA,VOLINT
WRITE(IOUT,6001) NAMSBR
ENDIF
RETURN
*----
* Output formats
*----
6000 FORMAT('(* Output from --',A6,'-- follows ')
6001 FORMAT(' Output from --',A6,'-- completed *)')
6010 FORMAT('POSANN={',2(F20.10,','),F20.10,'};')
6011 FORMAT('POSPIN={',2(F20.10,','),F20.10,'};')
6012 FORMAT(A6,'={',I5,',',F20.10,'};')
END
|
