c      
c     Computes the mean of array x
c    
      function mean(xc, n) bind(c)
      use ISO_C_BINDING
      implicit none
cGCC$ ATTRIBUTES STDCALL, DLLEXPORT :: mean
      
      real(kind=c_float)::mean
      
      type(c_ptr), value::xc
      integer(kind=c_int), value::n
      
      real(kind=c_float), dimension(:), pointer::x
      real::array_sum
          
          call C_F_POINTER(xc, x, [n])
          
          array_sum = sum(x,1)
          mean = array_sum / real(n)
          
      end function mean
      
c
c     A not-particularly-efficient sort algorithm
c
      subroutine sort(x, n)
      implicit none
          
      real, intent(inout), dimension(:)::x
      integer, intent(in)::n
      integer::i, minpos
      real::temp
      
          do i = 1, n-1
              minpos = MINLOC(x(i:n), 1) + i - 1
              if(minpos .NE. i) then
                  temp = x(i)
                  x(i) = x(minpos)
                  x(minpos) = temp
              end if
          end do
          
      end subroutine sort
c
c	  Computes and returns the median of array x      
c
      function median(xc, n) bind(c)
      use ISO_C_BINDING
      implicit none
cGCC$ ATTRIBUTES STDCALL, DLLEXPORT :: median
      
      real(kind=c_float)::median
      type(c_ptr), value::xc
      integer(kind=c_int), value::n
      
      real(kind=c_float), dimension(:), pointer::x
      real, dimension(:), allocatable::sorted
      
      interface
          subroutine sort(x, n)
          real, intent(inout), dimension(:)::x
          integer, intent(in)::n
          end subroutine
      end interface
          
          call C_F_POINTER(xc, x, [n])
          
          allocate(sorted(1:n))
          sorted = x
          call sort(sorted, n)
          median = sorted(n/2)
      
      end function median
      
c
c     Computes the mode of the array x and returns the mode or
c     modes in the "modes" argument.  The function returns the
c     number of modes located.
c
      function mode(xc, n, modesc) bind(c)
      use ISO_C_BINDING
      implicit none
cGCC$ ATTRIBUTES STDCALL, DLLEXPORT :: mode
      
      integer(kind=c_int)::mode
      type(c_ptr), value::xc, modesc
      integer(kind=c_int), value::n
      
      
      real(kind=c_float), dimension(:), pointer::x
      real(kind=c_float), dimension(:), pointer::modes
      
      real, dimension(:), allocatable::sorted
      real, dimension(:), allocatable::set
      integer, dimension(:), allocatable::counts
      
      real::working
      integer::i_sorted, i_set, i_modes
      integer::max_count
      
      interface
          subroutine sort(x, n)
          real, intent(inout), dimension(:)::x
          integer, intent(in)::n
          end subroutine
      end interface
          
          call C_F_POINTER(xc, x, [n])
          call C_F_POINTER(modesc, modes, [n])
          
          allocate(sorted(1:n))
          sorted = x
          call sort(sorted, n)
          
          allocate(set(1:n))
          allocate(counts(1:n))
          
          counts = 0
          set = 0
                    
          i_set = 0
          working = sorted(1)+1.0
          i_sorted = 1
          do i_sorted = 1, n
              if(working .NE. sorted(i_sorted)) then
                  i_set = i_set + 1
                  set(i_set) = sorted(i_sorted) 
                  counts(i_set) = 1
                  working = sorted(i_sorted)
              else
                  counts(i_set) = counts(i_set) + 1
              end if
          end do
          
          max_count = maxval(counts, 1)
          if(max_count .EQ. 1) then
c	          If all counts are 1, there is no mode          
              mode = 0
          else
c             There may be multiple modes...
              i_set = maxloc(counts, 1)
              i_modes = 0
              do while(counts(i_set) .GE. 1)
                  if(counts(i_set) .EQ. max_count) then
                      i_modes = i_modes + 1
                      modes(i_modes) = set(i_set)
                  end if
                  i_set = i_set + 1
              end do
              mode = i_modes
          end if
      
      end function mode
      