	PROGRAM D10R8
C	Driver for routine F1DIM
	PARAMETER(NDIM=3,NMAX=50)
	EXTERNAL F1DIM
	COMMON /F1COM/ NCOM,PCOM(NMAX),XICOM(NMAX)
	DIMENSION P(NDIM),XI(NDIM)
	DATA P/0.0,0.0,0.0/
	NCOM=NDIM
	WRITE(*,'(/1X,A)') 'Enter vector direction along which to'
	WRITE(*,'(1X,A)') 'plot the function. Minimum is in the'
	WRITE(*,'(1X,A)') 'direction 1.0,1.0,1.0 - Enter X,Y,Z:'
	READ(*,*) (XI(I),I=1,3)
	DO 11 J=1,NDIM
		PCOM(J)=P(J)
		XICOM(J)=XI(J)
11	CONTINUE
	CALL SCRSHO(F1DIM)
	END

	FUNCTION FUNC(X)
	DIMENSION X(3)
	FUNC=0.0
	DO 11 I=1,3
		FUNC=FUNC+(X(I)-1.0)**2
11	CONTINUE
	END
	PROGRAM D6R4
C	Driver for routine FACTLN
	CHARACTER TEXT*11
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'N-factorial') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) 'Log of N-factorial'
	WRITE(*,'(1X,T6,A1,T18,A6,T34,A9)')
     *		'N','Actual','FACTLN(N)'
	DO 11 I=1,NVAL
		READ(5,*) N,VALUE
		WRITE(*,'(I6,2F18.6)') N,LOG(VALUE),FACTLN(N)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R2
C	Driver for routine FACTRL
	CHARACTER TEXT*11
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'N-factorial') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T6,A1,T21,A6,T38,A9)')
     *		'N','Actual','FACTRL(N)'
	DO 11 I=1,NVAL
		READ(5,*) N,ACTUAL
		IF (ACTUAL.LT.(1.0E10)) THEN
			WRITE(*,'(I6,2F20.0)') N,ACTUAL,FACTRL(N)
		ELSE
			WRITE(*,'(I6,2E20.7)') N,ACTUAL,FACTRL(N)
		ENDIF
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D14R9
C	Driver for routine FGAUSS
	PARAMETER(NPT=3,NLIN=2,NA=3*NLIN)
	DIMENSION A(NA),DYDA(NA),DF(NA)
	DATA A/3.0,0.2,0.5,1.0,0.7,0.3/
	WRITE(*,'(/1X,T6,A,T14,A,T19,A,T27,A,T35,A,T43,A,T51,A,T59,A)')
     *		'X','Y','DYDA1','DYDA2','DYDA3','DYDA4','DYDA5','DYDA6'
	DO 11 I=1,NPT
		X=0.3*I
		CALL FGAUSS(X,A,Y,DYDA,NA)
		E1=EXP(-((X-A(2))/A(3))**2)
		E2=EXP(-((X-A(5))/A(6))**2)
		F=A(1)*E1+A(4)*E2
		DF(1)=E1
		DF(4)=E2
		DF(2)=A(1)*E1*2.0*(X-A(2))/(A(3)**2)
		DF(5)=A(4)*E2*2.0*(X-A(5))/(A(6)**2)
		DF(3)=A(1)*E1*2.0*((X-A(2))**2)/(A(3)**3)
		DF(6)=A(4)*E2*2.0*((X-A(5))**2)/(A(6)**3)
		WRITE(*,'(1X,A/,8F8.4)') 'from FGAUSS',X,Y,(DYDA(J),J=1,6)
		WRITE(*,'(1X,A/,8F8.4/)') 'independent calc.',X,F,(DF(J),J=1,6)
11	CONTINUE
	END
	PROGRAM D14R1
C	Driver for routine FIT
	PARAMETER(NPT=100,SPREAD=0.5)
	DIMENSION X(NPT),Y(NPT),SIG(NPT)
	IDUM=-117
	DO 11 I=1,NPT
		X(I)=0.1*I
		Y(I)=-2.0*X(I)+1.0+SPREAD*GASDEV(IDUM)
		SIG(I)=SPREAD
11	CONTINUE
	DO 12 MWT=0,1
		CALL FIT(X,Y,NPT,SIG,MWT,A,B,SIGA,SIGB,CHI2,Q)
		IF (MWT.EQ.0) THEN
			WRITE(*,'(//1X,A)') 'Ignoring standard deviation'
		ELSE
			WRITE(*,'(//1X,A)') 'Including standard deviation'
		ENDIF
		WRITE(*,'(1X,T5,A,F9.6,T24,A,F9.6)') 'A = ',A,'Uncertainty: ',
     *			SIGA
		WRITE(*,'(1X,T5,A,F9.6,T24,A,F9.6)') 'B = ',B,'Uncertainty: ',
     *			SIGB
		WRITE(*,'(1X,T5,A,4X,F10.6)') 'Chi-squared: ',CHI2
		WRITE(*,'(1X,T5,A,F10.6)') 'Goodness-of-fit: ',Q
12	CONTINUE
	END
	PROGRAM D12R11
C	Driver for routine FIXRTS
	PARAMETER(NPOLES=6,NPOL=NPOLES+1)
	DIMENSION D(NPOLES)
	COMPLEX ZCOEF(NPOL),ZEROS(NPOLES),Z
	LOGICAL POLISH
	DATA D/6.0,-15.0,20.0,-15.0,6.0,0.0/
C	Finding roots of (z-1.0)**6=1.0
C	First print roots
	ZCOEF(NPOLES+1)=CMPLX(1.0,0.0)
	DO 11 I=NPOLES,1,-1
		ZCOEF(I)=CMPLX(-D(NPOLES+1-I),0.0)
11	CONTINUE
	POLISH=.TRUE.
	CALL ZROOTS(ZCOEF,NPOLES,ZEROS,POLISH)
	WRITE(*,'(/1X,A)') 'Roots of (z-1.0)^6 = 1.0'
	WRITE(*,'(1X,T20,A,T42,A)') 'Root','(z-1.0)^6'
	DO 12 I=1,NPOLES
		Z=(ZEROS(I)-1.0)**6
		WRITE(*,'(1X,I6,4F12.6)') I,ZEROS(I),Z
12	CONTINUE
C	Now fix them to lie within unit circle
	CALL FIXRTS(D,NPOLES)
C	Check results
	ZCOEF(NPOLES+1)=CMPLX(1.0,0.0)
	DO 13 I=NPOLES,1,-1
		ZCOEF(I)=CMPLX(-D(NPOLES+1-I),0.0)
13	CONTINUE
	CALL ZROOTS(ZCOEF,NPOLES,ZEROS,POLISH)
	WRITE(*,'(/1X,A)') 'Roots reflected in unit circle'
	WRITE(*,'(1X,T20,A,T42,A)') 'Root','(z-1.0)^6'
	DO 14 I=1,NPOLES
		Z=(ZEROS(I)-1.0)**6
		WRITE(*,'(1X,I6,4F12.6)') I,ZEROS(I),Z
14	CONTINUE
	END
	PROGRAM D14R7
C	Driver for routine FLEG
	PARAMETER(NVAL=5,DX=0.2,NPOLY=5)
	DIMENSION AFUNC(NPOLY)
	WRITE(*,'(/1X,T25,A)') 'Legendre Polynomials'
	WRITE(*,'(/1X,T8,A,T18,A,T28,A,T38,A,T48,A)')
     *			'N=1','N=2','N=3','N=4','N=5'
	DO 11 I=1,NVAL
		X=I*DX
		CALL FLEG(X,AFUNC,NPOLY)
		WRITE(*,'(1X,A,F6.2)') 'X =',X
		WRITE(*,'(1X,5F10.4,A)') (AFUNC(J),J=1,NPOLY),'  routine FLEG'
		WRITE(*,'(1X,5F10.4,A/)') (PLGNDR(J-1,0,X),J=1,NPOLY),
     *			'  routine PLGNDR'
11	CONTINUE
	END
	PROGRAM D1R1
C	Driver for routine FLMOON
	PARAMETER(TZONE=5.0)
	CHARACTER PHASE(4)*15,TIMSTR(2)*3
	DATA PHASE/'new moon','first quarter',
     *			'full moon','last quarter'/
	DATA TIMSTR/' AM',' PM'/
	WRITE(*,*) 'Date of the next few phases of the moon'
	WRITE(*,*) 'Enter today''s date (e.g. 1,31,1982)'
	TIMZON=-TZONE/24.0
	READ(*,*) IM,ID,IY
C	Approximate number of full moons since January 1900
	N=12.37*(IY-1900+(IM-0.5)/12.0)
	NPH=2
	J1=JULDAY(IM,ID,IY)
	CALL FLMOON(N,NPH,J2,FRAC)
	N=N+(J1-J2)/28.0
	WRITE(*,'(/1X,T6,A,T19,A,T32,A)') 'Date','Time(EST)','Phase'
	DO 11 I=1,20
		CALL FLMOON(N,NPH,J2,FRAC)
		IFRAC=NINT(24.*(FRAC+TIMZON))
		IF (IFRAC.LT.0) THEN
			J2=J2-1
			IFRAC=IFRAC+24
		ENDIF
		IF (IFRAC.GE.12) THEN
			J2=J2+1
			IFRAC=IFRAC-12
		ELSE
			IFRAC=IFRAC+12
		ENDIF
		IF (IFRAC.GT.12) THEN
			IFRAC=IFRAC-12
			ISTR=2
		ELSE
			ISTR=1
		ENDIF
		CALL CALDAT(J2,IM,ID,IY)
		WRITE(*,'(1X,2I3,I5,T20,I2,A,5X,A)') IM,ID,IY,
     *			IFRAC,TIMSTR(ISTR),PHASE(NPH+1)
		IF (NPH.EQ.3) THEN
			NPH=0
			N=N+1
		ELSE
			NPH=NPH+1
		ENDIF
11	CONTINUE
	END
	PROGRAM D12R1
C	Driver for routine FOUR1
	PARAMETER (NN=32,NN2=2*NN)
	DIMENSION DATA(NN2),DCMP(NN2)
	WRITE(*,*) 'h(t)=real-valued even-function'
	WRITE(*,*) 'H(n)=H(N-n) and real?'
	DO 11 I=1,2*NN-1,2
		DATA(I)=1.0/(((I-NN-1.0)/NN)**2+1.0)
		DATA(I+1)=0.0
11	CONTINUE
	ISIGN=1
	CALL FOUR1(DATA,NN,ISIGN)
	CALL PRNTFT(DATA,NN2)
	WRITE(*,*) 'h(t)=imaginary-valued even-function'
	WRITE(*,*) 'H(n)=H(N-n) and imaginary?'
	DO 12 I=1,2*NN-1,2
		DATA(I+1)=1.0/(((I-NN-1.0)/NN)**2+1.0)
		DATA(I)=0.0
12	CONTINUE
	ISIGN=1
	CALL FOUR1(DATA,NN,ISIGN)
	CALL PRNTFT(DATA,NN2)
	WRITE(*,*) 'h(t)=real-valued odd-function'
	WRITE(*,*) 'H(n)=-H(N-n) and imaginary?'
	DO 13 I=1,2*NN-1,2
		DATA(I)=(I-NN-1.0)/NN/(((I-NN-1.0)/NN)**2+1.0)
		DATA(I+1)=0.0
13	CONTINUE
	DATA(1)=0.0
	ISIGN=1
	CALL FOUR1(DATA,NN,ISIGN)
	CALL PRNTFT(DATA,NN2)
	WRITE(*,*) 'h(t)=imaginary-valued odd-function'
	WRITE(*,*) 'H(n)=-H(N-n) and real?'
	DO 14 I=1,2*NN-1,2
		DATA(I+1)=(I-NN-1.0)/NN/(((I-NN-1.0)/NN)**2+1.0)
		DATA(I)=0.0
14	CONTINUE
	DATA(2)=0.0
	ISIGN=1
	CALL FOUR1(DATA,NN,ISIGN)
	CALL PRNTFT(DATA,NN2)
C	Transform, inverse-transform test
	DO 15 I=1,2*NN-1,2
		DATA(I)=1.0/((0.5*(I-NN-1)/NN)**2+1.0)
		DCMP(I)=DATA(I)
		DATA(I+1)=(0.25*(I-NN-1)/NN)*
     *			EXP(-(0.5*(I-NN-1.0)/NN)**2)
		DCMP(I+1)=DATA(I+1)
15	CONTINUE
	ISIGN=1
	CALL FOUR1(DATA,NN,ISIGN)
	ISIGN=-1
	CALL FOUR1(DATA,NN,ISIGN)
	WRITE(*,'(/1X,T10,A,T44,A)') 'Double Fourier Transform:',
     *		'Original Data:'
	WRITE(*,'(/1X,T5,A,T11,A,T24,A,T41,A,T53,A/)')
     *		'k','Real h(k)','Imag h(k)','Real h(k)','Imag h(k)'
	DO 16 I=1,NN,2
		J=(I+1)/2
		WRITE(*,'(1X,I4,2X,2F12.6,5X,2F12.6)') J,DCMP(I),
     *			DCMP(I+1),DATA(I)/NN,DATA(I+1)/NN
16	CONTINUE
	END

	SUBROUTINE PRNTFT(DATA,NN2)
	DIMENSION DATA(NN2)
	WRITE(*,'(/1X,T5,A,T11,A,T23,A,T39,A,T52,A)')
     *		'n','Real H(n)','Imag H(n)','Real H(N-n)','Imag H(N-n)'
	WRITE(*,'(1X,I4,2X,2F12.6,5X,2F12.6)') 0,DATA(1),DATA(2),
     *		DATA(1),DATA(2)
	DO 11 N=3,(NN2/2)+1,2
		M=(N-1)/2
		MM=NN2+2-N
		WRITE(*,'(1X,I4,2X,2F12.6,5X,2F12.6)') M,DATA(N),
     *			DATA(N+1),DATA(MM),DATA(MM+1)
11	CONTINUE
	WRITE(*,'(/1X,A)') ' press RETURN to continue ...'
	READ(*,*)
	RETURN
	END
	PROGRAM D12R13
C	Driver for routine FOURN
	PARAMETER(NDIM=3,NDAT=1024)
	DIMENSION NN(NDIM),DATA(NDAT)
	DO 11 I=1,NDIM
		NN(I)=2*(2**I)
11	CONTINUE
	DO 14 I=1,NN(3)
		DO 13 J=1,NN(2)
			DO 12 K=1,NN(1)
				L=K+(J-1)*NN(1)+(I-1)*NN(2)*NN(1)
 				LL=2*L-1
				DATA(LL)=FLOAT(LL)
				DATA(LL+1)=FLOAT(LL+1)
12			CONTINUE
13		CONTINUE
14	CONTINUE
	ISIGN=+1
	CALL FOURN(DATA,NN,NDIM,ISIGN)
	ISIGN=-1
	WRITE(*,'(1X,A)') 'Double 3-dimensional Transform'
	WRITE(*,'(/1X,T10,A,T35,A,T63,A)') 'Double Transf.',
     *		'Original Data','Ratio'
	WRITE(*,'(1X,T8,A,T20,A,T33,A,T45,A,T57,A,T69,A/)') 
     *		'Real','Imag.','Real','Imag.','Real','Imag.'
	CALL FOURN(DATA,NN,NDIM,ISIGN)
	DO 15 I=1,4
		J=2*I
		K=2*J
		L=K+(J-1)*NN(1)+(I-1)*NN(2)*NN(1)
		LL=2*L-1
		WRITE(*,'(1X,6F12.2)') DATA(LL),DATA(LL+1),FLOAT(LL),FLOAT(LL+1),
     *			DATA(LL)/LL,DATA(LL+1)/(LL+1)
15	CONTINUE
	WRITE(*,'(/1X,A,I4)') 'The product of transform lengths is:',
     *		NN(1)*NN(2)*NN(3)
	END
	PROGRAM D14R6
C	Driver for FPOLY
	PARAMETER(NVAL=15,DX=0.1,NPOLY=5)
	DIMENSION AFUNC(NPOLY)
	WRITE(*,'(/1X,T29,A)') 'Powers of X'
	WRITE(*,'(/1X,T9,A,T17,A,T27,A,T37,A,T47,A,T57,A)') 'X','X**0',
     *			'X**1','X**2','X**3','X**4'
	DO 11 I=1,NVAL
		X=I*DX
		CALL FPOLY(X,AFUNC,NPOLY)
		WRITE(*,'(1X,6F10.4)') X,(AFUNC(J),J=1,NPOLY)
11	CONTINUE
	END
	PROGRAM D10R9
C	Driver for routine FRPRMN
	PARAMETER(NDIM=3,FTOL=1.0E-6,PIO2=1.5707963)
	DIMENSION P(NDIM)
	WRITE(*,'(/1X,A)') 'Program finds the minimum of a function'
	WRITE(*,'(1X,A)') 'with different trial starting vectors.'
	WRITE(*,'(1X,A)') 'True minimum is (0.5,0.5,0.5)'
	DO 11 K=0,4
		ANGL=PIO2*K/4.0
		P(1)=2.0*COS(ANGL)
		P(2)=2.0*SIN(ANGL)
		P(3)=0.0
		WRITE(*,'(/1X,A,3(F6.4,A))') 'Starting vector: (',
     *			P(1),',',P(2),',',P(3),')'
		CALL FRPRMN(P,NDIM,FTOL,ITER,FRET)
		WRITE(*,'(1X,A,I3)') 'Iterations:',ITER
		WRITE(*,'(1X,A,3(F6.4,A))') 'Solution vector: (',
     *			P(1),',',P(2),',',P(3),')'
		WRITE(*,'(1X,A,E14.6)') 'Func. value at solution',FRET
11	CONTINUE
	END

	FUNCTION FUNC(X)
	DIMENSION X(3)
	FUNC=1.0-BESSJ0(X(1)-0.5)*BESSJ0(X(2)-0.5)*BESSJ0(X(3)-0.5)
	END

	SUBROUTINE DFUNC(X,DF)
	PARAMETER (NMAX=50)
	DIMENSION X(3),DF(NMAX)
	DF(1)=BESSJ1(X(1)-0.5)*BESSJ0(X(2)-0.5)*BESSJ0(X(3)-0.5)
	DF(2)=BESSJ0(X(1)-0.5)*BESSJ1(X(2)-0.5)*BESSJ0(X(3)-0.5)
	DF(3)=BESSJ0(X(1)-0.5)*BESSJ0(X(2)-0.5)*BESSJ1(X(3)-0.5)
	RETURN
	END
	PROGRAM D13R7
C	Driver for routine FTEST
	PARAMETER(NPTS=1000,MPTS=500,EPS=0.01,NVAL=10)
	DIMENSION DATA1(NPTS),DATA2(MPTS),DATA3(MPTS)
C	Generate two Gaussian distributions with
C	different variances
	IDUM=-13
	DO 11 J=1,NPTS
		DATA1(J)=GASDEV(IDUM)
11	CONTINUE
	DO 12 J=1,MPTS
		DATA2(J)=GASDEV(IDUM)
12	CONTINUE
	WRITE(*,'(1X,T5,A,F5.2)') 'Variance 1 = ',1.0
	WRITE(*,'(1X,T5,A,T21,A,T30,A)')
     *			'Variance 2','Ratio','Probability'
	DO 14 I=1,NVAL+1
		VAR=1.0+(I-1)*EPS
		FACTOR=SQRT(VAR)
		DO 13 J=1,MPTS
			DATA3(J)=FACTOR*DATA2(J)
13		CONTINUE
		CALL FTEST(DATA1,NPTS,DATA3,MPTS,F,PROB)
		WRITE(*,'(1X,F11.4,2X,2F12.4)') VAR,F,PROB
14	CONTINUE
	END
