C
C C Filenames: CHARACTER*80 FILE2 C C Logical unit numbers: INTEGER LU2 PARAMETER (LU2=2) C C Array dimensions: INTEGER MLE,MH,MSPINE PARAMETER (MLE=400,MH=MLE**2+1,MSPINE=MLE*(MLE+1)/2) C INTEGER ISP1(MSPINE),ISP2(MSPINE),ISP3(MSPINE) INTEGER IS1(MSPINE),IS2(MSPINE),IS3(MSPINE) C C....................................................................... C NH=100 FILE2='net.fs' WRITE(*,'(A)') *'+Enter maximum f.s. size (100), and f.s. filename (''net.fs''): ' READ(*,*) NH,FILE2 NH=NH*NH+1 IF(NH.GT.MH) THEN STOP 'ERROR: TOO LARGE FORWARD STAR.' END IF OPEN(LU2,FILE=FILE2) C DH=1./SQRT(12.*FLOAT(NH)) DH=DH/2. NFS=1 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 IF(I3I3.EQ.0) THEN I3=0 DO 10 J=2,I1 IF(MOD(I1,J).EQ.0.AND.MOD(I2,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 ISP1(NSPINE)=I3 ISP2(NSPINE)=I2 ISP3(NSPINE)=I1 C 56 CONTINUE END IF 57 CONTINUE 58 CONTINUE C IF(IH.EQ.NFS*NFS+1) THEN DO 66 ISPINE=NSPINE,3,-1 C Looking for two closest spines II=ISP2(ISPINE)*ISP2(ISPINE)+ISP3(ISPINE)*ISP3(ISPINE) JL=2 JLL=2 AL=999999. JR=1 JRR=1 AR=999999. DO 64 JSPINE=1,NSPINE IF(ISP3(JSPINE).GT.0) THEN IF(JSPINE.NE.ISPINE) THEN IJ=ISP2(ISPINE)*ISP3(JSPINE)-ISP3(ISPINE)*ISP2(JSPINE) JJ=ISP2(JSPINE)*ISP2(JSPINE)+ISP3(JSPINE)*ISP3(JSPINE) A=FLOAT(IJ*IJ)/FLOAT(II*JJ) IF(IJ.LT.0) THEN IF(A.LT.AL) THEN JL=JSPINE AL=A JLR=ISP2(JL)*ISP3(JR)-ISP3(JL)*ISP2(JR) JLL=ISP2(JL)*ISP2(JL)+ISP3(JL)*ISP3(JL) R=FLOAT(JLR*JLR*IH)/FLOAT(JLL*JRR) IF(R.LT.0.999999) THEN ISP3(ISPINE)=-ISP3(ISPINE) GO TO 65 END IF END IF ELSE IF(A.LT.AR) THEN JR=JSPINE AR=A JLR=ISP2(JL)*ISP3(JR)-ISP3(JL)*ISP2(JR) JRR=ISP2(JR)*ISP2(JR)+ISP3(JR)*ISP3(JR) R=FLOAT(JLR*JLR*IH)/FLOAT(JLL*JRR) IF(R.LT.0.999999) THEN ISP3(ISPINE)=-ISP3(ISPINE) GO TO 65 END IF END IF END IF END IF END IF 64 CONTINUE 65 CONTINUE 66 CONTINUE NS=0 DO 67 ISPINE=1,NSPINE IF(ISP3(ISPINE).GT.0) THEN NS=NS+1 IS1(NS)=ISP1(ISPINE) IS2(NS)=ISP2(ISPINE) IS3(NS)=ISP3(ISPINE) END IF 67 CONTINUE WRITE(*,'(A,I7,2I8)') '+',NFS,NS WRITE(LU2,'(I3,I6,I7,E13.6)') NFS,NS WRITE(LU2,'(8(3I3,1X))') (IS1(I),IS2(I),IS3(I),I=1,NS) NFS=NFS+1 DO 68 ISPINE=1,NSPINE ISP3(ISPINE)=IABS(ISP3(ISPINE)) 68 CONTINUE END IF 69 CONTINUE WRITE(LU2,'(I3,I6,I7,E13.6)') 0,0 C STOP END C C======================================================================= C C