C QUALT     SOURCE    PV090527  25/12/21    07:19:11     12434          
C   Calcul qualite d'un tetraedre par example le rapport du rayon de
C   la sphere inscrite au rayon de la sphere circonscrite
C
       FUNCTION QUALT(I1,I2,I3,I4)
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
-INC TDEMAIT

-INC PPARAM
-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
      if (abs(xden) .lt.xpetit) xden=xden+sign(xpetit,xden)
      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