% 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 = 1;

sigs = [26.245:.0005:26.265];
thetas = [19.24:.002:19.32];
sals = [36.7:.001:36.71];
presses = [75:2.5:115];
fls = [0:.5:20];

%!\rm hist_stat.txt
diary hist_stat.txt
disp('Histogram stats - total injection')
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));

for cast = [3:4],
  disp(int2str(cast))
  %cast = input('Enter cast number: ');
  if cast < 10
    fn = ['INJ_00',int2str(cast)];
  elseif cast < 100
    fn = ['sf010',int2str(cast)];
  else
    fn = ['sf01',int2str(cast)];
  end

  eval(['load ',mat_path,fn,'_fl'])
  if cast == 3
    eval(['load ',mat_path,fn,'_spd'])
    p1_spd = spd_1;
    p2_spd = spd_2;
  end

  psf6_1 = pump_1 + p1_sf6;
  psf6_2 = pump_2 + p2_sf6;
  ind_p1 = find(psf6_1 == 2);
  ind_p2 = find(psf6_2 == 2);

  for isig = 1:length(sigs)-1
    % determine which indicies have pump 1 on at this sigma0
    aa = find(sigma0(ind_p1) >= sigs(isig) & sigma0(ind_p1) < sigs(isig+1) );
    % weight by pump speed and add them to the count
    nsigs(isig) = nsigs(isig) + length(aa);
    spd_sum(isig) = spd_sum(isig) + sum(p1_spd(aa));
      
    % determine which indicies have pump 2 on at this sigma0
    aa = find(sigma0(ind_p2) >= sigs(isig) & sigma0(ind_p2) < sigs(isig+1) );
    % add them to the count
    nsigs(isig) = nsigs(isig) + length(aa);
    spd_sum(isig) = spd_sum(isig) + sum(p2_spd(aa));
  end  
  
  for itheta = 1:length(thetas)-1
    % determine which indicies have pump 1 on at this theta
    aa = find(theta(ind_p1) >= thetas(itheta) ...
	      & theta(ind_p1) < thetas(itheta+1) );    
    % add them to the count
    ntheta(itheta) = ntheta(itheta) + length(aa);
    % determine which indicies have pump 2 on at this theta
    aa = find(theta(ind_p2) >= thetas(itheta) ...
	      & theta(ind_p2) < 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_p1) >= sals(isal) ...
	      & sal0(ind_p1) < sals(isal+1) );    
    % add them to the count
    nsal(isal) = nsal(isal) + length(aa);
    % determine which indicies have pump 2 on at this sal
    aa = find(sal0(ind_p2) >= sals(isal) ...
	      & sal0(ind_p2) < 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_p1) >= presses(ip) ...
	      & press(ind_p1) < presses(ip+1) );    
    % add them to the count
    npress(ip) = npress(ip) + length(aa);
    % determine which indicies have pump 2 on at this press
    aa = find(press(ind_p2) >= presses(ip) ...
	      & press(ind_p2) < 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_p1) >= fls(ifl) ...
	      & fl(ind_p1) < fls(ifl+1) );    
    % add them to the count
    nfl(ifl) = nfl(ifl) + length(aa);
    % determine which indicies have pump 2 on at this fl
    aa = find(fl(ind_p2) >= fls(ifl) ...
	      & fl(ind_p2) < fls(ifl+1) );
    % add them to the count
    nfl(ifl) = nfl(ifl) + length(aa);
  end  
end

if plotit
  figure(1)
  bar(sigs,log10(nsigs))
  xlabel('Sigma Theta(kg/m^3)')
  ylabel('log_(10)(Number of Scans (each pump counted separatly))')
  title('EDDIES 07/05 Injection')
  ax = axis;
  axis([26.25 26.265 0 ax(4)])
  %print
  
  figure(2)
  bar(thetas,ntheta)
  xlabel(['Theta (^o C)'])
  ylabel('Number Scans (each pump counted separatly)')
  title('EDDIES 07/05 Injection')
  ax = axis;
  axis([19.23 19.3 0 ax(4)])
  %print
    
  figure(3)
  bar(sals,nsal)
  xlabel(['Salinity (PSU)'])
  ylabel('Number Scans (each pump counted separatly)')
  title('EDDIES 07/05 Injection')
  ax=axis;
  axis([36.699 36.71 ax(3:4)])
  %print
    
  figure(4)
  bar(presses,npress)
  xlabel(['Pressure (dbar)'])
  ylabel('Number Scans (each pump counted separatly)')
  title('EDDIES 07/05 Injection')
  ax=axis;
  %axis([36.699 36.71 ax(3:4)])
  %print
    
  figure(5)
  bar(fls,log10(nfl))
  xlabel(['Chlorophyll (/mug/l)'])
  ylabel('log_(10)(Number of Scans (each pump counted separatly))')
  title('EDDIES 07/05 Injection')
  ax=axis;
  %axis([36.699 36.71 ax(3:4)])
  %print

  figure(6)
  bar(fls,(nfl))
  xlabel(['Chlorophyll (/mug/l)'])
  ylabel('Number of Scans (each pump counted separatly)')
  title('EDDIES 07/05 Injection')
  ax=axis;
  %axis([36.699 36.71 ax(3:4)])
  %print

  figure(7)
  bar(sigs,nsigs)
  xlabel('Sigma Theta(kg/m^3)')
  ylabel('Number of Scans (each pump counted separatly)')
  title('EDDIES 07/05 Injection')
  ax = axis;
  %axis([26.25 26.265 0 ax(4)])
  %print

  figure(8)
  bar(sigs,spd_sum/sum(spd_sum))
  xlabel('Sigma Theta(kg/m^3)')
  ylabel('Nscans Weighted by pump speed ')
  title('EDDIES 07/05 Injection')
  ax = axis;
  %axis([26.25 26.265 0 ax(4)])
  %print

end

  mean_sig = sum(nsigs.*sigs)/sum(nsigs);
  rms_sig = sqrt(sum( ((sigs-mean_sig).^2).*nsigs)/sum(nsigs));

  mean_sigwt = sum(spd_sum.*sigs)/sum(spd_sum);
  rms_sigwt = sqrt(sum( ((sigs-mean_sigwt).^2).*spd_sum)/sum(spd_sum));

  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));


  diary hist_stat.txt
%  disp([' Mean sigma: ',num2str(mean_sig), ...
%	' RMS: ',num2str(rms_sig)])
  
  disp([' Mean weighted sigma: ',num2str(mean_sigwt), ...
	' RMS: ',num2str(rms_sigwt)])
if 0  
  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)])
end
  diary off




