	PROGRAM D8R10
C	Driver for routine ECLASS
	PARAMETER(N=15,M=11)
	DIMENSION LISTA(M),LISTB(M),NF(N),NFLAG(N),NSAV(N)
	DATA LISTA/1,1,5,2,6,2,7,11,3,4,12/
	DATA LISTB/5,9,13,6,10,14,3,7,15,8,4/
	CALL ECLASS(NF,N,LISTA,LISTB,M)
	DO 11 I=1,N
		NFLAG(I)=1
11	CONTINUE
	WRITE(*,'(/1X,A)') 'Numbers from 1-15 divided according to'
	WRITE(*,'(1X,A/)') 'their value modulo 4:'
	LCLAS=0
	DO 13 I=1,N
		NCLASS=NF(I)
		IF(NFLAG(NCLASS).NE.0) THEN
			NFLAG(NCLASS)=0
			LCLAS=LCLAS+1
			K=0
			DO 12 J=I,N
				IF (NF(J).EQ.NF(I)) THEN
					K=K+1
					NSAV(K)=J
				ENDIF
12			CONTINUE
			WRITE(*,'(1X,A,I2,A,3X,5I3)') 'Class',
     *				LCLAS,':',(NSAV(J),J=1,K)
		ENDIF
13	CONTINUE
	END
	PROGRAM D8R11
C	Driver for routine ECLAZZ
	EXTERNAL EQUIV
	LOGICAL EQUIV
	PARAMETER(N=15)
	DIMENSION NF(N),NFLAG(N),NSAV(N)
	CALL ECLAZZ(NF,N,EQUIV)
	DO 11 I=1,N
		NFLAG(I)=1
11	CONTINUE
	WRITE(*,'(/1X,A)') 'Numbers from 1-15 divided according to'
	WRITE(*,'(1X,A/)') 'their value modulo 4:'
	LCLAS=0
	DO 13 I=1,N
		NCLASS=NF(I)
		IF(NFLAG(NCLASS).NE.0) THEN
			NFLAG(NCLASS)=0
			LCLAS=LCLAS+1
			K=0
			DO 12 J=I,N
				IF (NF(J).EQ.NF(I)) THEN
					K=K+1
					NSAV(K)=J
				ENDIF
12			CONTINUE
			WRITE(*,'(1X,A,I2,A,3X,5I3)') 'Class',
     *				LCLAS,':',(NSAV(J),J=1,K)
		ENDIF
13	CONTINUE
	END

	LOGICAL FUNCTION EQUIV(I,J)
	EQUIV=.FALSE.
	IF (MOD(I,4).EQ.MOD(J,4)) EQUIV=.TRUE.
	RETURN
	END
	PROGRAM D11R2
C	Driver for routine EIGSRT
	PARAMETER(NP=10)
	DIMENSION D(NP),V(NP,NP),C(NP,NP)
	DATA C /5.0,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,-3.0,-4.0,
     *		4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,-3.0,
     *		3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,
     *		2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,
     *		1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,
     *		0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,
     *		-1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,
     *		-2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,
     *		-3.0,-2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,
     *		-4.0,-3.0,-2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0/
	CALL JACOBI(C,NP,NP,D,V,NROT)
	WRITE(*,*) 'Unsorted Eigenvectors:'
	DO 11 I=1,NP
		WRITE(*,'(/1X,A,I3,A,F12.6)') 'Eigenvalue',I,' =',D(I)
		WRITE(*,*) 'Eigenvector:'
		WRITE(*,'(10X,5F12.6)') (V(J,I),J=1,NP)
11	CONTINUE
	WRITE(*,'(//,A,//)') '****** sorting ******'
	CALL EIGSRT(D,V,NP,NP)
	WRITE(*,*) 'Sorted Eigenvectors:'
	DO 12 I=1,NP
		WRITE(*,'(/1X,A,I3,A,F12.6)') 'Eigenvalue',I,' =',D(I)
		WRITE(*,*) 'Eigenvector:'
		WRITE(*,'(10X,5F12.6)') (V(J,I),J=1,NP)
12	CONTINUE
	END
	PROGRAM D6R28
C	Driver for routine EL2
	COMMON AKC,A,B
	EXTERNAL FUNC
	WRITE(*,*) 'General Elliptic Integral of Second Kind'
	WRITE(*,'(1X,T8,A,T17,A,T28,A,T38,A,T47,A,T54,A)')
     *		'x','kc','a','b','EL2','Integral'
	IDUM=-55
	AGO=0.0
	DO 11 I=1,20
		AKC=5.0*RAN3(IDUM)	
		A=10.0*RAN3(IDUM)
		B=10.0*RAN3(IDUM)
		X=10.0*RAN3(IDUM)	
		ASTOP=ATAN(X)
		CALL QSIMP(FUNC,AGO,ASTOP,S)
		WRITE(*,'(1X,6F10.6)')
     *			X,AKC,A,B,EL2(X,AKC,A,B),S
11	CONTINUE
	END

	FUNCTION FUNC(PHI)
	COMMON AKC,A,B
	TN=TAN(PHI)
	TSQ=TN*TN
	FUNC=(A+B*TSQ)/SQRT((1.0+TSQ)*(1.0+AKC*AKC*TSQ))
	END
	PROGRAM D11R6
C	Driver for ELMHES
	PARAMETER(NP=5)
	DIMENSION A(NP,NP),R(NP),C(NP)
	DATA A/1.0,2.0,3.0,4.0,5.0,2.0,3.0,4.0,5.0,6.0,
     *		300.0,400.0,5.0,600.0,700.0,4.0,5.0,6.0,7.0,8.0,
     *		5.0,6.0,7.0,8.0,9.0/
	WRITE(*,'(/1X,A/)') '***** Original Matrix *****'
	DO 11 I=1,NP
		WRITE(*,'(1X,5F12.2)') (A(I,J),J=1,NP)
11	CONTINUE
	WRITE(*,'(/1X,A/)') '***** Balance Matrix *****'
	CALL BALANC(A,NP,NP)
	DO 12 I=1,NP
		WRITE(*,'(1X,5F12.2)') (A(I,J),J=1,NP)
12	CONTINUE
	WRITE(*,'(/1X,A/)') '***** Reduce to Hessenberg Form *****'
	CALL ELMHES(A,NP,NP)
	DO 14 J=1,NP-2
		DO 13 I=J+2,NP
			A(I,J)=0.0
13		CONTINUE
14	CONTINUE
	DO 15 I=1,NP
		WRITE(*,'(1X,5E12.4)') (A(I,J),J=1,NP)
15	CONTINUE
	END
	PROGRAM D6R10
C	Driver for routine ERF
	CHARACTER TEXT*14
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Error Function') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A1,T12,A6,T24,A6)')
     *		'X','Actual','ERF(X)'
	DO 11 I=1,NVAL
		READ(5,*) X,VALUE
		WRITE(*,'(F6.2,2F12.7)') X,VALUE,ERF(X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R11
C	Driver for routine ERFC
	CHARACTER TEXT*14
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Error Function') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) 'Complementary error function'
	WRITE(*,'(1X,T5,A1,T12,A6,T23,A7)')
     *		'X','Actual','ERFC(X)'
	DO 11 I=1,NVAL
		READ(5,*) X,VALUE
		VALUE=1.0-VALUE
		WRITE(*,'(F6.2,2F12.7)') X,VALUE,ERFC(X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R12
C	Driver for routine ERFCC
	CHARACTER TEXT*14
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Error Function') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) 'Complementary error function'
	WRITE(*,'(1X,T5,A1,T12,A6,T23,A8)')
     *		'X','Actual','ERFCC(X)'
	DO 11 I=1,NVAL
		READ(5,*) X,VALUE
		VALUE=1.0-VALUE
		WRITE(*,'(F6.2,2F12.7)') X,VALUE,ERFCC(X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D5R1
C	Driver for routine EULSUM
	PARAMETER (NVAL=40)
	DIMENSION WKSP(NVAL)
C	Evaluate ln(1+x)=x-x^2/2+x^3/3-x^4/4...  for -1<x<1
10	WRITE(*,*) 'How many terms in polynomial?'	
	WRITE(*,'(1X,A,I2,A)') 'Enter n between 1 and ',NVAL,
     *				'. Enter n=0 to end.'
	READ(*,*) MVAL
	IF ((MVAL.LE.0).OR.(MVAL.GT.NVAL)) STOP
	WRITE(*,'(1X,T9,A1,T18,A6,T28,A10)') 'X','Actual','Polynomial'
	DO 12 I=-8,8,1
		X=I/10.0
		SUM=0.0
		XPOWER=-1
		DO 11 J=1,MVAL
			XPOWER=-X*XPOWER
			TERM=XPOWER/J
			CALL EULSUM(SUM,TERM,J,WKSP)
11		CONTINUE
		WRITE(*,'(3F12.6)') X,LOG(1.0+X),SUM
12	CONTINUE
	GOTO 10
	END
	PROGRAM D12R10
C	Driver for routine EVLMEM
	PARAMETER(N=1000,M=10,NFDT=16)
	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(*,*) 'Power spectrum estimate of data in SPCTRL.DAT'
	WRITE(*,'(1X,T6,A,T20,A)') 'f*delta','power'
	DO 11 I=0,NFDT
		FDT=0.5*I/NFDT
		WRITE(*,'(1X,2F12.6)') FDT,EVLMEM(FDT,COF,M,PM)
11	CONTINUE
	END
	PROGRAM D7R5
C	Driver for routine EXPDEV
	PARAMETER(NPTS=10000,EE=2.718281828)
	DIMENSION TRIG(21),X(21)
	DO 11 I=1,21
		TRIG(I)=(I-1)/20.0
		X(I)=0.0
11	CONTINUE
	IDUM=-1
	DO 13 I=1,NPTS
		Y=EXPDEV(IDUM)
		DO 12 J=2,21
			IF((Y.LT.TRIG(J)).AND.(Y.GT.TRIG(J-1))) THEN
				X(J)=X(J)+1.0
			ENDIF
12		CONTINUE
13	CONTINUE
	TOTAL=0.0
	DO 14 I=2,21
		TOTAL=TOTAL+X(I)
14	CONTINUE
	WRITE(*,'(1X,A,I6,A)') 'Exponential distribution with',
     *		NPTS,' points:'
	WRITE(*,'(1X,T5,A,T19,A,T31,A)')
     *		'interval','observed','expected'
	DO 15 I=2,21
		X(I)=X(I)/TOTAL
		EXPECT=EXP(-(TRIG(I-1)+TRIG(I))/2.0)
		EXPECT=EXPECT*0.05*EE/(EE-1)
		WRITE(*,'(1X,2F6.2,2F12.4)')
     *			TRIG(I-1),TRIG(I),X(I),EXPECT
15	CONTINUE
	END
