	PROGRAM D5R9
C	Driver for routine PCSHFT
	PARAMETER(NVAL=40, PIO2=1.5707963)
	EXTERNAL FUNC
	DIMENSION C(NVAL),D(NVAL)
	A=-PIO2
	B=PIO2
	CALL CHEBFT(A,B,C,NVAL,FUNC)
10	WRITE(*,*) 'How many terms in Chebyshev evaluation?'
	WRITE(*,'(1X,A,I2,A)') 'Enter n between 6 and ',NVAL,
     *				'. Enter n=0 to end.'
	READ(*,*) MVAL
	IF ((MVAL.LE.0).OR.(MVAL.GT.NVAL)) GOTO 20
	CALL CHEBPC(C,D,MVAL)
	CALL PCSHFT(A,B,D,MVAL)
C	Test shifted polynomial
	WRITE(*,'(1X,T10,A,T19,A,T29,A)') 'X','Actual','Polynomial'
	DO 12 I=-8,8,1
		X=I*PIO2/10.0
		POLY=D(MVAL)
		DO 11 J=MVAL-1,1,-1
			POLY=POLY*X+D(J)
11		CONTINUE
		WRITE(*,'(1X,3F12.6)') X,FUNC(X),POLY
12	CONTINUE
	GOTO 10
20	END

	FUNCTION FUNC(X)
	FUNC=(X**2)*(X**2-2.0)*SIN(X)
	END
	PROGRAM D13R15
C	Driver for routine PEARSN
	DIMENSION DOSE(10),SPORE(10)
	DATA DOSE/56.1,64.1,70.0,66.6,82.,91.3,90.,99.7,115.3,110./
	DATA SPORE/0.11,0.4,0.37,0.48,0.75,0.66,0.71,1.2,1.01,0.95/
	WRITE(*,'(1X,A)') 
     *		'Effect of Gamma Rays on Man-in-the-Moon Marigolds'
	WRITE(*,'(1X,A,T29,A)') 'Count Rate (cpm)','Pollen Index'
	DO 11 I=1,10
		WRITE(*,'(1X,F10.2,F25.2)') DOSE(I),SPORE(I)
11	CONTINUE
	CALL PEARSN(DOSE,SPORE,10,R,PROB,Z)	
	WRITE(*,'(/1X,T24,A,T38,A)') 'PEARSN','Expected'
	WRITE(*,'(1X,A,T18,2E15.6)') 'Corr. Coeff.',R,0.906959
	WRITE(*,'(1X,A,T18,2E15.6)') 'Probability',PROB,0.292650E-3
	WRITE(*,'(1X,A,T18,2E15.6/)') 'Fisher''s Z',Z,1.51011
	END
	PROGRAM D8R2
C	Driver for routine PIKSR2
	DIMENSION A(100),B(100)
	OPEN(5,FILE='TARRAY.DAT',STATUS='OLD')
	READ(5,*) (A(I),I=1,100)
	CLOSE(5)
C	Generate B-array
	DO 11 I=1,100
		B(I)=I
11	CONTINUE
C	Sort A and mix B
	CALL PIKSR2(100,A,B)
	WRITE(*,*) 'After sorting A and mixing B, array A is:'
	DO 12 I=1,10
		WRITE(*,'(1X,10F6.2)') (A(10*(I-1)+J), J=1,10)
12	CONTINUE
	WRITE(*,*) '...and array B is:'
	DO 13 I=1,10
		WRITE(*,'(1X,10F6.2)') (B(10*(I-1)+J), J=1,10)
13	CONTINUE
	WRITE(*,*) 'press RETURN to continue...'
	READ(*,*)
C	Sort B and mix A
	CALL PIKSR2(100,B,A)
	WRITE(*,*) 'After sorting B and mixing A, array A is:'
	DO 14 I=1,10
		WRITE(*,'(1X,10F6.2)') (A(10*(I-1)+J), J=1,10)
14	CONTINUE
	WRITE(*,*) '...and array B is:'
	DO 15 I=1,10
		WRITE(*,'(1X,10F6.2)') (B(10*(I-1)+J), J=1,10)
15	CONTINUE
	END
	PROGRAM D8R1
C	Driver for routine PIKSRT
	DIMENSION A(100)
	OPEN(5,FILE='TARRAY.DAT',STATUS='OLD')
	READ(5,*) (A(I),I=1,100)
	CLOSE(5)
C	Print original array
	WRITE(*,*) 'Original array:'
	DO 11 I=1,10
		WRITE(*,'(1X,10F6.2)') (A(10*(I-1)+J),J=1,10)
11	CONTINUE
C	Sort array
	CALL PIKSRT(100,A)
C	Print sorted array
	WRITE(*,*) 'Sorted array:'
	DO 12 I=1,10
		WRITE(*,'(1X,10F6.2)') (A(10*(I-1)+J),J=1,10)
12	CONTINUE
	END
	PROGRAM D6R27
C	Driver for routine PLGNDR
	CHARACTER TEXT*27
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)',ERR=99) TEXT
	IF (TEXT.EQ.'Legendre Polynomials') THEN
		READ(5,*) NVAL
		WRITE(*,*) TEXT
		WRITE(*,'(1X,T5,A,T9,A,T20,A,T35,A,T49,A)')
     *			'N','M','X','Actual','PLGNDR(N,M,X)'
		DO 12 I=1,NVAL
			READ(5,*) N,M,X,VALUE
			FAC=1.0
			IF (M.GT.0) THEN
				DO 11 J=N-M+1,N+M
					FAC=FAC*J
11				CONTINUE
			ENDIF
			FAC=2.0*FAC/(2.0*N+1.0)
			VALUE=VALUE*SQRT(FAC)
			WRITE(*,'(1X,2I4,3E17.6)') N,M,X,VALUE,PLGNDR(N,M,X)
12		CONTINUE
	ENDIF
	GOTO 10
99	CLOSE(5)
	END
	PROGRAM D7R8
C	Driver for routine POIDEV
	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(*,*) 'Mean of Poisson distrib. (x=0 to 20); neg. to end'
	READ(*,*) XM
	IF (XM.LT.0.0) GOTO 99
	IF (IA.GT.20.0) GOTO 10		
	DO 12 I=1,NPTS
		J=INT(POIDEV(XM,IDUM))+1
		IF ((J.GE.1).AND.(J.LE.21)) DIST(J)=DIST(J)+1
12	CONTINUE
	WRITE(*,'(1X,A,F5.2,A,I6,A)')
     *		'Poisson-distributed deviate, mean ',
     *		XM,' 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 D3R7
C	Driver for routine POLCOE
	PARAMETER(NP=5,PI=3.1415926)
	DIMENSION XA(NP),YA(NP),COEFF(NP)
	DO 15 NFUNC=1,2
		IF (NFUNC.EQ.1) THEN
			WRITE(*,*) 'Sine function from 0 to PI'
			DO 11 I=1,NP
				XA(I)=I*PI/NP
				YA(I)=SIN(XA(I))
11			CONTINUE
		ELSE IF (NFUNC.EQ.2) THEN
			WRITE(*,*) 'Exponential function from 0 to 1'
			DO 12 I=1,NP
				XA(I)=1.0*I/NP
				YA(I)=EXP(XA(I))
12			CONTINUE
		ELSE
			STOP
		ENDIF
		CALL POLCOE(XA,YA,NP,COEFF)	
		WRITE(*,*) '    coefficients'
		WRITE(*,'(1X,6F12.6)') (COEFF(I),I=1,NP)
		WRITE(*,'(1X,T10,A1,T20,A4,T29,A10)')
     *			'x','f(x)','polynomial'
		DO 14 I=1,10
			IF (NFUNC.EQ.1) THEN
				X=(-0.05+I/10.0)*PI
				F=SIN(X)
			ELSE IF (NFUNC.EQ.2) THEN
				X=-0.05+I/10.0
				F=EXP(X)
			ENDIF
			SUM=COEFF(NP)
			DO 13 J=NP-1,1,-1
				SUM=COEFF(J)+SUM*X
13			CONTINUE
			WRITE(*,'(1X,3F12.6)') X,F,SUM
14		CONTINUE
		WRITE(*,*) '***********************************'
		WRITE(*,*) 'Press RETURN'
		READ(*,*)
15	CONTINUE	   
	END
	PROGRAM D3R8
C	Driver for routine POLCOF
	PARAMETER(NP=5,PI=3.141593)
	DIMENSION XA(NP),YA(NP),COEFF(NP)
	DO 15 NFUNC=1,2
		IF (NFUNC.EQ.1) THEN
			WRITE(*,*) 'Sine function from 0 to PI'
			DO 11 I=1,NP
				XA(I)=I*PI/NP
				YA(I)=SIN(XA(I))
11			CONTINUE
		ELSE IF (NFUNC.EQ.2) THEN
			WRITE(*,*) 'Exponential function from 0 to 1'
			DO 12 I=1,NP
				XA(I)=1.0*I/NP
				YA(I)=EXP(XA(I))
12			CONTINUE
		ELSE
			STOP
		ENDIF
		CALL POLCOF(XA,YA,NP,COEFF)	
		WRITE(*,*) '    coefficients'
		WRITE(*,'(1X,6F12.6)') (COEFF(I),I=1,NP)
		WRITE(*,'(1X,T10,A1,T20,A4,T29,A10)')
     *			'x','f(x)','polynomial'
		DO 14 I=1,10
			IF (NFUNC.EQ.1) THEN
				X=(-0.05+I/10.0)*PI
				F=SIN(X)
			ELSE IF (NFUNC.EQ.2) THEN
				X=-0.05+I/10.0
				F=EXP(X)
			ENDIF
			SUM=COEFF(NP)
			DO 13 J=NP-1,1,-1
				SUM=COEFF(J)+SUM*X
13			CONTINUE
			WRITE(*,'(1X,3F12.6)') X,F,SUM
14		CONTINUE
		WRITE(*,*) '***********************************'
		WRITE(*,*) 'Press RETURN'
		READ(*,*)
15	CONTINUE
	END
	PROGRAM D5R3
C	Driver for routine POLDIV
C	(X-1)**5/(X+1)**3
	PARAMETER(N=6,NV=4)
	DIMENSION U(N),V(NV),Q(N),R(N)
	DATA U/-1.0,5.0,-10.0,10.0,-5.0,1.0/
	DATA V/1.0,3.0,3.0,1.0/
	CALL POLDIV(U,N,V,NV,Q,R)
	WRITE(*,'(//1X,6(7X,A)/)')'X^0','X^1','X^2','X^3','X^4','X^5'
	WRITE(*,*) 'Quotient polynomial coefficients:'
	WRITE(*,'(1X,6F10.2/)') (Q(I),I=1,6)
	WRITE(*,*) 'Expected quotient coefficients:'
	WRITE(*,'(1X,6F10.2///)') 31.0,-8.0,1.0,0.0,0.0,0.0
	WRITE(*,*) 'Remainder polynomial coefficients:'
	WRITE(*,'(1X,4F10.2/)') (R(I),I=1,4)
	WRITE(*,*) 'Expected remainder coefficients:'
	WRITE(*,'(1X,4F10.2//)') -32.0,-80.0,-80.0,0.0
	END
	PROGRAM D3R9
C	Driver for routine POLIN2
	PARAMETER(N=5,PI=3.141593)
	DIMENSION X1A(N),X2A(N),YA(N,N)
	DO 12 I=1,N
		X1A(I)=I*PI/N
		DO 11 J=1,N
			X2A(J)=1.0*J/N
			YA(I,J)=SIN(X1A(I))*EXP(X2A(J))
11		CONTINUE
12  	CONTINUE
C	Test 2-dimensional interpolation
	WRITE(*,'(T9,A,T21,A,T32,A,T40,A,T58,A)')
     *  'x1','x2','f(x)','interpolated','error'
	DO 14 I=1,4
	X1=(-0.1+I/5.0)*PI
		DO 13 J=1,4
			X2=-0.1+J/5.0
			F=SIN(X1)*EXP(X2)
			CALL POLIN2(X1A,X2A,YA,N,N,X1,X2,Y,DY)
			WRITE(*,'(1X,4F12.6,F14.6)') X1,X2,F,Y,DY
13		CONTINUE
		WRITE(*,*) '***********************************'
14	CONTINUE
	END
	PROGRAM D3R1
C	Driver for routine POLINT
	PARAMETER(NP=10,PI=3.1415926)
	DIMENSION XA(NP),YA(NP)
	WRITE(*,*) 'Generation of interpolation tables'
	WRITE(*,*) ' ... sin(x)    0<x<pi'
	WRITE(*,*) ' ... exp(x)    0<x<1 '
	WRITE(*,*) 'How many entries go in these tables? (note: N<10)'
	READ(*,*) N
	DO 14 NFUNC=1,2
		IF (NFUNC.EQ.1) THEN
			WRITE(*,*) 'sine function from 0 to pi'
			DO 11 I=1,N
				XA(I)=I*PI/N
				YA(I)=SIN(XA(I))
11			CONTINUE
		ELSE IF (NFUNC.EQ.2) THEN
			WRITE(*,*) 'exponential function from 0 to 1'
			DO 12 I=1,N
				XA(I)=I*1.0/N
				YA(I)=EXP(XA(I))
12			CONTINUE
		ELSE
			STOP
		ENDIF
		WRITE(*,'(T10,A1,T20,A4,T28,A12,T46,A5)')
     *		      	'x','f(x)','interpolated','error'
		DO 13 I=1,10
			IF (NFUNC.EQ.1) THEN
				X=(-0.05+I/10.0)*PI
				F=SIN(X)
			ELSE IF (NFUNC.EQ.2) THEN
				X=(-0.05+I/10.0)
				F=EXP(X)
			ENDIF
			CALL POLINT(XA,YA,N,X,Y,DY)
			WRITE(*,'(1X,3F12.6,E15.4)') X,F,Y,DY
13		CONTINUE
		WRITE(*,*) '***********************************'
		WRITE(*,*) 'Press RETURN'
		READ(*,*)
14	CONTINUE
	END
	PROGRAM D10R6
C	Driver for routine POWELL
	PARAMETER(NDIM=3,FTOL=1.0E-6)
	DIMENSION P(NDIM),XI(NDIM,NDIM)
	NP=NDIM
	DATA XI/1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0/
	DATA P/1.5,1.5,2.5/
	CALL POWELL(P,XI,NDIM,NP,FTOL,ITER,FRET)
	WRITE(*,'(/1X,A,I3)') 'Iterations:',ITER
	WRITE(*,'(/1X,A/1X,3F12.6)') 'Minimum found at: ',(P(I),I=1,NDIM)
	WRITE(*,'(/1X,A,F12.6)') 'Minimum function value =',FRET
	WRITE(*,'(/1X,A)') 'True minimum of function is at:'
	WRITE(*,'(1X,3F12.6/)') 1.0,2.0,3.0
	END

	FUNCTION FUNC(X)
	DIMENSION X(3)
	FUNC=0.5-BESSJ0((X(1)-1.0)**2+(X(2)-2.0)**2+(X(3)-3.0)**2)
	END
	PROGRAM D12R12
C	Driver for routine PREDIC
	PARAMETER(NPTS=500,NPOLES=10,NFUT=20,PI=3.1415926)
	DIMENSION DATA(NPTS),D(NPOLES),WK1(NPTS),
     *		WK2(NPTS),WKM(NPOLES),FUTURE(NFUT)
	F(N)=EXP(-1.0*N/NPTS)*SIN(2.0*PI*N/50.0)
     *		+EXP(-2.0*N/NPTS)*SIN(2.2*PI*N/50.0)
	DO 11 I=1,NPTS
		DATA(I)=F(I)
11	CONTINUE
	CALL MEMCOF(DATA,NPTS,NPOLES,DUM,D,WK1,WK2,WKM)
	CALL FIXRTS(D,NPOLES)
	CALL PREDIC(DATA,NPTS,D,NPOLES,FUTURE,NFUT)
	WRITE(*,'(6X,A,T13,A,T25,A)') 'I','Actual','PREDIC'
	DO 12 I=1,NFUT
		WRITE(*,'(1X,I6,2F12.6)') I,F(I+NPTS),FUTURE(I)
12	CONTINUE
	END
	PROGRAM D13R12
C	Driver for routine PROBKS
	CHARACTER TEXT(50)*1
	WRITE(*,*) 'Probability func. for Kolmogorov-Smirnov statistic'
	WRITE(*,'(/1X,T3,A,T15,A,T27,A)') 'Lambda:','Value:','Graph:'
	NPTS=20
	EPS=0.1
	SCALE=40.0
	DO 12 I=1,NPTS
		ALAM=I*EPS
		VALUE=PROBKS(ALAM)
		TEXT(1)='*'
		DO 11 J=1,50
			IF(J.LE.NINT(SCALE*VALUE)) THEN
				TEXT(J)='*'
			ELSE
				TEXT(J)=' '
			ENDIF
11		CONTINUE
		WRITE(*,'(1X,F9.6,F12.6,4X,50A1)') ALAM,VALUE,
     *			(TEXT(J),J=1,50)
12	CONTINUE
	END
	PROGRAM D15R9
C	Driver for routine PZEXTR
C	Feed values from a rational function
C	Fn(x)=(1-x+x**3)/(x+1)**n
	PARAMETER(NV=4,NUSE=5)
	DIMENSION YEST(NV),YZ(NV),DY(NV)
	DO 12 I=1,10
		IEST=I
		XEST=1.0/FLOAT(I)
		DUM=1.0-XEST+XEST*XEST*XEST
		DO 11 J=1,NV
			DUM=DUM/(XEST+1.0)
			YEST(J)=DUM
11		CONTINUE
		CALL PZEXTR(IEST,XEST,YEST,YZ,DY,NV,NUSE)
		WRITE(*,'(/1X,A,I2)') 'I = ',I
		WRITE(*,'(1X,A,4F12.6)') 'Extrap. function:',(YZ(J),J=1,NV)
		WRITE(*,'(1X,A,4F12.6)') 'Estimated error: ',(DY(J),J=1,NV)
12	CONTINUE
	WRITE(*,'(/1X,A,4F12.6)') 'Actual values:   ',1.0,1.0,1.0,1.0
	END

