program stats

  !############################################################
  !
  !
  !     Calculates various stats on a simple ascii table of upto
  !     1000 values.
  !     
  !     Expects an ascii file calles 'test.dat' which contains
  !     on numerical value per row.
  !
  !     Output returned to screen.
  !
  !     Bugs: None known
  !
  !     Original version in f77 by Simon Driver (SPD) 23/04/06
  !
  !       DATE            MODIFICATION              PROGRAMMER
  !      25/03/06         Converted to f90              SPD
  !      10/04/06         Fixed f90 language            TR
  !
  !############################################################

  !     DECLARATIONS

  implicit none

  !     DECLARE VARIABLES

  integer :: i,n,fn,ofn,count
  integer :: ioerr
  real,dimension(1000) :: pixin
  real :: max,tmp,min,sum,sigma,nsum,mean,med

  !     DECLARE FUNCTIONS

  real :: getmax,getmin,getsum,getmed

  !############################################################

  !     INITIALISE VALUES

  ofn=0

  !     OPEN DATA FILE CALLED 'test.dat'

  open(unit=1,file='test.dat',status='old')

  !     READ IN DATA

  do i=1,1000
     read(1,*,iostat=ioerr) pixin(i)
  enddo

  !     CLOSE DATA FILE

  close(unit=1)

  !     DEFINE NUMBER OF ELEMENTS IN ARRAY

  n=i-1

  !     GET BASIC STATS

  max=getmax(pixin,n)
  min=getmin(pixin,n)
  sum=getsum(pixin,n)
  mean=sum/n
  call std(pixin,n,mean,0.,sigma,fn) ! 0. means no clipping

  !     WRITE OUT BASIC STATS

  write(*,*),' '
  write(*,*),' OPERATING ON: test.dat'
  write(*,100),' First   ',pixin(1)
  write(*,100),' Last    ',pixin(n)
  write(*,100),' Sum     ',sum
  write(*,100),' Highest ',max
  write(*,100),' Lowest  ',min
  write(*,100),' Mean    ',mean
  write(*,100),' Std     ',sigma

  !     CALCULATE MEDIAN

  med=getmed(pixin,n)

  write(*,100),' Median  ',med

  !     RECALCULATE MEAN AND SIGMA USING 3sigma CLIPPING
  !     DON"T CLIP MORE THAN 10 TIMES

  do i=1,10
     call std(pixin,n,mean,3.*sigma,sigma,fn)
     if(fn.eq.ofn) exit
     ofn=fn
  enddo

  !     WRITE CLIPPED OUTPUT

  write(*,110),' Mean(3sig)',mean,'(',fn,'/',n,')'
  write(*,110),' Std(3sig)',sigma
  write(*,*),' '

  !     FORMAT STATEMENTS

100 format(a12,f12.6)
110 format(a12,f12.6,3x,a1,i4.4,a1,i4.4,a1)

  !     EXIT PROGRAM

end program stats

!############################################################

subroutine std(x,elements,mean,limit,sigma,n)

  implicit none

  integer :: i,elements,n

  real,dimension(elements) :: x
  real :: mean,limit,sigma,sum

  n=0
  sigma=0.0
  sum=0.0

  do i=1,elements
     if(limit.eq.0.0)then
        sigma=sigma+(mean-x(i))**2
        sum=sum+x(i)
        n=n+1
     else
        if(abs(x(i)-mean).le.limit)then
           sigma=sigma+(mean-x(i))**2
           sum=sum+x(i)
           n=n+1
        endif
     endif
  enddo

  mean=sum/n
  sigma=(sigma/(n-1))**0.5

end subroutine std

!############################################################

function getmax(x,n)

  implicit none

  integer :: i,n
  real,dimension(n) :: x
  real :: max,getmax

  max=x(1)

  do i=1,n
     if(x(i).gt.max)max=x(i)
  enddo

  getmax=max

end function getmax

!############################################################

function getmin(x,n)

  implicit none

  integer :: i,n
  real,dimension(n) :: x
  real :: min,getmin

  min=x(1)

  do i=1,n
     if(x(i).lt.min)min=x(i)
  enddo

  getmin=min

end function getmin

!############################################################

function getsum(x,n)

  implicit none

  integer :: i,n
  real,dimension(n) :: x
  real :: sum,getsum

  sum=0.

  do i=1,n
     sum=sum+x(i)
  enddo

  getsum=sum

end function getsum

!############################################################

function getmed(x,n)

  implicit none

  integer :: n,i,cnt

  real,dimension(n) :: x
  real :: tmp,getmed

  do while(cnt==0)

     do i=1,n-1
        if(x(i).gt.x(i+1))then
           tmp=x(i)
           x(i)=x(i+1)
           x(i+1)=tmp
           cnt=1
        endif
     end do

  end do

  getmed=x(int(0.5*n))

end function getmed

!############################################################