Télécharger intpseg.eso

Retour à la liste

Numérotation des lignes :

  1. C INTPSEG SOURCE SH236661 13/11/25 21:15:08 7869
  2. C
  3. SUBROUTINE INTPSEG(NRA,NRB,NRC,NRD,X1,Y1,Z1,X2,Y2,Z2,
  4. & IINT,XI1,YI1,ZI1,XI2,YI2,ZI2,XZERO)
  5. C----------------------------------------------------------------------C
  6. C Determination de l'intersection entre un segment et C
  7. C un demi espace definit par un plan et un point dans le demi espace C
  8. C C
  9. C NRA, NRB, NRC = entiers, numeros des noeuds du plan (ABC) C
  10. C NRD = entier, numero du noeud dans le demi espace C
  11. C X1,Y1,Z1 C
  12. C X2,Y2,Z2 = coordonnees des points P1 et P2 du segment [P1,P2] C
  13. C l'ordre P1,P2 est important (voir ci dessous) C
  14. C C
  15. C IINT = indicateur sur l'intersection : C
  16. C - IINT = 0 P1P2 est hors du demi plan C
  17. C --> les coordonnees de I1 et I2 sont nulles C
  18. C - IINT = 1 P1P2 rentre dans le demi plan C
  19. C --> I1 est le point d'intersection C
  20. C --> I2 = P2 C
  21. C - IINT = 2 P1P2 est dans le demi plan C
  22. C --> I1 = P1 C
  23. C --> I2 = P2 C
  24. C - IINT = 3 P1P2 sort du demi plan, cad : C
  25. C --> I1 = P1 C
  26. C --> I2 est le point d'intersection C
  27. C XI1,YI1,ZI1 C
  28. C XI2,YI2,ZI2 = coordonnees des point du segment d'intersection C
  29. C oriente comme le segment [P1,P2] C
  30. C----------------------------------------------------------------------C
  31. IMPLICIT INTEGER(I-N)
  32. IMPLICIT REAL*8 (A-H,O-Z)
  33. -INC CCOPTIO
  34. -INC SMCOORD
  35. c PARAMETER (ZERO=1D-14)
  36. ZERO = XZERO
  37. C
  38. SEGACT,MCOORD
  39. IDIMP1=IDIM+1
  40. C Extraction des coordonnees des points A, B et C
  41. XA=XCOOR((NRA-1)*IDIMP1+1)
  42. YA=XCOOR((NRA-1)*IDIMP1+2)
  43. ZA=XCOOR((NRA-1)*IDIMP1+3)
  44. XB=XCOOR((NRB-1)*IDIMP1+1)
  45. YB=XCOOR((NRB-1)*IDIMP1+2)
  46. ZB=XCOOR((NRB-1)*IDIMP1+3)
  47. XC=XCOOR((NRC-1)*IDIMP1+1)
  48. YC=XCOOR((NRC-1)*IDIMP1+2)
  49. ZC=XCOOR((NRC-1)*IDIMP1+3)
  50. C Vecteur normal au plan ABC : n=AB^AC
  51. XN=((YB-YA)*(ZC-ZA))-((ZB-ZA)*(YC-YA))
  52. YN=((ZB-ZA)*(XC-XA))-((XB-XA)*(ZC-ZA))
  53. ZN=((XB-XA)*(YC-YA))-((YB-YA)*(XC-XA))
  54. C Signe du produit scalaire n.AD
  55. XD=XCOOR((NRD-1)*IDIMP1+1)
  56. YD=XCOOR((NRD-1)*IDIMP1+2)
  57. ZD=XCOOR((NRD-1)*IDIMP1+3)
  58. PSD=(XN*(XD-XA))+(YN*(YD-YA))+(ZN*(ZD-ZA))
  59. IF (ABS(PSD).LT.ZERO) THEN
  60. C Demi plan mal definit
  61. CALL ERREUR(363)
  62. ENDIF
  63. SIGD=SIGN(1D0,PSD)
  64. C Signe du produit scalaire n.AP1
  65. PS1=(XN*(X1-XA))+(YN*(Y1-YA))+(ZN*(Z1-ZA))
  66. SIG1=SIGN(1D0,PS1)
  67. IF (ABS(PS1).LT.ZERO) SIG1=SIGD
  68. C Signe du produit scalaire n.AP2
  69. PS2=(XN*(X2-XA))+(YN*(Y2-YA))+(ZN*(Z2-ZA))
  70. SIG2=SIGN(1D0,PS2)
  71. IF (ABS(PS2).LT.ZERO) SIG2=SIGD
  72. C Il y a 4 cas possibles
  73. IINT=-1
  74. XI1=0D0
  75. YI1=0D0
  76. ZI1=0D0
  77. XI2=0D0
  78. YI2=0D0
  79. ZI2=0D0
  80. C 1) [P1,P2] est hors du demi plan
  81. IF ((SIG1.NE.SIGD).AND.(SIG2.NE.SIGD)) THEN
  82. IINT=0
  83. GOTO 999
  84. ENDIF
  85. C 2) [P1,P2] rentre dans le demi plan
  86. IF ((SIG1.NE.SIGD).AND.(SIG2.EQ.SIGD)) THEN
  87. IINT=1
  88. XI2=X2
  89. YI2=Y2
  90. ZI2=Z2
  91. GOTO 100
  92. ENDIF
  93. C 3) [P1,P2] est entierment dans le demi plan
  94. IF ((SIG1.EQ.SIGD).AND.(SIG2.EQ.SIGD)) THEN
  95. IINT=2
  96. XI1=X1
  97. YI1=Y1
  98. ZI1=Z1
  99. XI2=X2
  100. YI2=Y2
  101. ZI2=Z2
  102. GOTO 999
  103. ENDIF
  104. C 4) [P1,P2] sort du demi plan
  105. IF ((SIG1.EQ.SIGD).AND.(SIG2.NE.SIGD)) THEN
  106. IINT=3
  107. XI1=X1
  108. YI1=Y1
  109. ZI1=Z1
  110. GOTO 100
  111. ENDIF
  112. C
  113. 100 CONTINUE
  114. C Calcul de l'intersection de la droite (P1,P2) avec le plan (ABC)
  115. C Le plan (ABC) est definit par l'equation catesienne
  116. C XN*x + YN*y + ZN*z + W = 0
  117. C calcul de W, sachant que le point A appartient au plan
  118. W=-1D0*((XN*XA)+(YN*YA)+(ZN*ZA))
  119. C La droite (P1,P2) est definie par le systeme d'equations de
  120. C parametre t suivant
  121. C x = F1*t + F10
  122. C y = F2*t + F20 (F1,F2,F3) est le vecteur directeur de la droite
  123. C z = F3*t + F30
  124. C calcul d'un vecteur directeur de (P1,P2)
  125. F1=X2-X1
  126. F2=Y2-Y1
  127. F3=Z2-Z1
  128. C on peut prendre les coordonnees du point P1 pour les termes
  129. C F10,F20 et F30 puisque ce point est sur la droite
  130. F10=X1
  131. F20=Y1
  132. F30=Z1
  133. C Le point d'intersection de la droite et du plan est le point J
  134. C dont les coordonnees verifient les 4 equations precedentes
  135. C Calcul de la valeur de t (parametre de la droite) pour
  136. C correspondre avec le plan :
  137. C XN*(F1*t+F10) + YN*(F2*t+F20) + ZN*(F3*t+F30) + W = 0
  138. XDENOM=(XN*F1)+(YN*F2)+(ZN*F3)
  139. XNOM=-1D0*((XN*F10)+(YN*F20)+(ZN*F30)+W)
  140. TJ=XNOM/XDENOM
  141. XJ=(F1*TJ)+F10
  142. YJ=(F2*TJ)+F20
  143. ZJ=(F3*TJ)+F30
  144. C Ce point correspond a I1 ou I2 selon la valeur de IINT
  145. IF (IINT.EQ.1) THEN
  146. XI1=XJ
  147. YI1=YJ
  148. ZI1=ZJ
  149. ELSEIF (IINT.EQ.3) THEN
  150. XI2=XJ
  151. YI2=YJ
  152. ZI2=ZJ
  153. ENDIF
  154. C Fin du programme
  155. 999 RETURN
  156. END
  157.  
  158.  
  159.  
  160.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales