% Code to recalculate chi using the code from the Slade and Boss 2006
% author: Wayne Slade
% date: July 16, 2017

% Ringarea file provided by Sequoia (Type B)
dcal = [1.0000000e+000  1.0038000e+000  9.9360000e-001  1.0027000e+000  9.9720000e-001  9.9570000e-001  9.9030000e-001  9.9430000e-001  9.9290000e-001  9.9000000e-001  9.9290000e-001  9.9300000e-001  9.9150000e-001  9.9300000e-001  9.9230000e-001  9.9090000e-001  1.1032000e+000  1.1123000e+000  1.2430000e+000  1.1562000e+000  1.3273000e+000  1.1999000e+000  1.0740000e+000  1.7489000e+000  1.5382000e+000  2.5109000e+000  2.5468000e+000  3.5504000e+000  3.9338000e+000  5.1747342e+000  7.5143548e+000  1.2528083e+001];

% VSF calibration as of March 2006
chi =   [3.0528348e-006  3.2074788e-006  3.1009823e-006  3.2282882e-006  3.0529240e-006  3.2936448e-006  2.7994703e-006  3.0813466e-006  3.0734541e-006  3.1125205e-006  3.1332614e-006  3.1428713e-006  3.1239683e-006  3.1443005e-006  3.3732488e-006  3.3661327e-006  3.2796991e-006  3.3096158e-006  3.3034023e-006  3.3257584e-006  3.3475003e-006  3.3049492e-006  3.3323031e-006  3.2894461e-006  3.2840309e-006  3.1988049e-006  3.1385122e-006  3.0374149e-006  2.9298281e-006  2.7713090e-006  2.2171952e-006  1.7589167e-006];

% Instrument pathlength [meters]
pathl = 0.05;
      
% phi factor for ring area calculation
phi = 1/6;
  
% used to determine ring angles (Type B)
theta0air = 0.1/180*pi;

% Both particle and background should be log files
rawdata = load('total_seawater.log')
% size of data file
[ns,col] = size(rawdata);

rawzscat = load('filtered_seawater.log')
zscat = mean(rawzscat);

% First calculate transmission and attenuation    
% rr is ratio of water-transmitted to reference laser counts in clean water
%   used to correct transmission data
rr = zscat(33)/zscat(36);
% transmission and beam attenuation
tau = rawdata(:,33)./rr./rawdata(:,36);
cpg = -log(tau)/pathl;

% Ring angle calculation
theta_i = (200.^((0:32)/32)).*theta0air;
theta_j = theta_i(2:33);
theta_i = theta_i(1:32);
theta_i = asin(sin(theta_i)./1.3308); 
theta_j = asin(sin(theta_j)./1.3308); 
% theta_i and theta_j are inner and outer ring angles IN WATER)
theta = geomean([theta_i; theta_j]);

% "calibrated scatter" implies that we've corrected the raw scatter for
% the detector area (determined based on geometry), zscat, dcal
% (responsivity correction, aka ringarea file), and attenuation 
% within the sample volume:
% (1) first correct for attenuation (using tau) and subtract zscat...
%     code in MatlabBundle.zip/nlia_func.m also has a z(36)/data(1,36)
%     factor to correct for drift in laser power; we make the 
%     correction later similar to vsfcode_new.m
cscat = rawdata(:,1:32)./repmat(tau,1,32) - repmat(zscat(1:32),ns,1);
cscat = cscat .* repmat(dcal,ns,1);
% (3) correct for laser power (counts) in water (relative to zscat)
cscat = cscat ./ repmat(rawdata(:,36)/zscat(36),1,32);
% (4) apply analytic formulation for ring area
pscat = cscat ./ repmat(pi*phi*pathl*(theta_j.^2-theta_i.^2),ns,1);
% (5) Apply VSF calibration factor for each ring (pscat has correct
% shape)
beta = repmat(chi,ns,1) .* pscat;

% Estimate of total particle scatter based on integration over
% LISST-B near-forward angles
bp = repmat(nan, size(beta,1),1);
for ki = 1:size(beta,1)
  bp(ki) = 2*pi*integral(theta, sin(theta).*beta(ki,:));
end