c     Some sample Fortran code for the binapprox algorithm.

      real function bapprox(n,x)
      integer n
      real x(n)
      integer i,bottom,counts(0:1000),bin,k,j,count,medbin
      real mu,sigma,scalefac,leftend,rghtend

c     Compute the mean and standard deviation
      mu = sum(x)/n
      sigma = sqrt(dot_product(x-mu,x-mu)/n)

c     Bin x across the interval [mu-sigma, mu+sigma]
      bottom = 0
      counts = 0
      scalefac = 1000/(2*sigma)
      leftend = mu-sigma
      rghtend = mu+sigma

      do 1 i = 1,n
         if (x(i).lt.leftend) then
            bottom = bottom+1
         else if (x(i).lt.rghtend) then
            bin = int((x(i)-leftend) * scalefac)
            counts(bin) = counts(bin)+1
         endif
 1    continue

c     If n is odd
      if (mod(n,2).eq.1) then
c        Find the bin that contains the median
         k = (n+1)/2
         count = bottom 
         do 2 i = 0,1000
            count = count+counts(i)
            if (count.ge.k) then
               bapprox = (i+0.5)/scalefac + leftend
               return
            endif
 2       continue
  
c     If n is even 
      else
c        Find the bins that contain the medians
         k = n/2
         count = bottom
         do 3 i = 0,1000
            count = count+counts(i)
            if (count.ge.k) then
               j = i
 4             continue
               if (count.eq.k) then
                  j = j+1
                  count = count+counts(j)
                  goto 4
               endif
               bapprox = (i+j+1)/(2*scalefac) + leftend
               return
            endif
 3       continue

      endif

      end
