1 : C INTERS SOURCE FANDEUR 13/01/16 21:15:22 7665 2 : SUBROUTINE INTERS 3 : C----------------------------------------------------------------------- 4 : C 5 : C INTERSECTIONS DE SURFACES (PLAN,SPHERE,CYLINDRE,CONE,TORE) 6 : C 7 : C----------------------------------------------------------------------- 8 : C 9 : C IOB1 EVENTUELLEMENT OBJET QUE L'ON PROLONGE 10 : C IOB2 EVENTUELLEMENT OBJET QUE L'ON AVANCE 11 : C IPOIN1 POINT INITIAL 12 : C IPOIN2 POINT FINAL 13 : C ICLE TYPE DE LA PREMIERE SURFACE 14 : C JCLE TYPE DE LA DEUXIEME SURFACE 15 : C XPAUX TABLEAU DE POINTS AUXILIAIRES SERVANT A RECTIFIER L'INTERSE 16 : C LA QUATRIEME COMPOSANTE EST L'ABCISSE CURVILIGNE 17 : C 18 : C----------------------------------------------------------------------- 19 : C 20 : C 17/09/93 : intersection de deux MELEMEs (sens ensembliste) 21 : C 22 : C PM 30/08/05 : Erreur si intersection vide, sauf si mot-clef 'NOVERIF' 23 : C 24 : C BP 28/08/2012 : Ajout de l'intersection geometrique de 2 maillages 25 : C 26 : C----------------------------------------------------------------------- 27 : c 28 : IMPLICIT INTEGER(I-N) 29 : IMPLICIT REAL*8 (A-H,O-Z) 30 : 31 : -INC CCOPTIO 32 : -INC CCGEOME 33 : 34 : -INC SMCOORD 35 : -INC SMELEME 36 : 37 : SEGMENT TABPAR(NPOIN) 38 : 39 : LOGICAL VERIF,ltelq 40 : CHARACTER*4 MCLE(8),MCLE2(1),MCLE3(1) 41 : C INITIALISATION A L'USAGE DE L'OPTIMISEUR 42 : PARAMETER ( NBT = 22 ) 43 : DIMENSION XPAUX(4,NBT) 44 : 45 : DATA RAYI/1./ 46 : DATA MCLE/'PLAN','SPHE','CYLI','CONI','TORI','DINI','DFIN','LONG'/ 47 : DATA MCLE2/'NOVE'/ 48 : DATA MCLE3/'GEOM'/ 49 : C 50 : C----------------------------------------------------------------------- 51 : C INTERSECTION DE DEUX MAILLAGES 52 : C 53 : C sens ensembliste ou geometrique? 53 : CALL LIRMOT(MCLE3,1,IGEO3,0) 54 : c write(ioimp,*) 'option GEOM ? ',IGEO3 55 : C 56 : c lecture des maillages 56 : CALL LIROBJ('MAILLAGE',IOB1,0,IRETOU) 57 : IF (IRETOU.EQ.1) THEN 57 : CALL LIROBJ('MAILLAGE',IOB2,0,IRETOV) 58 : c si pb, faut-il faire une erreur ou un maillage vide (NOVERIF)? 58 : CALL LIRMOT(MCLE2,1,ICLE,0) 59 : VERIF=(ICLE.EQ.0) 60 : 61 : C------ Intersection ensembliste ------------------ 62 : if(IGEO3.eq.0) then 63 : IF (IRETOV.EQ.1) THEN 64 : C Intersection de deux maillages 64 : CALL INTERB(IOB1,IOB2,IRETIB,IPT3) 65 : IF (IRETIB.EQ.1) THEN 66 : IF (VERIF) THEN 67 : C Impossible de calculer l'intersection 67 : CALL ERREUR(40) 68 : RETURN 69 : ELSE 70 : * On fournit un maillage vide 71 : NBNN=0 72 : NBELEM=0 73 : NBSOUS=0 74 : NBREF=0 75 : SEGINI,IPT3 76 : IPT3.ITYPEL=0 77 : SEGDES,IPT3 78 : ENDIF 79 : ENDIF 79 : CALL ECROBJ('MAILLAGE',IPT3) 80 : RETURN 81 : ELSE 81 : CALL REFUS 82 : ENDIF 83 : 84 : C------ Intersection geometrique ------------------ 85 : else 86 : IOB3=0 87 : IOB4=0 87 : CALL INTGEO(IOB1,IOB2,IOB3,IOB4,ICLE) 88 : C Cas ou n'a pas pu calculer l'intersection 89 : if(IOB3.eq.0) then 90 : IF (VERIF) THEN 91 : C Impossible de calculer l'intersection 91 : CALL ERREUR(40) 92 : RETURN 93 : ELSE 94 : * On fournit 2 maillage vide + les 2 initiaux 95 : NBNN=0 96 : NBELEM=0 97 : NBSOUS=0 98 : NBREF=0 99 : SEGINI,IPT3 100 : IPT3.ITYPEL=0 101 : IOB3=IPT3 102 : SEGDES,IPT3 103 : IF(IRETOV.EQ.1)THEN 104 : SEGINI,IPT3,IPT4 105 : IPT4.ITYPEL=0 106 : IOB4=IPT4 107 : SEGDES,IPT4 108 : ENDIF 109 : ENDIF 110 : endif 110 : IF(IRETOV.EQ.1) CALL ECROBJ('MAILLAGE',IOB4) 110 : CALL ECROBJ('MAILLAGE',IOB3) 110 : IF(IRETOV.EQ.1) CALL ECROBJ('MAILLAGE',IOB2) 110 : CALL ECROBJ('MAILLAGE',IOB1) 111 : RETURN 112 : endif 113 : 114 : ENDIF 115 : 116 : 117 : C----------------------------------------------------------------------- 118 : C CREATION D'UNE LIGNE RESULTAT D'UNE INTERSECTION DE SURFACES 119 : C DEFINIES ANALYTIQUEMENT (PLAN, SPHERE, ...) 120 : C 120 : IF (IDIM.NE.3) CALL ERREUR(14) 121 : IF (IERR.NE.0) RETURN 122 : 123 : C Y A T'IL UN FACTEUR DE DECOUPAGE 124 : INBR=0 124 : CALL LIRENT(INBR,0,IRETOU) 125 : IMPOI=0 126 : IMPOF=0 127 : IMPLON=0 128 : SEGACT MCOORD 129 : C LECTURE DU PREMIER POINT . SI C'EST UN ELEMENT ON EN EXTRAIT LE 130 : C DERNIER POINT . 131 : IOB1=0 131 : CALL LIROBJ('MAILLAGE',IOB1,0,IRETOU) 132 : IF (IRETOU.EQ.0) GOTO 1 133 : IPT6=IOB1 134 : SEGACT IPT6 134 : CALL EXTRPO(IOB1,1,IPOIN1) 135 : GOTO 2 135 : 1 CALL LIROBJ('POINT ',IPOIN1,1,IRETOU) 136 : 2 IF (IERR.NE.0) RETURN 137 : C PREMIERE SURFACE 137 : CALL LIRMOT(MCLE,8,ICLE,1) 138 : IF (IERR.NE.0) RETURN 139 : C Y A T-IL DES DENSITES IMPOSEES 140 : 82 IF (IMPOI.EQ.1.OR.ICLE.NE.6) GOTO 80 141 : IMPOI=1 141 : CALL LIRREE(TPOIN1,1,IRETOU) 141 : CALL LIRMOT(MCLE,8,ICLE,1) 142 : IF (IERR.NE.0) RETURN 143 : 80 IF (IMPOF.EQ.1.OR.ICLE.NE.7) GOTO 81 144 : IMPOF=1 144 : CALL LIRREE(TPOIN2,1,IRETOU) 144 : CALL LIRMOT(MCLE,8,ICLE,1) 145 : IF (IERR.NE.0) RETURN 146 : IF (IMPOI.EQ.0) GOTO 82 147 : 81 IF (IMPLON.EQ.1.OR.ICLE.NE.8) GOTO 800 148 : IMPLON=1 148 : CALL LIRREE(DLONG1,1,IRETOU) 148 : CALL LIRMOT(MCLE,8,ICLE,1) 149 : IF (IERR.NE.0) RETURN 150 : goto 82 151 : 800 continue 152 : C ACQUISITION DES DONNEES DE LA PREMIERE SURFACE 153 : GOTO (20,21,22,23,24),ICLE 153 : 20 CALL LIROBJ('POINT ',IPL,1,IRETOU) 154 : IF (IERR.NE.0) RETURN 155 : IREF=4*(IPL-1) 156 : XIPL=XCOOR(IREF+1) 157 : YIPL=XCOOR(IREF+2) 158 : ZIPL=XCOOR(IREF+3) 159 : GOTO 25 159 : 21 CALL LIROBJ('POINT ',IPL,1,IRETOU) 160 : IF (IERR.NE.0) RETURN 161 : IREF=4*(IPL-1) 162 : XIPL=XCOOR(IREF+1) 163 : YIPL=XCOOR(IREF+2) 164 : ZIPL=XCOOR(IREF+3) 165 : GOTO 25 165 : 22 CALL LIROBJ('POINT ',IPL,1,IRETOU) 166 : IF (IERR.NE.0) RETURN 167 : IREF=4*(IPL-1) 168 : XIPL=XCOOR(IREF+1) 169 : YIPL=XCOOR(IREF+2) 170 : ZIPL=XCOOR(IREF+3) 170 : CALL LIROBJ('POINT ',IPL2,1,IRETOU) 171 : IF (IERR.NE.0) RETURN 172 : IREF=4*(IPL2-1) 173 : XIV=XCOOR(IREF+1)-XIPL 174 : YIV=XCOOR(IREF+2)-YIPL 175 : ZIV=XCOOR(IREF+3)-ZIPL 176 : S=SQRT(XIV*XIV+YIV*YIV+ZIV*ZIV) 177 : IF (S.EQ.0.) THEN 177 : CALL ERREUR(21) 178 : RETURN 179 : ENDIF 180 : XIV=XIV/S 181 : YIV=YIV/S 182 : ZIV=ZIV/S 183 : GOTO 25 183 : 23 CALL LIROBJ('POINT ',IPL,1,IRETOU) 184 : IF (IERR.NE.0) RETURN 185 : IREF=4*(IPL-1) 186 : XIPL=XCOOR(IREF+1) 187 : YIPL=XCOOR(IREF+2) 188 : ZIPL=XCOOR(IREF+3) 188 : CALL LIROBJ('POINT ',IPL2,1,IRETOU) 189 : IF (IERR.NE.0) RETURN 190 : IREF=4*(IPL2-1) 191 : XIV=XCOOR(IREF+1)-XIPL 192 : YIV=XCOOR(IREF+2)-YIPL 193 : ZIV=XCOOR(IREF+3)-ZIPL 194 : S=SQRT(XIV**2+YIV**2+ZIV**2) 195 : IF (S.EQ.0.) THEN 195 : CALL ERREUR(21) 196 : RETURN 197 : ENDIF 198 : XIV=XIV/S 199 : YIV=YIV/S 200 : ZIV=ZIV/S 201 : GOTO 25 201 : 24 CALL LIROBJ('POINT ',IPL,1,IRETOU) 202 : IF (IERR.NE.0) RETURN 203 : IREF=4*(IPL-1) 204 : XIPL=XCOOR(IREF+1) 205 : YIPL=XCOOR(IREF+2) 206 : ZIPL=XCOOR(IREF+3) 206 : CALL LIROBJ('POINT ',IAXE,1,IRETOU) 207 : IF (IERR.NE.0) RETURN 208 : IREF=4*(IAXE-1) 209 : XIAXE=XCOOR(IREF+1)-XIPL 210 : YIAXE=XCOOR(IREF+2)-YIPL 211 : ZIAXE=XCOOR(IREF+3)-ZIPL 212 : SIAXE=SQRT(XIAXE**2+YIAXE**2+ZIAXE**2) 213 : IF (SIAXE.EQ.0.) THEN 213 : CALL ERREUR(21) 214 : RETURN 215 : ENDIF 216 : XIAXE=XIAXE/SIAXE 217 : YIAXE=YIAXE/SIAXE 218 : ZIAXE=ZIAXE/SIAXE 218 : CALL LIROBJ('POINT ',IP1,1,IRETOU) 219 : IF (IERR.NE.0) RETURN 220 : IREF=4*(IP1-1) 221 : XIP1=XCOOR(IREF+1) 222 : YIP1=XCOOR(IREF+2) 223 : ZIP1=XCOOR(IREF+3) 224 : GRAYI=SQRT((XIPL-XIP1)**2+(YIPL-YIP1)**2+(ZIPL-ZIP1)**2) 225 : IF (GRAYI.EQ.0.) THEN 225 : CALL ERREUR(21) 226 : RETURN 227 : ENDIF 228 : 25 CONTINUE 229 : C DEUXIEME SURFACE 229 : CALL LIRMOT(MCLE,8,JCLE,1) 230 : IF (IERR.NE.0) RETURN 231 : C Y A T-IL DES DENSITES IMPOSEES 232 : 85 IF (IMPOI.EQ.1.OR.JCLE.NE.6) GOTO 83 233 : IMPOI=1 233 : CALL LIRREE(TPOIN1,1,IRETOU) 233 : CALL LIRMOT(MCLE,8,JCLE,1) 234 : IF (IERR.NE.0) RETURN 235 : 83 IF (IMPOF.EQ.1.OR.JCLE.NE.7) GOTO 84 236 : IMPOF=1 236 : CALL LIRREE(TPOIN2,1,IRETOU) 236 : CALL LIRMOT(MCLE,8,JCLE,1) 237 : IF (IERR.NE.0) RETURN 238 : IF (IMPOI.EQ.0) GOTO 85 239 : 84 IF (IMPLON.EQ.1.OR.ICLE.NE.8) GOTO 801 240 : IMPLON=1 240 : CALL LIRREE(DLONG1,1,IRETOU) 240 : CALL LIRMOT(MCLE,8,ICLE,1) 241 : IF (IERR.NE.0) RETURN 242 : goto 85 243 : 801 continue 244 : C AQUISITION DES DONNEES DE LA DEUXIEME SURFACE 245 : GOTO (40,41,42,43,44),JCLE 245 : 40 CALL LIROBJ('POINT ',JPL,1,IRETOU) 246 : IF (IERR.NE.0) RETURN 247 : IREF=4*(JPL-1) 248 : XJPL=XCOOR(IREF+1) 249 : YJPL=XCOOR(IREF+2) 250 : ZJPL=XCOOR(IREF+3) 251 : GOTO 45 251 : 41 CALL LIROBJ('POINT ',JPL,1,IRETOU) 252 : IF (IERR.NE.0) RETURN 253 : IREF=4*(JPL-1) 254 : XJPL=XCOOR(IREF+1) 255 : YJPL=XCOOR(IREF+2) 256 : ZJPL=XCOOR(IREF+3) 257 : GOTO 45 257 : 42 CALL LIROBJ('POINT ',JPL,1,IRETOU) 258 : IF (IERR.NE.0) RETURN 259 : IREF=4*(JPL-1) 260 : XJPL=XCOOR(IREF+1) 261 : YJPL=XCOOR(IREF+2) 262 : ZJPL=XCOOR(IREF+3) 262 : CALL LIROBJ('POINT ',JPL2,1,IRETOU) 263 : IF (IERR.NE.0) RETURN 264 : IREF=4*(JPL2-1) 265 : XJV=XCOOR(IREF+1)-XJPL 266 : YJV=XCOOR(IREF+2)-YJPL 267 : ZJV=XCOOR(IREF+3)-ZJPL 268 : S=SQRT(XJV**2+YJV**2+ZJV**2) 269 : IF(S.EQ.0.) THEN 269 : CALL ERREUR(21) 270 : RETURN 271 : ENDIF 272 : XJV=XJV/S 273 : YJV=YJV/S 274 : ZJV=ZJV/S 275 : GOTO 45 275 : 43 CALL LIROBJ('POINT ',JPL,1,IRETOU) 276 : IF (IERR.NE.0) RETURN 277 : IREF=4*(JPL-1) 278 : XJPL=XCOOR(IREF+1) 279 : YJPL=XCOOR(IREF+2) 280 : ZJPL=XCOOR(IREF+3) 280 : CALL LIROBJ('POINT ',JPL2,1,IRETOU) 281 : IF (IERR.NE.0) RETURN 282 : IREF=4*(JPL2-1) 283 : XJV=XCOOR(IREF+1)-XJPL 284 : YJV=XCOOR(IREF+2)-YJPL 285 : ZJV=XCOOR(IREF+3)-ZJPL 286 : S=SQRT(XJV**2+YJV**2+ZJV**2) 287 : IF (S.EQ.0.) THEN 287 : CALL ERREUR(21) 288 : RETURN 289 : ENDIF 290 : XJV=XJV/S 291 : YJV=YJV/S 292 : ZJV=ZJV/S 293 : GOTO 45 293 : 44 CALL LIROBJ('POINT ',JPL,1,IRETOU) 294 : IF (IERR.NE.0) RETURN 295 : IREF=4*(JPL-1) 296 : XJPL=XCOOR(IREF+1) 297 : YJPL=XCOOR(IREF+2) 298 : ZJPL=XCOOR(IREF+3) 298 : CALL LIROBJ('POINT ',JAXE,1,IRETOU) 299 : IF (IERR.NE.0) RETURN 300 : IREF=4*(JAXE-1) 301 : XJAXE=XCOOR(IREF+1)-XJPL 302 : YJAXE=XCOOR(IREF+2)-YJPL 303 : ZJAXE=XCOOR(IREF+3)-ZJPL 304 : SJAXE=SQRT(XJAXE**2+YJAXE**2+ZJAXE**2) 305 : IF (SJAXE.EQ.0.) THEN 305 : CALL ERREUR(21) 306 : RETURN 307 : ENDIF 308 : XJAXE=XJAXE/SJAXE 309 : YJAXE=YJAXE/SJAXE 310 : ZJAXE=ZJAXE/SJAXE 310 : CALL LIROBJ('POINT ',JP1,1,IRETOU) 311 : IF (IERR.NE.0) RETURN 312 : IREF=4*(JP1-1) 313 : XJP1=XCOOR(IREF+1) 314 : YJP1=XCOOR(IREF+2) 315 : ZJP1=XCOOR(IREF+3) 316 : GRAYJ=SQRT((XJPL-XJP1)**2+(YJPL-YJP1)**2+(ZJPL-ZJP1)**2) 317 : IF (GRAYJ.EQ.0.) THEN 317 : CALL ERREUR(21) 318 : RETURN 319 : ENDIF 320 : 45 CONTINUE 321 : C Y A T-IL DES DENSITES IMPOSEES 322 : IRETOU=0 322 : 88 IF (IMPOI*IMPOF*implon.EQ.0) CALL LIRMOT(MCLE(6),3,MDPOS,0) 323 : IF (MDPOS.EQ.0) GOTO 89 324 : IF (IMPOI.EQ.1.OR.MDPOS.NE.1) GOTO 86 325 : IMPOI=1 325 : CALL LIRREE(TPOIN1,1,IRETOU) 326 : goto 88 327 : 86 IF (IMPOF.EQ.1.OR.MDPOS.NE.2) GOTO 87 328 : IMPOF=1 328 : CALL LIRREE(TPOIN2,1,IRETOU) 329 : GOTO 88 330 : 87 if (implon.EQ.1.OR.MDPOS.NE.3) GOTO 802 331 : implon=1 331 : CALL LIRREE(dlong1,1,IRETOU) 332 : goto 88 332 : 802 CALL REFUS 333 : 89 CONTINUE 334 : if (ierr.ne.0) return 335 : C LECTURE DU DEUXIEME POINT . SI C'EST UN ELEMENT ON EN EXTRAIT LE 336 : C PREMIER POINT . 337 : IOB2=0 338 : iretou=0 338 : if (implon.eq.0) CALL LIROBJ('MAILLAGE',IOB2,0,IRETOU) 339 : IF (IRETOU.EQ.0) GOTO 50 340 : IPT6=IOB2 341 : SEGACT IPT6 341 : CALL EXTRPO(IOB2,2,IPOIN2) 342 : GOTO 51 342 : 50 CALL LIROBJ('POINT ',IPOIN2,1,IRETOU) 343 : 51 IF (IERR.NE.0) RETURN 344 : NBP=NBT-2 345 : IREF1=(IPOIN1-1)*4 346 : IREF2=(IPOIN2-1)*4 347 : XPOIN1=XCOOR(IREF1+1) 348 : YPOIN1=XCOOR(IREF1+2) 349 : ZPOIN1=XCOOR(IREF1+3) 350 : IF (IMPOI.EQ.0) TPOIN1=XCOOR(IREF1+4) 351 : XPOIN2=XCOOR(IREF2+1) 352 : YPOIN2=XCOOR(IREF2+2) 353 : ZPOIN2=XCOOR(IREF2+3) 354 : IF (IMPOF.EQ.0) TPOIN2=XCOOR(IREF2+4) 355 : XVEC=XPOIN2-XPOIN1 356 : YVEC=YPOIN2-YPOIN1 357 : ZVEC=ZPOIN2-ZPOIN1 358 : SVEC=SQRT(XVEC**2+YVEC**2+ZVEC**2) 359 : IF (SVEC.EQ.0.) THEN 359 : CALL ERREUR(21) 360 : RETURN 361 : ENDIF 362 : EPS2=1E-10*SVEC**2 363 : C IL FAUT VERIFIER LA COMPATIBILITE DES DONNEES 364 : GOTO (60,61,62,63,64),ICLE 365 : 60 CONTINUE 366 : XPV1=XPOIN1-XIPL 367 : YPV1=YPOIN1-YIPL 368 : ZPV1=ZPOIN1-ZIPL 369 : XPV2=XPOIN2-XIPL 370 : YPV2=YPOIN2-YIPL 371 : ZPV2=ZPOIN2-ZIPL 372 : XIPN=YPV1*ZPV2-ZPV1*YPV2 373 : YIPN=ZPV1*XPV2-XPV1*ZPV2 374 : ZIPN=XPV1*YPV2-YPV1*XPV2 375 : SIPN=SQRT(XIPN**2+YIPN**2+ZIPN**2) 375 : IF (SIPN.EQ.0.) CALL ERREUR(21) 376 : IF (IERR.NE.0) RETURN 377 : XIPN=XIPN/SIPN 378 : YIPN=YIPN/SIPN 379 : ZIPN=ZIPN/SIPN 380 : GOTO 65 381 : 61 CONTINUE 382 : RAYI=SQRT((XPOIN1-XIPL)**2+(YPOIN1-YIPL)**2+(ZPOIN1-ZIPL)**2) 383 : RAYB=SQRT((XPOIN2-XIPL)**2+(YPOIN2-YIPL)**2+(ZPOIN2-ZIPL)**2) 384 : if (implon.eq.0) then 385 : IF (ABS((RAYI-RAYB)/RAYI).GE.1E-2) THEN 385 : CALL ERREUR(21) 386 : RETURN 387 : ENDIF 388 : RAYI=SQRT(RAYI*RAYB) 389 : endif 390 : GOTO 65 391 : 62 CONTINUE 392 : U=XPOIN1-XIPL 393 : V=YPOIN1-YIPL 394 : W=ZPOIN1-ZIPL 395 : SCA=U*XIV+V*YIV+W*ZIV 396 : U=U-SCA*XIV 397 : V=V-SCA*YIV 398 : W=W-SCA*ZIV 399 : RAYI=SQRT(U**2+V**2+W**2) 400 : U=XPOIN2-XIPL 401 : V=YPOIN2-YIPL 402 : W=ZPOIN2-ZIPL 403 : SCA=U*XIV+V*YIV+W*ZIV 404 : U=U-SCA*XIV 405 : V=V-SCA*YIV 406 : W=W-SCA*ZIV 407 : RAYB=SQRT(U**2+V**2+W**2) 408 : if (implon.eq.0) then 408 : IF (ABS(RAYI-RAYB)/RAYI.GT.1E-2) CALL ERREUR(21) 409 : RAYI=SQRT(RAYI*RAYB) 410 : endif 411 : GOTO 65 412 : 63 CONTINUE 413 : U=XPOIN1-XIPL 414 : V=YPOIN1-YIPL 415 : W=ZPOIN1-ZIPL 416 : SCA=ABS(U*XIV+V*YIV+W*ZIV) 417 : U1=V*ZIV-W*YIV 418 : V1=W*XIV-U*ZIV 419 : W1=U*YIV-V*XIV 420 : VEC=SQRT(U1**2+V1**2+W1**2) 420 : IF (SCA.EQ.0.) CALL ERREUR(21) 421 : IF (IERR.NE.0) RETURN 422 : TANGI=VEC/SCA 423 : U=XPOIN2-XIPL 424 : V=YPOIN2-YIPL 425 : W=ZPOIN2-ZIPL 426 : SCA=ABS(U*XIV+V*YIV+W*ZIV) 427 : U1=V*ZIV-W*YIV 428 : V1=W*XIV-U*ZIV 429 : W1=U*YIV-V*XIV 430 : VEC=SQRT(U1**2+V1**2+W1**2) 430 : IF (SCA.EQ.0.) CALL ERREUR(21) 431 : IF (IERR.NE.0) RETURN 432 : TANGB=VEC/SCA 433 : if (implon.eq.0) then 433 : IF (ABS(TANGI-TANGB)/TANGI.GE.1E-2) CALL ERREUR(21) 434 : TANGI=SQRT(TANGI*TANGB) 435 : endif 436 : GOTO 65 437 : 64 CONTINUE 438 : U=XPOIN1-XIPL 439 : V=YPOIN1-YIPL 440 : W=ZPOIN1-ZIPL 441 : SCA=U*XIAXE+V*YIAXE+W*ZIAXE 442 : U=U-SCA*XIAXE 443 : V=V-SCA*YIAXE 444 : W=W-SCA*ZIAXE 445 : S=SQRT(U**2+V**2+W**2) 445 : IF (S.EQ.0.) CALL ERREUR(21) 446 : IF (IERR.NE.0) RETURN 447 : XC=XIPL+U*GRAYI/S 448 : YC=YIPL+V*GRAYI/S 449 : ZC=ZIPL+W*GRAYI/S 450 : PRAYI=SQRT((XPOIN1-XC)**2+(YPOIN1-YC)**2+(ZPOIN1-ZC)**2) 451 : U=XPOIN2-XIPL 452 : V=YPOIN2-YIPL 453 : W=ZPOIN2-ZIPL 454 : SCA=U*XIAXE+V*YIAXE+W*ZIAXE 455 : U=U-SCA*XIAXE 456 : V=V-SCA*YIAXE 457 : W=W-SCA*ZIAXE 458 : S=SQRT(U**2+V**2+W**2) 458 : IF (S.EQ.0.) CALL ERREUR(21) 459 : IF (IERR.NE.0) RETURN 460 : XC=XIPL+U*GRAYI/S 461 : YC=YIPL+V*GRAYI/S 462 : ZC=ZIPL+W*GRAYI/S 463 : PRAYB=SQRT((XPOIN2-XC)**2+(YPOIN2-YC)**2+(ZPOIN2-ZC)**2) 464 : if (implon.eq.0) then 464 : IF (ABS(PRAYI-PRAYB)/PRAYI.GE.1E-2) CALL ERREUR(21) 465 : PRAYI=SQRT(PRAYI*PRAYB) 466 : endif 467 : 65 IF (IERR.NE.0) RETURN 468 : GOTO (70,71,72,73,74),JCLE 469 : 70 CONTINUE 470 : XPV1=XPOIN1-XJPL 471 : YPV1=YPOIN1-YJPL 472 : ZPV1=ZPOIN1-ZJPL 473 : XPV2=XPOIN2-XJPL 474 : YPV2=YPOIN2-YJPL 475 : ZPV2=ZPOIN2-ZJPL 476 : XJPN=YPV1*ZPV2-ZPV1*YPV2 477 : YJPN=ZPV1*XPV2-XPV1*ZPV2 478 : ZJPN=XPV1*YPV2-YPV1*XPV2 479 : SJPN=SQRT(XJPN**2+YJPN**2+ZJPN**2) 479 : IF (SJPN.EQ.0.) CALL ERREUR(21) 480 : IF (IERR.NE.0) RETURN 481 : XJPN=XJPN/SJPN 482 : YJPN=YJPN/SJPN 483 : ZJPN=ZJPN/SJPN 484 : GOTO 75 485 : 71 CONTINUE 486 : RAYJ=SQRT((XPOIN1-XJPL)**2+(YPOIN1-YJPL)**2+(ZPOIN1-ZJPL)**2) 487 : RAYB=SQRT((XPOIN2-XJPL)**2+(YPOIN2-YJPL)**2+(ZPOIN2-ZJPL)**2) 488 : if (implon.eq.0) then 488 : IF (ABS((RAYJ-RAYB)/RAYJ).GE.1E-2) CALL ERREUR(21) 489 : RAYJ=SQRT(RAYJ*RAYB) 490 : endif 491 : GOTO 75 492 : 72 CONTINUE 493 : U=XPOIN1-XJPL 494 : V=YPOIN1-YJPL 495 : W=ZPOIN1-ZJPL 496 : SCA=U*XJV+V*YJV+W*ZJV 497 : U=U-SCA*XJV 498 : V=V-SCA*YJV 499 : W=W-SCA*ZJV 500 : RAYJ=SQRT(U**2+V**2+W**2) 501 : U=XPOIN2-XJPL 502 : V=YPOIN2-YJPL 503 : W=ZPOIN2-ZJPL 504 : SCA=U*XJV+V*YJV+W*ZJV 505 : U=U-SCA*XJV 506 : V=V-SCA*YJV 507 : W=W-SCA*ZJV 508 : RAYB=SQRT(U**2+V**2+W**2) 509 : if (implon.eq.0) then 509 : IF (ABS(RAYJ-RAYB)/RAYJ.GT.1E-2) CALL ERREUR(21) 510 : RAYJ=SQRT(RAYJ*RAYB) 511 : endif 512 : GOTO 75 513 : 73 CONTINUE 514 : U=XPOIN1-XJPL 515 : V=YPOIN1-YJPL 516 : W=ZPOIN1-ZJPL 517 : SCA=ABS(U*XJV+V*YJV+W*ZJV) 518 : U1=V*ZJV-W*YJV 519 : V1=W*XJV-U*ZJV 520 : W1=U*YJV-V*XJV 521 : VEC=SQRT(U1**2+V1**2+W1**2) 521 : IF (SCA.EQ.0.) CALL ERREUR(21) 522 : IF (IERR.NE.0) RETURN 523 : TANGJ=VEC/SCA 524 : U=XPOIN2-XJPL 525 : V=YPOIN2-YJPL 526 : W=ZPOIN2-ZJPL 527 : SCA=ABS(U*XJV+V*YJV+W*ZJV) 528 : U1=V*ZJV-W*YJV 529 : V1=W*XJV-U*ZJV 530 : W1=U*YJV-V*XJV 531 : VEC=SQRT(U1**2+V1**2+W1**2) 531 : IF (SCA.EQ.0.) CALL ERREUR(21) 532 : IF (IERR.NE.0) RETURN 533 : TANGB=VEC/SCA 534 : if (implon.eq.0) then 534 : IF (ABS(TANGJ-TANGB)/TANGJ.GE.1E-2) CALL ERREUR(21) 535 : TANGJ=SQRT(TANGJ*TANGB) 536 : endif 537 : GOTO 75 538 : 74 CONTINUE 539 : U=XPOIN1-XJPL 540 : V=YPOIN1-YJPL 541 : W=ZPOIN1-ZJPL 542 : SCA=U*XJAXE+V*YJAXE+W*ZJAXE 543 : U=U-SCA*XJAXE 544 : V=V-SCA*YJAXE 545 : W=W-SCA*ZJAXE 546 : S=SQRT(U**2+V**2+W**2) 546 : IF (S.EQ.0.) CALL ERREUR(21) 547 : IF (IERR.NE.0) RETURN 548 : XC=XJPL+U*GRAYJ/S 549 : YC=YJPL+V*GRAYJ/S 550 : ZC=ZJPL+W*GRAYJ/S 551 : PRAYJ=SQRT((XPOIN1-XC)**2+(YPOIN1-YC)**2+(ZPOIN1-ZC)**2) 552 : U=XPOIN2-XJPL 553 : V=YPOIN2-YJPL 554 : W=ZPOIN2-ZJPL 555 : SCA=U*XJAXE+V*YJAXE+W*ZJAXE 556 : U=U-SCA*XJAXE 557 : V=V-SCA*YJAXE 558 : W=W-SCA*ZJAXE 559 : S=SQRT(U**2+V**2+W**2) 559 : IF (S.EQ.0.) CALL ERREUR(21) 560 : IF (IERR.NE.0) RETURN 561 : XC=XJPL+U*GRAYJ/S 562 : YC=YJPL+V*GRAYJ/S 563 : ZC=ZJPL+W*GRAYJ/S 564 : PRAYB=SQRT((XPOIN2-XC)**2+(YPOIN2-YC)**2+(ZPOIN2-ZC)**2) 565 : if (implon.eq.0) then 565 : IF (ABS(PRAYJ-PRAYB)/PRAYJ.GE.1E-2) CALL ERREUR(21) 566 : PRAYJ=SQRT(PRAYJ*PRAYB) 567 : endif 568 : 75 IF (IERR.NE.0) RETURN 569 : C ON PLACE NBP POINTS SUR L'INTERSECTION POUR CALCULER DES ABCISSES 570 : C CURVILIGNES 571 : XVECR=XVEC/(NBP+1) 572 : YVECR=YVEC/(NBP+1) 573 : ZVECR=ZVEC/(NBP+1) 574 : XVEC=XVEC/SVEC 575 : YVEC=YVEC/SVEC 576 : ZVEC=ZVEC/SVEC 577 : XPAUX(1,1)=XPOIN1 578 : XPAUX(2,1)=YPOIN1 579 : XPAUX(3,1)=ZPOIN1 580 : XPAUX(4,1)=0. 581 : IBP=0 582 : 102 IBP=IBP+1 583 : IF (IBP.GT.NBP) GOTO 100 584 : XP=XPAUX(1,IBP)+XVECR 585 : YP=XPAUX(2,IBP)+YVECR 586 : ZP=XPAUX(3,IBP)+ZVECR 587 : XPA=XP 588 : YPA=YP 589 : ZPA=ZP 590 : IRET = 101 591 : GOTO 400 592 : 101 CONTINUE 593 : IF (IIMPI.EQ.1) WRITE (IOIMP,8880) IBP,XPA,YPA,ZPA,XP1,YP1,ZP1 594 : 8880 FORMAT(I4,' POINT A PROJETER ',3G13.5,' POINT PROJETE ',3G13.5) 595 : C ON RANGE LE POINT D'INTERSECTION 596 : XPAUX(1,IBP+1)=XP1 597 : XPAUX(2,IBP+1)=YP1 598 : XPAUX(3,IBP+1)=ZP1 599 : XPAUX(4,IBP+1)=XPAUX(4,IBP)+SQRT((XP1-XPAUX(1,IBP))**2+ 600 : # (YP1-XPAUX(2,IBP))**2+(ZP1-XPAUX(3,IBP))**2) 601 : GOTO 102 602 : 100 CONTINUE 603 : XPAUX(1,NBT)=XPOIN2 604 : XPAUX(2,NBT)=YPOIN2 605 : XPAUX(3,NBT)=ZPOIN2 606 : XPAUX(4,NBT)=XPAUX(4,NBT-1)+SQRT((XP1-XPOIN2)**2+(YP1-YPOIN2)**2+ 607 : # (ZP1-ZPOIN2)**2) 608 : DLONG=XPAUX(4,NBT) 609 : if (implon.eq.1) then 610 : if ( iimpi.eq.1) write(ioimp,8888) dlong1 ,dlong 611 : 8888 format ( ' longueur demandee et longueur donnee',2g13.5) 611 : if (dlong1.gt.dlong) call erreur(21) 612 : if (ierr.ne.0) return 613 : dlong=dlong1 614 : endif 615 : DEN1=TPOIN1/DLONG 616 : DEN2=TPOIN2/DLONG 616 : CALL DECOUP(INBR,DEN1,DEN2,APROG,NBELEM,DENI,DECA,DLONG) 617 : NX=NBELEM-1 618 : DEN1=DEN1*DLONG 619 : IF (IIMPI.EQ.1) WRITE (IOIMP,9900) DLONG,NBELEM,APROG 620 : 9900 FORMAT (' LONGUEUR ',G12.5,' ELEMENT ',I5,' RAISON ',G12.5) 621 : C CREATION DU SEGMENT 622 : NBNN=NBNNE(KDEGRE(ILCOUR)) 622 : IF (ILCOUR.EQ.0) CALL ERREUR(16) 622 : IF (KDEGRE(ILCOUR).EQ.0) CALL ERREUR(16) 622 : IF (NBNN.EQ.0) CALL ERREUR(16) 623 : IF (IERR.NE.0) RETURN 624 : NBSOUS=0 625 : NBREF=0 626 : SEGINI MELEME 627 : ITYPEL=KDEGRE(ILCOUR) 628 : NPOIN=(ITYPEL-1)*NBELEM-1 629 : if (implon.eq.1) npoin=npoin+1 630 : SEGINI TABPAR 631 : C DEFINITION DE LA TOPOLOGIE 632 : IPOO=XCOOR(/1)/(IDIM+1) 633 : IADR=IPOO 634 : DIN=DEN1 635 : DL=0. 636 : NUM(1,1)=IPOIN1 637 : IF (NX.EQ.0) GOTO 510 638 : DO 512 I=1,NX 639 : DIN=DIN*APROG 640 : IF (ITYPEL.EQ.2) GOTO 513 641 : IPOO=IPOO+1 642 : NUM(2,I)=IPOO 643 : TABPAR(IPOO-IADR)=DIN/2.+DL 644 : 513 IPOO=IPOO+1 645 : NUM(ITYPEL,I)=IPOO 646 : NUM(1,I+1)=IPOO 647 : DL=DL+DIN 648 : TABPAR(IPOO-IADR)=DL 649 : 512 CONTINUE 650 : 510 CONTINUE 651 : IF (ITYPEL.NE.3) GOTO 514 652 : DIN=DIN*APROG 653 : IPOO=IPOO+1 654 : TABPAR(IPOO-IADR)=DL+DIN/2 655 : NUM(2,NBELEM)=IPOO 656 : 514 NUM(ITYPEL,NBELEM)=IPOIN2 657 : if (implon.eq.1) then 658 : NUM(ITYPEL,NBELEM)=ipoo+1 659 : TABPAR(IPOO+1-IADR)=DLong 660 : endif 661 : C CREATION DES POINTS LA LONGUEUR A RESPECTER EST DANS TABPAR 662 : C LA LONGUEUR DES POINTS CONNUS EST DANS XPAUX 663 : IP=1 664 : T=0. 665 : DLR=0. 666 : NBPTA=XCOOR(/1)/(IDIM+1) 667 : NBPTS=NBPTA+NPOIN 668 : SEGADJ MCOORD 669 : I=0 670 : 523 I=I+1 671 : IF (I.GT.NPOIN) GOTO 520 672 : DLV=TABPAR(I) 673 : 522 IF (DLV.LE.DLR) GOTO 521 674 : IP=IP+1 675 : DLR=XPAUX(4,IP) 676 : GOTO 522 677 : 521 CONTINUE 678 : DLA=XPAUX(4,IP-1) 679 : COF=(DLV-DLA)/(DLR-DLA) 680 : XP=XPAUX(1,IP-1)+COF*(XPAUX(1,IP)-XPAUX(1,IP-1)) 681 : YP=XPAUX(2,IP-1)+COF*(XPAUX(2,IP)-XPAUX(2,IP-1)) 682 : ZP=XPAUX(3,IP-1)+COF*(XPAUX(3,IP)-XPAUX(3,IP-1)) 683 : IF (IIMPI.NE.0) WRITE (IOIMP,8878) I,IP,XP,YP,ZP 684 : 8878 FORMAT(' AVANT ',I4,' INTERPOLE ',I4,3G13.5) 685 : IRET = 525 686 : GOTO 400 687 : 525 CONTINUE 688 : IF (IIMPI.NE.0) WRITE (IOIMP,8876) I,IP,XP1,YP1,ZP1 689 : 8876 FORMAT(' APRES ',I4,' INTERPOLE ',I4,3G13.5) 690 : XCOOR(NBPTA*4+1)=XP1 691 : XCOOR(NBPTA*4+2)=YP1 692 : XCOOR(NBPTA*4+3)=ZP1 693 : XCOOR(NBPTA*4+4)=TABPAR(I)-T 694 : NBPTA=NBPTA+1 695 : T=TABPAR(I) 696 : IADR=IADR+1 697 : GOTO 523 698 : 520 CONTINUE 699 : SEGSUP TABPAR 700 : IPT2=MELEME 701 : DO 527 I=1,NUM(/2) 702 : 527 ICOLOR(I)=IDCOUL 703 : IF (IOB1.NE.0) THEN 704 : ltelq=.false. 704 : CALL FUSE(IOB1,MELEME,IPT2,ltelq) 705 : IPT6=IOB1 706 : SEGDES IPT6 707 : ENDIF 708 : MELEME=IPT2 709 : IF (IOB2.NE.0) THEN 710 : ltelq=.false. 710 : CALL FUSE(MELEME,IOB2,IPT2,ltelq) 711 : IPT6=IOB2 712 : SEGDES IPT6 713 : ENDIF 714 : MELEME=IPT2 715 : SEGDES MELEME 715 : CALL ECROBJ('MAILLAGE',MELEME) 716 : RETURN 717 : C SOUS-PROGRAMME DE CALCUL ITERATIF D'INTERSECTION 718 : 400 CONTINUE 719 : ICYCL=0 720 : 402 CONTINUE 721 : C IL FAUT TROUVER L'INTERSECTION DU PLAN PERPENDICULAIRE A XVEC YVEC 722 : C ZVEC EN XPOIN YPOIN ZPOIN ET DES DEUX SURFACES 723 : C ON UTILISE UNE METHODE ITERATIVE 724 : C ON VA CHERCHER LE POINT DE PROJECTION DE XP SUR CHAQUE SURFACE DANS 725 : C LE PLAN ET UN VECTEUR NORMAL 726 : C XP YP ZP POINT A PROJETER 727 : C XP1 YP1 ZP1 INTERSECTION EN RETOUR 728 : GOTO (200,210,220,230,240),ICLE 729 : 250 CONTINUE 730 : GOTO (300,310,320,330,340),JCLE 731 : 350 CONTINUE 732 : C LES PTS ET LES VECTS SONT XPI... XPJ... XNI... XNJ... 733 : C ON CALCULE L'INTERSECTION DES DEUX VECTEURS (AFFINES) 734 : XDIF=XPJ-XPI 735 : YDIF=YPJ-YPI 736 : ZDIF=ZPJ-ZPI 737 : C COMPOSANTES SUR CHAQUE VECTEUR 738 : DET=(YNI*ZNJ-ZNI*YNJ)*XVEC 739 : # +(ZNI*XNJ-XNI*ZNJ)*YVEC 740 : # +(XNI*YNJ-YNI*XNJ)*ZVEC 740 : IF (DET**2.LT.EPS2) CALL ERREUR(40) 741 : A=((YDIF*ZNJ-ZDIF*YNJ)*XVEC 742 : # +(ZDIF*XNJ-XDIF*ZNJ)*YVEC 743 : # +(XDIF*YNJ-YDIF*XNJ)*ZVEC)/DET 744 : XP1=XPI+A*XNI 745 : YP1=YPI+A*YNI 746 : ZP1=ZPI+A*ZNI 747 : C TEST CONVERGENCE 748 : ECART=(XP-XP1)**2+(YP-YP1)**2+(ZP-ZP1)**2 749 : IF (ECART.LE.EPS2) THEN 750 : IF (IRET.EQ.101) GOTO 101 751 : IF (IRET.EQ.525) GOTO 525 752 : ENDIF 753 : XP=XP1 754 : YP=YP1 755 : ZP=ZP1 756 : ICYCL=ICYCL+1 757 : IF (IIMPI.EQ.1) WRITE (IOIMP,4478) ICYCL,ECART,EPS2,XPI,YPI,ZPI, 758 : # XNI,YNI,ZNI,XPJ,YPJ,ZPJ,XNJ,YNJ,ZNJ,XP,YP,ZP 759 : 4478 FORMAT(' ICYCL ',I4,' ECART ',G12.5,' EPS2 ',G12.5,/,' PI ', 760 : # 3G13.5,' NI ',3G13.5,/,' PJ ',3G13.5,' NJ ',3G13.5,' P ',3G13.5) 761 : IF (ICYCL.LE.25) GOTO 402 761 : CALL ERREUR(40) 762 : RETURN 763 : C PROJECTION SUR UN PLAN XIPL XIPN EN RESTANT DANS XVEC 764 : 200 CONTINUE 765 : C VECT PERPENDICULAIRE A XVEC ET XIPN 766 : U=YVEC*ZIPN-ZVEC*YIPN 767 : V=ZVEC*XIPN-XVEC*ZIPN 768 : W=XVEC*YIPN-YVEC*XIPN 769 : S=SQRT(U**2+V**2+W**2) 769 : IF (S.EQ.0.) CALL ERREUR(40) 770 : U1=U/S 771 : V1=V/S 772 : W1=W/S 773 : U=V1*ZVEC-W1*YVEC 774 : V=W1*XVEC-U1*ZVEC 775 : W=U1*YVEC-V1*XVEC 776 : C VECT SUIVANT LEQUEL ON DEPLACE LE POINT 777 : PAR=-((XP-XIPL)*XIPN+(YP-YIPL)*YIPN+(ZP-ZIPL)*ZIPN)/ 778 : # (U*XIPN+V*YIPN+W*ZIPN) 779 : XPI=XP+PAR*U 780 : YPI=YP+PAR*V 781 : ZPI=ZP+PAR*W 782 : XNI=U1 783 : YNI=V1 784 : ZNI=W1 785 : GOTO 250 786 : 300 CONTINUE 787 : U=YVEC*ZJPN-ZVEC*YJPN 788 : V=ZVEC*XJPN-XVEC*ZJPN 789 : W=XVEC*YJPN-YVEC*XJPN 790 : S=SQRT(U**2+V**2+W**2) 790 : IF (S.EQ.0.) CALL ERREUR(40) 791 : U1=U/S 792 : V1=V/S 793 : W1=W/S 794 : U=V1*ZVEC-W1*YVEC 795 : V=W1*XVEC-U1*ZVEC 796 : W=U1*YVEC-V1*XVEC 797 : PAR=-((XP-XJPL)*XJPN+(YP-YJPL)*YJPN+(ZP-ZJPL)*ZJPN)/ 798 : # (U*XJPN+V*YJPN+W*ZJPN) 799 : XPJ=XP+PAR*U 800 : YPJ=YP+PAR*V 801 : ZPJ=ZP+PAR*W 802 : XNJ=U1 803 : YNJ=V1 804 : ZNJ=W1 805 : GOTO 350 806 : C PROJECTION SUR SPHERE CENTRE XIPL RAYON RAYI EN RESTANT SUR LE PLAN X 807 : 210 CONTINUE 808 : U1=XIPL-XP 809 : V1=YIPL-YP 810 : W1=ZIPL-ZP 811 : SCA=U1*XVEC+V1*YVEC+W1*ZVEC 812 : XC=XIPL-SCA*XVEC 813 : YC=YIPL-SCA*YVEC 814 : ZC=ZIPL-SCA*ZVEC 814 : IF (SCA.GT.RAYI) CALL ERREUR(40) 815 : IF (IERR.NE.0) RETURN 816 : RAYC=SQRT(RAYI**2-SCA**2) 817 : U=XC-XP 818 : V=YC-YP 819 : W=ZC-ZP 820 : S=SQRT(U**2+V**2+W**2) 820 : IF (S.EQ.0.) CALL ERREUR(40) 821 : IF (IERR.NE.0) RETURN 822 : U=U/S 823 : V=V/S 824 : W=W/S 825 : XPI=XC-RAYC*U 826 : YPI=YC-RAYC*V 827 : ZPI=ZC-RAYC*W 828 : XNI=YVEC*W-ZVEC*V 829 : YNI=ZVEC*U-XVEC*W 830 : ZNI=XVEC*V-YVEC*U 831 : GOTO 250 832 : 310 CONTINUE 833 : U1=XJPL-XP 834 : V1=YJPL-YP 835 : W1=ZJPL-ZP 836 : SCA=U1*XVEC+V1*YVEC+W1*ZVEC 837 : XC=XJPL-SCA*XVEC 838 : YC=YJPL-SCA*YVEC 839 : ZC=ZJPL-SCA*ZVEC 839 : IF (SCA.GT.RAYJ) CALL ERREUR(40) 840 : IF (IERR.NE.0) RETURN 841 : RAYC=SQRT(RAYJ**2-SCA**2) 842 : U=XC-XP 843 : V=YC-YP 844 : W=ZC-ZP 845 : S=SQRT(U**2+V**2+W**2) 845 : IF (S.EQ.0.) CALL ERREUR(40) 846 : IF (IERR.NE.0) RETURN 847 : U=U/S 848 : V=V/S 849 : W=W/S 850 : XPJ=XC-RAYC*U 851 : YPJ=YC-RAYC*V 852 : ZPJ=ZC-RAYC*W 853 : XNJ=YVEC*W-ZVEC*V 854 : YNJ=ZVEC*U-XVEC*W 855 : ZNJ=XVEC*V-YVEC*U 856 : GOTO 350 857 : C PROJECTION SUR CYLINDRE AXE XIPL DIR XIV EN RESTANT DANS XVEC 858 : 220 CONTINUE 859 : U=XIPL-XP 860 : V=YIPL-YP 861 : W=ZIPL-ZP 862 : SCA=U*XIV+V*YIV+W*ZIV 863 : U=U-SCA*XIV 864 : V=V-SCA*YIV 865 : W=W-SCA*ZIV 866 : S=SQRT(U**2+V**2+W**2) 866 : IF (S.EQ.0.) CALL ERREUR(40) 867 : IF (IERR.NE.0) RETURN 868 : U1=U-U*RAYI/S 869 : V1=V-V*RAYI/S 870 : W1=W-W*RAYI/S 871 : DNUM=U1*XVEC+V1*YVEC+W1*ZVEC 872 : DNOM=XIV*XVEC+YIV*YVEC+ZIV*ZVEC 873 : IF (DNOM.EQ.0.) GOTO 221 874 : IF (IERR.NE.0) RETURN 875 : DPAR=-DNUM/DNOM 876 : XPI=XP+U1+DPAR*XIV 877 : YPI=YP+V1+DPAR*YIV 878 : ZPI=ZP+W1+DPAR*ZIV 879 : XNI=(V*ZVEC-W*YVEC)/S 880 : YNI=(W*XVEC-U*ZVEC)/S 881 : ZNI=(U*YVEC-V*XVEC)/S 882 : GOTO 250 883 : C L'AXE DU CYLINDRE EST PARALLELE AU PLAN XVEC 884 : 221 SCA=U*XVEC+V*YVEC+W*ZVEC 884 : IF (SCA.GT.RAYI) CALL ERREUR(40) 885 : IF (IERR.NE.0) RETURN 886 : RAYC=SQRT(RAYI**2-SCA**2) 887 : U=U-SCA*XVEC 888 : V=V-SCA*YVEC 889 : W=W-SCA*ZVEC 890 : S=SQRT(U**2+V**2+W**2) 890 : IF (S.EQ.0.) CALL ERREUR(40) 891 : IF (IERR.NE.0) RETURN 892 : XPI=XP+U-U*RAYC/S 893 : YPI=YP+V-V*RAYC/S 894 : ZPI=ZP+W-W*RAYC/S 895 : XNI=XIV 896 : YNI=YIV 897 : ZNI=ZIV 898 : GOTO 250 899 : 320 CONTINUE 900 : U=XJPL-XP 901 : V=YJPL-YP 902 : W=ZJPL-ZP 903 : SCA=U*XJV+V*YJV+W*ZJV 904 : U=U-SCA*XJV 905 : V=V-SCA*YJV 906 : W=W-SCA*ZJV 907 : S=SQRT(U**2+V**2+W**2) 907 : IF (S.EQ.0.) CALL ERREUR(40) 908 : IF (IERR.NE.0) RETURN 909 : U1=U-U*RAYJ/S 910 : V1=V-V*RAYJ/S 911 : W1=W-W*RAYJ/S 912 : DNUM=U1*XVEC+V1*YVEC+W1*ZVEC 913 : DNOM=XJV*XVEC+YJV*YVEC+ZJV*ZVEC 914 : IF (DNOM.EQ.0.) GOTO 321 915 : IF (IERR.NE.0) RETURN 916 : DPAR=-DNUM/DNOM 917 : XPJ=XP+U1+DPAR*XJV 918 : YPJ=YP+V1+DPAR*YJV 919 : ZPJ=ZP+W1+DPAR*ZJV 920 : XNJ=(V*ZVEC-W*YVEC)/S 921 : YNJ=(W*XVEC-U*ZVEC)/S 922 : ZNJ=(U*YVEC-V*XVEC)/S 923 : GOTO 350 924 : C L'AXE DU CYLINDRE EST PARALLELE AU PLAN XVEC 925 : 321 SCA=U*XVEC+V*YVEC+W*ZVEC 925 : IF (SCA.GT.RAYJ) CALL ERREUR(40) 926 : IF (IERR.NE.0) RETURN 927 : RAYC=SQRT(RAYJ**2-SCA**2) 928 : U=U-SCA*XVEC 929 : V=V-SCA*YVEC 930 : W=W-SCA*ZVEC 931 : S=SQRT(U**2+V**2+W**2) 931 : IF (S.EQ.0.) CALL ERREUR(40) 932 : IF (IERR.NE.0) RETURN 933 : XPJ=XP+U-U*RAYC/S 934 : YPJ=YP+V-V*RAYC/S 935 : ZPJ=ZP+W-W*RAYC/S 936 : XNJ=XJV 937 : YNJ=YJV 938 : ZNJ=ZJV 939 : GOTO 350 940 : C PROJECTION SUR CONE SOMMET XIPL DIR XIV TANG ANGLE SOMM TANGI EN 941 : C RESTANT DANS XVEC 942 : 230 CONTINUE 943 : U=XIPL-XP 944 : V=YIPL-YP 945 : W=ZIPL-ZP 946 : SCA=U*XIV+V*YIV+W*ZIV 947 : U=U-SCA*XIV 948 : V=V-SCA*YIV 949 : W=W-SCA*ZIV 950 : XC=XP+U 951 : YC=YP+U 952 : ZC=ZP+U 953 : RL2=(SCA*TANGI)**2 954 : SCA=U*XVEC+V*YVEC+W*ZVEC 955 : DL2=SCA**2 955 : IF (DL2.GT.RL2) CALL ERREUR(40) 956 : IF(IERR.NE.0) RETURN 957 : U=U-SCA*XVEC 958 : V=V-SCA*YVEC 959 : W=W-SCA*ZVEC 960 : S=SQRT(U**2+V**2+W**2) 960 : IF (S.EQ.0.) CALL ERREUR(40) 961 : IF (IERR.NE.0) RETURN 962 : U1=U/S 963 : V1=V/S 964 : W1=W/S 965 : VEC=SQRT(RL2-DL2) 966 : U=U-VEC*U1 967 : V=V-VEC*V1 968 : W=W-VEC*W1 969 : XPI=XP+U 970 : YPI=YP+V 971 : ZPI=ZP+W 972 : U=XIPL-XPI 973 : V=YIPL-YPI 974 : W=ZIPL-ZPI 975 : U1=V*ZIV-W*YIV 976 : V1=W*XIV-U*ZIV 977 : W1=U*YIV-V*XIV 978 : U2=V*W1-W*V1 979 : V2=W*U1-U*W1 980 : W2=U*V1-V*U1 981 : XNI=YVEC*W2-ZVEC*V2 982 : YNI=ZVEC*U2-XVEC*W2 983 : ZNI=XVEC*V2-YVEC*U2 984 : S=SQRT(XNI**2+YNI**2+ZNI**2) 984 : IF (S.EQ.0.) CALL ERREUR(40) 985 : IF (IERR.NE.0) RETURN 986 : XNI=XNI/S 987 : YNI=YNI/S 988 : ZNI=ZNI/S 989 : GOTO 250 990 : 330 CONTINUE 991 : U=XJPL-XP 992 : V=YJPL-YP 993 : W=ZJPL-ZP 994 : SCA=U*XJV+V*YJV+W*ZJV 995 : U=U-SCA*XJV 996 : V=V-SCA*YJV 997 : W=W-SCA*ZJV 998 : XC=XP+U 999 : YC=YP+U 1000 : ZC=ZP+U 1001 : RL2=(SCA*TANGJ)**2 1002 : SCA=U*XVEC+V*YVEC+W*ZVEC 1003 : DL2=SCA**2 1003 : IF (DL2.GT.RL2) CALL ERREUR(40) 1004 : IF(IERR.NE.0) RETURN 1005 : U=U-SCA*XVEC 1006 : V=V-SCA*YVEC 1007 : W=W-SCA*ZVEC 1008 : S=SQRT(U**2+V**2+W**2) 1008 : IF (S.EQ.0.) CALL ERREUR(40) 1009 : IF (IERR.NE.0) RETURN 1010 : U1=U/S 1011 : V1=V/S 1012 : W1=W/S 1013 : VEC=SQRT(RL2-DL2) 1014 : U=U-VEC*U1 1015 : V=V-VEC*V1 1016 : W=W-VEC*W1 1017 : XPJ=XP+U 1018 : YPJ=YP+V 1019 : ZPJ=ZP+W 1020 : U=XJPL-XPJ 1021 : V=YJPL-YPJ 1022 : W=ZJPL-ZPJ 1023 : U1=V*ZJV-W*YJV 1024 : V1=W*XJV-U*ZJV 1025 : W1=U*YJV-V*XJV 1026 : U2=V*W1-W*V1 1027 : V2=W*U1-U*W1 1028 : W2=U*V1-V*U1 1029 : XNJ=YVEC*W2-ZVEC*V2 1030 : YNJ=ZVEC*U2-XVEC*W2 1031 : ZNJ=XVEC*V2-YVEC*U2 1032 : S=SQRT(XNJ**2+YNJ**2+ZNJ**2) 1032 : IF (S.EQ.0.) CALL ERREUR(40) 1033 : IF (IERR.NE.0) RETURN 1034 : XNJ=XNJ/S 1035 : YNJ=YNJ/S 1036 : ZNJ=ZNJ/S 1037 : GOTO 350 1038 : C PROJECTION SUR UN TORE CENTRE XIPL AXE XIAXE CENTRE SECONDAIRE XIP1 1039 : C GRAND RAYON GRAYI PETIT RAYON PRAYI EN RESTANT DANS XVEC 1040 : 240 CONTINUE 1041 : U=XP-XIPL 1042 : V=YP-YIPL 1043 : W=ZP-ZIPL 1044 : SCA=U*XIAXE+V*YIAXE+W*ZIAXE 1045 : U=U-SCA*XIAXE 1046 : V=V-SCA*YIAXE 1047 : W=W-SCA*ZIAXE 1048 : S=SQRT(U**2+V**2+W**2) 1048 : IF (S.EQ.0.) CALL ERREUR(40) 1049 : IF (IERR.NE.0) RETURN 1050 : C CENTRE DU CERCLE INTERS 1051 : XC=XIPL+U*GRAYI/S 1052 : YC=YIPL+V*GRAYI/S 1053 : ZC=ZIPL+W*GRAYI/S 1054 : U=XP-XC 1055 : V=YP-YC 1056 : W=ZP-ZC 1057 : S=SQRT(U**2+V**2+W**2) 1057 : IF (S.EQ.0) CALL ERREUR(40) 1058 : IF (IERR.NE.0) RETURN 1059 : C PROJECTION SUR LE CERCLE 1060 : X1=XC+U*PRAYI/S 1061 : Y1=YC+V*PRAYI/S 1062 : Z1=ZC+W*PRAYI/S 1063 : C IL FAUT TOURNER SUR LE TORE JUSQU'AU PLAN XP XVEC 1064 : C CENTRE DE ROTATION 1065 : U=X1-XIPL 1066 : V=Y1-YIPL 1067 : W=Z1-ZIPL 1068 : SCA=U*XIAXE+V*YIAXE+W*ZIAXE 1069 : U=SCA*XIAXE 1070 : V=SCA*YIAXE 1071 : W=SCA*ZIAXE 1072 : XC1=XIPL+U 1073 : YC1=YIPL+V 1074 : ZC1=ZIPL+W 1075 : C RAYON DE CE CERCLE 1076 : RAY=SQRT((X1-XC1)**2+(Y1-YC1)**2+(Z1-ZC1)**2) 1077 : C DIRECTION DE L'INTERSECTION DU PLAN DU CERCLE ET DE XP XVEC 1078 : U=YIAXE*ZVEC-ZIAXE*YVEC 1079 : V=ZIAXE*XVEC-XIAXE*ZVEC 1080 : W=XIAXE*YVEC-YIAXE*XVEC 1081 : S=SQRT(U**2+V**2+W**2) 1081 : IF (S.EQ.0.) CALL ERREUR(40) 1082 : IF (IERR.NE.0) RETURN 1083 : U=U/S 1084 : V=V/S 1085 : W=W/S 1086 : U1=V*ZIAXE-W*YIAXE 1087 : V1=W*XIAXE-U*ZIAXE 1088 : W1=U*YIAXE-V*XIAXE 1089 : U2=XC1-XP 1090 : V2=YC1-YP 1091 : W2=ZC1-ZP 1092 : DNUM=U2*XVEC+V2*YVEC+W2*ZVEC 1093 : DNOM=U1*XVEC+V1*YVEC+W1*ZVEC 1093 : IF (DNOM.EQ.0.) CALL ERREUR(40) 1094 : IF (IERR.NE.0.) RETURN 1095 : PAR=-DNUM/DNOM 1096 : XC2=XC1+PAR*U1 1097 : YC2=YC1+PAR*V1 1098 : ZC2=ZC1+PAR*W1 1099 : S=SQRT((XC2-XC1)**2+(YC2-YC1)**2+(ZC2-ZC1)**2) 1099 : IF (S.GT.RAY) CALL ERREUR(40) 1100 : IF (IERR.NE.0) RETURN 1101 : D=SQRT(RAY**2-S**2) 1102 : IF ((X1-XC2)*U+(Y1-YC2)*V+(Z1-ZC2)*W.LE.0.) D=-D 1103 : XPI=XC2+D*U 1104 : YPI=YC2+D*V 1105 : ZPI=ZC2+D*W 1106 : C CENTRE DU PETIT CERCLE 1107 : U=XPI-XIPL 1108 : V=YPI-YIPL 1109 : W=ZPI-ZIPL 1110 : SCA=U*XIAXE+V*YIAXE+W*ZIAXE 1111 : U=U-SCA*XIAXE 1112 : V=V-SCA*YIAXE 1113 : W=W-SCA*ZIAXE 1114 : S=SQRT(U**2+V**2+W**2) 1115 : XC1=XIPL+U*GRAYI/S 1116 : YC1=YIPL+V*GRAYI/S 1117 : ZC1=ZIPL+W*GRAYI/S 1118 : U=XPI-XC1 1119 : V=YPI-YC1 1120 : W=ZPI-ZC1 1121 : XNI=V*ZVEC-W*YVEC 1122 : YNI=W*XVEC-U*ZVEC 1123 : ZNI=U*YVEC-V*XVEC 1124 : S=SQRT(XNI**2+YNI**2+ZNI**2) 1125 : IF (S.EQ.0.) GOTO 250 1126 : XNI=XNI/S 1127 : YNI=YNI/S 1128 : ZNI=ZNI/S 1129 : GOTO 250 1130 : 340 CONTINUE 1131 : U=XP-XJPL 1132 : V=YP-YJPL 1133 : W=ZP-ZJPL 1134 : SCA=U*XJAXE+V*YJAXE+W*ZJAXE 1135 : U=U-SCA*XJAXE 1136 : V=V-SCA*YJAXE 1137 : W=W-SCA*ZJAXE 1138 : S=SQRT(U**2+V**2+W**2) 1138 : IF (S.EQ.0.) CALL ERREUR(40) 1139 : IF (IERR.NE.0) RETURN 1140 : C CENTRE DU CERCLE INTERS 1141 : XC=XJPL+U*GRAYJ/S 1142 : YC=YJPL+V*GRAYJ/S 1143 : ZC=ZJPL+W*GRAYJ/S 1144 : U=XP-XC 1145 : V=YP-YC 1146 : W=ZP-ZC 1147 : S=SQRT(U**2+V**2+W**2) 1147 : IF (S.EQ.0) CALL ERREUR(40) 1148 : IF (IERR.NE.0) RETURN 1149 : C PROJECTION SUR LE CERCLE 1150 : X1=XC+U*PRAYJ/S 1151 : Y1=YC+V*PRAYJ/S 1152 : Z1=ZC+W*PRAYJ/S 1153 : C IL FAUT TOURNER SUR LE TORE JUSQU'AU PLAN XP XVEC 1154 : C CENTRE DE ROTATION 1155 : U=X1-XJPL 1156 : V=Y1-YJPL 1157 : W=Z1-ZJPL 1158 : SCA=U*XJAXE+V*YJAXE+W*ZJAXE 1159 : U=SCA*XJAXE 1160 : V=SCA*YJAXE 1161 : W=SCA*ZJAXE 1162 : XC1=XJPL+U 1163 : YC1=YJPL+V 1164 : ZC1=ZJPL+W 1165 : C RAYON DE CE CERCLE 1166 : RAY=SQRT((X1-XC1)**2+(Y1-YC1)**2+(Z1-ZC1)**2) 1167 : C DIRECTION DE L'INTERSECTION DU PLAN DU CERCLE ET DE XP XVEC 1168 : U=YJAXE*ZVEC-ZJAXE*YVEC 1169 : V=ZJAXE*XVEC-XJAXE*ZVEC 1170 : W=XJAXE*YVEC-YJAXE*XVEC 1171 : S=SQRT(U**2+V**2+W**2) 1171 : IF (S.EQ.0.) CALL ERREUR(40) 1172 : IF (IERR.NE.0) RETURN 1173 : U=U/S 1174 : V=V/S 1175 : W=W/S 1176 : U1=V*ZJAXE-W*YJAXE 1177 : V1=W*XJAXE-U*ZJAXE 1178 : W1=U*YJAXE-V*XJAXE 1179 : U2=XC1-XP 1180 : V2=YC1-YP 1181 : W2=ZC1-ZP 1182 : DNUM=U2*XVEC+V2*YVEC+W2*ZVEC 1183 : DNOM=U1*XVEC+V1*YVEC+W1*ZVEC 1183 : IF (DNOM.EQ.0.) CALL ERREUR(40) 1184 : IF (IERR.NE.0.) RETURN 1185 : PAR=-DNUM/DNOM 1186 : XC2=XC1+PAR*U1 1187 : YC2=YC1+PAR*V1 1188 : ZC2=ZC1+PAR*W1 1189 : S=SQRT((XC2-XC1)**2+(YC2-YC1)**2+(ZC2-ZC1)**2) 1189 : IF (S.GT.RAY) CALL ERREUR(40) 1190 : IF (IERR.NE.0) RETURN 1191 : D=SQRT(RAY**2-S**2) 1192 : IF ((X1-XC2)*U+(Y1-YC2)*V+(Z1-ZC2)*W.LE.0.) D=-D 1193 : XPJ=XC2+D*U 1194 : YPJ=YC2+D*V 1195 : ZPJ=ZC2+D*W 1196 : C CENTRE DU PETIT CERCLE 1197 : U=XPJ-XJPL 1198 : V=YPJ-YJPL 1199 : W=ZPJ-ZJPL 1200 : SCA=U*XJAXE+V*YJAXE+W*ZJAXE 1201 : U=U-SCA*XJAXE 1202 : V=V-SCA*YJAXE 1203 : W=W-SCA*ZJAXE 1204 : S=SQRT(U**2+V**2+W**2) 1205 : XC1=XJPL+U*GRAYJ/S 1206 : YC1=YJPL+V*GRAYJ/S 1207 : ZC1=ZJPL+W*GRAYJ/S 1208 : U=XPJ-XC1 1209 : V=YPJ-YC1 1210 : W=ZPJ-ZC1 1211 : XNJ=V*ZVEC-W*YVEC 1212 : YNJ=W*XVEC-U*ZVEC 1213 : ZNJ=U*YVEC-V*XVEC 1214 : S=SQRT(XNJ**2+YNJ**2+ZNJ**2) 1215 : IF (S.EQ.0.) GOTO 250 1216 : XNJ=XNJ/S 1217 : YNJ=YNJ/S 1218 : ZNJ=ZNJ/S 1219 : GOTO 350 1220 : END 1221 : 1222 : 1223 : 1224 : 1225 :
© Cast3M 2003 - All rights reserved.
Disclaimer