1 : C QUALT SOURCE PV 07/10/05 21:16:27 5890 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 :
© Cast3M 2003 - All rights reserved.
Disclaimer