% histfit.m: fits histograms from *.his files (see histg.c)
%            requires scale (*.scl) file used to obtain histogram range (see
%            scanpdg.c)

Flen = 16;   % 16 Hz (frequency range for image)
Ts = 1/64;   % 64 Hz sampling
P_CalVol = 25e-6;    % 50uV peak-to-peak

% locate histogram file
hisfil = input (' Histogram file ? ', 's');
if ~isempty(hisfil)
  [fb_nbins, fb_H, fb_x0, fb_length] = readhis (hisfil);
  Nhist = max(size(fb_nbins));
% fb_nbins = fb_nbins(1:Nhist);
% fb_H = fb_H(1:fb_nbins(1),1:Nhist);
% fb_x0 = fb_x0(1:fb_nbins(1),1:Nhist);
% fb_length=fb_length(1:fb_nbins(1),1:Nhist);

  % bug fix 03/22/94 BRK
  ind = find (hisfil=='.');
  inds = find (hisfil=='/');
  if ~isempty(ind)
    % get last '/'
    if isempty(inds), k = 0; else, k = max(inds); end;
    if k == 0
      ind_f = min(ind);     % the first dot
      ind_l = max(ind);     % the last dot
    else
      kk = find(ind>k);
      ind_f = min(ind(kk));    % the first dot after last '/'
      ind_l = max(ind(kk));    % the last dot after last '/'
    end;
    varnam = hisfil(k+1:ind_f-1);
    filnam = hisfil(1:ind_l-1);
  else
    if isempty(inds), k = 0; else, k = max(inds); end;
    varnam = hisfil(k+1:max(size(filnam)));
  end;

  % determine file name
  %ii = find(hisfil=='.');
  %if ~isempty(ii)
  %  %ii = ii(max(size(ii)));   % the last dot
  %  ii = ii(1);               % the first dot
  %  filnam = hisfil(1:ii-1);
  %else
  %  filnam = hisfil;
  %end;

  % load scale data (A matrix)
  comm = ['load ', filnam, '.scl']; eval(comm);
  comm = ['A = ', varnam, ';']; eval(comm);
  N = size(A,1)-1;
  dF = Flen/N; % N is number of pixels per image column

  % determine calibration constant
  Ac = input (' Calibration constant (calib2) ? ');
  if Ac <= 0.0, Ac = P_CalVol*sqrt(2e9*Ts/dF); end;
  alpha = P_CalVol / Ac;
  p_scale = 2 * alpha^2 * Ts / dF * 1e9;  % (V^2)s/Hz * 1e9

  % fetch statistics
  smin = A(1,1);
  smean = A(1,2);
  sstd = A(1,3);
  smax = A(1,4);

  bmin = A(2:N+1,1);
  bmean = A(2:N+1,2);
  bstd = A(2:N+1,3);
  bmax = A(2:N+1,4);

  % check number of histograms
  if Nhist ~= N, error('invalid number of frequency bins or histograms'); end;
  disp(sprintf(' Number of histograms: %d', Nhist));
end;

nnp = input (' Percentage of points to plot ? ');
while 1
  fb = input (sprintf (' Frequency bin ? [0:%d] ', N-1));
  if fb < 0 | fb >= N, break; end;
  fb = fb+1;

  y = fb_H(:,fb);
  x = fb_x0(:,fb);
  V = fb_length(:,fb);
  count = sum(y);
  pdf = y/count./V; cdf = cumsum(y/count);

  xfit = x + V/2; L = length(xfit);
  [a1, b1] = expfit (p_scale*xfit(2:L), pdf(2:L));
  pdf_hat = a1 * exp (b1 * p_scale * xfit);

  stx = x + V;
  stx = [0; stx];
  stpdf = [pdf; pdf(length(pdf))];
  stcdf = [cdf; cdf(length(cdf))];

  nn = round(nnp*length(x));
  nnst = round(nnp*length(stx));
  txt1 = sprintf (' Histogram for F = %gHz', (fb-1)*dF);
  txt2 = sprintf (' cdf for F = %gHz', (fb-1)*dF);
  clg;
  subplot(211);stairs(p_scale*stx(1:nnst),stpdf(1:nnst)); title (txt1);
     hold on;
       plot(p_scale*xfit(1:nn), pdf_hat(1:nn));
     hold off;
  subplot(212);stairs(p_scale*stx(1:nn),stcdf(1:nn)); title (txt2);
end;
