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