1 SUBROUTINE RAY(N,IND,XXX,YYY,ZZZ,RR,X0,Y0,ZC,IOPT) 2C 3 DIMENSION XXX(*),YYY(*),ZZZ(*),A(10),B(4) 4 INTEGER MDIAG(4) 5 DATA EPS /1.E-5/ 6 DATA MDIAG / 1 , 3 , 6 , 10 / 7 DATA NITMAX / 40 / 8C 9C IOPT=0 on calcule rayon 1 (cercle des moindres carres) 10C IOPT=1 on calcule rayon 1 (cercle des moindres carres - centre sur l'axe optique) 11C 12 IF (IOPT.EQ.0) THEN 13 IF (IND.GT.4) THEN 14 IND2 = 4 15 ELSE 16 IND2 = 3 17 ENDIF 18 CALL RAY0(IND2,XXX,YYY,ZZZ,XC,YC,ZC,RR0) 19 RR = RR0 20 ELSE 21 XC = X0 22 YC = Y0 23 RR0 = RR 24 ENDIF 25 DNI = REAL(IND) 26 DN4 = REAL(4*IND) 27 NITER = 0 28 D = 0. 29 SX = -DNI*XC 30 SY = -DNI*YC 31 SZ = -DNI*ZC 32 DO I=1,IND 33 D = D + (XXX(I)-XC)**2 + (YYY(I)-YC)**2 + (ZZZ(I)-ZC)**2 34 SX = SX + XXX(I) 35 SY = SY + YYY(I) 36 SZ = SZ + ZZZ(I) 37 ENDDO 38 USR = 1./RR 39 USR2 = USR**2 40 D = D*USR2 - DNI 41C 42 100 NITER = NITER+1 43ccc print*,niter,rr,d,' .. ',xc,yc,zc 44 SX = SX*4. 45 SY = SY*4. 46 SZ = SZ*4. 47 B(1) = -SX*D 48 B(2) = -SY*D 49 B(3) = -SZ*D 50 B(4) = -D*(2.*D+DN4)*RR 51C 52 A(1) = DN4*D + .5*SX**2 53 A(2) = .5*SY*SX 54 A(3) = DN4*D + .5*SY**2 55 A(4) = .5*SX*SZ 56 A(5) = .5*SZ*SY 57 A(6) = DN4*D + .5*SZ**2 58 A(7) = 2.*SX*(D+DNI)*USR 59 A(8) = 2.*SY*(D+DNI)*USR 60 A(9) = 2.*SZ*(D+DNI)*USR 61 A(10) = 6.*D*D + 12.*D+DNI + 8.*DNI*DNI 62 CALL SLVSPI(A,B,B,MDIAG,4,1,IERR) 63 IF (IERR.GT.0) THEN 64cc PRINT*,'*****ERREUR',N,' iter',NITER 65 RR = RR0 66 RETURN 67 ENDIF 68 DEL = MAX(ABS(B(1)),ABS(B(2)),ABS(B(3)),ABS(B(4))) 69 IF (IOPT.EQ.0) THEN 70 XC = XC - B(1) 71 YC = YC - B(2) 72 ENDIF 73 ZC = ZC - B(3) 74 RR = RR - B(4) 75C 76 D = 0. 77 SX = -DNI*XC 78 SY = -DNI*YC 79 SZ = -DNI*ZC 80 DO I=1,IND 81 D = D + (XXX(I)-XC)**2 + (YYY(I)-YC)**2 + (ZZZ(I)-ZC)**2 82 SX = SX + XXX(I) 83 SY = SY + YYY(I) 84 SZ = SZ + ZZZ(I) 85 ENDDO 86 USR = 1./RR 87 USR2 = USR**2 88 D = D*USR2 - DNI 89 IF (NITER.LT.NITMAX.AND.DEL.GE.EPS.AND.ABS(D).GE.EPS) GOTO 100 90C 91cc print*,'noeud',n,' ray',rr,rr0,' err',del,d,niter 92 IF (NITER.GE.NITMAX.OR.RR.LE.1.) THEN 93 print*, 94 $'noeud',n,' iopt',iopt,' ray',rr,rr0,' err',del,abs(d),niter 95 RR = RR0 96 ENDIF 97C 98 END 99C================================================================== 100 SUBROUTINE RAY0(IND,XXX,YYY,ZZZ,XC,YC,ZC,RR) 101C 102 DIMENSION XXX(*),YYY(*),ZZZ(*) 103 DIMENSION A(3,3),B(3,3),XX(4),YY(4),ZZ(4),BB(3) 104 INTEGER IPER(4,5) 105 DATA IPER / 1,2,3,4 , 1,2,3,5 , 1,2,4,5 , 1,3,4,5 , 2,3,4,5 / 106C 107 IF (IND.EQ.3) THEN 108 NR = 1 109 ELSE 110 NR = 5 111 ENDIF 112 RR = 0. 113 NBRAY = 0 114 XCC = 0. 115 YCC = 0. 116 ZCC = 0. 117 DO K=1,NR 118 DO I=1,4 119 II = IPER(I,K) 120 XX(I) = XXX(II) 121 YY(I) = YYY(II) 122 ZZ(I) = ZZZ(II) 123 ENDDO 124C 125 DO I=1,3 126 A(I,1) = XX(I+1)-XX(I) 127 A(I,2) = YY(I+1)-YY(I) 128 A(I,3) = ZZ(I+1)-ZZ(I) 129 BB(I) = XX(I+1)**2-XX(I)**2 130 & + YY(I+1)**2-YY(I)**2 + ZZ(I+1)**2-ZZ(I)**2 131 ENDDO 132C 133 CALL INV3X3(A,B,IERR) 134 IF (IERR.EQ.0) THEN 135 XC = .5*( B(1,1)*BB(1) + B(1,2)*BB(2) + B(1,3)*BB(3) ) 136 YC = .5*( B(2,1)*BB(1) + B(2,2)*BB(2) + B(2,3)*BB(3) ) 137 ZC = .5*( B(3,1)*BB(1) + B(3,2)*BB(2) + B(3,3)*BB(3) ) 138 RR = RR+SQRT((XX(1)-XC)**2 + (YY(1)-YC)**2 + (ZZ(1)-ZC)**2) 139 XCC = XCC+XC 140 YCC = YCC+YC 141 ZCC = ZCC+ZC 142 NBRAY = NBRAY+1 143 ENDIF 144 ENDDO 145C 146 IF (NBRAY.NE.0) THEN 147 RR = RR/REAL(NBRAY) 148 XC = XCC/REAL(NBRAY) 149 YC = YCC/REAL(NBRAY) 150 ZC = ZCC/REAL(NBRAY) 151 ELSE 152 XC = 0. 153 YC = 0. 154 ZC = 0. 155 RR = 8. 156 ENDIF 157C 158 END 159C===================================================================== 160 SUBROUTINE RAYTMS(X,Y,Z,ICOR,RAYON0,BIG,NUMNP,NBCORN) 161 DIMENSION X(*),Y(*),Z(*),ICOR(*),RAYON0(*) 162 REAL*8 BIG 163C 164 IF (NBCORN.GT.0) THEN 165 DISCEN = BIG 166 DO I=1,NUMNP 167 IF (ICOR(I).NE.0) THEN 168 IF ((X(I)**2+Y(I)**2).LT.DISCEN) THEN 169 DISCEN = X(I)**2+Y(I)**2 170 ICENT = I 171 ENDIF 172 ENDIF 173 ENDDO 174 IP1 = 0 175 IP2 = 0 176 IP3 = 0 177 DIS1 = BIG 178 DIS2 = BIG 179 DIS3 = BIG 180 DO I=1,NUMNP 181 IF (ICOR(I).NE.0.AND.I.NE.ICENT) THEN 182 DD = (X(I)-X(ICENT))**2 + (Y(I)-Y(ICENT))**2 183 IF (DD.LT.DIS3) THEN 184 IF (DD.LT.DIS2) THEN 185 IP3 = IP2 186 DIS3 = DIS2 187 IF (DD.LT.DIS1) THEN 188 IP2 = IP1 189 IP1 = I 190 DIS2 = DIS1 191 DIS1 = DD 192 ELSE 193 IP2 = I 194 DIS2 = DD 195 ENDIF 196 ELSE 197 IP3 = I 198 DIS3 = DD 199 ENDIF 200 ENDIF 201 ENDIF 202 ENDDO 203 ELSE 204 IP1 = 0 205 IP2 = 0 206 IP3 = 0 207 ICENT = 0 208 RETURN 209 ENDIF 210 DO I=1,NUMNP 211 IF (ICOR(I).NE.0.AND.I.NE.ICENT) THEN 212 IF (Z(I).NE.Z(ICENT)) THEN 213 RAYON0(I) = .5*( Z(ICENT)-Z(I) 214 & + (X(I)**2+Y(I)**2)/(Z(ICENT)-Z(I)) ) 215 ELSE 216 RAYON0(I) = 1000. 217 ENDIF 218 ELSE 219 RAYON0(I) = 0. 220 ENDIF 221 ENDDO 222 IF (NBCORN.NE.0) THEN 223 II = 0 224 RR = 0. 225 IF (IP1.NE.0) THEN 226 II = II+1 227 RR = RR+RAYON0(IP1) 228 ENDIF 229 IF (IP2.NE.0) THEN 230 II = II+1 231 RR = RR+RAYON0(IP2) 232 ENDIF 233 IF (IP3.NE.0) THEN 234 II = II+1 235 RR = RR+RAYON0(IP3) 236 ENDIF 237 IF (II.GT.0) THEN 238 RAYON0(ICENT) = RR/REAL(II) 239 ELSE 240 RAYON0(ICENT) = 1000. 241 ENDIF 242 ENDIF 243 END 244 245