initial commit

adding files from https://space.mit.edu/home//gsg/doc/Voyager/src/nrf/
with additional files removed
189 .f source code files
189 .dem demonstration files
10 .dat data files
2 .doc readme files

Author: jonlighthall

adi.dem

@@ -0,0 +1,61 @@
+	PROGRAM D17R2
+C	Driver for routine ADI
+	IMPLICIT REAL*8(A-H,O-Z)
+	PARAMETER(JMAX=11,PI=3.1415926)
+	DIMENSION A(JMAX,JMAX),B(JMAX,JMAX),C(JMAX,JMAX),D(JMAX,JMAX),
+     *		E(JMAX,JMAX),F(JMAX,JMAX),G(JMAX,JMAX),U(JMAX,JMAX)
+	DO 12 I=1,JMAX
+		DO 11 J=1,JMAX
+			A(I,J)=-1.0
+			B(I,J)=2.0
+			C(I,J)=-1.0
+			D(I,J)=-1.0
+			E(I,J)=2.0
+			F(I,J)=-1.0
+			G(I,J)=0.0
+			U(I,J)=0.0
+11		CONTINUE
+12	CONTINUE
+	MID=JMAX/2+1
+	G(MID,MID)=2.0
+	ALPHA=2.0*(1.0-COS(PI/JMAX))
+	BETA=2.0*(1.0-COS((JMAX-1)*PI/JMAX))
+	ALIM=LOG(4.0*JMAX/PI)
+	K=0
+1	K=K+1
+	IF (2**K .LT. ALIM) GOTO 1
+	EPS=1.0E-4
+	CALL ADI(A,B,C,D,E,F,G,U,JMAX,K,ALPHA,BETA,EPS)
+	WRITE(*, '(1X,A)') 'ADI Solution:'
+	DO 13 I=1,JMAX
+		WRITE(*,'(1X,11F7.2)') (U(I,J),J=1,JMAX)
+13	CONTINUE
+	WRITE(*,'(/1X,A)') 'Test that solution satisfies Difference Eqns:'
+	DO 15 I=2,JMAX-1
+		DO 14 J=2,JMAX-1
+			G(I,J)=-4.0*U(I,J)+U(I+1,J)
+     *				+U(I-1,J)+U(I,J-1)+U(I,J+1)
+14		CONTINUE
+		WRITE(*,'(8X,9F7.2)') (G(I,J),J=2,JMAX-1)
+15	CONTINUE
+	END
+
+	SUBROUTINE TRIDAG(A,B,C,R,U,N)
+C	This is a double precision version for use with ADI
+	IMPLICIT REAL*8(A-H,O-Z)
+	PARAMETER (NMAX=100)
+	DIMENSION GAM(NMAX),A(N),B(N),C(N),R(N),U(N)
+	IF(B(1).EQ.0.)PAUSE
+	BET=B(1)
+	U(1)=R(1)/BET
+	DO 11 J=2,N
+		GAM(J)=C(J-1)/BET
+		BET=B(J)-A(J)*GAM(J)
+		IF(BET.EQ.0.)PAUSE
+		U(J)=(R(J)-A(J)*U(J-1))/BET
+11	CONTINUE
+	DO 12 J=N-1,1,-1
+		U(J)=U(J)-GAM(J+1)*U(J+1)
+12	CONTINUE
+	RETURN
+	END

adi.f

@@ -0,0 +1,86 @@
+      SUBROUTINE ADI(A,B,C,D,E,F,G,U,JMAX,K,ALPHA,BETA,EPS)
+      IMPLICIT REAL*8(A-H,O-Z)
+      PARAMETER(JJ=100,KK=6,NRR=32,MAXITS=100,ZERO=0.D0,TWO=2.D0,HALF=.5
+     *D0)
+      DIMENSION A(JMAX,JMAX),B(JMAX,JMAX),C(JMAX,JMAX),D(JMAX,JMAX),
+     *    E(JMAX,JMAX),F(JMAX,JMAX),G(JMAX,JMAX),U(JMAX,JMAX),
+     *    AA(JJ),BB(JJ),CC(JJ),RR(JJ),UU(JJ),PSI(JJ,JJ),
+     *    ALPH(KK),BET(KK),R(NRR),S(NRR,KK)
+      IF(JMAX.GT.JJ)PAUSE 'Increase JJ'
+      IF(K.GT.KK-1)PAUSE 'Increase KK'
+      K1=K+1
+      NR=2**K
+      ALPH(1)=ALPHA
+      BET(1)=BETA
+      DO 11 J=1,K
+        ALPH(J+1)=SQRT(ALPH(J)*BET(J))
+        BET(J+1)=HALF*(ALPH(J)+BET(J))
+11    CONTINUE
+      S(1,1)=SQRT(ALPH(K1)*BET(K1))
+      DO 13 J=1,K
+        AB=ALPH(K1-J)*BET(K1-J)
+        DO 12 N=1,2**(J-1)
+          DISC=SQRT(S(N,J)**2-AB)
+          S(2*N,J+1)=S(N,J)+DISC
+          S(2*N-1,J+1)=AB/S(2*N,J+1)
+12      CONTINUE
+13    CONTINUE
+      DO 14 N=1,NR
+        R(N)=S(N,K1)
+14    CONTINUE
+      ANORMG=ZERO
+      DO 16 J=2,JMAX-1
+        DO 15 L=2,JMAX-1
+          ANORMG=ANORMG+ABS(G(J,L))
+          PSI(J,L)=-D(J,L)*U(J,L-1)+(R(1)-E(J,L))*U(J,L)
+     *        -F(J,L)*U(J,L+1)
+15      CONTINUE
+16    CONTINUE
+      NITS=MAXITS/NR
+      DO 27 KITS=1,NITS
+        DO 24 N=1,NR
+          IF(N.EQ.NR)THEN
+            NEXT=1
+          ELSE
+            NEXT=N+1
+          ENDIF
+          RFACT=R(N)+R(NEXT)
+          DO 19 L=2,JMAX-1
+            DO 17 J=2,JMAX-1
+              AA(J-1)=A(J,L)
+              BB(J-1)=B(J,L)+R(N)
+              CC(J-1)=C(J,L)
+              RR(J-1)=PSI(J,L)-G(J,L)
+17          CONTINUE
+            CALL TRIDAG(AA,BB,CC,RR,UU,JMAX-2)
+            DO 18 J=2,JMAX-1
+              PSI(J,L)=-PSI(J,L)+TWO*R(N)*UU(J-1)
+18          CONTINUE
+19        CONTINUE
+          DO 23 J=2,JMAX-1
+            DO 21 L=2,JMAX-1
+              AA(L-1)=D(J,L)
+              BB(L-1)=E(J,L)+R(N)
+              CC(L-1)=F(J,L)
+              RR(L-1)=PSI(J,L)
+21          CONTINUE
+            CALL TRIDAG(AA,BB,CC,RR,UU,JMAX-2)
+            DO 22 L=2,JMAX-1
+              U(J,L)=UU(L-1)
+              PSI(J,L)=-PSI(J,L)+RFACT*UU(L-1)
+22          CONTINUE
+23        CONTINUE
+24      CONTINUE
+        ANORM=ZERO
+        DO 26 J=2,JMAX-1
+          DO 25 L=2,JMAX-1
+            RESID=A(J,L)*U(J-1,L)+(B(J,L)+E(J,L))*U(J,L)
+     *          +C(J,L)*U(J+1,L)+D(J,L)*U(J,L-1)
+     *          +F(J,L)*U(J,L+1)+G(J,L)
+            ANORM=ANORM+ABS(RESID)
+25        CONTINUE
+26      CONTINUE
+        IF(ANORM.LT.EPS*ANORMG)RETURN
+27    CONTINUE
+      PAUSE 'MAXITS exceeded'
+      END

amoeba.dem

@@ -0,0 +1,29 @@
+	PROGRAM D10R5
+C	Driver for routine AMOEBA
+	EXTERNAL FAMOEB
+	PARAMETER(NP=3,MP=4,FTOL=1.0E-6)
+	DIMENSION P(MP,NP),X(NP),Y(MP)
+	DATA P/0.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,1.0/
+	NDIM=NP
+	DO 12 I=1,MP
+		DO 11 J=1,NP
+			X(J)=P(I,J)
+11		CONTINUE
+		Y(I)=FAMOEB(X)
+12	CONTINUE
+	CALL AMOEBA(P,Y,MP,NP,NDIM,FTOL,FAMOEB,ITER)
+	WRITE(*,'(/1X,A,I3)') 'Iterations: ',ITER
+	WRITE(*,'(/1X,A)') 'Vertices of final 3-D simplex and'
+	WRITE(*,'(1X,A)') 'function values at the vertices:'
+	WRITE(*,'(/3X,A,T11,A,T23,A,T35,A,T45,A/)') 'I',
+     *		'X(I)','Y(I)','Z(I)','FUNCTION'
+	DO 13 I=1,MP
+		WRITE(*,'(1X,I3,4F12.6)') I,(P(I,J),J=1,NP),Y(I)
+13	CONTINUE
+	WRITE(*,'(/1X,A)') 'True minimum is at (0.5,0.6,0.7)'
+	END
+
+	FUNCTION FAMOEB(X)
+	DIMENSION X(3)
+	FAMOEB=0.6-BESSJ0((X(1)-0.5)**2+(X(2)-0.6)**2+(X(3)-0.7)**2)
+	END

amoeba.f

@@ -0,0 +1,92 @@
+      SUBROUTINE AMOEBA(P,Y,MP,NP,NDIM,FTOL,FUNK,ITER)
+      PARAMETER (NMAX=20,ALPHA=1.0,BETA=0.5,GAMMA=2.0,ITMAX=500)
+      DIMENSION P(MP,NP),Y(MP),PR(NMAX),PRR(NMAX),PBAR(NMAX)
+      MPTS=NDIM+1
+      ITER=0
+1     ILO=1
+      IF(Y(1).GT.Y(2))THEN
+        IHI=1
+        INHI=2
+      ELSE
+        IHI=2
+        INHI=1
+      ENDIF
+      DO 11 I=1,MPTS
+        IF(Y(I).LT.Y(ILO)) ILO=I
+        IF(Y(I).GT.Y(IHI))THEN
+          INHI=IHI
+          IHI=I
+        ELSE IF(Y(I).GT.Y(INHI))THEN
+          IF(I.NE.IHI) INHI=I
+        ENDIF
+11    CONTINUE
+      RTOL=2.*ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))
+      IF(RTOL.LT.FTOL)RETURN
+      IF(ITER.EQ.ITMAX) PAUSE 'Amoeba exceeding maximum iterations.'
+      ITER=ITER+1
+      DO 12 J=1,NDIM
+        PBAR(J)=0.
+12    CONTINUE
+      DO 14 I=1,MPTS
+        IF(I.NE.IHI)THEN
+          DO 13 J=1,NDIM
+            PBAR(J)=PBAR(J)+P(I,J)
+13        CONTINUE
+        ENDIF
+14    CONTINUE
+      DO 15 J=1,NDIM
+        PBAR(J)=PBAR(J)/NDIM
+        PR(J)=(1.+ALPHA)*PBAR(J)-ALPHA*P(IHI,J)
+15    CONTINUE
+      YPR=FUNK(PR)
+      IF(YPR.LE.Y(ILO))THEN
+        DO 16 J=1,NDIM
+          PRR(J)=GAMMA*PR(J)+(1.-GAMMA)*PBAR(J)
+16      CONTINUE
+        YPRR=FUNK(PRR)
+        IF(YPRR.LT.Y(ILO))THEN
+          DO 17 J=1,NDIM
+            P(IHI,J)=PRR(J)
+17        CONTINUE
+          Y(IHI)=YPRR
+        ELSE
+          DO 18 J=1,NDIM
+            P(IHI,J)=PR(J)
+18        CONTINUE
+          Y(IHI)=YPR
+        ENDIF
+      ELSE IF(YPR.GE.Y(INHI))THEN
+        IF(YPR.LT.Y(IHI))THEN
+          DO 19 J=1,NDIM
+            P(IHI,J)=PR(J)
+19        CONTINUE
+          Y(IHI)=YPR
+        ENDIF
+        DO 21 J=1,NDIM
+          PRR(J)=BETA*P(IHI,J)+(1.-BETA)*PBAR(J)
+21      CONTINUE
+        YPRR=FUNK(PRR)
+        IF(YPRR.LT.Y(IHI))THEN
+          DO 22 J=1,NDIM
+            P(IHI,J)=PRR(J)
+22        CONTINUE
+          Y(IHI)=YPRR
+        ELSE
+          DO 24 I=1,MPTS
+            IF(I.NE.ILO)THEN
+              DO 23 J=1,NDIM
+                PR(J)=0.5*(P(I,J)+P(ILO,J))
+                P(I,J)=PR(J)
+23            CONTINUE
+              Y(I)=FUNK(PR)
+            ENDIF
+24        CONTINUE
+        ENDIF
+      ELSE
+        DO 25 J=1,NDIM
+          P(IHI,J)=PR(J)
+25      CONTINUE
+        Y(IHI)=YPR
+      ENDIF
+      GO TO 1
+      END

anneal.dem

@@ -0,0 +1,19 @@
+	PROGRAM D10R13
+	PARAMETER (NCITY=10)
+	DIMENSION X(NCITY),Y(NCITY),IORDER(NCITY)
+C	Create points of sale
+	IDUM=-111
+	DO 11 I=1,NCITY
+		X(I)=RAN3(IDUM)
+		Y(I)=RAN3(IDUM)
+		IORDER(I)=I
+11	CONTINUE
+	CALL ANNEAL(X,Y,IORDER,NCITY)
+	WRITE(*,*) '*** System Frozen ***'
+	WRITE(*,*) 'Final path:'
+	WRITE(*,'(1X,T3,A,T13,A,T23,A)') 'city','x','y'
+	DO 12 I=1,NCITY
+		II=IORDER(I)
+		WRITE(*,'(1X,I4,2F10.4)') II,X(II),Y(II)
+12	CONTINUE
+	END

anneal.f

@@ -0,0 +1,59 @@
+      SUBROUTINE ANNEAL(X,Y,IORDER,NCITY)
+C Rev. 12/13/85
+      DIMENSION X(NCITY),Y(NCITY),IORDER(NCITY),N(6)
+      LOGICAL ANS
+      ALEN(X1,X2,Y1,Y2)=SQRT((X2-X1)**2+(Y2-Y1)**2)
+      NOVER=100*NCITY
+      NLIMIT=10*NCITY
+      TFACTR=0.9
+      PATH=0.0
+      T=0.5
+      DO 11 I=1,NCITY-1
+          I1=IORDER(I)
+          I2=IORDER(I+1)
+          PATH=PATH+ALEN(X(I1),X(I2),Y(I1),Y(I2))
+11    CONTINUE
+      I1=IORDER(NCITY)
+      I2=IORDER(1)
+      PATH=PATH+ALEN(X(I1),X(I2),Y(I1),Y(I2))
+      IDUM=-1
+      ISEED=111
+      DO 15 J=1,100
+          NSUCC=0
+          DO 13 K=1,NOVER
+12            N(1)=1+INT(NCITY*RAN3(IDUM))
+              N(2)=1+INT((NCITY-1)*RAN3(IDUM))
+              IF (N(2).GE.N(1)) N(2)=N(2)+1
+              IDEC=IRBIT1(ISEED)
+              NN=1+MOD((N(1)-N(2)+NCITY-1),NCITY)
+              IF (NN.LT.3) GOTO 12
+              IF (IDEC.EQ.0) THEN
+                  N(3)=N(2)+INT(ABS(NN-2)*RAN3(IDUM))+1
+                  N(3)=1+MOD(N(3)-1,NCITY)
+                  CALL TRNCST(X,Y,IORDER,NCITY,N,DE)
+                  CALL METROP(DE,T,ANS)
+                  IF (ANS) THEN
+                      NSUCC=NSUCC+1
+                      PATH=PATH+DE
+                      CALL TRNSPT(IORDER,NCITY,N)
+                  ENDIF
+              ELSE
+                  CALL REVCST(X,Y,IORDER,NCITY,N,DE)
+                  CALL METROP(DE,T,ANS)
+                  IF (ANS) THEN
+                      NSUCC=NSUCC+1
+                      PATH=PATH+DE
+                      CALL REVERS(IORDER,NCITY,N)
+                  ENDIF
+              ENDIF
+              IF (NSUCC.GE.NLIMIT) GOTO 14
+13        CONTINUE
+14        WRITE(*,*)
+          WRITE(*,'(1X,A,F10.6,A,F10.6)') 'T =',T,
+     *            '   Path Length =',PATH
+          WRITE(*,'(1X,A,I6)') 'Successful Moves: ',NSUCC
+          T=T*TFACTR
+          IF (NSUCC.EQ.0) RETURN
+15    CONTINUE
+      RETURN
+      END

avevar.dem

@@ -0,0 +1,16 @@
+	PROGRAM D13R5
+C	Driver for routine AVEVAR
+	PARAMETER(NPTS=1000, EPS=0.1)
+	DIMENSION DATA(NPTS)
+C	Generate Gaussian distributed data
+	IDUM=-5
+	WRITE(*,'(1X,T4,A,T14,A,T26,A)') 'Shift','Average','Variance'
+	DO 12 I=1,11
+		SHIFT=(I-1)*EPS
+		DO 11 J=1,NPTS
+			DATA(J)=SHIFT+I*GASDEV(IDUM)
+11		CONTINUE
+		CALL AVEVAR(DATA,NPTS,AVE,VAR)
+		WRITE(*,'(1X,F6.2,2F12.2)') SHIFT,AVE,VAR
+12	CONTINUE
+	END

avevar.f

@@ -0,0 +1,15 @@
+      SUBROUTINE AVEVAR(DATA,N,AVE,VAR)
+      DIMENSION DATA(N)
+      AVE=0.0
+      VAR=0.0
+      DO 11 J=1,N
+        AVE=AVE+DATA(J)
+11    CONTINUE
+      AVE=AVE/N
+      DO 12 J=1,N
+        S=DATA(J)-AVE
+        VAR=VAR+S*S
+12    CONTINUE
+      VAR=VAR/(N-1)
+      RETURN
+      END