C PROGRAM TO DETERMINE, FOR A GIVEN RADIUS H=SQRT(IH) OF A FORWARD STAR,
C THE LARGEST ANGULAR CIRCLE NOT CONTAINING A SPINE.
C
C USAGE IF COMPILED WITH MS-FORTRAN V.5:
C     FS TABLE STARS
C
C INPUT FROM CONSOLE:
C     MAXIMUM VALUE OF H*H, SEE BELOW.
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 SPINE
C             OF THE FORWARD STAR.
C JUST THE FORWARD STARS WITH SMALLER R'S THAN THE PRECEDING ONES ARE
C LISTED.
C
C OUTPUT FILE STARS:
C (1) FOR EACH FORWARD STAR (1A) AND (1B):
C (1A) NLE, NSPINE, INT(H*H), R*R/2
C     NLE...  LEVEL OF A FORWARD STAR.
C     NSPINE... NUMBER OF SPINES.
C     H...    RADIUS OF A FORWARD STAR IN GRID INTERVALS.
C     R...    RADIUS OF THE LARGEST ANGULAR CIRCLE CONTAINING NO SPINE
C             OF THE FORWARD STAR.
C (1B) (I1(I),I2(I),I3(I),I=1,NSPINE)
C     LIST OF SPINES.  JUST SPINES WITH 0.LE.I1.LE.I2.LE.I3 ARE LISTED.
C (2) 0, 0, 0., 0
C
C DATE: 1994, FEBRUARY 9
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,MSPINE,MCIRC
      PARAMETER (MLE=80,MH=MLE**2+2,MSPINE=(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,NSPINE
      INTEGER KS1(MSPINE)
      INTEGER KS2(MSPINE)
      INTEGER KS3(MSPINE),KSPINE(MSPINE),LSPINE(MH),ISPINE,NS
      REAL SPINE1(MSPINE),S10,S11,S12,S13,A101,A112,A123,A130,C1
      REAL SPINE2(MSPINE),S20,S21,S22,S23,A201,A212,A223,A230,C2
      REAL SPINE3(MSPINE),S30,S31,S32,S33,A301,A312,A323,A330,C3
      REAL RCIRC(MCIRC),RADIUS(MH),DH,CS,C,RR,AUX
      REAL S1(MSPINE)
      REAL S2(MSPINE)
      REAL S3(MSPINE)
      EQUIVALENCE (KCIRC1,KS1),(KCIRC2,KS2),(KCIRC3,KS3)
      EQUIVALENCE (KCIRC1,S1),(KCIRC2,S2),(KCIRC3,S3),(RCIRC,KSPINE)
C
C     NCIRC.. NUMBER OF CIRCLES STORED.
C     ICIRC.. INDEX OF A CIRCLE.
C     KCIRC1(ICIRC),KCIRC2(ICIRC),KCIRC3(ICIRC)... SPINES 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 SPINES 2=(1,1,0) AND 3=(1,1,1).
C             KCIRC3(ICIRC)=-2: CIRCLE CENTRE AT THE SIDE OF THE REGION,
C               BETWEEN SPINES 3=(1,1,1) AND 2=(1,1,0).
C             KCIRC3(ICIRC)=-3: CIRCLE CENTRE AT THE SIDE OF THE REGION,
C               BETWEEN SPINES 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 A SPINE.
C     IH...   SQUARE OF THE LENGTH OF A SPINE.
C     I3I3... LIMIT FOR I3*I3.
C     I...    LIMIT FOR I2*I2+I3*I3.
C     J...    LOOP VARIABLE.
C     I1,I2,I3... COMPONENTS OF A SPINE.
C     NSPINE... NUMBER OF SPINES STORED.
C     KS1,KS2,KS3... SPINES OF THE RESULTING FORWARD STAR (VECTORS).
C     KSPINE... LIST OF SPINES OF THE RESULTING FORWARD STAR.
C     LSPINE(IH)... INDEX OF THE LAST SPINE OF SQUARE LENGTH IH.
C     ISPINE... INDEX OF A CURRENT SPINE WHEN DECIMATING A FORWARD STAR.
C     NS...   NUMBER OF SPINES IN VICINITY OF A CURRENT SPINE, WHEN
C             DECIMATING A FORWARD STAR.
C     S1,S2,S3... UNIT VECTORS OF SPINES IN VICINITY OF A CURRENT SPINE,
C             WHEN DECIMATING A FORWARD STAR.
C     SPINE1(IS),SPINE2(IS),SPINE3(IS)... UNIT VECTOR OF THE IS-TH
C             SPINE.
C     S10,S20,S30... UNIT VECTOR OF THE CURRENT SPINE.
C     S11,S21,S31, S12,S22,S32, S13,S23,S33... UNIT VECTORS OF THE
C             SPINES 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 SPINE 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  SPINES CIRCLES'
      WRITE(*,'(3I8,A)') NH,MSPINE,MCIRC,' MAXIMUM'
      WRITE(*,'(A)')
C
      DH=1./SQRT(12.*FLOAT(NH))
      DH=DH/2.
      NSPINE=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 SPINE
              NSPINE=NSPINE+1
              IF(NSPINE.GT.MSPINE) THEN
                STOP 'ERROR: TOO MANY SPINES.'
              END IF
              AUX=SQRT(FLOAT(IH))
              S10=FLOAT(I1)/AUX
              S20=FLOAT(I2)/AUX
              S30=FLOAT(I3)/AUX
              SPINE1(NSPINE)=S10
              SPINE2(NSPINE)=S20
              SPINE3(NSPINE)=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=SPINE1(K1)
                  S21=SPINE2(K1)
                  S31=SPINE3(K1)
                  S12=SPINE1(K2)
                  S22=SPINE2(K2)
                  S32=SPINE3(K2)
                  A112=S11-S12
                  A212=S21-S22
                  A312=S31-S32
                  IF(K3.GT.0) THEN
C                   CIRCLE CENTRE INSIDE THE REGION:
                    S13=SPINE1(K3)
                    S23=SPINE2(K3)
                    S33=SPINE3(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 SPINE 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,NSPINE,K1,K2,S10,S20,S30,
     *                              C1,C2,C3,SPINE1,SPINE2,SPINE3,
     *                           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,NSPINE,K2,K3,S10,S20,S30,
     *                              C1,C2,C3,SPINE1,SPINE2,SPINE3,
     *                           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,NSPINE,K3,K1,S10,S20,S30,
     *                              C1,C2,C3,SPINE1,SPINE2,SPINE3,
     *                           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,NSPINE,K1,K2,S10,S20,S30,
     *                              C1,C2,C3,SPINE1,SPINE2,SPINE3,
     *                           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,NSPINE,K1,K3,S10,S20,S30,
     *                              C1,C2,C3,SPINE1,SPINE2,SPINE3,
     *                           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,NSPINE,K2,K3,S10,S20,S30,
     *                              C1,C2,C3,SPINE1,SPINE2,SPINE3,
     *                           MCIRC,NCIRC,KCIRC1,KCIRC2,KCIRC3,RCIRC)
                    END IF
                    WRITE(*,'(A,I7,2I8)') '+',IH,NSPINE,NCIRC
C
                  END IF
                END IF
   49         CONTINUE
C
              WRITE(*,'(A,I7,2I8)') '+',IH,NSPINE,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
        LSPINE(IH)=NSPINE
        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 ISPINE=LSPINE(IH),1,-1
          KSPINE(ISPINE)=1
   81   CONTINUE
        DO 82 ISPINE=LSPINE(IH),4,-1
          WRITE(*,'(A,I3)') '+',ISPINE
          CALL PICK(ISPINE,LSPINE(IH),KSPINE,SPINE1,SPINE2,SPINE3,
     *                                           NS,S1,S2,S3,RADIUS(IH))
          CALL SPOT(NS,S1,S2,S3,RADIUS(IH))
          KSPINE(ISPINE)=NS
   82   CONTINUE
C
        NSPINE=0
        DO 95 ISPINE=1,LSPINE(IH)
          IF(KSPINE(ISPINE).NE.0) THEN
            NSPINE=NSPINE+1
            KSPINE(NSPINE)=ISPINE
          END IF
   95   CONTINUE
        I3I3=1
        DO 98 ISPINE=1,NSPINE
          I=KSPINE(ISPINE)
   96     CONTINUE
            IF(LSPINE(I3I3).GE.I) GO TO 97
            I3I3=I3I3+1
          GO TO 96
   97     CONTINUE
          AUX=SQRT(FLOAT(I3I3))
          KS1(ISPINE)=INT(SPINE1(I)*AUX+0.5)
          KS2(ISPINE)=INT(SPINE2(I)*AUX+0.5)
          KS3(ISPINE)=INT(SPINE3(I)*AUX+0.5)
   98   CONTINUE
        WRITE(*,'(A,I7,2I8)') '+',IH,NSPINE
        RR=ASIN(SQRT(RADIUS(IH)))**2
        WRITE(LU2,'(I3,I6,I7,E13.6)') J,NSPINE,IH,RR/2.
        WRITE(LU2,'(8(3I3,1X))') (KS3(I),KS2(I),KS1(I),I=1,NSPINE)
   99 CONTINUE
      WRITE(LU2,'(I3,I6,I7,E13.6)') 0,0,0,0.
C
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE CIRCLE(NEW,K1,K2,K3,S1,S2,S3,C1,C2,C3,
     *      SPINE1,SPINE2,SPINE3,MCIRC,NCIRC,KCIRC1,KCIRC2,KCIRC3,RCIRC)
      INTEGER NEW,K1,K2,K3,MCIRC,NCIRC,KCIRC1(*),KCIRC2(*),KCIRC3(*)
      REAL S1,S2,S3,C1,C2,C3
      REAL SPINE1(*),SPINE2(*),SPINE3(*),RCIRC(*)
C
C     NEW...  IF NON=ZERO, FREE MEMORY LOCATION FOR A CIRCLE.
C     K1,K2,K3... SPINES OF A NEW CIRCLE.
C     S1,S2,S3... UNIT VECTOR OF A SPINE AT THE PERIPHERY OF THE NEW
C             CIRCLE.
C     C1,C2,C3... AXIAL VECTOR OF THE NEW CIRCLE.
C     SPINE1(IS),SPINE2(IS),SPINE3(IS)... UNIT VECTOR OF THE IS-TH
C             SPINE.
C     MCIRC.. MAXIMUM NUMBER OF CIRCLES.
C     NCIRC.. NUMBER OF CIRCLES STORED.
C     KCIRC1(ICIRC),KCIRC2(ICIRC),KCIRC3(ICIRC)... SPINES AT THE
C             PERIPHERY OF A CIRCLE.
C
C-----------------------------------------------------------------------
C
      INTEGER IS,I
      REAL CC,CS,SS,C
C
C     IS...   INDEX OF A SPINE (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*SPINE1(IS)+C2*SPINE2(IS)+C3*SPINE3(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 SPINE S AND AXIAL SPINE 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(ISPINE,NSPINE,KSPINE,SPINE1,SPINE2,SPINE3,
     *                                                    NS,S1,S2,S3,R)
      INTEGER ISPINE,NSPINE,KSPINE(NSPINE),NS
      REAL SPINE1(NSPINE),SPINE2(NSPINE),SPINE3(NSPINE)
      REAL S1(NS),S2(NS),S3(NS),R
C
C SUBROUTINE DESIGNED TO PICK UP SPINES WITHIN A CIRCLE OF RADIUS
C 2*SQRT(R) AROUND THE GIVEN SPINE.
C
C     ISPINE... INDEX OF THE GIVEN SPINE.
C     SPINE1(IS),SPINE2(IS),SPINE3(IS)... UNIT VECTOR OF THE IS-TH
C             SPINE.
C     S1(IS),S2(IS),S3(IS)... UNIT VECTOR OF THE IS-TH PICKED UP SPINE.
C
C-----------------------------------------------------------------------
C
      INTEGER IS
      REAL C1,C2,C3,C,CR
C
C     IS1,IS2,IS3,IS...   INDICES OF SPINES (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)=SPINE1(ISPINE)
      S2(1)=SPINE2(ISPINE)
      S3(1)=SPINE3(ISPINE)
      C1=S1(1)
      C2=S2(1)
      C3=S3(1)
C
      DO 10 IS=1,NSPINE
        IF(KSPINE(IS).NE.0) THEN
          IF(IS.NE.ISPINE) THEN
            C=C1*SPINE1(IS)+C2*SPINE2(IS)+C3*SPINE3(IS)
            IF(ABS(C).GE.CR) THEN
              NS=NS+1
              S1(NS)=SPINE1(IS)
              S2(NS)=SPINE2(IS)
              S3(NS)=SPINE3(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 SPINE AND NO OTHER.
C
C     S1(IS),S2(IS),S3(IS)... UNIT VECTOR OF THE IS-TH
C             SPINE.
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 SPINES (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
