c     Some sample Fortran code for the binmedian algorithm.
c     Note: I haven't written the code for n even yet.

      real function bmedian(n,x)
      integer n
      real x(n)
      integer i,bottom,counts(0:1000),bin,k,count,medbin,r,oldbin,j
      real mu,sigma,scalefac,leftend,rghtend,a,oldscale,oldleft
      logical samepts;

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
         k = (n+1)/2
         r = 1
    
c        Beginning of recursive step:
c        Find the bin that contains the median, and the order
c        of the median within that bin
 2       count = bottom 
         do 3 i = 0,1000
            count = count+counts(i)
            if (count.ge.k) then
               medbin = i 
               k = k - (count-counts(i))
               goto 4
            endif
 3       continue
        
 4       bottom = 0
         counts = 0
         oldscale = scalefac
         oldleft = leftend
         scalefac = 1000*oldscale
         leftend = medbin/oldscale + oldleft
         rghtend = (medbin+1)/oldscale + oldleft

c        Determine which points map to medbin, and put
c        them in spots r,...n
         i = r
         r = n+1  
 5       if (i.lt.r) then 
            oldbin = int((x(i)-oldleft) * oldscale)
            
            if (oldbin.eq.medbin) then
               r = r-1
c              Swap x(i) and x(r)
               a = x(i)
               x(i) = x(r)
               x(r) = a

c              Re-bin on a finer scale
               if (x(r).lt.leftend) then
                  bottom = bottom+1
               else if (x(r).lt.rghtend) then
                  bin = int((x(r)-leftend) * scalefac)
                  counts(bin) = counts(bin)+1
               endif
            else 
               i = i+1
            endif
            goto 5
         endif
         
c        Stop if all points in medbin are the same
         samepts = .TRUE.
         do 6 i = r+1,n
            if (x(i).ne.x(r)) then
               samepts = .FALSE.
               goto 7
            endif
 6       continue
 7       if (samepts) then
            bmedian = x(r)
            return
         endif
         
c        Stop if there's <= 20 points left 
         if ((n+1-r).le.20) then
            goto 8
         endif
         goto 2

c        Perform insertion sort on the remaining points,
c        and then pick the kth smallest
 8       do 9 i = r+1,n
            a = x(i)
            do 10 j = i-1,r,-1
               if (x(j).lt.a) goto 11
               x(j+1) = x(j)
 10         continue
            j = r-1
 11         x(j+1) = a
 9       continue

         bmedian = x(r-1+k)
         return
         
c     If n is even 
      else 
         bmedian = 0
         return
      endif

      end
