Sources Esope | Cast3M

intpseg
C INTPSEG   SOURCE    SH236661  13/11/25    21:15:08     7869
C
      SUBROUTINE INTPSEG(NRA,NRB,NRC,NRD,X1,Y1,Z1,X2,Y2,Z2,
     &                   IINT,XI1,YI1,ZI1,XI2,YI2,ZI2,XZERO)
C----------------------------------------------------------------------C
C         Determination de l'intersection entre un segment et          C
C  un demi espace definit par un plan et un point dans le demi espace  C
C                                                                      C
C     NRA, NRB, NRC = entiers, numeros des noeuds du plan (ABC)        C
C     NRD = entier, numero du noeud dans le demi espace                C
C     X1,Y1,Z1                                                         C
C     X2,Y2,Z2 = coordonnees des points P1 et P2 du segment [P1,P2]    C
C                l'ordre P1,P2 est important (voir ci dessous)         C
C                                                                      C
C     IINT     = indicateur sur l'intersection :                       C
C                - IINT = 0 P1P2 est hors du demi plan                 C
C                         --> les coordonnees de I1 et I2 sont nulles  C
C                - IINT = 1 P1P2 rentre dans le demi plan              C
C                         --> I1 est le point d'intersection           C
C                         --> I2 = P2                                  C
C                - IINT = 2 P1P2 est dans le demi plan                 C
C                         --> I1 = P1                                  C
C                         --> I2 = P2                                  C
C                - IINT = 3 P1P2 sort du demi plan, cad :              C
C                         --> I1 = P1                                  C
C                         --> I2 est le point d'intersection           C
C     XI1,YI1,ZI1                                                      C
C     XI2,YI2,ZI2 = coordonnees des point du segment d'intersection    C
C                   oriente comme le segment [P1,P2]                   C
C----------------------------------------------------------------------C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
c      PARAMETER (ZERO=1D-14)
      ZERO = XZERO
C
      SEGACT,MCOORD
      IDIMP1=IDIM+1
C     Extraction des coordonnees des points A, B et C
      XA=XCOOR((NRA-1)*IDIMP1+1)
      YA=XCOOR((NRA-1)*IDIMP1+2)
      ZA=XCOOR((NRA-1)*IDIMP1+3)
      XB=XCOOR((NRB-1)*IDIMP1+1)
      YB=XCOOR((NRB-1)*IDIMP1+2)
      ZB=XCOOR((NRB-1)*IDIMP1+3)
      XC=XCOOR((NRC-1)*IDIMP1+1)
      YC=XCOOR((NRC-1)*IDIMP1+2)
      ZC=XCOOR((NRC-1)*IDIMP1+3)
C     Vecteur normal au plan ABC : n=AB^AC
      XN=((YB-YA)*(ZC-ZA))-((ZB-ZA)*(YC-YA))
      YN=((ZB-ZA)*(XC-XA))-((XB-XA)*(ZC-ZA))
      ZN=((XB-XA)*(YC-YA))-((YB-YA)*(XC-XA))
C     Signe du produit scalaire n.AD
      XD=XCOOR((NRD-1)*IDIMP1+1)
      YD=XCOOR((NRD-1)*IDIMP1+2)
      ZD=XCOOR((NRD-1)*IDIMP1+3)
      PSD=(XN*(XD-XA))+(YN*(YD-YA))+(ZN*(ZD-ZA))
      IF (ABS(PSD).LT.ZERO) THEN
C        Demi plan mal definit
         CALL ERREUR(363)
      ENDIF
      SIGD=SIGN(1D0,PSD)
C     Signe du produit scalaire n.AP1
      PS1=(XN*(X1-XA))+(YN*(Y1-YA))+(ZN*(Z1-ZA))
      SIG1=SIGN(1D0,PS1)
      IF (ABS(PS1).LT.ZERO) SIG1=SIGD
C     Signe du produit scalaire n.AP2
      PS2=(XN*(X2-XA))+(YN*(Y2-YA))+(ZN*(Z2-ZA))
      SIG2=SIGN(1D0,PS2)
      IF (ABS(PS2).LT.ZERO) SIG2=SIGD
C     Il y a 4 cas possibles
      IINT=-1
      XI1=0D0
      YI1=0D0
      ZI1=0D0
      XI2=0D0
      YI2=0D0
      ZI2=0D0
C     1) [P1,P2] est hors du demi plan
      IF ((SIG1.NE.SIGD).AND.(SIG2.NE.SIGD)) THEN
         IINT=0
         GOTO 999
      ENDIF
C     2) [P1,P2] rentre dans le demi plan
      IF ((SIG1.NE.SIGD).AND.(SIG2.EQ.SIGD)) THEN
         IINT=1
         XI2=X2
         YI2=Y2
         ZI2=Z2
         GOTO 100
      ENDIF
C     3) [P1,P2] est entierment dans le demi plan
      IF ((SIG1.EQ.SIGD).AND.(SIG2.EQ.SIGD)) THEN
         IINT=2
         XI1=X1
         YI1=Y1
         ZI1=Z1
         XI2=X2
         YI2=Y2
         ZI2=Z2
         GOTO 999
      ENDIF
C     4) [P1,P2] sort du demi plan
      IF ((SIG1.EQ.SIGD).AND.(SIG2.NE.SIGD)) THEN
         IINT=3
         XI1=X1
         YI1=Y1
         ZI1=Z1
         GOTO 100
      ENDIF
C
100   CONTINUE
C     Calcul de l'intersection de la droite (P1,P2) avec le plan (ABC)
C     Le plan (ABC) est definit par l'equation catesienne
C     XN*x + YN*y + ZN*z + W = 0
C     calcul de W, sachant que le point A appartient au plan
      W=-1D0*((XN*XA)+(YN*YA)+(ZN*ZA))
C     La droite (P1,P2) est definie par le systeme d'equations de
C     parametre t suivant
C     x = F1*t + F10
C     y = F2*t + F20    (F1,F2,F3) est le vecteur directeur de la droite
C     z = F3*t + F30
C     calcul d'un vecteur directeur de (P1,P2)
      F1=X2-X1
      F2=Y2-Y1
      F3=Z2-Z1
C     on peut prendre les coordonnees du point P1 pour les termes
C     F10,F20 et F30 puisque ce point est sur la droite
      F10=X1
      F20=Y1
      F30=Z1
C     Le point d'intersection de la droite et du plan est le point J
C     dont les coordonnees verifient les 4 equations precedentes
C     Calcul de la valeur de t (parametre de la droite) pour
C     correspondre avec le plan :
C     XN*(F1*t+F10) + YN*(F2*t+F20) + ZN*(F3*t+F30) + W = 0
      XDENOM=(XN*F1)+(YN*F2)+(ZN*F3)
      XNOM=-1D0*((XN*F10)+(YN*F20)+(ZN*F30)+W)
      TJ=XNOM/XDENOM
      XJ=(F1*TJ)+F10
      YJ=(F2*TJ)+F20
      ZJ=(F3*TJ)+F30
C     Ce point correspond a I1 ou I2 selon la valeur de IINT
      IF (IINT.EQ.1) THEN
         XI1=XJ
         YI1=YJ
         ZI1=ZJ
      ELSEIF (IINT.EQ.3) THEN
         XI2=XJ
         YI2=YJ
         ZI2=ZJ
      ENDIF
C     Fin du programme
999   RETURN
      END