Numérotation des lignes :

C QUALT     SOURCE    JC220346  16/11/29    21:15:32     9221           C   Calcul qualite d'un tetraedre par example le rapport du rayon deC   la sphere inscrite au rayon de la sphere circonscriteC       FUNCTION QUALT(I1,I2,I3,I4)C      IMPLICIT INTEGER(I-N)      IMPLICIT REAL*8(A-H,O-Z)-INC TDEMAIT-INC CCOPTIO-INC CCREEL       VTOT=VOL(I1,I2,I3,I4)*      x1=xyz(1,i1)      y1=xyz(2,i1)      z1=xyz(3,i1)      x2=xyz(1,i2)      y2=xyz(2,i2)      z2=xyz(3,i2)      x3=xyz(1,i3)      y3=xyz(2,i3)      z3=xyz(3,i3)      x4=xyz(1,i4)      y4=xyz(2,i4)      z4=xyz(3,i4)*  sphere circonscrite le centre est l'intersection des plan mediateurs*  de 12 13 et 14      x21=x2-x1      y21=y2-y1      z21=z2-z1      s21=sqrt(x21**2+y21**2+z21**2)      x21n=x21/s21      y21n=y21/s21      z21n=z21/s21      x12m=(x1+x2)/2      y12m=(y1+y2)/2      z12m=(z1+z2)/2*      x31=x3-x1      y31=y3-y1      z31=z3-z1      s31=sqrt(x31**2+y31**2+z31**2)      x31n=x31/s31      y31n=y31/s31      z31n=z31/s31      x13m=(x1+x3)/2      y13m=(y1+y3)/2      z13m=(z1+z3)/2*      x41=x4-x1      y41=y4-y1      z41=z4-z1      s41=sqrt(x41**2+y41**2+z41**2)      x41n=x41/s41      y41n=y41/s41      z41n=z41/s41      x14m=(x1+x4)/2      y14m=(y1+y4)/2      z14m=(z1+z4)/2* mediateur de 12      xs=x12m      ys=y12m      zs=z12m*  direction perpendiculaire a 12 dans le plan 123      scal=x31*x21n+y31*y21n+z31*z21n      xd1=x31-scal*x21n      yd1=y31-scal*y21n      zd1=z31-scal*z21n*  distance au plan mediateur de 13      scal=(x13m-xs)*x31n+(y13m-ys)*y31n+(z13m-zs)*z31n      scald=xd1*x31n+yd1*y31n+zd1*z31n+xspeti      coef=scal/scald* mediateur de 123      xs=xs+coef*xd1      ys=ys+coef*yd1      zs=zs+coef*zd1*  direction normale a 123      xd1=y21*z31-z21*y31      yd1=z21*x31-x21*z31      zd1=x21*y31-y21*x31*  distance au plan mediateur de 14      scal=(x14m-xs)*x41n+(y14m-ys)*y41n+(z14m-zs)*z41n      scald=xd1*x41n+yd1*y41n+zd1*z41n+xspeti      coef=scal/scald* mediateur de 1234      xs=xs+coef*xd1      ys=ys+coef*yd1      zs=zs+coef*zd1* verification      rc=sqrt((xs-x1)**2+(ys-y1)**2+(zs-z1)**2)*     r2=sqrt((xs-x2)**2+(ys-y2)**2+(zs-z2)**2)*     r3=sqrt((xs-x3)**2+(ys-y3)**2+(zs-z3)**2)*     r4=sqrt((xs-x4)**2+(ys-y4)**2+(zs-z4)**2)*     write (6,*) ' rayons ',r1,r2,r3,r4*     write (6,*) ' rayon sphere circonscrite ',rc**  sphere inscrite** normale a 123      xn123=(y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)      yn123=(z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)      zn123=(x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)      sn123=sqrt(xn123**2+yn123**2+zn123**2)      xn123=xn123/sn123      yn123=yn123/sn123      zn123=zn123/sn123* normale a 134      xn134=(y3-y1)*(z4-z1)-(z3-z1)*(y4-y1)      yn134=(z3-z1)*(x4-x1)-(x3-x1)*(z4-z1)      zn134=(x3-x1)*(y4-y1)-(y3-y1)*(x4-x1)      sn134=sqrt(xn134**2+yn134**2+zn134**2)      xn134=xn134/sn134      yn134=yn134/sn134      zn134=zn134/sn134* normale a 142      xn142=(y4-y1)*(z2-z1)-(z4-z1)*(y2-y1)      yn142=(z4-z1)*(x2-x1)-(x4-x1)*(z2-z1)      zn142=(x4-x1)*(y2-y1)-(y4-y1)*(x2-x1)      sn142=sqrt(xn142**2+yn142**2+zn142**2)      xn142=xn142/sn142      yn142=yn142/sn142      zn142=zn142/sn142*  direction du centre a partir de 1      xd1=(yn123-yn134)*(zn134-zn142)-(zn123-zn134)*(yn134-yn142)      yd1=(zn123-zn134)*(xn134-xn142)-(xn123-xn134)*(zn134-zn142)      zd1=(xn123-xn134)*(yn134-yn142)-(yn123-yn134)*(xn134-xn142)*  normale a 432      xn432=(y3-y4)*(z2-z4)-(z3-z4)*(y2-y4)      yn432=(z3-z4)*(x2-x4)-(x3-x4)*(z2-z4)      zn432=(x3-x4)*(y2-y4)-(y3-y4)*(x2-x4)      sn432=sqrt(xn432**2+yn432**2+zn432**2)      xn432=xn432/sn432      yn432=yn432/sn432      zn432=zn432/sn432*  petit calcul      xnum=(x2-x1)*xn432+(y2-y1)*yn432+(z2-z1)*zn432      xden=xd1*(xn123-xn432)+yd1*(yn123-yn432)+zd1*(zn123-zn432)+xspeti      rap=-xnum/xden      xi=x1+rap*xd1      yi=y1+rap*yd1      zi=z1+rap*zd1*  verif*     ri123=(xi-x1)*xn123+(yi-y1)*yn123+(zi-z1)*zn123*     ri134=(xi-x1)*xn134+(yi-y1)*yn134+(zi-z1)*zn134*     ri142=(xi-x1)*xn142+(yi-y1)*yn142+(zi-z1)*zn142*     ri234=(xi-x2)*xn432+(yi-y2)*yn432+(zi-z2)*zn432      ri=abs((xi-x1)*xn123+(yi-y1)*yn123+(zi-z1)*zn123)*     write (6,*) ' rayon sphere inscrite ',ri*  arete maxi      x32=x3-x2      y32=y3-y2      z32=z3-z2      s32=sqrt(x32**2+y32**2+z32**2)      x42=x4-x2      y42=y4-y2      z42=z4-z2      s42=sqrt(x42**2+y42**2+z42**2)      x43=x4-x3      y43=y4-y3      z43=z4-z3      s43=sqrt(x43**2+y43**2+z43**2)      s=max(s21,s31,s41,s32,s42,s43)*     qualt=ri/rc      qualt=ri/s      if (vtot.le.0.) qualt=-1e8*     write (6,*) ' qualt ',qualt      END

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