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