% hist_stat2.m - to generate histogram statistics of injection pump data
% this one does each cast separatly
% sigma, theta, salinity, chlor, pressure

mat_path = './';
plotit = 0;

%sigs = [26.245:.001:26.28];
sigs = [26.205:.001:26.28];
%thetas = [19.12:.002:19.32];
thetas = [19.12:.002:19.5];
sals = [36.68:.001:36.72];
presses = [70:2.5:120];
fls = [-1:.5:10];


for cast = [1:13],
  eval(['!\rm hist_stat',int2str(cast),'.txt'])
  eval(['diary hist_stat',int2str(cast),'.txt'])
  disp(['     Histogram stats - Sample Cast'])
  disp('      ------------------------------')
  diary off
  nsigs = zeros(size(sigs));
  ntheta = zeros(size(thetas));
  nsal = zeros(size(sals));
  npress = zeros(size(presses));
  nfl = zeros(size(fls));
  spd_sum = zeros(size(sigs));

  %diary on
  eval(['diary hist_stat',int2str(cast),'.txt'])
  disp(['     Cast ',int2str(cast)])
  diary off
  %cast = input('Enter cast number: ');
  if cast < 10
    fn = ['SAM_00',int2str(cast)];
  elseif cast < 100
    fn = ['SAM_0',int2str(cast)];
  end

  eval(['load ',mat_path,fn,'_fl'])

  ind = find(sampl==1);
  
  for isig = 1:length(sigs)-1
    % determine which indicies have pump 1 on at this sigma0
    aa = find(sigma0(ind) >= sigs(isig) & sigma0(ind) < sigs(isig+1) );
    % weight by pump speed and add them to the count
    nsigs(isig) = nsigs(isig) + length(aa);
  end  
  
  for itheta = 1:length(thetas)-1
    % determine which indicies have pump 1 on at this theta
    aa = find(theta(ind) >= thetas(itheta) ...
	      & theta(ind) < thetas(itheta+1) );    
    % add them to the count
    ntheta(itheta) = ntheta(itheta) + length(aa);
  end  

  for isal = 1:length(sals)-1
    % determine which indicies have pump 1 on at this sal
    aa = find(sal0(ind) >= sals(isal) ...
	      & sal0(ind) < sals(isal+1) );    
    % add them to the count
    nsal(isal) = nsal(isal) + length(aa);
  end  

  for ip = 1:length(presses)-1
    % determine which indicies have pump 1 on at this press
    aa = find(press(ind) >= presses(ip) ...
	      & press(ind) < presses(ip+1) );    
    % add them to the count
    npress(ip) = npress(ip) + length(aa);
  end  

  for ifl = 1:length(fls)-1    
    % determine which indicies have pump 1 on at this fl
    aa = find(fl(ind) >= fls(ifl) ...
	      & fl(ind) < fls(ifl+1) );    
    % add them to the count
    nfl(ifl) = nfl(ifl) + length(aa);
  end  


if plotit
  if 0
  figure(1)
  bar(sigs,log10(nsigs))
  xlabel('Sigma Theta(kg/m^3)')
  ylabel('log_(10)(Number of Scans (each pump counted separatly))')
  title(['EDDIES oc415-4 Sample Cast - ',fn])
  ax = axis;
  axis([26.24 26.285 0 ax(4)])
  %print
  end
  
  figure(2)
  bar(thetas,ntheta)
  xlabel(['Theta (^o C)'])
  ylabel('Number Scans (each pump counted separatly)')
  title(['EDDIES oc415-4 Sample Cast - ',fn],'interpreter','none')
  ax = axis;
  axis([19.1 19.7 0 ax(4)])
  orient portrait
  eval(['print -dpng ',fn,'_hist_theta.png'])
  %print
    
  figure(3)
  bar(sals,nsal)
  xlabel(['Salinity (PSU)'])
  ylabel('Number Scans (each pump counted separatly)')
  title(['EDDIES oc415-4 Sample Cast - ',fn],'interpreter','none')
  ax=axis;
  axis([36.67 36.72 ax(3:4)])
  %print
  orient portrait
  eval(['print -dpng ',fn,'_hist_sal.png'])
    
  figure(4)
  bar(presses,npress)
  xlabel(['Pressure (dbar)'])
  ylabel('Number Scans (each pump counted separatly)')
  title(['EDDIES oc415-4 Sample Cast - ',fn],'interpreter','none')
  ax=axis;
  %axis([36.699 36.71 ax(3:4)])
  %print
  orient portrait
  eval(['print -dpng ',fn,'_hist_press.png'])

  if 0
  figure(5)
  bar(fls,log10(nfl))
  xlabel(['Chlorophyll (/mug/l)'])
  ylabel('log_(10)(Number of Scans (each pump counted separatly))')
  title(['EDDIES oc415-4 Sample Cast - ',fn],'interpreter','none')
  ax=axis;
  %axis([36.699 36.71 ax(3:4)])
  %print
  end
  
  figure(6)
  bar(fls,(nfl))
  xlabel(['Chlorophyll (/mug/l)'])
  ylabel('Number of Scans (each pump counted separatly)')
  title(['EDDIES oc415-4 Sample Cast - ',fn],'interpreter','none')
  ax=axis;
  %axis([36.699 36.71 ax(3:4)])
  %print
  orient portrait
  eval(['print -dpng ',fn,'_hist_chl.png'])

  figure(7)
  bar(sigs,nsigs)
  xlabel('Sigma Theta(kg/m^3)')
  ylabel('Number of Scans (each pump counted separatly)')
  title(['EDDIES oc415-4 Sample Cast - ',fn],'interpreter','none')
  ax = axis;
  %axis([26.25 26.265 0 ax(4)])
  %print
  orient portrait
  eval(['print -dpng ',fn,'_hist_sig.png'])

end

mean_sig = sum(nsigs.*sigs)/sum(nsigs);
rms_sig = sqrt(sum( ((sigs-mean_sig).^2).*nsigs)/sum(nsigs));
if 0
  mean_sigwt = sum(spd_sum.*sigs)/sum(spd_sum);
  rms_sigwt = sqrt(sum( ((sigs-mean_sigwt).^2).*spd_sum)/sum(spd_sum));
end
mean_th = sum(ntheta.*thetas)/sum(ntheta);
rms_th = sqrt(sum( ((thetas-mean_th).^2).*ntheta)/sum(ntheta));

mean_sal = sum(nsal.*sals)/sum(nsal);
rms_sal = sqrt(sum( ((sals-mean_sal).^2).*nsal)/sum(nsal));

mean_pr = sum(npress.*presses)/sum(npress);
rms_pr = sqrt(sum( ((presses-mean_pr).^2).*npress)/sum(npress));

mean_fl = sum(nfl.*fls)/sum(nfl);
rms_fl = sqrt(sum( ((fls-mean_fl).^2).*nfl)/sum(nfl));


eval(['diary hist_stat',int2str(cast),'.txt'])
%diary hist_stat.txt
disp(['       Mean sigma:       ',num2str(mean_sig), ...
      '       RMS: ',num2str(rms_sig)])


if 0
  disp(['Mean weighted sigma:   ',num2str(mean_sigwt), ...
	' RMS: ',num2str(rms_sigwt)])
end
disp(['       Mean theta:       ',num2str(mean_th), ...
      '       RMS: ',num2str(rms_th)])

disp(['       Mean salinity:    ',num2str(mean_sal), ...
      '       RMS: ',num2str(rms_sal)])

disp(['       Mean Pressure:    ',num2str(mean_pr), ...
      '       RMS: ',num2str(rms_pr)])

disp(['       Mean Chlorophyll: ',num2str(mean_fl), ...
      '       RMS: ',num2str(rms_fl)])

diary off
end

return

% to print all
for ii = [2:4 6 7]
  figure(ii)
  orient landscape
  print
  orient portrait
end


