C FACTR1 SOURCE CHAT 05/01/12 23:55:53 5004 SUBROUTINE FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2 $ ,XN1,XN2,XN3,XINT,ITEST) ************************************************************************ *** SP 'FACTR1' : A partir d'un plan et d'une droite definies par des *** pts, donne les normales (repere local lie au plan) et le pt *** d'intersection droite-plan. *** *** APPELES 1 = 'PROVEC', 'SCVECT' *** APPELES 2 = aucun *** *** E = 'NDIM' dimension de l'espace *** 'PT1', PT2', PT3' pts definissant le plan considere *** 'XDEP2', 'XARI2' pts definissant la droite consideree *** *** S = 'XN1' vecteur unitaire normal au plan *** 'XN2' vecteur unitaire normal à droite consideree et à 'XN1' *** 'XN3' vecteur unitaire norma à 'XN1' et 'XN2' *** 'XINT' coordonnees du pt d'intersection *** 'ITEST' vaut 1 si intersection, 0 si parallellisme droite-plan *** *** Rq : 'XZPREC' (-INC CCREEL), erreur precision calcul machine *** *** Auteur Cyril Nou ************************************************************************ IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC CCREEL DIMENSION XDEP2(3),XARI2(3) DIMENSION XINT(3),XN1(3),XN2(3),XN3(3) DIMENSION PT1(3),PT2(3),PT3(3) DIMENSION X1(3),X2(3),X3(3),X4(3) *** construction des vecteurs directeurs droite,plan, etc... DCARAC=0.D0 DO 10 I=1,NDIM X1(I)=XARI2(I)-XDEP2(I) X2(I)=XDEP2(I)-PT1(I) X3(I)=PT2(I)-PT1(I) X4(I)=PT3(I)-PT1(I) *** distance caracteristique pour adimensionner DCARAC=DCARAC+X1(I)**2 10 CONTINUE DCARAC=SQRT(DCARAC) *** 'ITEST'=1 si non parallellisme trajectoire-plan, 0 sinon ITEST=1 ************** *** CAS 2D *** ************** IF (NDIM.EQ.2) THEN *** 'XN1' normale au plan IF (ABS(PT1(2)-PT2(2)).GT.(XZPREC*DCARAC)) THEN XN1(1)=1.D0 XN1(2)=-(PT2(1)-PT1(1))/(PT2(2)-PT1(2)) XNORM=SQRT(XN1(1)**2+XN1(2)**2) XN1(1)=XN1(1)/XNORM XN1(2)=XN1(2)/XNORM ELSE XN1(1)=0.D0 XN1(2)=1.D0 ENDIF *** 'XN2' directeur du plan IF (ABS(PT1(1)-PT2(1)).GT.(XZPREC*DCARAC)) THEN XN2(1)=1.D0 XN2(2)=(PT2(2)-PT1(2))/(PT2(1)-PT1(1)) XNORM=SQRT(XN2(1)**2+XN2(2)**2) XN2(1)=XN2(1)/XNORM XN2(2)=XN2(2)/XNORM ELSE XN2(1)=0.D0 XN2(2)=1.D0 ENDIF *** 'XINT' intersection trajectoire/plan A=SCVECT(X2,XN2,NDIM) B=SCVECT(X2,XN1,NDIM) C=SCVECT(X1,XN2,NDIM) D=SCVECT(X1,XN1,NDIM) *** verification non parallelisme trajectoire-plan IF (ABS(D).GT.(XZPREC*DCARAC)) THEN COEFF=A-B*C/D XINT(1)=PT1(1)+COEFF*XN2(1) XINT(2)=PT1(2)+COEFF*XN2(2) ELSE ITEST=0 ENDIF ************** *** CAS 3D *** ************** ELSEIF (NDIM.EQ.3) THEN *** 'XN1' normale au plan CALL PROVEC(X3,X4,XN1) XNORM=SQRT(XN1(1)**2+XN1(2)**2+XN1(3)**2) XN1(1)=XN1(1)/XNORM XN1(2)=XN1(2)/XNORM XN1(3)=XN1(3)/XNORM *** 'XN2' normale à trajectoire et à 'XN1' CALL PROVEC(XN1,X1,XN2) XNORM=SQRT(XN2(1)**2+XN2(2)**2+XN2(3)**2) IF (ABS(XNORM).GT.(XZPREC*DCARAC)) THEN XN2(1)=XN2(1)/XNORM XN2(2)=XN2(2)/XNORM XN2(3)=XN2(3)/XNORM ELSE XNORM=SQRT(X3(1)**2+X3(2)**2+X3(3)**2) XN2(1)=X3(1)/XNORM XN2(2)=X3(2)/XNORM XN2(3)=X3(3)/XNORM ENDIF *** 'XN3' normale restante du repere local 'XN1','XN2',XN3' CALL PROVEC(XN1,XN2,XN3) *** 'XINT' intersection trajectoire/plan A=SCVECT(X2,XN2,NDIM) B=SCVECT(X2,XN3,NDIM) C=SCVECT(X2,XN1,NDIM) D=SCVECT(X1,XN3,NDIM) E=SCVECT(X1,XN1,NDIM) *** verification du non parallelisme trajectoire-plan IF (ABS(E).GT.(XZPREC*DCARAC)) THEN COEFF=B-C*D/E XINT(1)=PT1(1)+A*XN2(1)+COEFF*XN3(1) XINT(2)=PT1(2)+A*XN2(2)+COEFF*XN3(2) XINT(3)=PT1(3)+A*XN2(3)+COEFF*XN3(3) ELSE ITEST=0 ENDIF ENDIF RETURN END