	PROGRAM D8R8
C	Driver for routine QCKSRT
	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 QCKSRT(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 D4R8
C	Driver for routine QGAUS
	EXTERNAL FUNC
	PARAMETER(X1=0.0,X2=5.0,NVAL=10)
	DX=(X2-X1)/NVAL
	WRITE(*,'(/1X,A,T12,A,T23,A/)') '0.0 to','QGAUS','Expected'
	DO 11 I=1,NVAL
		X=X1+I*DX
		CALL QGAUS(FUNC,X1,X,SS)
		WRITE(*,'(1X,F5.2,2F12.6)') X,SS,
     *			-(1.0+X)*EXP(-X)+(1.0+X1)*EXP(-X1)
11	CONTINUE
	END

	FUNCTION FUNC(X)
	FUNC=X*EXP(-X)
	END
	PROGRAM D4R4
C	Driver for routine QROMB
	PARAMETER(PIO2=1.5707963)
	EXTERNAL FUNC,FINT
	A=0.0
	B=PIO2
	WRITE(*,'(1X,A)') 'Integral of FUNC computed with QROMB'
	WRITE(*,'(1X,A,F10.6)') 'Actual value of integral is',
     *		FINT(B)-FINT(A)
	CALL QROMB(FUNC,A,B,S)
	WRITE(*,'(1X,A,F10.6)') 'Result from routine QROMB is',S
	END

	FUNCTION FUNC(X)
	FUNC=(X**2)*(X**2-2.0)*SIN(X)
	END

	FUNCTION FINT(X)
C	Integral of FUNC
	FINT=4.0*X*((X**2)-7.0)*SIN(X)
     *		-((X**4)-14.0*(X**2)+28.0)*COS(X)
	END
	PROGRAM D4R6
C	Driver for QROMO
	PARAMETER(X1=0.0,X2=1.5707963,X3=3.1415926,AINF=1.0E20)
	EXTERNAL FUNCL,MIDSQL,FUNCU,MIDSQU,FNCINF,MIDINF,FNCEND
	WRITE(*,'(/1X,A)') 'Improper integrals:'
	CALL QROMO(FUNCL,X1,X2,RESULT,MIDSQL)
	WRITE(*,'(/1X,A)')
     *		'Function: SQRT(x)/SIN(x)       Interval: (0,pi/2)'
	WRITE(*,'(1X,A,F8.4)')
     *		'Using: MIDSQL                  Result:',RESULT
	CALL QROMO(FUNCU,X2,X3,RESULT,MIDSQU)
	WRITE(*,'(/1X,A)')
     *		'Function: SQRT(pi-x)/SIN(x)    Interval: (pi/2,pi)'
	WRITE(*,'(1X,A,F8.4)')
     *		'Using: MIDSQU                  Result:',RESULT
	CALL QROMO(FNCINF,X2,AINF,RESULT,MIDINF)
	WRITE(*,'(/1X,A)')
     *		'Function: SIN(x)/x**2          Interval: (pi/2,infty)'
	WRITE(*,'(1X,A,F8.4)')
     *		'Using: MIDINF                  Result:',RESULT
	CALL QROMO(FNCINF,-AINF,-X2,RESULT,MIDINF)
	WRITE(*,'(/1X,A)')
     *		'Function: SIN(x)/x**2          Interval: (-infty,-pi/2)'
	WRITE(*,'(1X,A,F8.4)')
     *		'Using: MIDINF                  Result:',RESULT
	CALL QROMO(FNCEND,X1,X2,RES1,MIDSQL)
	CALL QROMO(FNCEND,X2,AINF,RES2,MIDINF)
	WRITE(*,'(/1X,A)')
     *		'Function: EXP(-x)/SQRT(x)      Interval: (0.0,infty)'
	WRITE(*,'(1X,A,F8.4/)')
     *		'Using: MIDSQL,MIDINF           Result:',RES1+RES2
	END

	SUBROUTINE MIDSQL(FUNK,AA,BB,S,N)
	FUNC(X)=2.*X*FUNK(AA+X**2)
	B=SQRT(BB-AA)
	A=0.0
	IF (N.EQ.1) THEN
		S=(B-A)*FUNC(0.5*(A+B))
		IT=1
	ELSE
		TNM=IT
		DEL=(B-A)/(3.*TNM)
		DDEL=DEL+DEL
		X=A+0.5*DEL
		SUM=0.
		DO 11 J=1,IT
			SUM=SUM+FUNC(X)
			X=X+DDEL
			SUM=SUM+FUNC(X)
			X=X+DEL
11		CONTINUE
		S=(S+(B-A)*SUM/TNM)/3.
		IT=3*IT
	ENDIF
	RETURN
	END

	SUBROUTINE MIDSQU(FUNK,AA,BB,S,N)
	FUNC(X)=2.0*X*FUNK(BB-X**2)
	B=SQRT(BB-AA)
	A=0.0
	IF (N.EQ.1) THEN
		S=(B-A)*FUNC(0.5*(A+B))
		IT=1
	ELSE
		TNM=IT
		DEL=(B-A)/(3.*TNM)
		DDEL=DEL+DEL
		X=A+0.5*DEL
		SUM=0.
		DO 11 J=1,IT
			SUM=SUM+FUNC(X)
			X=X+DDEL
			SUM=SUM+FUNC(X)
			X=X+DEL
11		CONTINUE
		S=(S+(B-A)*SUM/TNM)/3.
		IT=3*IT
	ENDIF
	RETURN
	END

	FUNCTION FUNCL(X)
	FUNCL=SQRT(X)/SIN(X)
	END

	FUNCTION FUNCU(X)
	PARAMETER(PI=3.1415926)
	FUNCU=SQRT(PI-X)/SIN(X)
	END

	FUNCTION FNCINF(X)
	FNCINF=SIN(X)/(X**2)
	END

	FUNCTION FNCEND(X)
	FNCEND=EXP(-X)/SQRT(X)
	END
	PROGRAM D9R12
C	Driver for routine QROOT
	PARAMETER(N=7,EPS=1.0E-6,NTRY=10,TINY=1.0E-5)
	DIMENSION P(N),B(NTRY),C(NTRY)
	DATA P/10.0,-18.0,25.0,-24.0,16.0,-6.0,1.0/
	WRITE(*,'(/1X,A)') 'P(x)=x^6-6x^5+16x^4-24x^3+25x^2-18x+10'
	WRITE(*,'(1X,A)') 'Quadratic factors x^2+Bx+C'
	WRITE(*,'(/1X,A,T15,A,T27,A/)') 'Factor','B','C'
	NROOT=0
	DO 12 I=1,NTRY
		C(I)=0.5*I
		B(I)=-0.5*I
		CALL QROOT(P,N,B(I),C(I),EPS)
		IF (NROOT.EQ.0) THEN
			WRITE(*,'(1X,I3,2X,2F12.6)') NROOT,B(I),C(I)
			NROOT=1
		ELSE
			NFLAG=0
			DO 11 J=1,NROOT		     
				IF (ABS(B(I)-B(J)).LT.TINY
     *				.AND. ABS(C(I)-C(J)).LT.TINY) NFLAG=1
11			CONTINUE
			IF (NFLAG.EQ.0) THEN
				WRITE(*,'(1X,I3,2X,2F12.6)') NROOT,B(I),C(I)
				NROOT=NROOT+1
			ENDIF
		ENDIF
12	CONTINUE
	END
	PROGRAM D4R3
C	Driver for routine QSIMP
	PARAMETER(PIO2=1.5707963)
	EXTERNAL FUNC,FINT
	A=0.0
	B=PIO2
	WRITE(*,'(1X,A)') 'Integral of FUNC computed with QSIMP'
	WRITE(*,'(1X,A,F10.6)') 'Actual value of integral is',
     *		FINT(B)-FINT(A)
	CALL QSIMP(FUNC,A,B,S)
	WRITE(*,'(1X,A,F10.6)') 'Result from routine QSIMP is',S
	END

	FUNCTION FUNC(X)
	FUNC=(X**2)*(X**2-2.0)*SIN(X)
	END

	FUNCTION FINT(X)
C	Integral of FUNC
	FINT=4.0*X*((X**2)-7.0)*SIN(X)
     *		-((X**4)-14.0*(X**2)+28.0)*COS(X)
	END
	PROGRAM D4R2
C	Driver for routine QTRAP
	PARAMETER(PIO2=1.5707963)
	EXTERNAL FUNC,FINT
	A=0.0
	B=PIO2
	WRITE(*,'(1X,A)') 'Integral of FUNC computed with QTRAP'
	WRITE(*,'(1X,A,F10.6)') 'Actual value of integral is',
     *		FINT(B)-FINT(A)
	CALL QTRAP(FUNC,A,B,S)
	WRITE(*,'(1X,A,F10.6)') 'Result from routine QTRAP is',S
	END

	FUNCTION FUNC(X)
	FUNC=(X**2)*(X**2-2.0)*SIN(X)
	END

	FUNCTION FINT(X)
C	Integral of FUNC
	FINT=4.0*X*((X**2)-7.0)*SIN(X)
     *		-((X**4)-14.0*(X**2)+28.0)*COS(X)
	END
	PROGRAM D4R10
C	Driver for routine QUAD3D
	COMMON XMAX
	PARAMETER(PI=3.1415926,NVAL=10)
	WRITE(*,*) 'Integral of r^2 over a spherical volume'
	WRITE(*,'(/4X,A,T14,A,T24,A)') 'Radius','QUAD3D','Actual'
	DO 11 I=1,NVAL
		XMAX=0.1*I
		XMIN=-XMAX
		CALL QUAD3D(XMIN,XMAX,S)
		WRITE(*,'(1X,F8.2,2F10.4)') XMAX,S,4.0*PI*(XMAX**5)/5.0
11	CONTINUE
	END

	FUNCTION FUNC(X,Y,Z)
	FUNC=X**2+Y**2+Z**2
	END

	FUNCTION Z1(X,Y)
	COMMON XMAX
	Z1=-SQRT(ABS(XMAX**2-X**2-Y**2))
	END

	FUNCTION Z2(X,Y)
	COMMON XMAX
	Z2=SQRT(ABS(XMAX**2-X**2-Y**2))
	END

	FUNCTION Y1(X)
	COMMON XMAX
	Y1=-SQRT(ABS(XMAX**2-X**2))
	END

	FUNCTION Y2(X)
	COMMON XMAX
	Y2=SQRT(ABS(XMAX**2-X**2))
	END

	SUBROUTINE QGAUSX(FUNC,A,B,SS)
	DIMENSION X(5),W(5)
	DATA X/.1488743389,.4333953941,.6794095682,
     *		.8650633666,.9739065285/
	DATA W/.2955242247,.2692667193,.2190863625,
     *		.1494513491,.0666713443/
	XM=0.5*(B+A)
	XR=0.5*(B-A)
	SS=0
	DO 11 J=1,5
		DX=XR*X(J)
		SS=SS+W(J)*(FUNC(XM+DX)+FUNC(XM-DX))
11	CONTINUE
	SS=XR*SS
	RETURN
	END

	SUBROUTINE QGAUSY(FUNC,A,B,SS)
	DIMENSION X(5),W(5)
	DATA X/.1488743389,.4333953941,.6794095682,
     *		.8650633666,.9739065285/
	DATA W/.2955242247,.2692667193,.2190863625,
     *		.1494513491,.0666713443/
	XM=0.5*(B+A)
	XR=0.5*(B-A)
	SS=0
	DO 11 J=1,5
		DX=XR*X(J)
		SS=SS+W(J)*(FUNC(XM+DX)+FUNC(XM-DX))
11	CONTINUE
	SS=XR*SS
	RETURN
	END

	SUBROUTINE QGAUSZ(FUNC,A,B,SS)
	DIMENSION X(5),W(5)
	DATA X/.1488743389,.4333953941,.6794095682,
     *		.8650633666,.9739065285/
	DATA W/.2955242247,.2692667193,.2190863625,
     *		.1494513491,.0666713443/
	XM=0.5*(B+A)
	XR=0.5*(B-A)
	SS=0
	DO 11 J=1,5
		DX=XR*X(J)
		SS=SS+W(J)*(FUNC(XM+DX)+FUNC(XM-DX))
11	CONTINUE
	SS=XR*SS
	RETURN
	END
