      SUBROUTINE RAY(N,IND,XXX,YYY,ZZZ,RR,X0,Y0,ZC,IOPT)
C
      DIMENSION XXX(*),YYY(*),ZZZ(*),A(10),B(4)
      INTEGER MDIAG(4)
      DATA EPS /1.E-5/
      DATA MDIAG / 1 , 3 , 6 , 10 /
      DATA NITMAX / 40 /
C
C IOPT=0 on calcule rayon 1 (cercle des moindres carres)
C IOPT=1 on calcule rayon 1 (cercle des moindres carres - centre sur l'axe optique)
C
      IF (IOPT.EQ.0) THEN
        IF (IND.GT.4) THEN
          IND2 = 4
        ELSE
          IND2 = 3
        ENDIF
        CALL RAY0(IND2,XXX,YYY,ZZZ,XC,YC,ZC,RR0)
        RR = RR0
      ELSE
        XC = X0
        YC = Y0
        RR0 = RR
      ENDIF
      DNI = REAL(IND)
      DN4 = REAL(4*IND)
      NITER = 0
      D  = 0.
      SX = -DNI*XC
      SY = -DNI*YC
      SZ = -DNI*ZC
      DO I=1,IND
        D = D + (XXX(I)-XC)**2 + (YYY(I)-YC)**2 + (ZZZ(I)-ZC)**2
        SX = SX + XXX(I)
        SY = SY + YYY(I)
        SZ = SZ + ZZZ(I)
      ENDDO
      USR  = 1./RR
      USR2 = USR**2
      D = D*USR2 - DNI
C
  100 NITER = NITER+1
ccc      print*,niter,rr,d,' .. ',xc,yc,zc
        SX = SX*4.
        SY = SY*4.
        SZ = SZ*4.
        B(1) = -SX*D
        B(2) = -SY*D
        B(3) = -SZ*D
        B(4) = -D*(2.*D+DN4)*RR
C
        A(1) = DN4*D + .5*SX**2        
        A(2) = .5*SY*SX                
        A(3) = DN4*D + .5*SY**2        
        A(4) = .5*SX*SZ                
        A(5) = .5*SZ*SY                
        A(6) = DN4*D + .5*SZ**2        
        A(7) = 2.*SX*(D+DNI)*USR
        A(8) = 2.*SY*(D+DNI)*USR
        A(9) = 2.*SZ*(D+DNI)*USR
        A(10) = 6.*D*D + 12.*D+DNI + 8.*DNI*DNI
        CALL SLVSPI(A,B,B,MDIAG,4,1,IERR)
        IF (IERR.GT.0) THEN
cc          PRINT*,'*****ERREUR',N,' iter',NITER
          RR = RR0
          RETURN
        ENDIF
        DEL = MAX(ABS(B(1)),ABS(B(2)),ABS(B(3)),ABS(B(4)))
        IF (IOPT.EQ.0) THEN
          XC = XC - B(1)
          YC = YC - B(2)
        ENDIF
        ZC = ZC - B(3)
        RR = RR - B(4)
C
        D  = 0.
        SX = -DNI*XC
        SY = -DNI*YC
        SZ = -DNI*ZC
        DO I=1,IND
          D = D + (XXX(I)-XC)**2 + (YYY(I)-YC)**2 + (ZZZ(I)-ZC)**2
          SX = SX + XXX(I)
          SY = SY + YYY(I)
          SZ = SZ + ZZZ(I)
        ENDDO
        USR  = 1./RR
        USR2 = USR**2
        D = D*USR2 - DNI
        IF (NITER.LT.NITMAX.AND.DEL.GE.EPS.AND.ABS(D).GE.EPS) GOTO 100
C
cc      print*,'noeud',n,' ray',rr,rr0,' err',del,d,niter
      IF (NITER.GE.NITMAX.OR.RR.LE.1.) THEN
        print*,
     $'noeud',n,' iopt',iopt,' ray',rr,rr0,' err',del,abs(d),niter
        RR = RR0
      ENDIF
C     
      END
C==================================================================
      SUBROUTINE RAY0(IND,XXX,YYY,ZZZ,XC,YC,ZC,RR)
C
      DIMENSION XXX(*),YYY(*),ZZZ(*)
      DIMENSION A(3,3),B(3,3),XX(4),YY(4),ZZ(4),BB(3)
      INTEGER IPER(4,5)
      DATA IPER / 1,2,3,4 , 1,2,3,5 , 1,2,4,5 , 1,3,4,5 , 2,3,4,5 /
C
      IF (IND.EQ.3) THEN
        NR = 1
      ELSE
        NR = 5
      ENDIF
      RR = 0.
      NBRAY = 0
      XCC = 0.
      YCC = 0.
      ZCC = 0.
      DO K=1,NR
        DO I=1,4
          II = IPER(I,K)
          XX(I) = XXX(II)
          YY(I) = YYY(II)
          ZZ(I) = ZZZ(II)
        ENDDO
C
        DO I=1,3
          A(I,1) = XX(I+1)-XX(I)
          A(I,2) = YY(I+1)-YY(I)
          A(I,3) = ZZ(I+1)-ZZ(I)
          BB(I) = XX(I+1)**2-XX(I)**2 
     &          + YY(I+1)**2-YY(I)**2 + ZZ(I+1)**2-ZZ(I)**2
        ENDDO
C
        CALL INV3X3(A,B,IERR)
        IF (IERR.EQ.0) THEN
          XC = .5*( B(1,1)*BB(1) + B(1,2)*BB(2) + B(1,3)*BB(3) )
          YC = .5*( B(2,1)*BB(1) + B(2,2)*BB(2) + B(2,3)*BB(3) )
          ZC = .5*( B(3,1)*BB(1) + B(3,2)*BB(2) + B(3,3)*BB(3) )
          RR = RR+SQRT((XX(1)-XC)**2 + (YY(1)-YC)**2 + (ZZ(1)-ZC)**2)
          XCC = XCC+XC
          YCC = YCC+YC
          ZCC = ZCC+ZC
          NBRAY = NBRAY+1
        ENDIF
      ENDDO
C
      IF (NBRAY.NE.0) THEN
        RR = RR/REAL(NBRAY)
        XC = XCC/REAL(NBRAY) 
        YC = YCC/REAL(NBRAY) 
        ZC = ZCC/REAL(NBRAY) 
      ELSE
        XC = 0.
        YC = 0.
        ZC = 0.
        RR = 8.
      ENDIF
C     
      END
C=====================================================================
      SUBROUTINE RAYTMS(X,Y,Z,ICOR,RAYON0,BIG,NUMNP,NBCORN)
      DIMENSION X(*),Y(*),Z(*),ICOR(*),RAYON0(*)
      REAL*8 BIG
C
      IF (NBCORN.GT.0) THEN
        DISCEN = BIG
        DO I=1,NUMNP
          IF (ICOR(I).NE.0) THEN
            IF ((X(I)**2+Y(I)**2).LT.DISCEN) THEN
              DISCEN = X(I)**2+Y(I)**2
              ICENT = I
            ENDIF
          ENDIF
        ENDDO
        IP1 = 0
        IP2 = 0
        IP3 = 0
        DIS1 = BIG
        DIS2 = BIG
        DIS3 = BIG
        DO I=1,NUMNP
          IF (ICOR(I).NE.0.AND.I.NE.ICENT) THEN
            DD = (X(I)-X(ICENT))**2 + (Y(I)-Y(ICENT))**2
            IF (DD.LT.DIS3) THEN
              IF (DD.LT.DIS2) THEN
                IP3 = IP2
                DIS3 = DIS2
                IF (DD.LT.DIS1) THEN
                  IP2 = IP1
                  IP1 = I
                  DIS2 = DIS1
                  DIS1 = DD
                ELSE
                  IP2 = I
                  DIS2 = DD
                ENDIF
              ELSE
                IP3 = I
                DIS3 = DD
              ENDIF
            ENDIF
          ENDIF
        ENDDO
      ELSE
        IP1 = 0
        IP2 = 0
        IP3 = 0
        ICENT = 0
        RETURN
      ENDIF
      DO I=1,NUMNP
        IF (ICOR(I).NE.0.AND.I.NE.ICENT) THEN
          IF (Z(I).NE.Z(ICENT)) THEN
            RAYON0(I) = .5*( Z(ICENT)-Z(I)
     &                     + (X(I)**2+Y(I)**2)/(Z(ICENT)-Z(I)) )
          ELSE
            RAYON0(I) = 1000.
          ENDIF
        ELSE
          RAYON0(I) = 0.
        ENDIF
      ENDDO
      IF (NBCORN.NE.0) THEN
        II = 0
        RR = 0.
        IF (IP1.NE.0) THEN
          II = II+1
          RR = RR+RAYON0(IP1)
        ENDIF
        IF (IP2.NE.0) THEN
          II = II+1
          RR = RR+RAYON0(IP2)
        ENDIF
        IF (IP3.NE.0) THEN
          II = II+1
          RR = RR+RAYON0(IP3)
        ENDIF
        IF (II.GT.0) THEN
          RAYON0(ICENT) = RR/REAL(II)
        ELSE
          RAYON0(ICENT) = 1000.
        ENDIF
      ENDIF
      END

