Télécharger factr2.eso

Retour à la liste

Numérotation des lignes :

  1. C FACTR2 SOURCE CHAT 05/01/12 23:55:58 5004
  2. SUBROUTINE FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4,XDEP2,XARI2
  3. $ ,XN1,XN2,XN3,XINT,ITEST)
  4. ***********************************************************************
  5. *** SP 'FACTR2' : a partir d'une face definie par des pts, permet de
  6. *** dire si un segment (pt depart-pt arrivee) rencontre cette face.
  7. ***
  8. *** APPELES 1 = 'TRPOSE', 'MATMAT', 'SCVECT' (fonction)
  9. *** APPELES 2 = aucun
  10. ***
  11. *** E = 'NDIM' dimension de l'espace
  12. *** 'ITYG' entier caracterisant la geometrie de l'elemnt considere
  13. *** 'PT1', 'PT2', 'PT3', 'PT4' noeuds de la face
  14. *** 'XDEP2' coordonnees reelles de depart de la particule
  15. *** 'XARI2' coordonnees reelles d'arrivee de la particule
  16. *** 'XN1' vecteur unitaire normal à la face
  17. *** 'XN2' vecteur unitaire normal à trajectoire et 'XN1'
  18. *** 'XN3' vecteur unitaire normal à 'XN1' et 'XN2'
  19. *** 'XINT' pt intersection plan associé à face/droite associé à trajectoire
  20. ***
  21. *** S = 'ITEST' vaut 1 si face effectivement traversee, 0 sinon
  22. ***
  23. *** Rq : 'XZPREC' (-INC CCREEL) erreur precision calcul machine
  24. ***
  25. *** ORIGINE = CYRIL NOU
  26. ***********************************************************************
  27. IMPLICIT INTEGER(I-N)
  28. IMPLICIT REAL*8 (A-H,O-Z)
  29. -INC CCREEL
  30. DIMENSION XDEP2(3),XARI2(3)
  31. DIMENSION XINT(3),XN1(3),XN2(3),XN3(3)
  32. DIMENSION PT1(3),PT2(3),PT3(3),PT4(3)
  33. DIMENSION X1(3),X2(3),X3(3),X4(4),XPASS(3,3)
  34. DIMENSION PPT1(3),PPT2(3),PPT3(3),PPT4(3),XXINT(3)
  35. DIMENSION PPPT1(3),PPPT2(3),PPPT3(3),PPPT4(3),XXXINT(3)
  36. *** construction des vecteurs directeurs droite,plan, etc...
  37. DO 10 I=1,NDIM
  38. X1(I)=XARI2(I)-XDEP2(I)
  39. X2(I)=XINT(I)-XDEP2(I)
  40. X3(I)=PT2(I)-PT1(I)
  41. X4(I)=XINT(I)-PT1(I)
  42. 10 CONTINUE
  43. ITEST=0
  44. *** 'ITEST1' teste appartenance à trajectoire, 'ITEST2' appartenance à face
  45. ITEST1=0
  46. ITEST2=0
  47.  
  48. **************************************************
  49. *** PT INTERSECTION APPARTIENT À TRAJECTOIRE ? ***
  50. **************************************************
  51.  
  52. *** calcul de la distance entre les pts départ et arrivee
  53. X1NORM=0.D0
  54. DO 20 I=1,NDIM
  55. X1NORM=X1NORM+X1(I)**2
  56. 20 CONTINUE
  57. *** verification si 'XINT' appartient au segment depart-arrivee
  58. SCAL=SCVECT(X2,X1,NDIM)
  59. IF ((SCAL.GT.(-XZPREC*X1NORM))
  60. $ .AND.(SCAL.LT.(X1NORM*(1+XZPREC)))) ITEST1=1
  61.  
  62. *******************************************
  63. *** PT INTERSECTION APPARTIENT À FACE ? ***
  64. *******************************************
  65.  
  66. ********** CAS 2D **********
  67. IF (NDIM.EQ.2) THEN
  68. *** calcul de la norme de la face
  69. X3NORM=0.D0
  70. DO 30 I=1,NDIM
  71. X3NORM=X3NORM+X3(I)**2
  72. 30 CONTINUE
  73. *** verification si 'XINT' appartient au segment face
  74. SCAL=SCVECT(X4,X3,NDIM)
  75. IF ((SCAL.GT.(-XZPREC*X3NORM))
  76. $ .AND.(SCAL.LT.(X3NORM*(1+XZPREC)))) ITEST2=1
  77.  
  78. ********** CAS 3D **********
  79. ELSEIF (NDIM.EQ.3) THEN
  80. *** matrice de passage du repere lié à la face vers repere absolu
  81. DO 40 I=1,NDIM
  82. XPASS(I,1)=XN1(I)
  83. XPASS(I,2)=XN2(I)
  84. XPASS(I,3)=XN3(I)
  85. 40 CONTINUE
  86. *** transpose de 'XPASS' (repere absolu vers repere face)
  87. CALL TRPOSE(XPASS)
  88. *** passage des coordonnees des pts du repere absolu au repere face
  89. CALL MATMAT(XPASS,PT1,NDIM,NDIM,1,PPT1)
  90. CALL MATMAT(XPASS,PT2,NDIM,NDIM,1,PPT2)
  91. CALL MATMAT(XPASS,PT3,NDIM,NDIM,1,PPT3)
  92. CALL MATMAT(XPASS,PT4,NDIM,NDIM,1,PPT4)
  93. CALL MATMAT(XPASS,XINT,NDIM,NDIM,1,XXINT)
  94. *** test par rapport à droite ppt1-ppt2, origine repere face=ppt3
  95. DO 50 I=1,NDIM
  96. PPPT1(I)=PPT1(I)-PPT3(I)
  97. PPPT2(I)=PPT2(I)-PPT3(I)
  98. XXXINT(I)=XXINT(I)-PPT3(I)
  99. 50 CONTINUE
  100. A=(XXXINT(2)-PPPT1(2))*(PPPT2(3)-PPPT1(3))
  101. B=(XXXINT(3)-PPPT1(3))*(PPPT2(2)-PPPT1(2))
  102. C=PPPT1(3)*(PPPT2(2)-PPPT1(2))
  103. D=PPPT1(2)*(PPPT2(3)-PPPT1(3))
  104. TEST=(A-B)/(C-D)
  105. IF (TEST.GT.(-XZPREC)) THEN
  106. ** test par rapport à droite ppt2-ppt3, origine repere face=ppt1
  107. DO 60 I=1,NDIM
  108. PPPT2(I)=PPT2(I)-PPT1(I)
  109. PPPT3(I)=PPT3(I)-PPT1(I)
  110. XXXINT(I)=XXINT(I)-PPT1(I)
  111. 60 CONTINUE
  112. A=(XXXINT(2)-PPPT2(2))*(PPPT3(3)-PPPT2(3))
  113. B=(XXXINT(3)-PPPT2(3))*(PPPT3(2)-PPPT2(2))
  114. C=PPPT2(3)*(PPPT3(2)-PPPT2(2))
  115. D=PPPT2(2)*(PPPT3(3)-PPPT2(3))
  116. TEST=(A-B)/(C-D)
  117. IF (TEST.GT.(-XZPREC)) THEN
  118. *** test par rapport à droite ppt3-ppt4, origine repere face=ppt2
  119. IF ((ABS(PPT3(2)-PPT4(2)).GT.(XZPREC*SQRT(X1NORM)))
  120. $ .OR.(ABS(PPT3(3)-PPT4(3)).GT.(XZPREC*SQRT(X1NORM)))) THEN
  121. DO 70 I=1,NDIM
  122. PPPT3(I)=PPT3(I)-PPT2(I)
  123. PPPT4(I)=PPT4(I)-PPT2(I)
  124. XXXINT(I)=XXINT(I)-PPT2(I)
  125. 70 CONTINUE
  126. A=(XXXINT(2)-PPPT3(2))*(PPPT4(3)-PPPT3(3))
  127. B=(XXXINT(3)-PPPT3(3))*(PPPT4(2)-PPPT3(2))
  128. C=PPPT3(3)*(PPPT4(2)-PPPT3(2))
  129. D=PPPT3(2)*(PPPT4(3)-PPPT3(3))
  130. TEST=(A-B)/(C-D)
  131. IF (TEST.GT.(-XZPREC)) THEN
  132. *** test par rapport à droite ppt4-ppt1, origine repere face=ppt3
  133. DO 80 I=1,NDIM
  134. PPPT4(I)=PPT4(I)-PPT3(I)
  135. PPPT1(I)=PPT1(I)-PPT3(I)
  136. XXXINT(I)=XXINT(I)-PPT3(I)
  137. 80 CONTINUE
  138. A=(XXXINT(2)-PPPT4(2))*(PPPT1(3)-PPPT4(3))
  139. B=(XXXINT(3)-PPPT4(3))*(PPPT1(2)-PPPT4(2))
  140. C=PPPT4(3)*(PPPT1(2)-PPPT4(2))
  141. D=PPPT4(2)*(PPPT1(3)-PPPT4(3))
  142. TEST=(A-B)/(C-D)
  143. IF (TEST.GT.(-XZPREC)) ITEST2=1
  144. ENDIF
  145. ELSE
  146. *** cas ou ppt3=ppt4 (3 pts seulement definissent la face)
  147. *** test par rapport à droite ppt3-ppt1, origine repere face=ppt2
  148. DO 90 I=1,NDIM
  149. PPPT3(I)=PPT3(I)-PPT2(I)
  150. PPPT1(I)=PPT1(I)-PPT2(I)
  151. XXXINT(I)=XXINT(I)-PPT2(I)
  152. 90 CONTINUE
  153. A=(XXXINT(2)-PPPT3(2))*(PPPT1(3)-PPPT3(3))
  154. B=(XXXINT(3)-PPPT3(3))*(PPPT1(2)-PPPT3(2))
  155. C=PPPT3(3)*(PPPT1(2)-PPPT3(2))
  156. D=PPPT3(2)*(PPPT1(3)-PPPT3(3))
  157. TEST=(A-B)/(C-D)
  158. IF (TEST.GT.(-XZPREC)) ITEST2=1
  159. ENDIF
  160. ENDIF
  161. ENDIF
  162. ENDIF
  163. *** verification des deux conditions à la fois
  164. IF ((ITEST1.EQ.1).AND.(ITEST2.EQ.1)) ITEST=1
  165.  
  166. RETURN
  167. END
  168.  
  169.  
  170.  
  171.  
  172.  
  173.  
  174.  
  175.  
  176.  

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