1 *
2 * \$Id\$
3 *
4 * \$Log\$
5 * Revision 1.1.1.1  1995/10/24 10:20:19  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.25  by  S.Giani
11 *-- Author :
12       SUBROUTINE GD3D3D(XIN,NPOINT,XOUT,INV)
13 C.
14 C.    ******************************************************************
15 C.    *                                                                *
16 C.    *       Routine to transform the 3D points given by array XIN    *
17 C.    *       into the 3D points in XOUT:                              *
18 C.    *       if INV=0 the Z-axis for XOUT is the line of sight for XIN*
19 C.    *       if INV=1 the Z-axis for XIN is the line of sight for XOUT*
20 C.    *                                                                *
21 C.    *    ==>Called by : <USER>, GDCUT                                *
22 C.    *       Author : P.Zanarini   *********                          *
23 C.    *                                                                *
24 C.    ******************************************************************
25 C.
26 #include "geant321/gcvolu.inc"
27 #include "geant321/gcdraw.inc"
28       DIMENSION XC(3),XOUT(3,1),VL(3),VM(3),VN(3),XIN(3,1)
29       SAVE VL,VM,VN
30 C.
31 C.    ------------------------------------------------------------------
32 C.
33       IF(NPOINT.EQ.0)GO TO 999
34       N=NPOINT
35       IF(NPOINT.LT.0)N=-NPOINT
36       IF(NGVIEW.NE.0) GO TO 40
37 C
38 C             First call compute the rotation matrix
39 C
40       PH     = ABS(MOD(GPHI,360.))
41       THET   = ABS(MOD(GTHETA,360.))
42       IF(THET.LE.180.)GO TO 10
43       PH     = PH + 180.
44       THET   = 360. - THET
45 C
46    10 ST     = SIN(THET * 0.017453)
47       CT     = COS(THET * 0.017453)
48       SP     = SIN(PH * 0.017453)
49       CP     = COS(PH * 0.017453)
50 C
51 C             VN is new nu axis
52 C
53       VN(1)  = ST * CP
54       VN(2)  = ST * SP
55       VN(3)  = CT
56 C
57       IF(ABS(VN(2)).GT.0.99999)GO TO 20
58 C
59       VM(1)  = 0.
60       VM(2)  = 1.
61       VM(3)  = 0.
62 C
63 C             Define new lambda axis
64 C
65       CALL CROSS(VM,VN,VL)
66       CALL VUNIT(VL,VL,3)
67 C
68 C             Define new mu axis
69 C
70       CALL CROSS(VN,VL,VM)
71       GO TO 30
72 C
73 C             Special case when observer line of sight is along mu:
74 C             in this case one chooses arbitrarily the vertical axis of
75 C             plane of projection as the lambda axis and the horizontal
76 C             as the nu axis
77 C
78    20 VL(1)  = 0.
79       VL(2)  = 0.
80       VL(3)  = 1.
81       VM(1)  = 1.
82       VM(2)  = 0.
83       VM(3)  = 0.
84    30 CONTINUE
85 C
86       NGVIEW=1
87 C
88    40 CONTINUE
89 C
90 C             Begin of a normal call (i.e. with NGVIEW=1)
91 C
92       DO 70  I=1,N
93 C
94 C             Loop over the points
95 C
96          IF (NPOINT.LT.0) THEN
97 C
98 C             NPOINT < 0 : XIN is in general reference system
99 C
100             IF (INV.EQ.1) THEN
101                XC(1)=VL(1)*XIN(1,I)+VM(1)*XIN(2,I)+VN(1)*XIN(3,I)
102                XC(2)=VL(2)*XIN(1,I)+VM(2)*XIN(2,I)+VN(2)*XIN(3,I)
103                XC(3)=VL(3)*XIN(1,I)+VM(3)*XIN(2,I)+VN(3)*XIN(3,I)
104             ELSE
105                XC(1)=XIN(1,I)*VL(1)+XIN(2,I)*VL(2)+XIN(3,I)*VL(3)
106                XC(2)=XIN(1,I)*VM(1)+XIN(2,I)*VM(2)+XIN(3,I)*VM(3)
107                XC(3)=XIN(1,I)*VN(1)+XIN(2,I)*VN(2)+XIN(3,I)*VN(3)
108             ENDIF
109 C
110             DO 50  J=1,3
111                XOUT(J,I)=XC(J)
112    50       CONTINUE
113 C
114          ELSE
115 C
116 C             NPOINT > 0 : XIN is in volume reference system
117 C
118             CALL GINROT(XIN(1,I),GRMAT(1,NLEVEL),XC)
119             DO 60  J=1,3
120                XOUT(J,I)=XC(J)+GTRAN(J,NLEVEL)
121    60       CONTINUE
122 C
123          ENDIF
124 C
125    70 CONTINUE
126 C
127   999 RETURN
128       END