	PROGRAM D7R7
C	Driver for routine GAMDEV
	PARAMETER(N=20,NPTS=10000,ISCAL=200,LLEN=50)
	DIMENSION DIST(21)
	CHARACTER TEXT(50)*1
	IDUM=-13
10	DO 11 J=1,21
		DIST(J)=0.0
11	CONTINUE
	WRITE(*,*) 'Order of Gamma distribution (n=1..20); -1 to end.'
	READ(*,*) IA
	IF (IA.LE.0) GOTO 99
	IF (IA.GT.20) GOTO 10		
	DO 12 I=1,NPTS
		J=INT(GAMDEV(IA,IDUM))+1
		IF ((J.GE.1).AND.(J.LE.21)) DIST(J)=DIST(J)+1
12	CONTINUE
	WRITE(*,'(1X,A,I2,A,I6,A)') 'Gamma-distribution deviate, order ',
     *		IA,' of ',NPTS,' points'
	WRITE(*,'(1X,T6,A,T14,A,T23,A)') 'x','p(x)','graph:'
	DO 15 J=1,20
		DIST(J)=DIST(J)/NPTS
		DO 13 K=1,50
			TEXT(K)=' '
13		CONTINUE
		KLIM=INT(ISCAL*DIST(J))
		IF (KLIM.GT.LLEN) KLIM=LLEN
		DO 14 K=1,KLIM
			TEXT(K)='*'
14		CONTINUE
		WRITE(*,'(1X,F7.2,F10.4,4X,50A1)')
     *			FLOAT(J),DIST(J),(TEXT(K),K=1,50)
15	CONTINUE
	GOTO 10
99	END
	PROGRAM D6R1
C	Driver for routine GAMMLN
	CHARACTER TEXT*14
	PARAMETER(PI=3.1415926)
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Gamma Function') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) 'Log of gamma function:'
	WRITE(*,'(1X,T11,A1,T24,A6,T40,A10)')
     *		'X','Actual','GAMMLN(X)'
	DO 11 I=1,NVAL
		READ(5,*) X,ACTUAL
		IF (X.GT.0.0) THEN
			IF (X.GE.1.0) THEN
				CALC=GAMMLN(X)
			ELSE
				CALC=GAMMLN(X+1.0)-LOG(X)
			ENDIF
			WRITE(*,'(F12.2,2F18.6)') X,LOG(ACTUAL),CALC
		ENDIF
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R6
C	Driver for routine GAMMP
	CHARACTER TEXT*25
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Incomplete Gamma Function') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A,T16,A,T25,A,T35,A)')
     *		'A','X','Actual','GAMMP(A,X)'
	DO 11 I=1,NVAL
		READ(5,*) A,X,VALUE
		WRITE(*,'(1X,F6.2,3F12.6)') A,X,VALUE,GAMMP(A,X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R7
C	Driver for routine GAMMQ
	CHARACTER TEXT*25
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Incomplete Gamma Function') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A,T16,A,T25,A,T35,A)')
     *		'A','X','Actual','GAMMQ(A,X)'
	DO 11 I=1,NVAL
		READ(5,*) A,X,VALUE
		WRITE(*,'(1X,F6.2,3F12.6)') A,X,1.0-VALUE,GAMMQ(A,X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D7R6
C	Driver for routine GASDEV
	PARAMETER(N=20,NP1=N+1,NOVER2=N/2,NPTS=10000,ISCAL=400,LLEN=50)
	DIMENSION DIST(NP1)
	CHARACTER TEXT(50)*1
	IDUM=-13
	DO 11 J=1,NP1
		DIST(J)=0.0
11	CONTINUE
	DO 12 I=1,NPTS
		J=NINT(0.25*N*GASDEV(IDUM))+NOVER2+1
		IF ((J.GE.1).AND.(J.LE.NP1)) DIST(J)=DIST(J)+1
12	CONTINUE
	WRITE(*,'(1X,A,I6,A)')
     *		'Normally distributed deviate of ',NPTS,' points'
	WRITE(*,'(1X,T6,A,T14,A,T23,A)') 'x','p(x)','graph:'
	DO 15 J=1,NP1
		DIST(J)=DIST(J)/NPTS
		DO 13 K=1,50
			TEXT(K)=' '
13		CONTINUE
		KLIM=INT(ISCAL*DIST(J))
		IF (KLIM.GT.LLEN) KLIM=LLEN
		DO 14 K=1,KLIM
			TEXT(K)='*'
14		CONTINUE
		X=FLOAT(J)/(0.25*N)
		WRITE(*,'(1X,F7.2,F10.4,4X,50A1)')
     *			X,DIST(J),(TEXT(K),K=1,50)
15	CONTINUE
	END
	PROGRAM D4R9
C	Driver for routine GAULEG
	PARAMETER(NPOINT=10,X1=0.0,X2=1.0,X3=10.0)
	DIMENSION X(NPOINT),W(NPOINT)
	CALL GAULEG(X1,X2,X,W,NPOINT)
	WRITE(*,'(/1X,T3,A,T10,A,T22,A/)') '#','X(I)','W(I)'
	DO 11 I=1,NPOINT
		WRITE(*,'(1X,I2,2F12.6)') I,X(I),W(I)
11	CONTINUE
C	Demonstrate the use of GAULEG for an integral
	CALL GAULEG(X1,X3,X,W,NPOINT)
	XX=0.0
	DO 12 I=1,NPOINT
		XX=XX+W(I)*FUNC(X(I))
12	CONTINUE
	WRITE(*,'(/1X,A,F12.6)') 'Integral from GAULEG:',XX
	WRITE(*,'(1X,A,F12.6)') 'Actual value:',1.0-(1.0+X3)*EXP(-X3)
	END

	FUNCTION FUNC(X)
	FUNC=X*EXP(-X)
	END
	PROGRAM D2R1
C	Driver program for subroutine GAUSSJ
C	Reads Matrices from file MATRIX.DAT and feeds them to GAUSSJ
	PARAMETER(NP=20)
	DIMENSION A(NP,NP),B(NP,NP),AI(NP,NP),X(NP,NP)
	DIMENSION U(NP,NP),T(NP,NP)
	CHARACTER DUMMY*3
	OPEN(5,FILE='MATRX1.DAT',STATUS='OLD')
10	READ(5,'(A)') DUMMY
	IF (DUMMY.EQ.'END') GOTO 99
	READ(5,*)
	READ(5,*) N,M
	READ(5,*)
	READ(5,*) ((A(K,L), L=1,N), K=1,N)
	READ(5,*)
	READ(5,*) ((B(K,L), K=1,N), L=1,M)
C	Save Matrices for later testing of results
	DO 13 L=1,N
		DO 11 K=1,N
			AI(K,L)=A(K,L)
11		CONTINUE
		DO 12 K=1,M
			X(L,K)=B(L,K)
12		CONTINUE
13	CONTINUE
C	Invert Matrix
	CALL GAUSSJ(AI,N,NP,X,M,NP)
	WRITE(*,*) 'Inverse of Matrix A : '
	DO 14 K=1,N
		WRITE(*,'(1H,(6F12.6))') (AI(K,L), L=1,N)
14	CONTINUE
C	Test Results
C	Check Inverse
	WRITE(*,*) 'A times A-inverse (compare with unit matrix)'
	DO 17 K=1,N
		DO 16 L=1,N
			U(K,L)=0.0
			DO 15 J=1,N
				U(K,L)=U(K,L)+A(K,J)*AI(J,L)
15			CONTINUE
16		CONTINUE
		WRITE(*,'(1H,(6F12.6))') (U(K,L), L=1,N)
17	CONTINUE
C	Check Vector Solutions
	WRITE(*,*) 'Check the following vectors for equality:'
	WRITE(*,'(T12,A8,T23,A12)') 'Original','Matrix*Sol''n'
	DO 21 L=1,M
	        WRITE(*,'(1X,A,I2,A)') 'Vector ',L,':'
		DO 19 K=1,N
			T(K,L)=0.0
			DO 18 J=1,N
				T(K,L)=T(K,L)+A(K,J)*X(J,L)
18			CONTINUE
		WRITE(*,'(8X,2F12.6)') B(K,L),T(K,L)
19		CONTINUE
21	CONTINUE
	WRITE(*,*) '***********************************'
	WRITE(*,*) 'Press RETURN for next problem:'
	READ(*,*)
	GOTO 10
99	CLOSE(5)
	END	
	PROGRAM D6R9
C	Driver for routine GCF
	CHARACTER TEXT*25
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Incomplete Gamma Function') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A,T16,A,T25,A,T36,A,T47,A,T62,A)')
     *		'A','X','Actual','GCF(A,X)','GAMMLN(A)','GLN'
	DO 11 I=1,NVAL
		READ(5,*) A,X,VALUE
		IF (X.GE.A+1.0) THEN
			CALL GCF(GAMMCF,A,X,GLN)
			WRITE(*,'(1X,F6.2,5F12.6)') A,X,1.0-VALUE,
     *				GAMMCF,GAMMLN(A),GLN
		ENDIF
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D10R2
C	Driver for routine GOLDEN
	PARAMETER(TOL=1.0E-6,EQL=1.E-3)
	EXTERNAL BESSJ0
	DIMENSION AMIN(20)
	NMIN=0
	WRITE(*,'(/1X,A)') 'Minima of the function BESSJ0'
	WRITE(*,'(/1X,T6,A,T19,A,T27,A,T40,A/)') 'Min. #','X',
     *		'BESSJ0(X)','BESSJ1(X)'
	DO 12 I=1,100
		AX=I
		BX=I+1.0
		CALL MNBRAK(AX,BX,CX,FA,FB,FC,BESSJ0)
		G=GOLDEN(AX,BX,CX,BESSJ0,TOL,XMIN)
		IF (NMIN.EQ.0) THEN
			AMIN(1)=XMIN
			NMIN=1
			WRITE(*,'(1X,5X,I2,3X,3F12.6)') NMIN,XMIN,
     *				BESSJ0(XMIN),BESSJ1(XMIN)
		ELSE
			IFLAG=0
			DO 11 J=1,NMIN
				IF (ABS(XMIN-AMIN(J)).LE.EQL*XMIN)
     *					IFLAG=1
11			CONTINUE
			IF (IFLAG.EQ.0) THEN
				NMIN=NMIN+1
				AMIN(NMIN)=XMIN
				WRITE(*,'(1X,5X,I2,3X,3F12.6)') NMIN,
     *					XMIN,BESSJ0(XMIN),BESSJ1(XMIN)
			ENDIF
		ENDIF
12	CONTINUE
	END
	PROGRAM D6R8
C	Driver for routine GSER
	CHARACTER TEXT*25
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Incomplete Gamma Function') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A,T16,A,T25,A,T36,A,T47,A,T62,A)')
     *		'A','X','Actual','GSER(A,X)','GAMMLN(A)','GLN'
	DO 11 I=1,NVAL
		READ(5,*) A,X,VALUE
		CALL GSER(GAMSER,A,X,GLN)
		WRITE(*,'(1X,F6.2,5F12.6)') A,X,VALUE,GAMSER,
     *			GAMMLN(A),GLN
11	CONTINUE
	CLOSE(5)
	END
