function [H, X, GaussFit, means, stds ] = djcmhisto(Y,bins)
  ## [H, X, GaussFit, means, stds ] = djcmhisto(Y,bins)
  ##
  ##  Makes a histogram of each _column_ of the matrix Y
  ##  and returns a Gaussian fit to each histogram.
  ##  The integer "bins" defines the number of bins that should be 
  ##  used for each histogram. 
  ##
  ##  Each column is processed separately. H(b,c) is the NUMBER of 
  ##  occurrences in bin b. X(b,c) is the midpoint of bin b (for column c).
  ##  The bins are determined by finding the range (min,max) of each 
  ##  column, and chopping it up into (bins) equal pieces.
  ##
  ##  djcmhisto also computes the mean and standard deviation of each
  ##  column (the vectors "means" and "stds"). It also creates a Gaussian
  ##  function, GaussFit, scaled such that it is comparable to H. 
  ##
  ## Example: make random numbers, find histogram, and display Gaussian fit.
  ## 
  ##   a = rand(2000,3) + rand(2000,3) + rand(2000,3);
  ##   [H,X,G] = djcmhisto(a,50);
  ##   f=1;   
  ##   hold off; plot( X(:,f), H(:,f) ); hold on; plot( X(:,f), G(:,f));
  ##   
  ##  or alternatively using the gplot function (which is faster than plot)
  ##
  ##   XH =  [ X(:,f) , H(:,f) , G(:,f) ] ;
  ##   gset yrange [0:]
  ##   gplot XH u 1:2 t 'Data' w histeps,  XH u 1:3 t 'Gaussian fit' w l 4 ;
  verbose = 1 ; 
  [rows,cols] = size(Y) ; 
  maxs = max( Y );
  mins = min( Y ) ;
  means = mean( Y ) ;
  stds  = std( Y ) ; 
  vars  = stds.**2 ;
  range = maxs - mins ; 

  for c = 1 : cols
    if( range(1,c) <= 0.0 ) 
      disp("warning no range in column"),disp(c) ;
      range(1,c) = 1.0 ;
    endif
  endfor
  if(verbose)
    range,   bins
  endif
  dx = range / bins ; 

  if(verbose>1)
    maxs
    mins
    range
    means
    vars
    dx
    if(verbose>3)
      ans=input("ready");
    endif
  endif

  H = zeros( bins , cols ) ;
  X = zeros( bins , cols ) ;
  GaussFit = zeros( bins , cols ) ;

  for c = 1:cols
##    mins(1,c),    dx(1,c)
    X(:,c) = ones(bins,1) * mins(1,c) + ( 0.5 + [0:bins-1]' )  * dx(1,c) ; 
    ## X(b,c) is the midpoints of the bins
  endfor
  
  for c = 1:cols
    ## subtract the 0.5 in order to get rounding DOWN to the floor integer
    Y(:,c) = round( ( Y(:,c) - mins(1,c)  ) / dx(1,c) - 0.5 ) + 1 ;
##    Y(:,c) = myround( ( Y(:,c) - mins(1,c)  ) / dx(1,c) ) + 1 ;
    for r =1:rows
##      [c,r,Y(r,c)]
      if (( Y(r,c) <= bins ) &&  ( Y(r,c) >= 1 ) ) 
	H(Y(r,c),c) ++ ; 
      else
	if (( Y(r,c) >= bins+2 ) ||  ( Y(r,c) <= -1 ) ) 
	  disp ("unexpected overflow") ;
	  [r,c,mins(1,c),dx(1,c),Y(r,c), bins ]
	endif
	if (( Y(r,c) > bins ))
	  H(bins,c) ++ ; 
	else
	  H(1,c) ++ ;
	endif
      endif
    endfor
  endfor

  for c = 1:cols
    for b = 1:bins
      x = X(b,c) ; ##  This is the midpoint of the bin
      GaussFit(b,c) = rows * exp( -( x - means(1,c) )**2/(2*vars(1,c))) * dx(1,c) / sqrt( 2.0 * pi * vars(1,c) ) ; 
    endfor
  endfor
endfunction
    

  ##
  ##
  ## (c) David MacKay Sept 2005. Free software, Distributed under the GPL.
  ##
  ##
## This library is free software; you can redistribute it and/or
## modify it under the terms of the GNU Library General Public
## License as published by the Free Software Foundation; either
## Version 2 of the License, or (at your option) any later version.
