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