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