REAL FUNCTION SIMP( DX, Y, NPTS )
C
C	This function will integrate a set of equally spaced points Y.
C	If NPTS is odd, Simpson's 1/3 rule for integration will be
C	employed; if NPTS equals 2, trapezodial rule; if NPTS equals 4,
C	the integral of a cubic through the points; if NPTS is even
C	and greater than 4, Simpson's 1/3 rule is applied to the first
C	NPTS-1 points with a parabolic end correction.  FCTR: SI
C
C	DX	Spacing between points.
C	Y	Array of values.
C	NPTS	Number of points.
C
	DIMENSION Y( NPTS )
C
C - CHECK TO SEE IF THE NUMBER OF POINTS IS ODD.
	IF( MOD( NPTS, 2 ) .EQ. 0 ) GO TO 1100
C
C - PERFORM SIMPSON'S 1/3 RULE INTEGRATION.
	SIMP = Y(1) - Y( NPTS )
	NEND = NPTS - 1
	DO 1000 I=2,NEND,2
1000	SIMP = SIMP + 4.*Y( I ) + 2.*Y( I+1 )
	SIMP = 0.33333333*DX*SIMP
	RETURN
C
C - EVEN NUMBER OF POINTS, GREATER THAN TWO?
1100	IF( NPTS .GT. 2 ) GO TO 1200
C
C - PERFORM TRAPEZODIAL INTEGRATION.
	SIMP = 0.5*DX*( Y(1) + Y(2) )
	RETURN
C
C - SUFFICIENT POINTS TO APPLY SIMPSON'S RULE WITH CORRECTION?
1200	IF( NPTS .GE. 6 ) GO TO 1300
C
C - FIT A CUBIC AND INTEGRATE.
	Y2Y1 = Y( 2 ) - Y( 1 )
	Y3Y1 = Y( 3 ) - Y( 1 )
	Y4Y1 = Y( 4 ) - Y( 1 )
	A = 0.5*Y2Y1 - 0.5*Y3Y1 + 0.16666667*Y4Y1
	B = -2.5*Y2Y1 + 2.*Y3Y1 - 0.5*Y4Y1
	C = 3.*Y2Y1 - 1.5*Y3Y1 +0.33333333*Y4Y1
	D = Y( 1 )
	DX3 = DX*3.
	SIMP = DX3*( D + DX3*( 0.5*C + DX3*
	1  ( 0.33333333*B + 0.25*DX3*A ) ) )
	RETURN
C
C - SIMPSON'S 1/3 RULE INTEGRATION.
1300	SIMP = Y( 1 ) - Y( NPTS-1 )
	NEND = NPTS - 2
	DO 2000 I=2,NEND,2
2000	SIMP = SIMP + 4.*Y( I ) + 2.*Y( I+1 )
	SIMP = 0.33333333*DX*SIMP
C
C - INCLUDE PARABOLIC CORRECTION.
	YM0YM1 = Y( NPTS ) - Y( NPTS-1 )
	YM2YM1 = Y( NPTS-2 ) - Y( NPTS-1 )
	A = 0.5*( YM2YM1 + YM0YM1 )
	B = 0.5*( YM0YM1 - YM2YM1 )
	C = Y( NPTS-1 )
	SIMP = SIMP + DX*( C + DX*( 0.5*B + 0.33333333*DX*A ) )
	RETURN
	END