subroutine statistics ( NPoints, NBlocks, X, XBavg, Xavg, Std, Err ) 

implicit none

integer,          intent(in)	:: NPoints       ! Number of points to be averaged
integer,          intent(in)    :: NBlocks       ! Number of blocks to be used

real, dimension(NPoints), intent(in)	:: X     ! Array of points to average

real, dimension(NBlocks), intent(out)	:: XBavg ! Block Averages
real, intent(out)						:: Xavg  ! Average
real, intent(out)						:: Std   ! Standard Deviation
real, intent(out)						:: Err   ! Error

real									:: Sign_Level = 0.95
! Level of significance for evaluating student-t for
! Nblock-1 degrees of freedom.

! Local variables.
integer                     ::  iblock, ipoint
integer                     ::  FirstPoint, Pointsperblock
real                        ::  tt
real, external              ::  student_t   ! Function


Pointsperblock = NPoints / NBlocks	   ! Integer division

if (Pointsperblock.eq.0) then
	write (*,*) 'The number of Points per block is too large for the specified number of blocks'
	write (*,*) 'Please specify a smaller number of Points per block, or increase the number of Points'
	stop
endif

FirstPoint = MOD( NPoints, NBlocks )
tt = student_t ( NBlocks, Sign_Level )
	
! Initializing
XBavg = 0.0
Xavg = 0.0
Std = 0.0
Err = 0.0

! Calculate the block averages.
do iblock = 1, NBlocks
	do ipoint = FirstPoint + (iblock - 1) * Pointsperblock + 1, &
	       		FirstPoint + iblock * Pointsperblock
		XBavg(iblock) = XBavg(iblock) + X(ipoint)
	end do
	XBavg(iblock) = XBavg(iblock) / real(Pointsperblock)
end do						  

! Averege over the simulation.	
Xavg = SUM( XBavg( 1:NBlocks ) ) / real( NBlocks )

Do iblock = 1, NBlocks
	Std = Std + ( XBavg(iblock)	- Xavg ) ** 2
end do

Std = sqrt( Std / ( NBlocks -1 ) )

Err = Std * tt / sqrt( real( NBlocks ) )

return

end	subroutine statistics