]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1995/10/24 10:20:50 cernlib | |
6 | * Geant | |
7 | * | |
8 | * | |
9 | #include "geant321/pilot.h" | |
10 | *FCA : 05/01/99 09:58:02 by Federico Carminati | |
11 | * Effectively print the message when a shape is | |
12 | * not implemented | |
13 | *CMZ : 3.21/02 29/03/94 15.41.28 by S.Giani | |
14 | *-- Author : | |
15 | SUBROUTINE GGPERP (X,U,IERR) | |
16 | C. | |
17 | C. **************************************************************** | |
18 | C. * * | |
19 | C. * This routine solves the general problem of calculating the * | |
20 | C. * unit vector normal to the surface of the current volume at * | |
21 | C. * the point X. The result is returned in the array U. X is * | |
22 | C. * assumed to be on or near a boundary of the current volume. * | |
23 | C. * The current volume is indicated by the common /GCVOLU/. * | |
24 | C. * U points from inside to outside in that neighbourhood. * | |
25 | C. * If X is equidistant to more than one boundary (in a corner) * | |
26 | C. * an arbitrary choice is made based upon the order of * | |
27 | C. * precedence implied by the IF statements below. If the * | |
28 | C. * routine fails to find the unit normal, it returns with * | |
29 | C. * IERR=1, otherwise IERR=0. * | |
30 | C. * * | |
31 | C. * Called by : GSURFP, GDSTEP * | |
32 | C. * Authors : F.Carminati, R.Jones, F.Ohlsson-Malek * | |
33 | C. * * | |
34 | C. **************************************************************** | |
35 | #include "geant321/gcvolu.inc" | |
36 | #include "geant321/gconsp.inc" | |
37 | #include "geant321/gcbank.inc" | |
38 | #include "geant321/gcshno.inc" | |
39 | #include "geant321/gctmed.inc" | |
40 | #include "geant321/gcunit.inc" | |
41 | DIMENSION X(3),U(3),XL(3),UL(3),DXL(3),PAR(50),SPAR(50),ATT(20) | |
42 | DIMENSION PERP(10) | |
43 | #if !defined(CERNLIB_SINGLE) | |
44 | DOUBLE PRECISION PERP,PMIN0 | |
45 | DOUBLE PRECISION PAR,DXL,RHO,R,RINV,PHI,THE | |
46 | DOUBLE PRECISION PHI1,PHI2,THE1,THE2,XWID | |
47 | DOUBLE PRECISION GUARD,DPHI,PHI0,SPHI0,CPHI0 | |
48 | DOUBLE PRECISION FACT,CALPH,SALPH,TALPH | |
49 | DOUBLE PRECISION RAT,RATL,RATH,H,BL,TL,DX,DY,DZ,DU | |
50 | DOUBLE PRECISION UU0,VV0,UU,W1,W2,W3,W4 | |
51 | DOUBLE PRECISION SEW1,SEW2,SEW3,SEW4 | |
52 | DOUBLE PRECISION TAN1,TAN2,TAN3,TAN4 | |
53 | DOUBLE PRECISION SEC1,SEC2,SEC3,SEC4 | |
54 | DOUBLE PRECISION U0,V0,U1,U1L,U2,U2L | |
55 | DOUBLE PRECISION ONE,TWO | |
56 | DOUBLE PRECISION DSECT,ZERO,FULL,FULL10,DBY2 | |
57 | #endif | |
58 | LOGICAL LNTFOU | |
59 | PARAMETER (ONE=1,TWO=2) | |
60 | PARAMETER (ZERO=0.,DBY2=0.5,FULL=360.,FULL10=3600.) | |
61 | C. | |
62 | C. ------------------------------------------------------------------ | |
63 | C. | |
64 | LNTFOU = .FALSE. | |
65 | * | |
66 | * *** Transform current point into local reference system | |
67 | CALL GMTOD(X,XL,1) | |
68 | DXL(1) = XL(1) | |
69 | DXL(2) = XL(2) | |
70 | DXL(3) = XL(3) | |
71 | * | |
72 | * *** Fetch the parameters of the current volume | |
73 | JVO = LQ(JVOLUM-LVOLUM(NLEVEL)) | |
74 | IN = LINDEX(NLEVEL) | |
75 | IF (NLEVEL.GT.1) THEN | |
76 | JVOM = LQ(JVOLUM-LVOLUM(NLEVEL-1)) | |
77 | JIN = LQ(JVOM-IN) | |
78 | ENDIF | |
79 | ISH = Q(JVO+2) | |
80 | NIN = Q(JVO+3) | |
81 | IF (NLEVEL.LT.NLDEV(NLEVEL)) THEN | |
82 | JPAR = 0 | |
83 | ELSE | |
84 | * (case with structure JVOLUM locally developed) | |
85 | JPAR = LQ(LQ(JVOLUM-LVOLUM(NLDEV(NLEVEL)))) | |
86 | IF (NLEVEL.EQ.NLDEV(NLEVEL)) GO TO 20 | |
87 | DO 10 ILEV = NLDEV(NLEVEL), NLEVEL-1 | |
88 | IF (IQ(JPAR+1).EQ.0) THEN | |
89 | JPAR = LQ(JPAR-LINDEX(ILEV+1)) | |
90 | IF (JPAR.EQ.0) GO TO 20 | |
91 | ELSE IF (IQ(JPAR-3).GT.1) THEN | |
92 | JPAR = LQ(JPAR-LINDEX(ILEV+1)) | |
93 | ELSE | |
94 | JPAR = LQ(JPAR-1) | |
95 | ENDIF | |
96 | IF (ILEV.EQ.NLEVEL-1) THEN | |
97 | JPAR = JPAR + 5 | |
98 | NPAR = IQ(JPAR) | |
99 | CALL UCOPY (Q(JPAR+1), SPAR, NPAR) | |
100 | DO 100 I=1,NPAR | |
101 | PAR(I)=SPAR(I) | |
102 | 100 CONTINUE | |
103 | ENDIF | |
104 | 10 CONTINUE | |
105 | GO TO 30 | |
106 | ENDIF | |
107 | * (normal case) | |
108 | 20 CONTINUE | |
109 | CALL GFIPAR(JVO,JIN,IN,NPAR,NATT,SPAR,ATT) | |
110 | DO 101 I=1,NPAR | |
111 | PAR(I)=SPAR(I) | |
112 | 101 CONTINUE | |
113 | 30 CONTINUE | |
114 | * | |
115 | * *** Case of the BOX: | |
116 | IF (ISH.EQ.NSBOX) THEN | |
117 | PERP(1) = ABS(ABS(DXL(1))-PAR(1)) | |
118 | PERP(2) = ABS(ABS(DXL(2))-PAR(2)) | |
119 | PERP(3) = ABS(ABS(DXL(3))-PAR(3)) | |
120 | PMIN0 = MIN(PERP(1),PERP(2),PERP(3)) | |
121 | IF (PERP(1).EQ.PMIN0) THEN | |
122 | UL(1) = SIGN(ONE,DXL(1)) | |
123 | UL(2) = 0. | |
124 | UL(3) = 0. | |
125 | ELSE IF (PERP(2).EQ.PMIN0) THEN | |
126 | UL(1) = 0. | |
127 | UL(2) = SIGN(ONE,DXL(2)) | |
128 | UL(3) = 0. | |
129 | ELSE IF (PERP(3).EQ.PMIN0) THEN | |
130 | UL(1) = 0. | |
131 | UL(2) = 0. | |
132 | UL(3) = SIGN(ONE,DXL(3)) | |
133 | ELSE | |
134 | LNTFOU=.TRUE. | |
135 | ENDIF | |
136 | * | |
137 | * *** Case of the TUBE, TUBeSection: | |
138 | ELSE IF (ISH.EQ.NSTUBE.OR.ISH.EQ.NSTUBS) THEN | |
139 | RHO = SQRT(DXL(1)**2 + DXL(2)**2) | |
140 | PERP(1) = ABS(RHO-PAR(1)) | |
141 | PERP(2) = ABS(RHO-PAR(2)) | |
142 | PERP(3) = ABS(ABS(DXL(3))-PAR(3)) | |
143 | IF (ISH.EQ.NSTUBE) THEN | |
144 | PMIN0 = MIN(PERP(1),PERP(2),PERP(3)) | |
145 | ELSE | |
146 | PHI = ATAN2(DXL(2),DXL(1)) | |
147 | IF (PHI.LT.0.) PHI = PHI+TWOPI | |
148 | PHI1 = MOD(PAR(4)+FULL10,FULL)*DEGRAD | |
149 | PERP(4) = ABS(PHI-PHI1) | |
150 | IF (PERP(4).GT.PI) PERP(4) = TWOPI-PERP(4) | |
151 | PHI2 = MOD(PAR(5)+FULL10,FULL)*DEGRAD | |
152 | PERP(5) = ABS(PHI-PHI2) | |
153 | IF (PERP(5).GT.PI) PERP(5) = TWOPI-PERP(5) | |
154 | PERP(4) = PERP(4)*RHO | |
155 | PERP(5) = PERP(5)*RHO | |
156 | PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),PERP(5)) | |
157 | ENDIF | |
158 | IF (PERP(1).EQ.PMIN0) THEN | |
159 | UL(1) = -DXL(1)/RHO | |
160 | UL(2) = -DXL(2)/RHO | |
161 | UL(3) = 0. | |
162 | ELSE IF (PERP(2).EQ.PMIN0) THEN | |
163 | UL(1) = DXL(1)/RHO | |
164 | UL(2) = DXL(2)/RHO | |
165 | UL(3) = 0. | |
166 | ELSE IF (PERP(3).EQ.PMIN0) THEN | |
167 | UL(1) = 0. | |
168 | UL(2) = 0. | |
169 | UL(3) = SIGN(ONE,DXL(3)) | |
170 | ELSE IF (PERP(4).EQ.PMIN0) THEN | |
171 | UL(1) = SIN(PHI1) | |
172 | UL(2) = -COS(PHI1) | |
173 | UL(3) = 0. | |
174 | ELSE IF (PERP(5).EQ.PMIN0) THEN | |
175 | UL(1) = -SIN(PHI2) | |
176 | UL(2) = COS(PHI2) | |
177 | UL(3) = 0. | |
178 | ELSE | |
179 | LNTFOU=.TRUE. | |
180 | ENDIF | |
181 | * | |
182 | * *** Case of the CONE, CONeSection: | |
183 | ELSE IF (ISH.EQ.NSCONE.OR.ISH.EQ.NSCONS) THEN | |
184 | RHO = SQRT(DXL(1)**2 + DXL(2)**2) | |
185 | TAN1 = (PAR(4)-PAR(2))/(TWO*PAR(1)) | |
186 | SEC1 = SQRT(ONE+TAN1**2) | |
187 | U1 = RHO-DXL(3)*TAN1 | |
188 | U1L = PAR(4)-PAR(1)*TAN1 | |
189 | TAN2 = (PAR(5)-PAR(3))/(TWO*PAR(1)) | |
190 | SEC2 = SQRT(ONE+TAN2**2) | |
191 | U2 = RHO-DXL(3)*TAN2 | |
192 | U2L = PAR(5)-PAR(1)*TAN2 | |
193 | PERP(1) = ABS(ABS(DXL(3))-PAR(1)) | |
194 | PERP(2) = ABS(U1-U1L)/SEC1 | |
195 | PERP(3) = ABS(U2-U2L)/SEC2 | |
196 | IF (ISH.EQ.NSCONE) THEN | |
197 | PMIN0 = MIN(PERP(1),PERP(2),PERP(3)) | |
198 | ELSE | |
199 | PHI = ATAN2(DXL(2),DXL(1)) | |
200 | IF (PHI.LT.0.) PHI = PHI+TWOPI | |
201 | PHI1 = MOD(PAR(6)+FULL10,FULL)*DEGRAD | |
202 | PERP(4) = ABS(PHI-PHI1) | |
203 | IF (PERP(4).GT.PI) PERP(4) = TWOPI-PERP(4) | |
204 | PHI2 = MOD(PAR(7)+FULL10,FULL)*DEGRAD | |
205 | PERP(5) = ABS(PHI-PHI2) | |
206 | IF (PERP(5).GT.PI) PERP(5) = TWOPI-PERP(5) | |
207 | PERP(4) = PERP(4)*RHO | |
208 | PERP(5) = PERP(5)*RHO | |
209 | PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),PERP(5)) | |
210 | ENDIF | |
211 | IF (PERP(1).EQ.PMIN0) THEN | |
212 | UL(1) = 0. | |
213 | UL(2) = 0. | |
214 | UL(3) = SIGN(ONE,DXL(3)) | |
215 | ELSE IF (PERP(2).EQ.PMIN0) THEN | |
216 | RHO = RHO*SEC1 | |
217 | UL(1) = -DXL(1)/RHO | |
218 | UL(2) = -DXL(2)/RHO | |
219 | UL(3) = TAN1/SEC1 | |
220 | ELSE IF (PERP(3).EQ.PMIN0) THEN | |
221 | RHO = RHO*SEC2 | |
222 | UL(1) = DXL(1)/RHO | |
223 | UL(2) = DXL(2)/RHO | |
224 | UL(3) = -TAN2/SEC2 | |
225 | ELSE IF (PERP(4).EQ.PMIN0) THEN | |
226 | UL(1) = SIN(PHI1) | |
227 | UL(2) = -COS(PHI1) | |
228 | UL(3) = 0. | |
229 | ELSE IF (PERP(5).EQ.PMIN0) THEN | |
230 | UL(1) = -SIN(PHI2) | |
231 | UL(2) = COS(PHI2) | |
232 | UL(3) = 0. | |
233 | ELSE | |
234 | LNTFOU=.TRUE. | |
235 | ENDIF | |
236 | * | |
237 | * *** Case of the PolyCONe: | |
238 | ELSE IF (ISH.EQ.NSPCON) THEN | |
239 | PERP(1) = ABS(DXL(3)-PAR(4)) | |
240 | DO 400 I=7,NPAR,3 | |
241 | PERP(2) = ABS(DXL(3)-PAR(I)) | |
242 | IF (PERP(2).GT.PERP(1)) GOTO 401 | |
243 | PERP(1) = PERP(2) | |
244 | 400 CONTINUE | |
245 | 401 I = I-3 | |
246 | IF (I.GT.4) THEN | |
247 | PERP(1) = 100. | |
248 | RHO = SQRT(DXL(1)**2 + DXL(2)**2) | |
249 | DZ = PAR(I)-PAR(I-3)+1.e-10 | |
250 | TAN1 = (PAR(I+1)-PAR(I-2))/DZ | |
251 | SEC1 = SQRT(ONE+TAN1**2) | |
252 | U1 = RHO-DXL(3)*TAN1 | |
253 | U1L = PAR(I+1)-PAR(I)*TAN1 | |
254 | TAN2 = (PAR(I+2)-PAR(I-1))/DZ | |
255 | SEC2 = SQRT(ONE+TAN2**2) | |
256 | U2 = RHO-DXL(3)*TAN2 | |
257 | U2L = PAR(I+2)-PAR(I)*TAN2 | |
258 | GUARD = MAX(DXL(3)-PAR(I),ZERO) | |
259 | PERP(3) = ABS(U1-U1L)/SEC1 + GUARD*SEC1 | |
260 | PERP(4) = ABS(U2-U2L)/SEC2 + GUARD*SEC2 | |
261 | ELSE | |
262 | PERP(3) = 100. | |
263 | PERP(4) = 100. | |
264 | ENDIF | |
265 | IF (I.LT.NPAR-2) THEN | |
266 | PERP(2) = 100. | |
267 | RHO = SQRT(DXL(1)**2 + DXL(2)**2) | |
268 | DZ = PAR(I+3)-PAR(I)+1.e-10 | |
269 | TAN3 = (PAR(I+4)-PAR(I+1))/DZ | |
270 | SEC3 = SQRT(ONE+TAN3**2) | |
271 | U1 = RHO-DXL(3)*TAN3 | |
272 | U1L = PAR(I+1)-PAR(I)*TAN3 | |
273 | TAN4 = (PAR(I+5)-PAR(I+2))/DZ | |
274 | SEC4 = SQRT(ONE+TAN4**2) | |
275 | U2 = RHO-DXL(3)*TAN4 | |
276 | U2L = PAR(I+2)-PAR(I)*TAN4 | |
277 | GUARD = MAX(PAR(I)-DXL(3),ZERO) | |
278 | PERP(5) = ABS(U1-U1L)/SEC3 + GUARD*SEC3 | |
279 | PERP(6) = ABS(U2-U2L)/SEC4 + GUARD*SEC4 | |
280 | ELSE | |
281 | PERP(5) = 100. | |
282 | PERP(6) = 100. | |
283 | ENDIF | |
284 | PHI = ATAN2(DXL(2),DXL(1)) | |
285 | IF (PHI.LT.0.) PHI = PHI+TWOPI | |
286 | PHI1 = MOD(PAR(1)+FULL10,FULL)*DEGRAD | |
287 | PERP(7) = ABS(PHI-PHI1) | |
288 | IF (PERP(7).GT.PI) PERP(7) = TWOPI-PERP(7) | |
289 | PHI2 = MOD(PAR(1)+PAR(2)+FULL10,FULL)*DEGRAD | |
290 | PERP(8) = ABS(PHI-PHI2) | |
291 | IF (PERP(8).GT.PI) PERP(8) = TWOPI-PERP(8) | |
292 | PERP(7) = PERP(7)*RHO | |
293 | PERP(8) = PERP(8)*RHO | |
294 | PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4), | |
295 | + PERP(5),PERP(6),PERP(7),PERP(8)) | |
296 | IF (PERP(1).EQ.PMIN0) THEN | |
297 | UL(1) = 0. | |
298 | UL(2) = 0. | |
299 | UL(3) = -1. | |
300 | ELSE IF (PERP(2).EQ.PMIN0) THEN | |
301 | UL(1) = 0. | |
302 | UL(2) = 0. | |
303 | UL(3) = 1. | |
304 | ELSE IF (PERP(3).EQ.PMIN0) THEN | |
305 | RHO = RHO*SEC1 | |
306 | UL(1) = -DXL(1)/RHO | |
307 | UL(2) = -DXL(2)/RHO | |
308 | UL(3) = TAN1/SEC1 | |
309 | ELSE IF (PERP(4).EQ.PMIN0) THEN | |
310 | RHO = RHO*SEC2 | |
311 | UL(1) = DXL(1)/RHO | |
312 | UL(2) = DXL(2)/RHO | |
313 | UL(3) = -TAN2/SEC2 | |
314 | ELSE IF (PERP(5).EQ.PMIN0) THEN | |
315 | RHO = RHO*SEC3 | |
316 | UL(1) = -DXL(1)/RHO | |
317 | UL(2) = -DXL(2)/RHO | |
318 | UL(3) = TAN3/SEC3 | |
319 | ELSE IF (PERP(6).EQ.PMIN0) THEN | |
320 | RHO = RHO*SEC4 | |
321 | UL(1) = DXL(1)/RHO | |
322 | UL(2) = DXL(2)/RHO | |
323 | UL(3) = -TAN4/SEC4 | |
324 | ELSE IF (PERP(7).EQ.PMIN0) THEN | |
325 | UL(1) = SIN(PHI1) | |
326 | UL(2) = -COS(PHI1) | |
327 | UL(3) = 0. | |
328 | ELSE IF (PERP(8).EQ.PMIN0) THEN | |
329 | UL(1) = -SIN(PHI2) | |
330 | UL(2) = COS(PHI2) | |
331 | UL(3) = 0. | |
332 | ELSE | |
333 | LNTFOU=.TRUE. | |
334 | ENDIF | |
335 | * | |
336 | * *** Case of the PolyGON: | |
337 | ELSE IF (ISH.EQ.NSPGON) THEN | |
338 | RHO = SQRT(DXL(1)**2+DXL(2)**2) | |
339 | PHI = ATAN2(DXL(2),DXL(1)) | |
340 | IF (PHI.LT.0.) PHI = PHI+TWOPI | |
341 | DPHI = MOD(PHI*RADDEG-PAR(1)+FULL10,FULL) | |
342 | PDIV = PAR(2)/PAR(3) | |
343 | DSECT = INT(DPHI/PDIV + ONE) | |
344 | IF (DSECT.GT.PAR(3)) THEN | |
345 | IF (DPHI.GT.(180.+PAR(2)*DBY2)) THEN | |
346 | DSECT = ONE | |
347 | ELSE | |
348 | DSECT = PAR(3) | |
349 | ENDIF | |
350 | ENDIF | |
351 | PHI0 = MOD(PAR(1)+(DSECT-DBY2)*PDIV+FULL10,FULL)*DEGRAD | |
352 | SPHI0 = SIN(PHI0) | |
353 | CPHI0 = COS(PHI0) | |
354 | U0 = DXL(1)*CPHI0 + DXL(2)*SPHI0 | |
355 | V0 = DXL(2)*CPHI0 - DXL(1)*SPHI0 | |
356 | PERP(1) = ABS(DXL(3)-PAR(5)) | |
357 | DO 500 I=8,NPAR,3 | |
358 | PERP(2) = ABS(DXL(3)-PAR(I)) | |
359 | IF (PERP(2).GT.PERP(1)) GOTO 501 | |
360 | PERP(1) = PERP(2) | |
361 | 500 CONTINUE | |
362 | 501 I = I-3 | |
363 | IF (I.GT.5) THEN | |
364 | PERP(1) = 100. | |
365 | DZ = PAR(I)-PAR(I-3)+1.e-10 | |
366 | TAN1 = (PAR(I+1)-PAR(I-2))/DZ | |
367 | SEC1 = SQRT(ONE+TAN1**2) | |
368 | U1 = U0-DXL(3)*TAN1 | |
369 | U1L = PAR(I+1)-PAR(I)*TAN1 | |
370 | TAN2 = (PAR(I+2)-PAR(I-1))/DZ | |
371 | SEC2 = SQRT(ONE+TAN2**2) | |
372 | U2 = U0-DXL(3)*TAN2 | |
373 | U2L = PAR(I+2)-PAR(I)*TAN2 | |
374 | GUARD = MAX(DXL(3)-PAR(I),ZERO) | |
375 | PERP(3) = ABS(U1-U1L)/SEC1 + GUARD*SEC1 | |
376 | PERP(4) = ABS(U2-U2L)/SEC2 + GUARD*SEC2 | |
377 | ELSE | |
378 | PERP(3) = 100. | |
379 | PERP(4) = 100. | |
380 | ENDIF | |
381 | IF (I.LT.NPAR-2) THEN | |
382 | PERP(2) = 100. | |
383 | DZ = PAR(I+3)-PAR(I)+1.e-10 | |
384 | TAN3 = (PAR(I+4)-PAR(I+1))/DZ | |
385 | SEC3 = SQRT(ONE+TAN3**2) | |
386 | U1 = U0-DXL(3)*TAN3 | |
387 | U1L = PAR(I+1)-PAR(I)*TAN3 | |
388 | TAN4 = (PAR(I+5)-PAR(I+2))/DZ | |
389 | SEC4 = SQRT(ONE+TAN4**2) | |
390 | U2 = U0-DXL(3)*TAN4 | |
391 | U2L = PAR(I+2)-PAR(I)*TAN4 | |
392 | GUARD = MAX(PAR(I)-DXL(3),ZERO) | |
393 | PERP(5) = ABS(U1-U1L)/SEC3 + GUARD*SEC3 | |
394 | PERP(6) = ABS(U2-U2L)/SEC4 + GUARD*SEC4 | |
395 | ELSE | |
396 | PERP(5) = 100. | |
397 | PERP(6) = 100. | |
398 | ENDIF | |
399 | PHI1 = MOD(PAR(1)+FULL10,FULL)*DEGRAD | |
400 | PERP(7) = ABS(PHI-PHI1) | |
401 | IF (PERP(7).GT.PI) PERP(7) = TWOPI-PERP(7) | |
402 | PHI2 = MOD(PAR(1)+PAR(2)+FULL10,FULL)*DEGRAD | |
403 | PERP(8) = ABS(PHI-PHI2) | |
404 | IF (PERP(8).GT.PI) PERP(8) = TWOPI-PERP(8) | |
405 | PERP(7) = PERP(7)*RHO | |
406 | PERP(8) = PERP(8)*RHO | |
407 | PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4), | |
408 | + PERP(5),PERP(6),PERP(7),PERP(8)) | |
409 | IF (PERP(1).EQ.PMIN0) THEN | |
410 | UL(1) = 0. | |
411 | UL(2) = 0. | |
412 | UL(3) = -1. | |
413 | ELSE IF (PERP(2).EQ.PMIN0) THEN | |
414 | UL(1) = 0. | |
415 | UL(2) = 0. | |
416 | UL(3) = 1. | |
417 | ELSE IF (PERP(3).EQ.PMIN0) THEN | |
418 | FACT = ONE/SEC1 | |
419 | UL(1) = -CPHI0*FACT | |
420 | UL(2) = -SPHI0*FACT | |
421 | UL(3) = TAN1*FACT | |
422 | ELSE IF (PERP(4).EQ.PMIN0) THEN | |
423 | FACT = ONE/SEC2 | |
424 | UL(1) = CPHI0*FACT | |
425 | UL(2) = SPHI0*FACT | |
426 | UL(3) = -TAN2*FACT | |
427 | ELSE IF (PERP(5).EQ.PMIN0) THEN | |
428 | FACT = ONE/SEC3 | |
429 | UL(1) = -CPHI0*FACT | |
430 | UL(2) = -SPHI0*FACT | |
431 | UL(3) = TAN3*FACT | |
432 | ELSE IF (PERP(6).EQ.PMIN0) THEN | |
433 | FACT = ONE/SEC4 | |
434 | UL(1) = CPHI0*FACT | |
435 | UL(2) = SPHI0*FACT | |
436 | UL(3) = -TAN4*FACT | |
437 | ELSE IF (PERP(7).EQ.PMIN0) THEN | |
438 | UL(1) = SIN(PHI1) | |
439 | UL(2) = -COS(PHI1) | |
440 | UL(3) = 0. | |
441 | ELSE IF (PERP(8).EQ.PMIN0) THEN | |
442 | UL(1) = -SIN(PHI2) | |
443 | UL(2) = COS(PHI2) | |
444 | UL(3) = 0. | |
445 | ELSE | |
446 | LNTFOU=.TRUE. | |
447 | ENDIF | |
448 | * | |
449 | * *** Case of the SPHEre: | |
450 | ELSE IF (ISH.EQ.NSSPHE) THEN | |
451 | R = SQRT(DXL(1)**2+DXL(2)**2+DXL(3)**2) | |
452 | RHO = SQRT(DXL(1)**2+DXL(2)**2) | |
453 | THE = ATAN2(RHO,DXL(3)) | |
454 | PHI = ATAN2(DXL(2),DXL(1)) | |
455 | IF (PHI.LT.0.) PHI = PHI+TWOPI | |
456 | THE1 = MOD(PAR(3)+FULL10,FULL)*DEGRAD | |
457 | THE2 = MOD(PAR(4)+FULL10,FULL)*DEGRAD | |
458 | PHI1 = MOD(PAR(5)+FULL10,FULL)*DEGRAD | |
459 | PHI2 = MOD(PAR(6)+FULL10,FULL)*DEGRAD | |
460 | PERP(1) = ABS(R-PAR(1)) | |
461 | PERP(2) = ABS(R-PAR(2)) | |
462 | PERP(3) = ABS(THE-THE1)*R | |
463 | PERP(4) = ABS(THE-THE2)*R | |
464 | PERP(5) = ABS(PHI-PHI1) | |
465 | IF (PERP(5).GT.PI) PERP(5) = TWOPI-PERP(5) | |
466 | PERP(5) = PERP(5)*RHO | |
467 | PERP(6) = ABS(PHI-PHI2) | |
468 | IF (PERP(6).GT.PI) PERP(6) = TWOPI-PERP(6) | |
469 | PERP(6) = PERP(6)*RHO | |
470 | PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),PERP(5),PERP(6)) | |
471 | IF (PERP(1).EQ.PMIN0) THEN | |
472 | RINV = ONE/R | |
473 | UL(1) = -DXL(1)*RINV | |
474 | UL(2) = -DXL(2)*RINV | |
475 | UL(3) = -DXL(3)*RINV | |
476 | ELSE IF (PERP(2).EQ.PMIN0) THEN | |
477 | RINV = ONE/R | |
478 | UL(1) = DXL(1)*RINV | |
479 | UL(2) = DXL(2)*RINV | |
480 | UL(3) = DXL(3)*RINV | |
481 | ELSE IF (PERP(3).EQ.PMIN0) THEN | |
482 | UL(1) = -COS(THE1)*COS(PHI) | |
483 | UL(2) = -COS(THE1)*SIN(PHI) | |
484 | UL(3) = +SIN(THE1) | |
485 | ELSE IF (PERP(4).EQ.PMIN0) THEN | |
486 | UL(1) = +COS(THE2)*COS(PHI) | |
487 | UL(2) = +COS(THE2)*SIN(PHI) | |
488 | UL(3) = -SIN(THE2) | |
489 | ELSE IF (PERP(5).EQ.PMIN0) THEN | |
490 | UL(1) = +SIN(PHI1) | |
491 | UL(2) = -COS(PHI1) | |
492 | UL(3) = 0 | |
493 | ELSE IF (PERP(6).EQ.PMIN0) THEN | |
494 | UL(1) = -SIN(PHI2) | |
495 | UL(2) = +COS(PHI2) | |
496 | UL(3) = 0 | |
497 | ELSE | |
498 | LNTFOU=.TRUE. | |
499 | ENDIF | |
500 | * | |
501 | * *** Case of the PARAllelpiped: | |
502 | *************************************************************** | |
503 | * Warning: the parameters for this shape are NOT stored in * | |
504 | * the data structure as the user supplies them. Rather, the * | |
505 | * user supplies PAR(4)=alph, PAR(5)=the, PAR(6)=phi, and the * | |
506 | * data structure contains PAR(4)=Tan(alph), PAR(5)=Tan(the)* * | |
507 | * Cos(phi), PAR(6)=Tan(the)*Sin(phi). * | |
508 | *************************************************************** | |
509 | ELSE IF (ISH.EQ.NSPARA) THEN | |
510 | DX = PAR(5) | |
511 | DY = PAR(6) | |
512 | U0 = DXL(1)-DX*DXL(3) | |
513 | V0 = DXL(2)-DY*DXL(3) | |
514 | CALPH = ONE/SQRT(ONE+PAR(4)**2) | |
515 | SALPH = -CALPH*PAR(4) | |
516 | U1 = U0*CALPH+V0*SALPH | |
517 | U1L = PAR(1)*CALPH | |
518 | PERP(1) = ABS(ABS(U1)-U1L) | |
519 | PERP(2) = ABS(ABS(V0)-PAR(2)) | |
520 | PERP(3) = ABS(ABS(DXL(3))-PAR(3)) | |
521 | PMIN0 = MIN(PERP(1),PERP(2),PERP(3)) | |
522 | IF (PERP(1).EQ.PMIN0) THEN | |
523 | DU = DX*CALPH+DY*SALPH | |
524 | FACT = SIGN(ONE/SQRT(ONE+DU**2),U1) | |
525 | UL(1) = CALPH*FACT | |
526 | UL(2) = SALPH*FACT | |
527 | UL(3) = -DU*FACT | |
528 | ELSE IF (PERP(2).EQ.PMIN0) THEN | |
529 | FACT = SIGN(ONE/SQRT(ONE+DY**2),V0) | |
530 | UL(1) = 0. | |
531 | UL(2) = FACT | |
532 | UL(3) = -DY*FACT | |
533 | ELSE IF (PERP(3).EQ.PMIN0) THEN | |
534 | UL(1) = 0. | |
535 | UL(2) = 0. | |
536 | UL(3) = SIGN(ONE,DXL(3)) | |
537 | ELSE | |
538 | LNTFOU=.TRUE. | |
539 | ENDIF | |
540 | * | |
541 | * *** Case of the trapezoid TRD1 | |
542 | ELSE IF (ISH.EQ.NSTRD1) THEN | |
543 | DZ = TWO*PAR(4)+1.e-10 | |
544 | TAN1 = (PAR(2)-PAR(1))/DZ | |
545 | SEC1 = SQRT(ONE+TAN1**2) | |
546 | U1 = ABS(DXL(1))-DXL(3)*TAN1 | |
547 | U1L = PAR(2)-PAR(4)*TAN1 | |
548 | PERP(1) = ABS(U1-U1L)/SEC1 | |
549 | PERP(2) = ABS(ABS(DXL(2))-PAR(3)) | |
550 | PERP(3) = ABS(ABS(DXL(3))-PAR(4)) | |
551 | PMIN0 = MIN(PERP(1),PERP(2),PERP(3)) | |
552 | IF (PERP(1).EQ.PMIN0) THEN | |
553 | FACT = ONE/SEC1 | |
554 | UL(1) = SIGN(FACT,DXL(1)) | |
555 | UL(2) = 0. | |
556 | UL(3) = -TAN1*FACT | |
557 | ELSE IF (PERP(2).EQ.PMIN0) THEN | |
558 | UL(1) = 0. | |
559 | UL(2) = SIGN(ONE,DXL(2)) | |
560 | UL(3) = 0. | |
561 | ELSE IF (PERP(3).EQ.PMIN0) THEN | |
562 | UL(1) = 0. | |
563 | UL(2) = 0. | |
564 | UL(3) = SIGN(ONE,DXL(3)) | |
565 | ELSE | |
566 | LNTFOU=.TRUE. | |
567 | ENDIF | |
568 | * | |
569 | * *** Case of the trapezoid TRD2 | |
570 | ELSE IF (ISH.EQ.NSTRD2) THEN | |
571 | DZ = TWO*PAR(5)+1.e-10 | |
572 | TAN1 = (PAR(2)-PAR(1))/DZ | |
573 | SEC1 = SQRT(ONE+TAN1**2) | |
574 | U1 = ABS(DXL(1))-DXL(3)*TAN1 | |
575 | U1L = PAR(2)-PAR(5)*TAN1 | |
576 | TAN2 = (PAR(4)-PAR(3))/DZ | |
577 | SEC2 = SQRT(ONE+TAN2**2) | |
578 | U2 = ABS(DXL(2))-DXL(3)*TAN2 | |
579 | U2L = PAR(4)-PAR(5)*TAN2 | |
580 | PERP(1) = ABS(U1-U1L)/SEC1 | |
581 | PERP(2) = ABS(U2-U2L)/SEC2 | |
582 | PERP(3) = ABS(ABS(DXL(3))-PAR(5)) | |
583 | PMIN0 = MIN(PERP(1),PERP(2),PERP(3)) | |
584 | IF (PERP(1).EQ.PMIN0) THEN | |
585 | FACT = ONE/SEC1 | |
586 | UL(1) = SIGN(FACT,DXL(1)) | |
587 | UL(2) = 0. | |
588 | UL(3) = -TAN1*FACT | |
589 | ELSE IF (PERP(2).EQ.PMIN0) THEN | |
590 | FACT = ONE/SEC2 | |
591 | UL(1) = 0. | |
592 | UL(2) = SIGN(FACT,DXL(2)) | |
593 | UL(3) = -TAN2*FACT | |
594 | ELSE IF (PERP(3).EQ.PMIN0) THEN | |
595 | UL(1) = 0. | |
596 | UL(2) = 0. | |
597 | UL(3) = SIGN(ONE,DXL(3)) | |
598 | ELSE | |
599 | LNTFOU=.TRUE. | |
600 | ENDIF | |
601 | * | |
602 | * *** Case of the TRAPezoid | |
603 | *************************************************************** | |
604 | * Warning: the parameters for this shape are NOT stored in * | |
605 | * the data structure as the user supplies them. Rather, the * | |
606 | * user supplies PAR(2)=thet, PAR(3)=phi, PAR(7)=alp1, and * | |
607 | * PAR(11)=alp2, while the data structure contains PAR(2)= * | |
608 | * Tan(thet)*Cos(phi), PAR(3)=Tan(thet)*Sin(phi), PAR(7)= * | |
609 | * Tan(alp1), and PAR(11)=Tan(alp2). * | |
610 | *************************************************************** | |
611 | ELSE IF (ISH.EQ.NSTRAP) THEN | |
612 | PERP(1) = ABS(ABS(DXL(3))-PAR(1)) | |
613 | DX = PAR(2) | |
614 | DY = PAR(3) | |
615 | U0 = DX*DXL(3) | |
616 | V0 = DY*DXL(3) | |
617 | UU0 = DX*PAR(1) | |
618 | VV0 = DY*PAR(1) | |
619 | RAT = DXL(3)/PAR(1) | |
620 | RATL = (ONE-RAT)/TWO | |
621 | RATH = (ONE+RAT)/TWO | |
622 | H = PAR(4)*RATL+PAR(8)*RATH | |
623 | BL = PAR(5)*RATL+PAR(9)*RATH | |
624 | TL = PAR(6)*RATL+PAR(10)*RATH | |
625 | TALPH = PAR(7)*RATL+PAR(11)*RATH | |
626 | XWID = (TL+BL)/TWO | |
627 | TAN1 = TALPH+(TL-BL)/(TWO*H) | |
628 | SEC1 = SQRT(ONE+TAN1**2) | |
629 | U1 = DXL(1)-DXL(2)*TAN1 | |
630 | U1L = U0+XWID-V0*TAN1 | |
631 | TAN2 = TAN1-TWO*TALPH | |
632 | SEC2 = SQRT(ONE+TAN2**2) | |
633 | U2 = DXL(1)+DXL(2)*TAN2 | |
634 | U2L = U0-XWID+V0*TAN2 | |
635 | IF (DXL(3).LT.0) THEN | |
636 | DZ = PAR(1)-DXL(3)+1.e-10 | |
637 | UU = UU0+(PAR(9)+PAR(10))/TWO | |
638 | W1 = (UU-VV0*TAN1-U1L)/DZ | |
639 | UU = TWO*UU0-UU | |
640 | W2 = (UU+VV0*TAN2-U2L)/DZ | |
641 | ELSE | |
642 | DZ = -PAR(1)-DXL(3)+1.e-10 | |
643 | UU = -UU0+(PAR(5)+PAR(6))/TWO | |
644 | W1 = (UU+VV0*TAN1-U1L)/DZ | |
645 | UU = -TWO*UU0-UU | |
646 | W2 = (UU-VV0*TAN2-U2L)/DZ | |
647 | ENDIF | |
648 | W3 = DY+(PAR(8)-PAR(4))/(TWO*PAR(1)) | |
649 | W4 = TWO*DY-W3 | |
650 | SEW1 = SQRT(ONE+W1**2) | |
651 | SEW2 = SQRT(ONE+W2**2) | |
652 | SEW3 = SQRT(ONE+W3**2) | |
653 | SEW4 = SQRT(ONE+W4**2) | |
654 | PERP(2) = ABS(U1-U1L)/(SEC1*SEW1) | |
655 | PERP(3) = ABS(U2-U2L)/(SEC2*SEW2) | |
656 | PERP(4) = ABS(DXL(2)-V0-H)/SEW3 | |
657 | PERP(5) = ABS(DXL(2)-V0+H)/SEW4 | |
658 | PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),PERP(5)) | |
659 | IF (PERP(1).EQ.PMIN0) THEN | |
660 | UL(1) = 0. | |
661 | UL(2) = 0. | |
662 | UL(3) = SIGN(ONE,DXL(3)) | |
663 | ELSE IF (PERP(2).EQ.PMIN0) THEN | |
664 | FACT = ONE/(SEC1*SEW1) | |
665 | UL(1) = FACT | |
666 | UL(2) = -TAN1*FACT | |
667 | UL(3) = -W1/SEW1 | |
668 | ELSE IF (PERP(3).EQ.PMIN0) THEN | |
669 | FACT = ONE/(SEC2*SEW2) | |
670 | UL(1) = -FACT | |
671 | UL(2) = -TAN2*FACT | |
672 | UL(3) = W2/SEW2 | |
673 | ELSE IF (PERP(4).EQ.PMIN0) THEN | |
674 | FACT = ONE/SEW3 | |
675 | UL(1) = 0. | |
676 | UL(2) = FACT | |
677 | UL(3) = -W3*FACT | |
678 | ELSE IF (PERP(5).EQ.PMIN0) THEN | |
679 | FACT = ONE/SEW4 | |
680 | UL(1) = 0. | |
681 | UL(2) = -FACT | |
682 | UL(3) = W4*FACT | |
683 | ELSE | |
684 | LNTFOU=.TRUE. | |
685 | ENDIF | |
686 | * | |
687 | * *** everything else (currently NOT IMPLEMENTED) | |
688 | ELSE | |
689 | WRITE(CHMAIL,10100) ISH | |
690 | CALL GMAIL(0,0) | |
691 | IERR = 1 | |
692 | GOTO 999 | |
693 | ENDIF | |
694 | ||
695 | IF(LNTFOU) THEN | |
696 | WRITE(CHMAIL,10000) ISH | |
697 | CALL GMAIL(0,0) | |
698 | IERR = 1 | |
699 | ELSE | |
700 | * | |
701 | * *** Transform back into the MCS | |
702 | CALL GDTOM(UL,U,2) | |
703 | IERR = 0 | |
704 | ENDIF | |
705 | ||
706 | 10000 FORMAT(' GGPERP - geometry check error for shape #',I2,'!') | |
707 | 10100 FORMAT(' GGPERP - non implemented for shape #',I2) | |
708 | 999 END |