Télécharger qualt.eso

Retour à la liste

Numérotation des lignes :

  1. C QUALT SOURCE JC220346 16/11/29 21:15:32 9221
  2. C Calcul qualite d'un tetraedre par example le rapport du rayon de
  3. C la sphere inscrite au rayon de la sphere circonscrite
  4. C
  5. FUNCTION QUALT(I1,I2,I3,I4)
  6. C
  7. IMPLICIT INTEGER(I-N)
  8. IMPLICIT REAL*8(A-H,O-Z)
  9. -INC TDEMAIT
  10. -INC CCOPTIO
  11. -INC CCREEL
  12. VTOT=VOL(I1,I2,I3,I4)
  13. *
  14. x1=xyz(1,i1)
  15. y1=xyz(2,i1)
  16. z1=xyz(3,i1)
  17. x2=xyz(1,i2)
  18. y2=xyz(2,i2)
  19. z2=xyz(3,i2)
  20. x3=xyz(1,i3)
  21. y3=xyz(2,i3)
  22. z3=xyz(3,i3)
  23. x4=xyz(1,i4)
  24. y4=xyz(2,i4)
  25. z4=xyz(3,i4)
  26. * sphere circonscrite le centre est l'intersection des plan mediateurs
  27. * de 12 13 et 14
  28. x21=x2-x1
  29. y21=y2-y1
  30. z21=z2-z1
  31. s21=sqrt(x21**2+y21**2+z21**2)
  32. x21n=x21/s21
  33. y21n=y21/s21
  34. z21n=z21/s21
  35. x12m=(x1+x2)/2
  36. y12m=(y1+y2)/2
  37. z12m=(z1+z2)/2
  38. *
  39. x31=x3-x1
  40. y31=y3-y1
  41. z31=z3-z1
  42. s31=sqrt(x31**2+y31**2+z31**2)
  43. x31n=x31/s31
  44. y31n=y31/s31
  45. z31n=z31/s31
  46. x13m=(x1+x3)/2
  47. y13m=(y1+y3)/2
  48. z13m=(z1+z3)/2
  49. *
  50. x41=x4-x1
  51. y41=y4-y1
  52. z41=z4-z1
  53. s41=sqrt(x41**2+y41**2+z41**2)
  54. x41n=x41/s41
  55. y41n=y41/s41
  56. z41n=z41/s41
  57. x14m=(x1+x4)/2
  58. y14m=(y1+y4)/2
  59. z14m=(z1+z4)/2
  60. * mediateur de 12
  61. xs=x12m
  62. ys=y12m
  63. zs=z12m
  64. * direction perpendiculaire a 12 dans le plan 123
  65. scal=x31*x21n+y31*y21n+z31*z21n
  66. xd1=x31-scal*x21n
  67. yd1=y31-scal*y21n
  68. zd1=z31-scal*z21n
  69. * distance au plan mediateur de 13
  70. scal=(x13m-xs)*x31n+(y13m-ys)*y31n+(z13m-zs)*z31n
  71. scald=xd1*x31n+yd1*y31n+zd1*z31n+xspeti
  72. coef=scal/scald
  73. * mediateur de 123
  74. xs=xs+coef*xd1
  75. ys=ys+coef*yd1
  76. zs=zs+coef*zd1
  77. * direction normale a 123
  78. xd1=y21*z31-z21*y31
  79. yd1=z21*x31-x21*z31
  80. zd1=x21*y31-y21*x31
  81. * distance au plan mediateur de 14
  82. scal=(x14m-xs)*x41n+(y14m-ys)*y41n+(z14m-zs)*z41n
  83. scald=xd1*x41n+yd1*y41n+zd1*z41n+xspeti
  84. coef=scal/scald
  85. * mediateur de 1234
  86. xs=xs+coef*xd1
  87. ys=ys+coef*yd1
  88. zs=zs+coef*zd1
  89. * verification
  90. rc=sqrt((xs-x1)**2+(ys-y1)**2+(zs-z1)**2)
  91. * r2=sqrt((xs-x2)**2+(ys-y2)**2+(zs-z2)**2)
  92. * r3=sqrt((xs-x3)**2+(ys-y3)**2+(zs-z3)**2)
  93. * r4=sqrt((xs-x4)**2+(ys-y4)**2+(zs-z4)**2)
  94. * write (6,*) ' rayons ',r1,r2,r3,r4
  95. * write (6,*) ' rayon sphere circonscrite ',rc
  96. *
  97. * sphere inscrite
  98. *
  99. * normale a 123
  100. xn123=(y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)
  101. yn123=(z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)
  102. zn123=(x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)
  103. sn123=sqrt(xn123**2+yn123**2+zn123**2)
  104. xn123=xn123/sn123
  105. yn123=yn123/sn123
  106. zn123=zn123/sn123
  107. * normale a 134
  108. xn134=(y3-y1)*(z4-z1)-(z3-z1)*(y4-y1)
  109. yn134=(z3-z1)*(x4-x1)-(x3-x1)*(z4-z1)
  110. zn134=(x3-x1)*(y4-y1)-(y3-y1)*(x4-x1)
  111. sn134=sqrt(xn134**2+yn134**2+zn134**2)
  112. xn134=xn134/sn134
  113. yn134=yn134/sn134
  114. zn134=zn134/sn134
  115. * normale a 142
  116. xn142=(y4-y1)*(z2-z1)-(z4-z1)*(y2-y1)
  117. yn142=(z4-z1)*(x2-x1)-(x4-x1)*(z2-z1)
  118. zn142=(x4-x1)*(y2-y1)-(y4-y1)*(x2-x1)
  119. sn142=sqrt(xn142**2+yn142**2+zn142**2)
  120. xn142=xn142/sn142
  121. yn142=yn142/sn142
  122. zn142=zn142/sn142
  123. * direction du centre a partir de 1
  124. xd1=(yn123-yn134)*(zn134-zn142)-(zn123-zn134)*(yn134-yn142)
  125. yd1=(zn123-zn134)*(xn134-xn142)-(xn123-xn134)*(zn134-zn142)
  126. zd1=(xn123-xn134)*(yn134-yn142)-(yn123-yn134)*(xn134-xn142)
  127. * normale a 432
  128. xn432=(y3-y4)*(z2-z4)-(z3-z4)*(y2-y4)
  129. yn432=(z3-z4)*(x2-x4)-(x3-x4)*(z2-z4)
  130. zn432=(x3-x4)*(y2-y4)-(y3-y4)*(x2-x4)
  131. sn432=sqrt(xn432**2+yn432**2+zn432**2)
  132. xn432=xn432/sn432
  133. yn432=yn432/sn432
  134. zn432=zn432/sn432
  135. * petit calcul
  136. xnum=(x2-x1)*xn432+(y2-y1)*yn432+(z2-z1)*zn432
  137. xden=xd1*(xn123-xn432)+yd1*(yn123-yn432)+zd1*(zn123-zn432)+xspeti
  138. rap=-xnum/xden
  139. xi=x1+rap*xd1
  140. yi=y1+rap*yd1
  141. zi=z1+rap*zd1
  142. * verif
  143. * ri123=(xi-x1)*xn123+(yi-y1)*yn123+(zi-z1)*zn123
  144. * ri134=(xi-x1)*xn134+(yi-y1)*yn134+(zi-z1)*zn134
  145. * ri142=(xi-x1)*xn142+(yi-y1)*yn142+(zi-z1)*zn142
  146. * ri234=(xi-x2)*xn432+(yi-y2)*yn432+(zi-z2)*zn432
  147. ri=abs((xi-x1)*xn123+(yi-y1)*yn123+(zi-z1)*zn123)
  148. * write (6,*) ' rayon sphere inscrite ',ri
  149. * arete maxi
  150. x32=x3-x2
  151. y32=y3-y2
  152. z32=z3-z2
  153. s32=sqrt(x32**2+y32**2+z32**2)
  154. x42=x4-x2
  155. y42=y4-y2
  156. z42=z4-z2
  157. s42=sqrt(x42**2+y42**2+z42**2)
  158. x43=x4-x3
  159. y43=y4-y3
  160. z43=z4-z3
  161. s43=sqrt(x43**2+y43**2+z43**2)
  162. s=max(s21,s31,s41,s32,s42,s43)
  163. * qualt=ri/rc
  164. qualt=ri/s
  165. if (vtot.le.0.) qualt=-1e8
  166. * write (6,*) ' qualt ',qualt
  167. END
  168.  
  169.  
  170.  
  171.  
  172.  

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