	PROGRAM D13R2
C	Driver for routine MDIAN1
	PARAMETER(NPTS=50)
	DIMENSION DATA(NPTS)
	IDUM=-5
	DO 11 I=1,NPTS
		DATA(I)=GASDEV(IDUM)
11	CONTINUE
	CALL MDIAN1(DATA,NPTS,XMED)
	WRITE(*,'(1X,A/)') 'Gaussian distrib., zero mean, unit variance'
	WRITE(*,'(1X,A,F10.6/)') 'Median of data set is',XMED
	WRITE(*,'(1X,A/,(5F12.6))') 'Sorted data',
     *			(DATA(I),I=1,50)
	END
	PROGRAM D13R3
C	Driver for routine MDIAN2
	PARAMETER(NPTS=50)
	DIMENSION DATA(NPTS)
	IDUM=-5
	DO 11 I=1,NPTS
		DATA(I)=GASDEV(IDUM)
11	CONTINUE
	CALL MDIAN2(DATA,NPTS,XMED)
	WRITE(*,'(1X,A/)') 'Gaussian distrib., zero mean, unit variance'
	WRITE(*,'(1X,A,F12.6)') 'Median according to MDIAN2 is',XMED
	CALL MDIAN1(DATA,NPTS,XMED)
	WRITE(*,'(1X,A,F12.6/)') 'Median according to MDIAN1 is',XMED
	END
	PROGRAM D14R10
C	Driver for routine MEDFIT
	PARAMETER(NPT=100,SPREAD=0.1)
	DIMENSION X(NPT),Y(NPT),SIG(NPT)
	IDUM=-1984
	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
	MWT=1
	CALL FIT(X,Y,NPT,SIG,MWT,A,B,SIGA,SIGB,CHI2,Q)
	WRITE(*,'(/1X,A)') 'According to routine FIT the result is:'
	WRITE(*,'(1X,T5,A,F8.4,T20,A,F8.4)') 'A = ',A,'Uncertainty: ',
     *		SIGA
	WRITE(*,'(1X,T5,A,F8.4,T20,A,F8.4)') 'B = ',B,'Uncertainty: ',
     *		SIGB
	WRITE(*,'(1X,T5,A,F8.4,A,I4,A)') 'Chi-squared: ',CHI2,
     *		' for ',NPT,' points'
	WRITE(*,'(1X,T5,A,F8.4)') 'Goodness-of-fit: ',Q
	WRITE(*,'(/1X,A)') 'According to routine MEDFIT the result is:'
	CALL MEDFIT(X,Y,NPT,A,B,ABDEV)
	WRITE(*,'(1X,T5,A,F8.4)') 'A = ',A
	WRITE(*,'(1X,T5,A,F8.4)') 'B = ',B
	WRITE(*,'(1X,T5,A,F8.4)') 'Absolute deviation (per data point): ',ABDEV
	WRITE(*,'(1X,T5,A,F8.4,A)') '(note: Gaussian spread is',SPREAD,')'
	END
	PROGRAM D12R9
C	Driver for routine MEMCOF
	PARAMETER(N=1000,M=10)
	DIMENSION DATA(N),COF(M),WK1(N),WK2(N),WKM(M)
	OPEN(5,FILE='SPCTRL.DAT',STATUS='OLD')
	READ(5,*) (DATA(I),I=1,N)
	CLOSE(5)
	CALL MEMCOF(DATA,N,M,PM,COF,WK1,WK2,WKM)
	WRITE(*,'(/1X,A/)') 'Coeff. for spectral estim. of SPCTRL.DAT'
	DO 11 I=1,M
		WRITE(*,'(1X,A,I2,A,F12.6)') 'a[',I,'] =',COF(I)
11	CONTINUE
	WRITE(*,'(/1X,A,F12.6/)') 'b0 =',PM
	END
	PROGRAM D4R5
C	Driver for routine MIDPNT
	PARAMETER(NMAX=10)
	EXTERNAL FUNC2,FINT2
	A=0.0
	B=1.0
	WRITE(*,*) 'Integral of FUNC2 computed with MIDPNT'
	WRITE(*,*) 'Actual value of integral is',FINT2(B)-FINT2(A)
	WRITE(*,'(1X,T7,A,T20,A)') 'n','Approx. Integral'
	DO 11 I=1,NMAX
		CALL MIDPNT(FUNC2,A,B,S,I)
		WRITE(*,'(1X,I6,F24.6)') I,S
11	CONTINUE
	END

	FUNCTION FUNC2(X)
	FUNC2=1.0/SQRT(X)
	END

	FUNCTION FINT2(X)
C	Integral of FUNC2
	FINT2=2.0*SQRT(X)
	END
	PROGRAM D15R5
C	Driver for routine MMID
	EXTERNAL DERIVS
	PARAMETER(NVAR=4,X1=1.0,HTOT=0.5)
	DIMENSION Y(NVAR),YOUT(NVAR),DYDX(NVAR)
	Y(1)=BESSJ0(X1)
	Y(2)=BESSJ1(X1)
	Y(3)=BESSJ(2,X1)
	Y(4)=BESSJ(3,X1)
	DYDX(1)=-Y(2)
	DYDX(2)=Y(1)-Y(2)
	DYDX(3)=Y(2)-2.0*Y(3)
	DYDX(4)=Y(3)-3.0*Y(4)
	XF=X1+HTOT
	B1=BESSJ0(XF)
	B2=BESSJ1(XF)
	B3=BESSJ(2,XF)
	B4=BESSJ(3,XF)
	WRITE(*,'(1X,A/)') 'First four Bessel functions'
	DO 11 I=5,50,5
		CALL MMID(Y,DYDX,NVAR,X1,HTOT,I,YOUT,DERIVS)
		WRITE(*,'(1X,A,F6.4,A,F6.4,A,I2,A)') 'X = ',X1,
     *			' to ',X1+HTOT,' in ',I,' steps'
		WRITE(*,'(1X,T5,A,T20,A)') 'Integration','BESSJ'
		WRITE(*,'(1X,2F12.6)') YOUT(1),B1
		WRITE(*,'(1X,2F12.6)') YOUT(2),B2
		WRITE(*,'(1X,2F12.6)') YOUT(3),B3
		WRITE(*,'(1X,2F12.6)') YOUT(4),B4
		WRITE(*,'(/1X,A)') 'press RETURN to continue...'
		READ(*,*)
11	CONTINUE
	END

	SUBROUTINE DERIVS(X,Y,DYDX)
	DIMENSION Y(1),DYDX(1)
	DYDX(1)=-Y(2)
	DYDX(2)=Y(1)-(1.0/X)*Y(2)
	DYDX(3)=Y(2)-(2.0/X)*Y(3)
	DYDX(4)=Y(3)-(3.0/X)*Y(4)
	END

	PROGRAM D10R1
C	Driver for routine MNBRAK
	EXTERNAL BESSJ0
	DO 11 I=1,10
		AX=I*0.5
		BX=(I+1.0)*0.5
		CALL MNBRAK(AX,BX,CX,FA,FB,FC,BESSJ0)
		WRITE(*,'(1X,T13,A,T25,A,T37,A)') 'A','B','C'
		WRITE(*,'(1X,A3,T5,3F12.6)') 'X',AX,BX,CX
		WRITE(*,'(1X,A3,T5,3F12.6)') 'F',FA,FB,FC
11	CONTINUE
	END
	PROGRAM D9R13
C	Driver for routine MNEWT
	PARAMETER(NTRIAL=5,TOLX=1.0E-6,N=4,TOLF=1.0E-6,NP=15)
	DIMENSION X(NP),ALPHA(NP,NP),BETA(NP)
	DO 15 KK=-1,1,2
		DO 14 K=1,3
			XX=0.2*K*KK
			WRITE(*,'(/1X,A,I2)') 'Starting vector number',K
			DO 11 I=1,4
				X(I)=XX+0.2*I
				WRITE(*,'(1X,T5,A,I1,A,F5.2)')'X(',I,') = ',X(I)
11			CONTINUE
			DO 13 J=1,NTRIAL
				CALL MNEWT(1,X,N,TOLX,TOLF)
				CALL USRFUN(X,ALPHA,BETA)
				WRITE(*,'(/1X,T5,A,T14,A,T29,A/)')'I','X(I)','F'
				DO 12 I=1,N
					WRITE(*,'(1X,I4,2E15.6)')I,X(I),-BETA(I)
12				CONTINUE
				WRITE(*,'(/1X,A)') 'press RETURN to continue...'
				READ(*,*)
13			CONTINUE
14		CONTINUE
15	CONTINUE
	END

	SUBROUTINE USRFUN(X,ALPHA,BETA)
	PARAMETER(NP=15)
	DIMENSION ALPHA(NP,NP),BETA(NP),X(NP)
	ALPHA(1,1)=-2.0*X(1)
	ALPHA(1,2)=-2.0*X(2)
	ALPHA(1,3)=-2.0*X(3)
	ALPHA(1,4)=1.0
	ALPHA(2,1)=2.0*X(1)
	ALPHA(2,2)=2.0*X(2)
	ALPHA(2,3)=2.0*X(3)
	ALPHA(2,4)=2.0*X(4)
	ALPHA(3,1)=1.0
	ALPHA(3,2)=-1.0
	ALPHA(3,3)=0.0
	ALPHA(3,4)=0.0
	ALPHA(4,1)=0.0
	ALPHA(4,2)=1.0
	ALPHA(4,3)=-1.0
	ALPHA(4,4)=0.0
	BETA(1)=X(1)**2+X(2)**2+X(3)**2-X(4)
	BETA(2)=-X(1)**2-X(2)**2-X(3)**2-X(4)**2+1.0
	BETA(3)=-X(1)+X(2)
	BETA(4)=-X(2)+X(3)
	END
	PROGRAM D13R1
C	Driver for routine MOMENT
	PARAMETER(PI=3.14159265,NPTS=10000,NBIN=100,NDAT=NPTS+NBIN)
	DIMENSION DATA(NDAT)
	I=1
	DO 12 J=1,NBIN
		X=PI*J/NBIN
		NLIM=NINT(SIN(X)*PI/2.0*NPTS/NBIN)
		DO 11 K=1,NLIM
			DATA(I)=X
			I=I+1
11		CONTINUE
12	CONTINUE
	WRITE(*,'(1X,A/)') 'Moments of a sinusoidal distribution'
	CALL MOMENT(DATA,NPTS,AVE,ADEV,SDEV,VAR,SKEW,CURT)
	WRITE(*,'(1X,T29,A,T42,A/)') 'Calculated','Expected'
	WRITE(*,'(1X,A,T25,2F12.4)') 'Mean :',AVE,PI/2.0
	WRITE(*,'(1X,A,T25,2F12.4)') 'Average Deviation :',ADEV,0.570796
	WRITE(*,'(1X,A,T25,2F12.4)') 'Standard Deviation :',SDEV,0.683667
	WRITE(*,'(1X,A,T25,2F12.4)') 'Variance :',VAR,0.467401
	WRITE(*,'(1X,A,T25,2F12.4)') 'Skewness :',SKEW,0.0
	WRITE(*,'(1X,A,T25,2F12.4)') 'Kurtosis :',CURT,-0.806249
	END
	PROGRAM D2R5
C	Driver for routine MPROVE
	PARAMETER(N=5,NP=5)
	DIMENSION A(NP,NP),INDX(N),B(N),X(N),AA(NP,NP)
	DATA A/1.0,2.0,1.0,4.0,5.0,2.0,3.0,1.0,5.0,1.0,
     *		3.0,4.0,1.0,1.0,2.0,4.0,5.0,1.0,2.0,3.0,
     *		5.0,1.0,1.0,3.0,4.0/
	DATA B/1.0,1.0,1.0,1.0,1.0/
	DO 12 I=1,N
		X(I)=B(I)
		DO 11 J=1,N
			AA(I,J)=A(I,J)
11		CONTINUE
12	CONTINUE
	CALL LUDCMP(AA,N,NP,INDX,D)
	CALL LUBKSB(AA,N,NP,INDX,X)
	WRITE(*,'(/1X,A)') 'Solution vector for the equations:'
	WRITE(*,'(1X,5F12.6)') (X(I),I=1,N)
C	Now phoney up X and let MPROVE fit it
	IDUM=-13
	DO 13 I=1,N
		X(I)=X(I)*(1.0+0.2*RAN3(IDUM))
13	CONTINUE
	WRITE(*,'(/1X,A)') 'Solution vector with noise added:'
	WRITE(*,'(1X,5F12.6)') (X(I),I=1,N)
	CALL MPROVE(A,AA,N,NP,INDX,B,X)
	WRITE(*,'(/1X,A)') 'Solution vector recovered by MPROVE:'
	WRITE(*,'(1X,5F12.6)') (X(I),I=1,N)
	END
	PROGRAM D14R8B
C	Driver for routine MRQCOF
	EXTERNAL FGAUSS
	PARAMETER(NPT=100,MA=6,SPREAD=0.1)
	DIMENSION X(NPT),Y(NPT),SIG(NPT),A(MA),LISTA(MA),
     *		COVAR(MA,MA),ALPHA(MA,MA),BETA(MA),GUES(MA)
	DATA A/5.0,2.0,3.0,2.0,5.0,3.0/
	DATA GUES/4.9,2.1,2.9,2.1,4.9,3.1/
	IDUM=-911
C	First try sum of two gaussians
	DO 12 I=1,100
		X(I)=0.1*I
		Y(I)=0.0
		DO 11 J=1,4,3
			Y(I)=Y(I)+A(J)*EXP(-((X(I)-A(J+1))/A(J+2))**2)
11		CONTINUE
		Y(I)=Y(I)*(1.0+SPREAD*GASDEV(IDUM))
		SIG(I)=SPREAD*Y(I)
12	CONTINUE
	MFIT=MA
	DO 13 I=1,MFIT
		LISTA(I)=I
13	CONTINUE
	DO 14 I=1,MA
		A(I)=GUES(I)
14	CONTINUE
	CALL MRQCOF(X,Y,SIG,NPT,A,MA,LISTA,MFIT,ALPHA,
     *		BETA,MA,CHISQ,FGAUSS)
	WRITE(*,'(/1X,A)') 'matrix alpha'
	DO 15 I=1,MA
		WRITE(*,'(1X,6F12.4)') (ALPHA(I,J),J=1,MA)
15	CONTINUE
	WRITE(*,'(1X,A)') 'vector beta'
	WRITE(*,'(1X,6F12.4)') (BETA(I),I=1,MA)
	WRITE(*,'(1X,A,F12.4/)') 'Chi-squared:',CHISQ
C	Next fix one line and improve the other
	DO 16 I=1,3
		LISTA(I)=I+3
16	CONTINUE
	MFIT=3
	DO 17 I=1,MA
		A(I)=GUES(I)
17	CONTINUE
	CALL MRQCOF(X,Y,SIG,NPT,A,MA,LISTA,MFIT,
     *		ALPHA,BETA,MA,CHISQ,FGAUSS)
	WRITE(*,'(1X,A)') 'matrix alpha'
	DO 18 I=1,MFIT
		WRITE(*,'(1X,6F12.4)') (ALPHA(I,J),J=1,MFIT)
18 	CONTINUE
	WRITE(*,'(1X,A)') 'vector beta'
	WRITE(*,'(1X,6F12.4)') (BETA(I),I=1,MFIT)
	WRITE(*,'(1X,A,F12.4/)') 'Chi-squared:',CHISQ
	END
	PROGRAM D14R8A
C	Driver for routine MRQMIN
	EXTERNAL FGAUSS
	PARAMETER(NPT=100,MA=6,SPREAD=0.001)
	DIMENSION X(NPT),Y(NPT),SIG(NPT),A(MA),LISTA(MA),
     *			COVAR(MA,MA),ALPHA(MA,MA),GUES(MA)
	DATA A/5.0,2.0,3.0,2.0,5.0,3.0/
	DATA GUES/4.5,2.2,2.8,2.5,4.9,2.8/
	IDUM=-911
C	First try a sum of two Gaussians
	DO 12 I=1,100
		X(I)=0.1*I
		Y(I)=0.0
		DO 11 J=1,4,3
			Y(I)=Y(I)+A(J)*EXP(-((X(I)-A(J+1))/A(J+2))**2)
11		CONTINUE
		Y(I)=Y(I)*(1.0+SPREAD*GASDEV(IDUM))
		SIG(I)=SPREAD*Y(I)
12	CONTINUE
	MFIT=MA
	DO 13 I=1,MFIT
		LISTA(I)=I
13	CONTINUE
	ALAMDA=-1
	DO 14 I=1,MA
		A(I)=GUES(I)
14	CONTINUE
	CALL MRQMIN(X,Y,SIG,NPT,A,MA,LISTA,MFIT,COVAR,ALPHA,
     *			MA,CHISQ,FGAUSS,ALAMDA)		
	K=1
	ITST=0
1	WRITE(*,'(/1X,A,I2,T18,A,F10.4,T43,A,E9.2)') 'Iteration #',K,
     *		'Chi-squared:',CHISQ,'ALAMDA:',ALAMDA
	WRITE(*,'(1X,T5,A,T13,A,T21,A,T29,A,T37,A,T45,A)') 'A(1)',
     *		'A(2)','A(3)','A(4)','A(5)','A(6)'
	WRITE(*,'(1X,6F8.4)') (A(I),I=1,6)
	K=K+1
	OCHISQ=CHISQ
	CALL MRQMIN(X,Y,SIG,NPT,A,MA,LISTA,MFIT,COVAR,ALPHA,
     *			MA,CHISQ,FGAUSS,ALAMDA)		
	IF (CHISQ.GT.OCHISQ) THEN
		ITST=0
	ELSE IF (ABS(OCHISQ-CHISQ).LT.0.1) THEN
		ITST=ITST+1
	ENDIF
	IF (ITST.LT.2) THEN
		GOTO 1
	ENDIF
	ALAMDA=0.0
	CALL MRQMIN(X,Y,SIG,NPT,A,MA,LISTA,MFIT,COVAR,ALPHA,
     *			MA,CHISQ,FGAUSS,ALAMDA)		
	WRITE(*,*) 'Uncertainties:'
	WRITE(*,'(1X,6F8.4/)') (SQRT(COVAR(I,I)),I=1,6)
	WRITE(*,'(1X,A)') 'Expected results:'
	WRITE(*,'(1X,F7.2,5F8.2/)') 5.0,2.0,3.0,2.0,5.0,3.0
	END
