1      SUBROUTINE TTCBOXES(TTC)
2      IMPLICIT REAL*8 (A-H,O-Z)
3      REAL*8 TTCVEC(3,6,8),TTC(3,13),TTC12(3,2,8),TTCPOINT(3,7,8)
4      INTEGER TTCBOX(5,8),NTTCVEC(8),TTCPAIR(2,6,8)
5      REAL*8 MINTOL,NULL
6
7      COMMON/TTCBOXS/ TTCVEC,NTTCVEC,TTCPAIR,TTCPOINT
8
9      PARAMETER (CON13=1.0d0/3.0d0, CON23=2.0d0/3.0d0, MINTOL=1.d-6,
10     1     NULL=0.0d0, TWO=2.0d0, TWON=-2.0d0)
11
12C     ************************************************************************
13C     TTCVEC --- vectors for all 8 boxes
14C     TTC    --- coordinates for RTTC/HTTC points
15C     TTC12  --- to TTC points we add one point in front and one on the back
16C     TTCBOX --- indexes of TTC points, which are in box
17C     NTTCVEC--- number of vectors in box
18C     TTCPAIR--- which pair of points are basis for vector
19C     ************************************************************************
20C     ***** to each "ttcbox box" we add first and last point,
21C     here are coordinates *****
22      DATA TTC12/
23     1     NULL, TWO, NULL,
24     1     TWO, NULL, NULL,
25     2     TWON,TWON, NULL,
26     2     TWO, NULL, NULL,
27     3     TWON,TWON, NULL,
28     3     TWO, TWO, NULL,
29     4     NULL,TWON, NULL,
30     4     TWO, TWO, NULL,
31     5     NULL,TWON, NULL,
32     5     TWON, NULL, NULL,
33     6     TWO, TWO, NULL,
34     6     TWON, NULL, NULL,
35     7     TWO, TWO, NULL,
36     7     TWON,TWON, NULL,
37     8     NULL, TWO, NULL,
38     8     TWON,TWON, NULL/
39
40C     *** indexes of "BOXs" POINTS ***
41      DATA TTCBOX/
42     1     3,4,5,6,1,
43     2     5,6,1,0,0,
44     3     5,6,1,2,0,
45     4     6,1,2,0,0,
46     5     6,1,2,3,4,
47     6     2,3,4,0,0,
48     7     2,3,4,5,0,
49     8     3,4,5,0,0/
50C     *** how many points are in each box ***
51      DATA NTTCVEC/6,4,5,4,6,4,5,4/
52
53C     *** from TTC & TTC12 make TTCPOINTs ***
54C     *** insert first point ***
55      DO I=1,8
56         DO K=1,3
57            TTCPOINT(K,1,I)=TTC12(K,1,I)
58         ENDDO
59      ENDDO
60C     *** TTC -> TTCPOINT ***
61      DO I=1,8   !four boxes
62         DO J=2,NTTCVEC(I)       !vectors per boxes
63            DO K=1,3  !xyz coor of point
64               TTCPOINT(K,J,I)=TTC(K,TTCBOX(J-1,I))
65            ENDDO
66         ENDDO
67      ENDDO
68C     *** insert last point ***
69      DO I=1,8
70         DO K=1,3
71            TTCPOINT(K,NTTCVEC(I)+1,I)=TTC12(K,2,I)
72         ENDDO
73      ENDDO
74
75C     FROM POINTS -> VECTORS
76c      print *,'TTCVEC>'
77      DO I=1,8   !four boxes
78         DO J=1,NTTCVEC(I)       !vectors per boxes
79            DO K=1,3  !xyz coor of point
80               TTCVEC(K,J,I)=TTCPOINT(K,J+1,I)-TTCPOINT(K,J,I)
81            ENDDO
82            write(6,*) (TTCVEC(k,j,i),k=1,3)
83         ENDDO
84c         print *,'------------------------------'
85      ENDDO
86
87c      print *,'TTCPOINT'
88c      do i=1,8
89c         do j=1,7
90c            print *,(ttcpoint(k,j,i),k=1,3)
91c         enddo
92c         print *,'-------------------------------'
93c      enddo
94
95C     **** VECPAIR --- prvi vector je vector, ki skupaj z C tvori ravnino,
96C     ****             drugi pa pove katera stran je "prava"
97      DATA TTCPAIR/
98C     *** "BOX1" PAIRS ***
99     1     1,-2,
100     2     2,3,
101     3     3,4,
102     4     4,5,
103     5     5,-6,
104     6     6,5,
105C     *** "BOX2" PAIRS ***
106     1     1,-2,
107     2     2,3,
108     3     3,-4,
109     4     4,3,
110     5     0,0,
111     6     0,0,
112C     *** "BOX3" PAIRS ***
113     1     1,-2,
114     2     2,3,
115     3     3,4,
116     4     4,-5,
117     5     5,4,
118     6     0,0,
119C     *** "BOX4" PAIRS ***
120     1     1,-2,
121     2     2,3,
122     3     3,-4,
123     4     4,3,
124     5     0,0,
125     6     0,0,
126C     *** "BOX5" PAIRS
127     1     1,-2,
128     2     2,3,
129     3     3,4,
130     4     4,5,
131     5     5,-6,
132     6     6,5,
133C     *** "BOX6" PAIRS
134     1     1,-2,
135     2     2,3,
136     3     3,-4,
137     4     4,3,
138     5     0,0,
139     6     0,0,
140C     *** "BOX7" PAIRS
141     1     1,-2,
142     2     2,3,
143     3     3,4,
144     4     4,-5,
145     5     5,4,
146     6     0,0,
147C     *** "BOX8" PAIRS
148     1     1,-2,
149     2     2,3,
150     3     3,-4,
151     4     4,3,
152     5     0,0,
153     6     0,0/
154
155      RETURN
156      END
157
158
159C     tocka1(p1,p2,p3) pove kje je skatla, tocka2(r1,r2,r3) pa katero tocko
160C     ocenjujemo
161      LOGICAL FUNCTION TTCBOX(NBOX,P1,P2,P3,R1,R2,R3,DVEC)
162      IMPLICIT REAL*8 (A-H,O-Z)
163      REAL*8 TTCPOINT(3,7,8),TTCVEC(3,6,8),DVEC(3,3),IDVEC(3,3)
164      INTEGER NTTCVEC(8),TTCPAIR(2,6,8),TMPVECP1,TMPVECP2
165      LOGICAL ISINTTC,LTTC(2,8),LDONE
166
167      COMMON/TTCBOXS/ TTCVEC,NTTCVEC,TTCPAIR,TTCPOINT
168      COMMON/IVEC/ IDVEC  !inverse matrix of DVEC(3,3)
169
170      PARAMETER(CON13=1.0d0/3.0d0, CON23=2.0d0/3.0d0)
171
172      TTCBOX=.TRUE.
173      LDONE=.FALSE.
174
175C     x,y,z are "relative coord." of point in box
176      x=r1-p1
177      y=r2-p2
178      z=r3-p3
179C     transform (x,y,z) to DVEC basis
180      xx=x*IDVEC(1,1)+y*IDVEC(2,1)+z*IDVEC(3,1)
181      yy=x*IDVEC(1,2)+y*IDVEC(2,2)+z*IDVEC(3,2)
182      zz=x*IDVEC(1,3)+y*IDVEC(2,3)+z*IDVEC(3,3)
183c      print *,'xx,yy,zz>',xx,yy,zz
184
185      LTTC(1,1)=YY.GT.1.0d0
186      LTTC(2,1)=XX.GT.1.0d0
187      LTTC(1,2)=XX.LT.-1.0d0.AND.YY.LT.-1.0d0
188      LTTC(2,2)=XX.GT.1.0d0
189      LTTC(1,3)=XX.LT.-1.0d0.AND.YY.LT.-1.0d0
190      LTTC(2,3)=XX.GT.1.0d0.AND.YY.GT.1.0d0
191      LTTC(1,4)=YY.LT.-1.0d0
192      LTTC(2,4)=XX.GT.1.0d0.AND.YY.GT.1.0d0
193      LTTC(1,5)=YY.LT.-1.0d0
194      LTTC(2,5)=XX.LT.-1.0d0
195      LTTC(1,6)=XX.GT.1.0d0.AND.YY.GT.1.0d0
196      LTTC(2,6)=XX.LT.-1.0d0
197      LTTC(1,7)=XX.GT.1.0d0.AND.YY.GT.1.0d0
198      LTTC(2,7)=XX.LT.-1.0d0.AND.YY.LT.-1.0d0
199      LTTC(1,8)=YY.GT.1.0d0
200      LTTC(2,8)=XX.LT.-1.0d0.AND.YY.LT.-1.0d0
201
202C     *** first check 3th part of the box **
203      I=NTTCVEC(NBOX)
204      PX=P1+TTCPOINT(1,I,NBOX)*DVEC(1,1)+
205     1     TTCPOINT(2,I,NBOX)*DVEC(2,1)+TTCPOINT(3,I,NBOX)*DVEC(3,1)
206      PY=P2+TTCPOINT(1,I,NBOX)*DVEC(1,2)+
207     1     TTCPOINT(2,I,NBOX)*DVEC(2,2)+TTCPOINT(3,I,NBOX)*DVEC(3,2)
208      PZ=P3+TTCPOINT(1,I,NBOX)*DVEC(1,3)+
209     1     TTCPOINT(2,I,NBOX)*DVEC(2,3)+TTCPOINT(3,I,NBOX)*DVEC(3,3)
210C     is (R1,R2,R3) in the 3TH part of the box
211      if(lttc(2,nbox))then
212c         print *,'PART> 3th',xx,yy
213         ldone=.true.
214         if(.not.isinttc(nbox,ttcpair(1,i,nbox),ttcpair(2,i,nbox),
215     1        px,py,pz,r1,r2,r3,dvec)) ttcbox=.false.
216      endif
217
218      DO I=1,NTTCVEC(NBOX)-1
219c         print *,'NTTCVEC> i=',i,'; xx,yy,zz>',xx,yy,zz
220         TMPVECP1=TTCPAIR(1,I,NBOX)
221         TMPVECP2=TTCPAIR(2,I,NBOX)
222c         print *,'vec1=',tmpvecp1
223c         print *,'vec2=',tmpvecp2
224C     IS POINT(R1,R2,R3) ON THE RIGHT SIDE OF THE PLANE?????
225         PX=P1+TTCPOINT(1,I,NBOX)*DVEC(1,1)+
226     1        TTCPOINT(2,I,NBOX)*DVEC(2,1)+TTCPOINT(3,I,NBOX)*DVEC(3,1)
227         PY=P2+TTCPOINT(1,I,NBOX)*DVEC(1,2)+
228     1        TTCPOINT(2,I,NBOX)*DVEC(2,2)+TTCPOINT(3,I,NBOX)*DVEC(3,2)
229         PZ=P3+TTCPOINT(1,I,NBOX)*DVEC(1,3)+
230     1        TTCPOINT(2,I,NBOX)*DVEC(2,3)+TTCPOINT(3,I,NBOX)*DVEC(3,3)
231c         print *,'PXYZ> ',px,py,pz
232C     is (R1,R2,R3) in the 1ST part of the box;
233C     for 1ST part i must be 1 -> this means: first point and first vector
234         if(i.eq.1.and.lttc(1,nbox))then
235c            print *,'PART> 1st',xx,yy
236            ldone=.true.
237            if(.not.isinttc(nbox,tmpvecp1,tmpvecp2,
238     1           px,py,pz,r1,r2,r3,dvec)) ttcbox=.false.
239         endif
240C     is (R1,R2,R3) in the 2ND part of the box:
241         if(i.gt.1.and.(.not.ldone))then
242c            print *,'PART> 2nd',xx,yy
243            if(.not.isinttc(nbox,tmpvecp1,tmpvecp2,
244     1           px,py,pz,r1,r2,r3,dvec)) ttcbox=.false.
245         endif
246      ENDDO                     !DO I=1,NVEC(NBOX)
247
248      RETURN
249      END
250
251
252      LOGICAL FUNCTION ISINTTC(NBOX,NVC1,NVC2,P1,P2,P3,R1,R2,R3,dvec)
253      REAL*8 P1,P2,P3,R1,R2,R3
254      REAL*8 DVEC(3,3),MINTOL
255      REAL*8 A,B,C    !components of normal vector
256      REAL*8 D,KFC   ! Ax + By + Cz + D = 0
257      REAL*8 AVEC,BVEC,CVEC,AVEC2,BVEC2,CVEC2
258C     if nvec2<0 negative=.true.; third vector is turned upside down
259      LOGICAL NEGATIVE
260      REAL*8 TTCPOINT(3,7,8),TTCVEC(3,6,8)
261      INTEGER NTTCVEC(8),TTCPAIR(2,6,8)
262
263      COMMON/TTCBOXS/ TTCVEC,NTTCVEC,TTCPAIR,TTCPOINT
264
265      PARAMETER (MINTOL=1.0d-6)
266
267      NVEC1=NVC1
268      NVEC2=NVC2
269c      print *,'nvec1=',nvec1,', nvec2=',nvec2
270      NEGATIVE=.FALSE.
271      ISINTTC=.FALSE.
272      IF(NVEC2.LT.0.)THEN
273         NVEC2=-NVEC2
274         NEGATIVE=.TRUE.
275      ENDIF
276
277C     plvec is in coloumn major mode, but dvec is in row-major mode
278C     *** transform plvec1 in XYZ basis ***
279      AVEC=DVEC(1,1)*TTCVEC(1,NVEC1,NBOX)+
280     1     DVEC(2,1)*TTCVEC(2,NVEC1,NBOX)+DVEC(3,1)*TTCVEC(3,NVEC1,NBOX)
281      BVEC=DVEC(1,2)*TTCVEC(1,NVEC1,NBOX)+
282     1     DVEC(2,2)*TTCVEC(2,NVEC1,NBOX)+DVEC(3,2)*TTCVEC(3,NVEC1,NBOX)
283      CVEC=DVEC(1,3)*TTCVEC(1,NVEC1,NBOX)+
284     1     DVEC(2,3)*TTCVEC(2,NVEC1,NBOX)+DVEC(3,3)*TTCVEC(3,NVEC1,NBOX)
285      A=BVEC*DVEC(3,3)-DVEC(3,2)*CVEC
286      B=CVEC*DVEC(3,1)-DVEC(3,3)*AVEC
287      C=AVEC*DVEC(3,2)-DVEC(3,1)*BVEC
288
289C     *** transform ttcvec2 in XYZ basis ***
290      AVEC2=DVEC(1,1)*TTCVEC(1,NVEC2,NBOX)+
291     1     DVEC(2,1)*TTCVEC(2,NVEC2,NBOX)+DVEC(3,1)*TTCVEC(3,NVEC2,NBOX)
292      BVEC2=DVEC(1,2)*TTCVEC(1,NVEC2,NBOX)+
293     1     DVEC(2,2)*TTCVEC(2,NVEC2,NBOX)+DVEC(3,2)*TTCVEC(3,NVEC2,NBOX)
294      CVEC2=DVEC(1,3)*TTCVEC(1,NVEC2,NBOX)+
295     1     DVEC(2,3)*TTCVEC(2,NVEC2,NBOX)+DVEC(3,3)*TTCVEC(3,NVEC2,NBOX)
296      IF(A*AVEC2+B*BVEC2+C*CVEC2.LT.0.0d0)THEN
297         A=-A
298         B=-B
299         C=-C
300      ENDIF
301
302c      print *,'TTCVEC:',TTCVEC(1,NVEC1,NBOX),TTCVEC(2,NVEC1,NBOX),
303c     1     TTCVEC(3,NVEC1,NBOX)*DVEC(3,1)
304c      print *,'DVEC',DVEC(3,1),DVEC(3,2),DVEC(3,3)
305c      print *,'p1,p2,p3=',p1,p2,p3
306c      print *,'r1,r2,r3=',r1,r2,r3
307      D=-A*P1-B*P2-C*P3
308c      print *,'a,b,c,d::',a,b,c,d
309      KFC=A*R1+B*R2+C*R3+D
310c      print *,'kfc1=',kfc,', kfc2=',A*R1+B*R2+C*R3-(-A*P1-B*P2-C*P3)
311      IF(KFC.GE.(0.0d0-MINTOL).AND.(.NOT.NEGATIVE)) ISINTTC=.TRUE.
312      IF(KFC.LE.(0.0d0+MINTOL).AND.NEGATIVE) ISINTTC=.TRUE.
313c      print *,'ISINTTC=',isinttc
314
315      RETURN
316      END
317
318
319
320