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

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