C
C Program to determine optimized spherical 3-D forward stars
C
C Input from console:
C (1) NH,'STARS','TABLE',/
C     NH...   Maximum forward star size.
C     'STARS'... Name of the output file with 3-D forward stars.
C             Renamed to 'net.fs3', it may serve as the data file for
C             the NET program.
C     'TABLE'... Table filename.  Just informative output.
C     Default: NH=27, 'STARS'='net.fs', 'TABLE'='table.fs'.
C
C Output file 'STARS':
C (1) For each forward star (1A) and (1B):
C (1A) NLE, NEDGE
C     NLE...  Level of a forward star.
C     NEDGE...Number of edges.
C (1B) (I1(I),I2(I),I3(I),I=1,NEDGE)
C     List of edges.  Just edges with 0.LE.I1.LE.I2.LE.I3 are listed.
C (2) 0, 0
C
C Output file 'TABLE':
C For each forward star (1):
C (1) INT(H*H), H, R*R/2, H*H*R*R
C     H...    Radius of a forward star in grid intervals.
C     R...    Radius of the largest angular circle containing no edge
C             of the forward star.
C Just the forward stars with smaller R's than the preceding ones are
C listed.
C
C Date: 1997, January 19
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Filenames:
      CHARACTER*80 FILE1,FILE2
C
C     Logical unit numbers:
      INTEGER LU1,LU2
      PARAMETER (LU1=1)
      PARAMETER (LU2=2)
C
C     Array dimensions:
      INTEGER MLE,MH,MEDGE,MCIRC
      PARAMETER (MLE=80,MH=MLE**2+2,MEDGE=(MLE+1)**3/13)
      PARAMETER (MCIRC=(MLE+1)**3/5)
C     Memory: more than 4.2*(MLE+1)**3 B.
C     Time 386/20MHZ: (NH/56)**3 s = (NH/219)**3 min =(NH/860)**3 h.
C
      INTEGER NCIRC,ICIRC,KCIRC1(MCIRC),KCIRC2(MCIRC),KCIRC3(MCIRC)
      INTEGER K1,K2,K3,NH,IH,I1,I2,I3,I3I3,I,J,NEDGE
      INTEGER KS1(MEDGE)
      INTEGER KS2(MEDGE)
      INTEGER KS3(MEDGE),KEDGE(MEDGE),LEDGE(MH),IEDGE,NS
      REAL EDGE1(MEDGE),S10,S11,S12,S13,A101,A112,A123,A130,C1
      REAL EDGE2(MEDGE),S20,S21,S22,S23,A201,A212,A223,A230,C2
      REAL EDGE3(MEDGE),S30,S31,S32,S33,A301,A312,A323,A330,C3
      REAL RCIRC(MCIRC),RADIUS(MH),DH,CS,C,RR,AUX
      REAL S1(MEDGE)
      REAL S2(MEDGE)
      REAL S3(MEDGE)
      EQUIVALENCE (KCIRC1,KS1),(KCIRC2,KS2),(KCIRC3,KS3)
      EQUIVALENCE (KCIRC1,S1),(KCIRC2,S2),(KCIRC3,S3),(RCIRC,KEDGE)
C
C     NCIRC.. Number of circles stored.
C     ICIRC.. Index of a circle.
C     KCIRC1(ICIRC),KCIRC2(ICIRC),KCIRC3(ICIRC)... Edges at the
C             periphery of a circle.
C             KCIRC3(ICIRC)= 0: The circle does not exist.
C             KCIRC3(ICIRC)=-1: Circle centre at the side of the region,
C               between edges 2=(1,1,0) and 3=(1,1,1).
C             KCIRC3(ICIRC)=-2: Circle centre at the side of the region,
C               between edges 3=(1,1,1) and 2=(1,1,0).
C             KCIRC3(ICIRC)=-3: Circle centre at the side of the region,
C               between edges 1=(1,0,0) and 2=(1,1,0).
C     K1,K2,K3... Above indices of a current circle.
C     NH...   Square of the maximum length of an edge.
C     IH...   Square of the length of an edge.
C     I3I3... Limit for I3*I3.
C     I...    Limit for I2*I2+I3*I3.
C     J...    Loop variable.
C     I1,I2,I3... Components of an edge.
C     NEDGE...Number of edges stored.
C     KS1,KS2,KS3... Edges of the resulting forward star (vectors).
C     KEDGE...List of edges of the resulting forward star.
C     LEDGE(IH)... Index of the last edge of square length IH.
C     IEDGE...Index of a current edge when decimating a forward star.
C     NS...   Number of edges in vicinity of the current edge, when
C             decimating a forward star.
C     S1,S2,S3... Unit vectors of edges in vicinity of the current edge,
C             when decimating a forward star.
C     EDGE1(IS),EDGE2(IS),EDGE3(IS)... Unit vector of the IS-th edge.
C     S10,S20,S30... Unit vector of the current edge.
C     S11,S21,S31, S12,S22,S32, S13,S23,S33... Unit vectors of the
C             edges at the periphery of the current circle.
C     AIJK... I-th conponent of the vector connecting points J and K of
C             the unit sphere.
C     C1,C2,C3... Axial vector of the current circle.
C     RCIRC(ICIRC)... Radius**2 of a ICIRC-th circle.
C     RADIUS(IH)... Radius**2 of the maximum circle.
C     DH...   A small real to remove rounding errors.
C     CS,C... Cosines of angular radii of circles.
C     RR...   Square of radius of the largest angular circle containing
C             no edge of the forward star.
C     AUX...  Auxiliary starage location.
C
C.......................................................................
C
C     Opening data files:
      NH=27
      FILE2='net.fs'
      FILE1='table.fs'
      WRITE(*,'(A)')
     * '+Enter maximum f.s. size (27), f.s. filename (''net.fs'')',
     * ' and table filename (''table.fs''): '
      READ(*,*) NH,FILE2,FILE1
      NH=NH*NH+2
      IF(NH.GT.MH) THEN
        STOP 'ERROR: TOO LARGE FORWARD STAR.'
      END IF
      OPEN(LU2,FILE=FILE2)
      OPEN(LU1,FILE=FILE1)
C
      WRITE(*,'(A)') '+    H*H   edges circles'
      WRITE(*,'(3I8,A)') NH,MEDGE,MCIRC,' maximum'
      WRITE(*,'(A)')
C
      DH=1./SQRT(12.*FLOAT(NH))
      DH=DH/2.
      NEDGE=0
      DO 69 IH=1,NH
        DO 58 I1=INT(SQRT(FLOAT(IH)/3.)-DH+1.),INT(SQRT(FLOAT(IH))+DH)
          I=IH-I1*I1
          DO 57 I2=INT(SQRT(FLOAT(I)/2.)-DH+1.),
     *                                   MIN0(INT(SQRT(FLOAT(I))+DH),I1)
            I3I3=I-I2*I2
            I3=INT(SQRT(FLOAT(I3I3))+0.500)
            IF(I3*I3.EQ.I3I3.AND.I3.LE.I2) THEN
              DO 10 J=2,I1
                IF(MOD(I1,J).EQ.0.AND.MOD(I2,J).EQ.0.AND.MOD(I3,J).EQ.0)
     *                                                              THEN
                  GO TO 56
                END IF
   10         CONTINUE
C
C             New edge
              NEDGE=NEDGE+1
              IF(NEDGE.GT.MEDGE) THEN
                STOP 'ERROR: TOO MANY EDGES.'
              END IF
              AUX=SQRT(FLOAT(IH))
              S10=FLOAT(I1)/AUX
              S20=FLOAT(I2)/AUX
              S30=FLOAT(I3)/AUX
              EDGE1(NEDGE)=S10
              EDGE2(NEDGE)=S20
              EDGE3(NEDGE)=S30
C
C             Loop over circles
              DO 49 ICIRC=1,NCIRC
                K3=KCIRC3(ICIRC)
                IF(K3.NE.0) THEN
                  K1=KCIRC1(ICIRC)
                  K2=KCIRC2(ICIRC)
                  S11=EDGE1(K1)
                  S21=EDGE2(K1)
                  S31=EDGE3(K1)
                  S12=EDGE1(K2)
                  S22=EDGE2(K2)
                  S32=EDGE3(K2)
                  A112=S11-S12
                  A212=S21-S22
                  A312=S31-S32
                  IF(K3.GT.0) THEN
C                   Circle centre inside the region:
                    S13=EDGE1(K3)
                    S23=EDGE2(K3)
                    S33=EDGE3(K3)
                    A123=S12-S13
                    A223=S22-S23
                    A323=S32-S33
                    C1=A212*A323-A312*A223
                    C2=A312*A123-A112*A323
                    C3=A112*A223-A212*A123
                  ELSE
C                   Circle centre at a side of the region:
                    IF(K3.EQ.-1) THEN
                      C1= A312
                      C2= A312
                      C3=-A112-A212
                    ELSE IF(K3.EQ.-2) THEN
                      C1=-A212-A312
                      C2= A112
                      C3= A112
                    ELSE
                      C1= A212
                      C2=-A112
                      C3= 0.
                    END IF
                  END IF
                  CS=C1*S11+C2*S21+C3*S31
                  C =C1*S10+C2*S20+C3*S30
                  IF(ABS(C).GT.ABS(CS)) THEN
C
C                   Current edge inside the current circle.
                    NEW=ICIRC
C
                    A101=S10-S11
                    A201=S20-S21
                    A301=S30-S31
                    IF(K3.GT.0) THEN
C                     Circle centre inside the region:
                      A130=S13-S10
                      A230=S23-S20
                      A330=S33-S30
C                     New circle 0-1-2:
                      C1=A201*A312-A301*A212
                      C2=A301*A112-A101*A312
                      C3=A101*A212-A201*A112
                      CALL CIRCLE(NEW,NEDGE,K1,K2,S10,S20,S30,
     *                              C1,C2,C3,EDGE1,EDGE2,EDGE3,
     *                           MCIRC,NCIRC,KCIRC1,KCIRC2,KCIRC3,RCIRC)
C                     New circle 2-3-0:
                      C1=A223*A330-A323*A230
                      C2=A323*A130-A123*A330
                      C3=A123*A230-A223*A130
                      CALL CIRCLE(NEW,NEDGE,K2,K3,S10,S20,S30,
     *                              C1,C2,C3,EDGE1,EDGE2,EDGE3,
     *                           MCIRC,NCIRC,KCIRC1,KCIRC2,KCIRC3,RCIRC)
C                     New circle 3-0-1:
                      C1=A230*A301-A330*A201
                      C2=A330*A101-A130*A301
                      C3=A130*A201-A230*A101
                      CALL CIRCLE(NEW,NEDGE,K3,K1,S10,S20,S30,
     *                              C1,C2,C3,EDGE1,EDGE2,EDGE3,
     *                           MCIRC,NCIRC,KCIRC1,KCIRC2,KCIRC3,RCIRC)
                    ELSE
C                     Circle centre at a side of the region:
                      A120=S12-S10
                      A220=S22-S20
                      A320=S32-S30
C                     New circle 2-0-1:
                      C1=A220*A301-A320*A201
                      C2=A320*A101-A120*A301
                      C3=A120*A201-A220*A101
                      CALL CIRCLE(NEW,NEDGE,K1,K2,S10,S20,S30,
     *                              C1,C2,C3,EDGE1,EDGE2,EDGE3,
     *                           MCIRC,NCIRC,KCIRC1,KCIRC2,KCIRC3,RCIRC)
C                     New circle 0-1:
                      IF(K3.EQ.-1) THEN
                        C1= A301
                        C2= A301
                        C3=-A101-A201
                      ELSE IF(K3.EQ.-2) THEN
                        C1=-A201-A301
                        C2= A101
                        C3= A101
                      ELSE IF(K3.EQ.-3) THEN
                        C1= A201
                        C2=-A101
                        C3= 0.
                      END IF
                      CALL CIRCLE(NEW,NEDGE,K1,K3,S10,S20,S30,
     *                              C1,C2,C3,EDGE1,EDGE2,EDGE3,
     *                           MCIRC,NCIRC,KCIRC1,KCIRC2,KCIRC3,RCIRC)
C                     New circle 2-0:
                      IF(K3.EQ.-1) THEN
                        C1= A320
                        C2= A320
                        C3=-A120-A220
                      ELSE IF(K3.EQ.-2) THEN
                        C1=-A220-A320
                        C2= A120
                        C3= A120
                      ELSE
                        C1= A220
                        C2=-A120
                        C3= 0.
                      END IF
                      CALL CIRCLE(NEW,NEDGE,K2,K3,S10,S20,S30,
     *                              C1,C2,C3,EDGE1,EDGE2,EDGE3,
     *                           MCIRC,NCIRC,KCIRC1,KCIRC2,KCIRC3,RCIRC)
                    END IF
                    WRITE(*,'(A,I7,2I8)') '+',IH,NEDGE,NCIRC
C
                  END IF
                END IF
   49         CONTINUE
C
              WRITE(*,'(A,I7,2I8)') '+',IH,NEDGE,NCIRC
   56         CONTINUE
            END IF
   57     CONTINUE
   58   CONTINUE
        IF(IH.EQ.3) THEN
          NCIRC=4
          KCIRC1(1)=2
          KCIRC2(1)=3
          KCIRC3(1)=-1
          RCIRC(1)=0.0917517
          KCIRC1(2)=3
          KCIRC2(2)=1
          KCIRC3(2)=-2
          RCIRC(2)=0.2113249
          KCIRC1(3)=1
          KCIRC2(3)=2
          KCIRC3(3)=-3
          RCIRC(3)=0.1464466
          KCIRC1(4)=1
          KCIRC2(4)=2
          KCIRC3(4)=3
          RCIRC(4)=0.2142030
          RADIUS(1)=0.6666667
          RADIUS(2)=0.3333333
        END IF
        AUX=0.
        DO 68 ICIRC=1,NCIRC
          IF(KCIRC3(ICIRC).NE.0) THEN
            AUX=AMAX1(RCIRC(ICIRC),AUX)
          END IF
   68   CONTINUE
        LEDGE(IH)=NEDGE
        RADIUS(IH)=AUX
   69 CONTINUE
C
C.......................................................................
C
      DO 71 IH=NH,2,-1
        IF(RADIUS(IH).EQ.RADIUS(IH-1)) THEN
          RADIUS(IH)=0.
        END IF
   71 CONTINUE
C
      DO 72 IH=1,NH
        IF(RADIUS(IH).NE.0.) THEN
          AUX=SQRT(FLOAT(IH))
          I3=INT(AUX+0.00001)
          I2=IH-I3*I3
          RR=ASIN(SQRT(RADIUS(IH)))**2
          WRITE(LU1,'(I4,A,I2,A,I2,A,I2,F12.6,E13.6,F9.6)')
     *       IH,'=',I3,'*',I3,'+',I2,AUX,RR/2.,FLOAT(IH)*RR
        END IF
   72 CONTINUE
C
C.......................................................................
C
      WRITE(*,'(A,I7,2I8)') ' ',0,0
      J=INT(SQRT(FLOAT(NH-2)+DH))
      DO 99 J=1,J
        IH=J*J+2
        DO 81 IEDGE=LEDGE(IH),1,-1
          KEDGE(IEDGE)=1
   81   CONTINUE
        DO 82 IEDGE=LEDGE(IH),4,-1
          WRITE(*,'(A,I7)') '+',IEDGE
          CALL PICK(IEDGE,LEDGE(IH),KEDGE,EDGE1,EDGE2,EDGE3,
     *                                           NS,S1,S2,S3,RADIUS(IH))
          CALL SPOT(NS,S1,S2,S3,RADIUS(IH))
          KEDGE(IEDGE)=NS
   82   CONTINUE
C
        NEDGE=0
        DO 95 IEDGE=1,LEDGE(IH)
          IF(KEDGE(IEDGE).NE.0) THEN
            NEDGE=NEDGE+1
            KEDGE(NEDGE)=IEDGE
          END IF
   95   CONTINUE
        I3I3=1
        DO 98 IEDGE=1,NEDGE
          I=KEDGE(IEDGE)
   96     CONTINUE
            IF(LEDGE(I3I3).GE.I) GO TO 97
            I3I3=I3I3+1
          GO TO 96
   97     CONTINUE
          AUX=SQRT(FLOAT(I3I3))
          KS1(IEDGE)=INT(EDGE1(I)*AUX+0.5)
          KS2(IEDGE)=INT(EDGE2(I)*AUX+0.5)
          KS3(IEDGE)=INT(EDGE3(I)*AUX+0.5)
   98   CONTINUE
        WRITE(*,'(A,I7,2I8)') '+',IH,NEDGE
        RR=ASIN(SQRT(RADIUS(IH)))**2
C-old   WRITE(LU2,'(I3,I6,I7,E13.6)') J,NEDGE,IH,RR/2.
C-old   RR... SQUARE OF THE RADIUS OF THE LARGEST ANGULAR CIRCLE
C-old         CONTAINING NO EDGE OF THE FORWARD STAR.
C-old   IH... SQUARE OF THE RADIUS OF A FORWARD STAR IN GRID INTERVALS.
        WRITE(LU2,'(I3,I6,I7,E13.6)') J,NEDGE
        WRITE(LU2,'(8(3I3,1X))') (KS3(I),KS2(I),KS1(I),I=1,NEDGE)
   99 CONTINUE
C-old WRITE(LU2,'(I3,I6,I7,E13.6)') 0,0,0,0.
      WRITE(LU2,'(I3,I6,I7,E13.6)') 0,0
C
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE CIRCLE(NEW,K1,K2,K3,S1,S2,S3,C1,C2,C3,
     *      EDGE1,EDGE2,EDGE3,MCIRC,NCIRC,KCIRC1,KCIRC2,KCIRC3,RCIRC)
      INTEGER NEW,K1,K2,K3,MCIRC,NCIRC,KCIRC1(*),KCIRC2(*),KCIRC3(*)
      REAL S1,S2,S3,C1,C2,C3
      REAL EDGE1(*),EDGE2(*),EDGE3(*),RCIRC(*)
C
C     NEW...  If non-zero, free memory location for a circle.
C     K1,K2,K3... Edges of a new circle.
C     S1,S2,S3... Unit vector of an edge at the periphery of the new
C             circle.
C     C1,C2,C3... Axial vector of the new circle.
C     EDGE1(IS),EDGE2(IS),EDGE3(IS)... Unit vector of the IS-th
C             edge.
C     MCIRC.. Maximum number of circles.
C     NCIRC.. Number of circles stored.
C     KCIRC1(ICIRC),KCIRC2(ICIRC),KCIRC3(ICIRC)... Edges at the
C             periphery of a circle.
C
C-----------------------------------------------------------------------
C
      INTEGER IS,I
      REAL CC,CS,SS,C
C
C     IS...   Index of an edge (loop variable).
C     I...    Auxiliary starage location.
C     CS,C... Cosines of angular radii of circles.
C     CC,CS,SS... Scalar products C*C, C*S, S*S of the input vectors.
C
C.......................................................................
C
      IF(NEW.GT.0) THEN
        KCIRC3(NEW)=0
      END IF
      IF((C1+0.00001.GE.C2.AND.C2+0.00001.GE.C3.AND.C3+0.00001.GE.0.).OR
     *  .(C1-0.00001.LE.C2.AND.C2-0.00001.LE.C3.AND.C3-0.00001.LE.0.))
     *                                                              THEN
        CS=ABS(C1*S1+C2*S2+C3*S3)
        DO 21 IS=1,K1-1
          IF(IS.NE.K2.AND.IS.NE.K3) THEN
            C=C1*EDGE1(IS)+C2*EDGE2(IS)+C3*EDGE3(IS)
            IF(ABS(C).GT.CS) THEN
              GO TO 22
            END IF
          END IF
   21   CONTINUE
C         New circle is to be stored:
          IF(NEW.GT.0) THEN
            I=NEW
            NEW=0
          ELSE
            I=NCIRC+1
            NCIRC=I
            IF(NCIRC.GT.MCIRC) THEN
              STOP 'ERROR: TOO MANY CIRCLES.'
            END IF
          END IF
          KCIRC1(I)=K1
          KCIRC2(I)=K2
          KCIRC3(I)=K3
C         Given: periphery edge S and axial edge C.
          SS=S1*S1+S2*S2+S3*S3
          CS=C1*S1+C2*S2+C3*S3
          CC=C1*C1+C2*C2+C3*C3
          RCIRC(I)=1.-CS*CS/(CC*SS)
   22   CONTINUE
      END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE PICK(IEDGE,NEDGE,KEDGE,EDGE1,EDGE2,EDGE3,
     *                                                    NS,S1,S2,S3,R)
      INTEGER IEDGE,NEDGE,KEDGE(NEDGE),NS
      REAL EDGE1(NEDGE),EDGE2(NEDGE),EDGE3(NEDGE)
      REAL S1(NEDGE),S2(NEDGE),S3(NEDGE),R
C
C Subroutine designed to pick up edges within a circle of radius
C 2*SQRT(R) around the given edge.
C
C     IEDGE... Index of the given edge.
C     EDGE1(IS),EDGE2(IS),EDGE3(IS)... Unit vector of the IS-th
C             edge.
C     S1(IS),S2(IS),S3(IS)... Unit vector of the IS-th picked up edge.
C
C-----------------------------------------------------------------------
C
      INTEGER IS
      REAL C1,C2,C3,C,CR
C
C     IS1,IS2,IS3,IS...   Indices of edges (loop variables).
C     C... Cosine of angular radus of a circle.
C
C.......................................................................
C
      CR=SQRT(1.-4.*R)-0.00001
      NS=1
      S1(1)=EDGE1(IEDGE)
      S2(1)=EDGE2(IEDGE)
      S3(1)=EDGE3(IEDGE)
      C1=S1(1)
      C2=S2(1)
      C3=S3(1)
C
      DO 10 IS=1,NEDGE
        IF(KEDGE(IS).NE.0) THEN
          IF(IS.NE.IEDGE) THEN
            C=C1*EDGE1(IS)+C2*EDGE2(IS)+C3*EDGE3(IS)
            IF(ABS(C).GE.CR) THEN
              NS=NS+1
              S1(NS)=EDGE1(IS)
              S2(NS)=EDGE2(IS)
              S3(NS)=EDGE3(IS)
            END IF
          END IF
        END IF
   10 CONTINUE
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE SPOT(NS,S1,S2,S3,R)
      INTEGER NS
      REAL S1(NS),S2(NS),S3(NS),R
C
C Subroutine designed to look for a circle of radius .GT. SQRT(R),
C containing the first given edge and no other.
C
C     S1(IS),S2(IS),S3(IS)... Unit vector of the IS-th
C             edge.
C
C Output:
C     NS=0 if a circle has been found.
C
C-----------------------------------------------------------------------
C
      INTEGER IS1,IS2,IS3,IS
      REAL S11,S12,S13,A112,A123,C1
      REAL S21,S22,S23,A212,A223,C2
      REAL S31,S32,S33,A312,A323,C3,CS,C,CR
C
C     IS1,IS2,IS3,IS...   Indices of edges (loop variables).
C     CS,C... Cosines of angular radii of circles.
C
C.......................................................................
C
      CR=SQRT(1.-R)-0.00001
      S10=S1(1)
      S20=S2(1)
      S30=S3(1)
C
      DO 33 IS3=2,NS
        S13=S1(IS3)
        S23=S2(IS3)
        S33=S3(IS3)
        DO 32 IS2=2,IS3-1
          S12=S1(IS2)
          S22=S2(IS2)
          S32=S3(IS2)
          A123=S12-S13
          A223=S22-S23
          A323=S32-S33
          DO 31 IS1=-1,IS2-1
            IF(IS1.GE.2) THEN
              S11=S1(IS1)
              S21=S2(IS1)
              S31=S3(IS1)
              A112=S11-S12
              A212=S21-S22
              A312=S31-S32
C             Circle axis:
              C1=A212*A323-A312*A223
              C2=A312*A123-A112*A323
              C3=A112*A223-A212*A123
            ELSE IF(IS1.EQ.1) THEN
              C1= A323
              C2= A323
              C3=-A123-A223
            ELSE IF(IS1.EQ.0) THEN
              C1=-A223-A323
              C2= A123
              C3= A123
            ELSE
              C1= A223
              C2=-A123
              C3= 0.
            END IF
            C =SQRT(C1*C1+C2*C2+C3*C3)
            C1=C1/C
            C2=C2/C
            C3=C3/C
            CS=ABS(C1*S13+C2*S23+C3*S33)
            IF(CS.LT.CR) THEN
              C =C1*S10+C2*S20+C3*S30
              IF(ABS(C).GT.CS) THEN
                DO 21 IS=2,NS
                  IF(IS.NE.IS1.AND.IS.NE.IS2.AND.IS.NE.IS3) THEN
                    C=C1*S1(IS)+C2*S2(IS)+C3*S3(IS)
                    IF(ABS(C).GT.CS) THEN
                      GO TO 22
                    END IF
                  END IF
   21           CONTINUE
C                 There is a large circle:
                  RETURN
   22           CONTINUE
              END IF
            END IF
   31     CONTINUE
   32   CONTINUE
   33 CONTINUE
C
      NS=0
      RETURN
      END
C
C=======================================================================
C
C