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