	PROGRAM D7R1
C	Driver for routine RAN0
C	Calculates pi statistically using volume of unit n-sphere
	PARAMETER(PI=3.1415926)
	DIMENSION IY(3),YPROB(3)
	FNC(X1,X2,X3,X4)=SQRT(X1**2+X2**2+X3**2+X4**2)
	IDUM=-1
	DO 11 I=1,3
		IY(I)=0
11	CONTINUE
	WRITE(*,'(1X,/,T15,A)') 'Volume of unit n-sphere, n=2,3,4'
	WRITE(*,'(1X,/,T3,A,T17,A,T26,A,T37,A)')
     *		'# points','pi','(4/3)*pi','(1/2)*pi^2'
	DO 14 J=1,15
		DO 12 K=2**(J-1),2**J
			X1=RAN0(IDUM)
			X2=RAN0(IDUM)
			X3=RAN0(IDUM)
			X4=RAN0(IDUM)
			IF(FNC(X1,X2,0.0,0.0).LT.1.0) IY(1)=IY(1)+1
			IF(FNC(X1,X2,X3,0.0).LT.1.0) IY(2)=IY(2)+1
			IF(FNC(X1,X2,X3,X4).LT.1.0) IY(3)=IY(3)+1
12		CONTINUE
		DO 13 I=1,3
			YPROB(I)=1.0*(2**(I+1))*IY(I)/(2**J)
13		CONTINUE
		WRITE(*,'(1X,I8,3F12.6)') 2**J,(YPROB(I),I=1,3)
14	CONTINUE
	WRITE(*,'(1X,/,T4,A,3F12.6,/)') 'actual',PI,4.0*PI/3.0,0.5*(PI**2)
	END
	PROGRAM D7R2
C	Driver for routine RAN1
C	Calculates pi statistically using volume of unit n-sphere
	PARAMETER(PI=3.1415926)
	DIMENSION IY(3),YPROB(3)
	FNC(X1,X2,X3,X4)=SQRT(X1**2+X2**2+X3**2+X4**2)
	IDUM=-1
	DO 11 I=1,3
		IY(I)=0
11	CONTINUE
	WRITE(*,'(1X,/,T15,A)') 'Volume of unit n-sphere, n=2,3,4'
	WRITE(*,'(1X,/,T3,A,T17,A,T26,A,T37,A)')
     *		'# points','pi','(4/3)*pi','(1/2)*pi^2'
	DO 14 J=1,15
		DO 12 K=2**(J-1),2**J
			X1=RAN1(IDUM)
			X2=RAN1(IDUM)
			X3=RAN1(IDUM)
			X4=RAN1(IDUM)
			IF(FNC(X1,X2,0.0,0.0).LT.1.0) IY(1)=IY(1)+1
			IF(FNC(X1,X2,X3,0.0).LT.1.0) IY(2)=IY(2)+1
			IF(FNC(X1,X2,X3,X4).LT.1.0) IY(3)=IY(3)+1
12		CONTINUE
		DO 13 I=1,3
			YPROB(I)=1.0*(2**(I+1))*IY(I)/(2**J)
13		CONTINUE
		WRITE(*,'(1X,I8,3F12.6)') 2**J,(YPROB(I),I=1,3)
14	CONTINUE
	WRITE(*,'(1X,/,T4,A,3F12.6,/)') 'actual',PI,4.0*PI/3.0,0.5*(PI**2)
	END
	PROGRAM D7R3
C	Driver for routine RAN2
C	Calculates pi statistically using volume of unit n-sphere
	PARAMETER(PI=3.1415926)
	DIMENSION IY(3),YPROB(3)
	FNC(X1,X2,X3,X4)=SQRT(X1**2+X2**2+X3**2+X4**2)
	IDUM=-1
	DO 11 I=1,3
		IY(I)=0
11	CONTINUE
	WRITE(*,'(1X,/,T15,A)') 'Volume of unit n-sphere, n=2,3,4'
	WRITE(*,'(1X,/,T3,A,T17,A,T26,A,T37,A)')
     *		'# points','pi','(4/3)*pi','(1/2)*pi^2'
	DO 14 J=1,15
		DO 12 K=2**(J-1),2**J
			X1=RAN2(IDUM)
			X2=RAN2(IDUM)
			X3=RAN2(IDUM)
			X4=RAN2(IDUM)
			IF(FNC(X1,X2,0.0,0.0).LT.1.0) IY(1)=IY(1)+1
			IF(FNC(X1,X2,X3,0.0).LT.1.0) IY(2)=IY(2)+1
			IF(FNC(X1,X2,X3,X4).LT.1.0) IY(3)=IY(3)+1
12		CONTINUE
		DO 13 I=1,3
			YPROB(I)=1.0*(2**(I+1))*IY(I)/(2**J)
13		CONTINUE
		WRITE(*,'(1X,I8,3F12.6)') 2**J,(YPROB(I),I=1,3)
14	CONTINUE
	WRITE(*,'(1X,/,T4,A,3F12.6,/)') 'actual',PI,4.0*PI/3.0,0.5*(PI**2)
	END
	PROGRAM D7R4
C	Driver for routine RAN3
C	Calculates pi statistically using volume of unit n-sphere
	PARAMETER(PI=3.1415926)
	DIMENSION IY(3),YPROB(3)
	FNC(X1,X2,X3,X4)=SQRT(X1**2+X2**2+X3**2+X4**2)
	IDUM=-1
	DO 11 I=1,3
		IY(I)=0
11	CONTINUE
	WRITE(*,'(1X,/,T15,A)') 'Volume of unit n-sphere, n=2,3,4'
	WRITE(*,'(1X,/,T3,A,T17,A,T26,A,T37,A)')
     *		'# points','pi','(4/3)*pi','(1/2)*pi^2'
	DO 14 J=1,15
		DO 12 K=2**(J-1),2**J
			X1=RAN3(IDUM)
			X2=RAN3(IDUM)
			X3=RAN3(IDUM)
			X4=RAN3(IDUM)
			IF(FNC(X1,X2,0.0,0.0).LT.1.0) IY(1)=IY(1)+1
			IF(FNC(X1,X2,X3,0.0).LT.1.0) IY(2)=IY(2)+1
			IF(FNC(X1,X2,X3,X4).LT.1.0) IY(3)=IY(3)+1
12		CONTINUE
		DO 13 I=1,3
			YPROB(I)=1.0*(2**(I+1))*IY(I)/(2**J)
13		CONTINUE
		WRITE(*,'(1X,I8,3F12.6)') 2**J,(YPROB(I),I=1,3)
14	CONTINUE
	WRITE(*,'(1X,/,T4,A,3F12.6,/)') 'actual',PI,4.0*PI/3.0,0.5*(PI**2)
	END
	PROGRAM D7R12
C	Driver for routine RAN4
	PARAMETER(NPT=50)
	DIMENSION Y(NPT)
	IDUM=-123
	AVE=0.0
	WRITE(*,'(1X,A,I5/)') 'First 10 random numbers with IDUM =',IDUM
	WRITE(*,'(1X,A4,A11)') '#','RAN4'
	DO 11 J=1,10
		WRITE(*,'(1X,I4,F12.6)') J,RAN4(IDUM)
11	CONTINUE
	WRITE(*,'(/1X,A,I3)') 'Average and Variance of next ',NPT
	DO 12 J=1,NPT
		Y(J)=RAN4(IDUM)
12	CONTINUE
	CALL AVEVAR(Y,NPT,AVE,VAR)
	WRITE(*,'(/1X,A,F10.4)') 'Average: ',AVE
	WRITE(*,'(1X,A,F10.4/)') 'Variance:',VAR
	WRITE(*,'(/1X,A)') 'Expected Result for an Infinite Sample:'
	WRITE(*,'(/1X,A,F10.4)') 'Average: ',0.5
	WRITE(*,'(1X,A,F10.4/)') 'Variance:',1./12.
	END
	PROGRAM D8R7
C	Driver for routine RANK
	DIMENSION A(100),B(10),INDX(100),IRANK(100)
	OPEN(5,FILE='TARRAY.DAT',STATUS='OLD')
	READ(5,*) (A(I),I=1,100)
	CLOSE(5)
	CALL INDEXX(100,A,INDX)
	CALL RANK(100,INDX,IRANK)
	WRITE(*,*) 'Original array is:'
	DO 11 I=1,10
		WRITE(*,'(1X,10F6.2)') (A(10*(I-1)+J), J=1,10)
11	CONTINUE
	WRITE(*,*) 'Table of ranks is:'
	DO 12 I=1,10
		WRITE(*,'(1X,10I6)') (IRANK(10*(I-1)+J), J=1,10)
12	CONTINUE
	WRITE(*,*) 'press RETURN to continue...'
	READ(*,*)
	WRITE(*,*) 'Array sorted according to rank table:'
	DO 15 I=1,10
		DO 14 J=1,10
			K=10*(I-1)+J
			DO 13 L=1,100
				IF (IRANK(L).EQ.K) B(J)=A(L)
13			CONTINUE
14		CONTINUE
		WRITE(*,'(1X,10F6.2)') (B(J),J=1,10)
15	CONTINUE
	END
	PROGRAM D3R2
C	Driver for routine RATINT
	PARAMETER(NPT=6,EPSSQ=1.0)
	DIMENSION X(NPT),Y(NPT)
	F(X)=X*EXP(-X)/((X-1.0)**2+EPSSQ)
	DO 11 I=1,NPT
		X(I)=I*2.0/NPT
		Y(I)=F(X(I))
11	CONTINUE
	WRITE(*,'(/1X,A/)') 'Diagonal rational function interpolation'
	WRITE(*,'(1X,T6,A,T13,A,T26,A,T40,A)') 
     *		'x','interp.','accuracy','actual'
	DO 12 I=1,10
		XX=0.2*I
		CALL RATINT(X,Y,NPT,XX,YY,DYY)
		YEXP=F(XX)
		WRITE(*,'(1X,F6.2,F12.6,E15.4,F12.6)') XX,YY,DYY,YEXP
12	CONTINUE
	END
	PROGRAM D12R3
C	Driver for routine REALFT
	PARAMETER(EPS=1.0E-3,NP=32,NPP2=NP+2,WIDTH=50.0,PI=3.14159)
	DIMENSION DATA(NPP2),SIZE(NP)
	N=NP/2
1	WRITE(*,'(1X,A,I2,A)') 'Period of sinusoid in channels (2-',
     *		NP,', OR 0 TO STOP)'
	READ(*,*) PER
	IF(PER.LE.0.)STOP
	DO 11 I=1,NP
		DATA(I)=COS(2.0*PI*(I-1)/PER)
11	CONTINUE
	CALL REALFT(DATA,N,+1)
	BIG=-1.0E10
	DO 12 I=1,N
		SIZE(I)=SQRT(DATA(2*I-1)**2+DATA(2*I)**2)
		IF (I.EQ.1) SIZE(I)=DATA(I)
		IF(SIZE(I).GT.BIG) BIG=SIZE(I)
12	CONTINUE
	SCAL=WIDTH/BIG
	DO 13 I=1,N
		NLIM=SCAL*SIZE(I)+EPS
		WRITE(*,'(1X,I4,1X,60A1)') I,('*',J=1,NLIM+1)
13	CONTINUE
	WRITE(*,*) 'press continue ...'
	READ(*,*)
	CALL REALFT(DATA,N,-1)
	BIG=-1.0E10
	SMALL=1.0E10
	DO 14 I=1,NP
		IF(DATA(I).LT.SMALL) SMALL=DATA(I)
		IF(DATA(I).GT.BIG) BIG=DATA(I)
14	CONTINUE
	SCAL=WIDTH/(BIG-SMALL)
	DO 15 I=1,NP
		NLIM=SCAL*(DATA(I)-SMALL)+EPS
		WRITE(*,'(1X,I4,1X,60A1)') I,
     *		  ('*',J=1,NLIM+1)
15	CONTINUE
	GOTO 1
	END
	PROGRAM D15R1
C	Driver for routine RK4
	EXTERNAL DERIVS
	PARAMETER(N=4)
	DIMENSION Y(N),DYDX(N),YOUT(N)
	X=1.0
	Y(1)=BESSJ0(X)
	Y(2)=BESSJ1(X)
	Y(3)=BESSJ(2,X)
	Y(4)=BESSJ(3,X)
	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)
	WRITE(*,'(/1X,A,T19,A,T31,A,T43,A,T55,A)') 
     *		'Bessel Function:','J0','J1','J3','J4'
	DO 11 I=1,5
		H=0.2*I
		CALL RK4(Y,DYDX,N,X,H,YOUT,DERIVS)
		WRITE(*,'(/1X,A,F6.2)') 'For a step size of:',H
		WRITE(*,'(1X,A10,4F12.6)') 'RK4:',(YOUT(J),J=1,4)
		WRITE(*,'(1X,A10,4F12.6)') 'Actual:',BESSJ0(X+H),
     *				BESSJ1(X+H),BESSJ(2,X+H),BESSJ(3,X+H)
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)
	RETURN
	END
	PROGRAM D15R2
C	Driver for routine RKDUMB
	PARAMETER(NVAR=4,NSTEP=150)
	DIMENSION VSTART(NVAR)
	COMMON /PATH/X(200),Y(10,200)
	EXTERNAL DERIVS
	X1=1.0
	VSTART(1)=BESSJ0(X1)
	VSTART(2)=BESSJ1(X1)
	VSTART(3)=BESSJ(2,X1)
	VSTART(4)=BESSJ(3,X1)
	X2=20.0
	CALL RKDUMB(VSTART,NVAR,X1,X2,NSTEP,DERIVS)
	WRITE(*,'(/1X,T9,A,T17,A,T31,A/)') 'X','Integrated','BESSJ3'
	DO 11 I=1,(NSTEP/10)
		J=10*I
		WRITE(*,'(1X,F10.4,2X,2F12.6)') X(J),Y(4,J),BESSJ(3,X(J))
11	CONTINUE
	END

	SUBROUTINE DERIVS(X,Y,DYDX)
	DIMENSION Y(*),DYDX(*)
	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)
	RETURN
	END
	PROGRAM D15R3
C	Driver for routine RKQC
	EXTERNAL DERIVS
	PARAMETER(N=4)
	DIMENSION Y(N),DYDX(N),YSCAL(N)
	X=1.0
	Y(1)=BESSJ0(X)
	Y(2)=BESSJ1(X)
	Y(3)=BESSJ(2,X)
	Y(4)=BESSJ(3,X)
	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)
	DO 11 I=1,N
		YSCAL(I)=1.0
11	CONTINUE
	HTRY=0.1
	WRITE(*,'(/1X,T8,A,T19,A,T31,A,T43,A)')
     *		'eps','htry','hdid','hnext'
	DO 12 I=1,15
		EPS=EXP(-FLOAT(I))
		CALL RKQC(Y,DYDX,N,X,HTRY,EPS,YSCAL,HDID,HNEXT,DERIVS)
		WRITE(*,'(2X,E12.4,F8.2,2X,2F12.6)') EPS,HTRY,HDID,HNEXT
12	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)
	RETURN
	END
	PROGRAM D14R11
C	Driver for routine ROFUNC
	PARAMETER(NMAX=1000,SPREAD=0.05)
	COMMON /ARRAYS/ NPT,X(NMAX),Y(NMAX),ARR(NMAX),AA,ABDEV
	IDUM=-11
	NPT=100
	DO 11 I=1,NPT
		X(I)=0.1*I
		Y(I)=-2.0*X(I)+1.0+SPREAD*GASDEV(IDUM)
11	CONTINUE
	WRITE(*,'(/1X,T10,A,T20,A,T26,A,T37,A/)') 'B','A','ROFUNC','ABDEV'
	DO 12 I=-5,5
		B=-2.0+0.02*I
		RF=ROFUNC(B)
		WRITE(*,'(1X,4F10.2)') B,AA,ROFUNC(B),ABDEV
12	CONTINUE	
	END
	PROGRAM D9R4
C	Driver for routine RTBIS
	EXTERNAL BESSJ0
	PARAMETER(N=100,NBMAX=20,X1=1.0,X2=50.0)
	DIMENSION XB1(NBMAX),XB2(NBMAX)
	NB=NBMAX
	CALL ZBRAK(BESSJ0,X1,X2,N,XB1,XB2,NB)
	WRITE(*,'(/1X,A)') 'Roots of BESSJ0:'
	WRITE(*,'(/1X,T19,A,T31,A/)') 'x','F(x)'
	DO 11 I=1,NB
		XACC=(1.0E-6)*(XB1(I)+XB2(I))/2.0
		ROOT=RTBIS(BESSJ0,XB1(I),XB2(I),XACC)
		WRITE(*,'(1X,A,I2,2X,F12.6,E16.4)') 'Root ',I,ROOT,BESSJ0(ROOT)
11	CONTINUE
	END
	PROGRAM D9R5
C	Driver for routine RTFLSP
	EXTERNAL BESSJ0
	PARAMETER(N=100,NBMAX=20,X1=1.0,X2=50.0)
	DIMENSION XB1(NBMAX),XB2(NBMAX)
	NB=NBMAX
	CALL ZBRAK(BESSJ0,X1,X2,N,XB1,XB2,NB)
	WRITE(*,'(/1X,A)') 'Roots of BESSJ0:'
	WRITE(*,'(/1X,T19,A,T31,A/)') 'x','F(x)'
	DO 11 I=1,NB
		XACC=(1.0E-6)*(XB1(I)+XB2(I))/2.0
		ROOT=RTFLSP(BESSJ0,XB1(I),XB2(I),XACC)
		WRITE(*,'(1X,A,I2,2X,F12.6,E16.4)') 'Root ',I,ROOT,BESSJ0(ROOT)
11	CONTINUE
	END
	PROGRAM D9R8
C	Driver for routine RTNEWT
	EXTERNAL FUNCD,BESSJ0
	PARAMETER(N=100,NBMAX=20,X1=1.0,X2=50.0)
	DIMENSION XB1(NBMAX),XB2(NBMAX)
	NB=NBMAX
	CALL ZBRAK(BESSJ0,X1,X2,N,XB1,XB2,NB)
	WRITE(*,'(/1X,A)') 'Roots of BESSJ0:'
	WRITE(*,'(/1X,T19,A,T31,A/)') 'x','F(x)'
	DO 11 I=1,NB
		XACC=(1.0E-6)*(XB1(I)+XB2(I))/2.0
		ROOT=RTNEWT(FUNCD,XB1(I),XB2(I),XACC)
		WRITE(*,'(1X,A,I2,2X,F12.6,E16.4)') 'Root ',I,ROOT,BESSJ0(ROOT)
11	CONTINUE
	END

	SUBROUTINE FUNCD(X,FN,DF)
	FN=BESSJ0(X)
	DF=-BESSJ1(X)
	RETURN
	END
	PROGRAM D9R9
C	Driver for routine RTSAFE
	EXTERNAL FUNCD,BESSJ0
	PARAMETER(N=100,NBMAX=20,X1=1.0,X2=50.0)
	DIMENSION XB1(NBMAX),XB2(NBMAX)
	NB=NBMAX
	CALL ZBRAK(BESSJ0,X1,X2,N,XB1,XB2,NB)
	WRITE(*,'(/1X,A)') 'Roots of BESSJ0:'
	WRITE(*,'(/1X,T19,A,T31,A/)') 'x','F(x)'
	DO 11 I=1,NB
		XACC=(1.0E-6)*(XB1(I)+XB2(I))/2.0
		ROOT=RTSAFE(FUNCD,XB1(I),XB2(I),XACC)
		WRITE(*,'(1X,A,I2,2X,F12.6,E16.4)') 'Root ',I,ROOT,BESSJ0(ROOT)
11	CONTINUE
	END

	SUBROUTINE FUNCD(X,FN,DF)
	FN=BESSJ0(X)
	DF=-BESSJ1(X)
	RETURN
	END
	PROGRAM D9R6
C	Driver for routine RTSEC
	EXTERNAL BESSJ0
	PARAMETER(N=100,NBMAX=20,X1=1.0,X2=50.0)
	DIMENSION XB1(NBMAX),XB2(NBMAX)
	NB=NBMAX
	CALL ZBRAK(BESSJ0,X1,X2,N,XB1,XB2,NB)
	WRITE(*,'(/1X,A)') 'Roots of BESSJ0:'
	WRITE(*,'(/1X,T19,A,T31,A/)') 'x','F(x)'
	DO 11 I=1,NB
		XACC=(1.0E-6)*(XB1(I)+XB2(I))/2.0
		ROOT=RTSEC(BESSJ0,XB1(I),XB2(I),XACC)
		WRITE(*,'(1X,A,I2,2X,F12.6,E16.4)') 'Root ',I,ROOT,BESSJ0(ROOT)
11	CONTINUE
	END
	PROGRAM D15R8
C	Driver for routine RZEXTR
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**3
		DO 11 J=1,NV
			DUM=DUM/(XEST+1.0)
			YEST(J)=DUM
11		CONTINUE
		CALL RZEXTR(IEST,XEST,YEST,YZ,DY,NV,NUSE)
		WRITE(*,'(/1X,A,I2,A,F8.4)') 'IEST = ',I,'   XEST =',XEST
		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

