% function program to convert data into VSF, see also Application Note L012 at our website
%
%	Usage:
%
%		vsf=makevsf(datafile,zscfile,RingArea_file,laser_ref_cal,i_type)
%
%	where scat is the scattering data, instrument_type defines the angles and gains
%	 use value 1 for 5-500 microns, 2 for 1.2 to 250 microns
%   laser_ref_cal is counts in channel 36 per mW laser power in water 
%	

%		function  [vsf,b,beam-c]=makevsf(datafile,zscfile,dcal_file,laser_ref_cal,i_type)
      
clear all
datafile='pseudonana_102003.log';
zscfile='SWF_102003.asc';
dcal_file='ringarea_1102.asc';
laser_ref_cal=.66;
i_type=2;
%C=3125;
C=6250;

      
      % establish instrument type and related parameters
      if (i_type ==1)
         rho=100^(1/32); f=12.0; theta0=0.85e-3;l=5;
      elseif(i_type ==2)
         rho=200^(1/32); f=6.0; theta0=1.7e-3;l=5;
      elseif(i_type ==3)
         rho=200^(1/32); f=12.0; theta0=0.85e-3;l=5;
      end
%_____

%zsc=mean(tt2mat(zscfile,40));
%s=num2str(datafile);
%if(s(1:4)=='cool')
%  data=readfile(datafile);
%else
 %  data=tt2mat(datafile,40);
%end

zsc=load(zscfile);
data=load(datafile);


% The value of laser_ref_cal is = laser_power_in_water/zsc(36); to get laser_power_in_water in mW, call factory.

eval(['load ',dcal_file])
s=num2str(dcal_file);
i=find(s=='.');
eval(['dcal= ',s(1:i-1),';']);
ra=zsc(33)/zsc(36);
tau=data(:,33)./ra./data(:,36);
beam_c=-log(tau)/0.05;
%beam_c=-log(tau)/0.025;
mean_c=mean(beam_c)
laser_power_in_water=data(:,36)/zsc(36)*laser_ref_cal; %CHANGED, DIVIDE BY ZSC(36) this is P0
%laser_power_in_water=data(:,36)*laser_ref_cal*1e-3; %units Watts;   

%abc=(data(i,1:32)/(exp(-(beam_c*0.05))*0.05*laser_power_in_water(i)*pi*6*((rho^2-1).*rho.^([0:2:62])*(theta0^2))))*5/4095/5e5/0.4

const=(rho^2-1).*rho.^([0:2:62]);const=1./const;
const=flipud(const);
const=const/(theta0^2)*6*20/pi; % per eq(8) in memo "Measuring the Vol Scat Func.doc"
[r,c]=size(data);

for i=1:r
   scat(i,:)=data(i,1:32)/tau(i)-zsc(1:32)';
   %scat(i,:)=data(i,1:32)-zsc(1:32)';
   cscat(i,:)=dcal.*scat(i,:);
   %vsf(i,:)=const.*cscat(i,:)/tau(i).*5/4095/5e5/0.4;  % in Watt; 5V/4095counts/0.5Mohm/0.4(A/Watt)
   vsf(i,:)=const.*cscat(i,:)*5/4095/5e5/0.4;
   vsf(i,:)=vsf(i,:)/laser_power_in_water(i);
   pf(i,:)=vsf(i,:)/mean(beam_c);
   j=find(vsf(i,:)<=0);  % at times for very fine particles, small rings may show negative net scattering due to noise
   vsf(i,j)=-1*vsf(i,j); % **** fix negative numbers to avoid messing up logplots*****
end   
for i=1:r
  
   theta=[0.106:0.1:17.906]*pi/180;
   %ll=(rho^2-1).*rho.^([0:2:62]);const=1./const;
	%ll=flipud(const);
	%ll=const/(theta0^2)
   thetam=[0.106 0.125 0.148 0.174 0.206 0.243 0.287 0.338 0.40 0.47 0.56 0.66 0.77 0.91 1.08 1.27 1.5 1.77 2.09 2.46 2.91 3.43 4.05 4.78 5.64 6.65 7.85 9.26 10.93 12.90 15.22 17.96];
   b_new(i,:)=interp1(thetam*(pi/180),vsf(i,:),theta,'linear');
   b_new_factor(i,:)=b_new(i,:).*sin(theta);
   b(i)=2*pi*trapz(b_new_factor(i,:)).*(0.1*(pi/180));
   %bb(i)=b(i)/2*pi;
   mean_b=mean(b(i))
   pf2(i,:)=vsf(i,:)/(mean_b);
   %b_yogi(i)=4.3e-4*sum(scat(i,:).*dcal(1:32))/laser_power_in_water; %here is Yogi's new formula

 end  
 plot(theta*180/pi,mean(b_new), '-',thetam,mean(vsf),'*')
 %phase=vsf(i,:)/b(i)
 
% figure out denominator for this instrument
meanvsf=mean(vsf);
stdvsf=std(vsf);
meanpf=mean(pf2);
stdpf=std(pf2);
mean_c=mean(beam_c);

angles = [0.106 0.125 0.148 0.174 0.206 0.243 0.287 0.338 0.40 0.47 0.56 0.66 0.77 0.91 1.08 1.27 1.5 1.77 2.09 2.46 2.91 3.43 4.05 4.78 5.64 6.65 7.85 9.26 10.93 12.90 15.22 17.96];
              
 % VSF inversion to PSD
[r,c]=size(cscat);
load m2kb_v.asc;     % matrix produced by Mie Thoery for polystyrene beads, provided by Sequoia
load Dias32_b.asc;  %  diameters at which the vd is computed, from Sequoia
load H2kb_v.asc;
for i=1:r
vd(i,:)=nlia(cscat(i,:),m2kb_v,H2kb_v);
end
% Volume dist (vd) divided by "volume conversion" factor C = volume concentration = "particle size distribution"
psd=median(vd)/C;


%   Do some plotting to visualize the data
figure(1);
subplot (2,1,1); loglog(angles,meanvsf,'k.-');
hold on
loglog(angles,meanvsf+stdvsf,'--');
loglog(angles,meanvsf-stdvsf,'--');
hold off
title('Volume Scattering Function at small angles measured with the LISST\bf');
xlabel('angle [degrees]'); ylabel('VSF [m^-^1 sr^-^1]');
subplot (2,1,2); loglog(angles,meanpf','r.-');
hold on
loglog(angles,meanpf+stdpf,'--');
loglog(angles,meanpf-stdpf,'--');
hold off
title('Phase function at small angles measured with the LISST\bf');
xlabel('angle [degrees]'); ylabel('phase function [sr-1]');

figure(2);
subplot(2,1,1);semilogx(Dias32_b,psd,'k.-'); 
title('Particle size distribution measured with the LISST');
xlabel('Mean particle diamteter [microns]'); ylabel('volume concentration [ul/l]')
subplot(2,1,2);
plot(mean(cscat),'b.-');
title('Mean scattering in each ring');
xlabel('ring number'); ylabel('counts');
 
if(0)
figure(3);
% subplot(2,1,1);
errorbar(angles,meanvsf,stdev,'k');
xlabel('angle [degrees]'); ylabel('VSF and standard deviation [m^-^1 sr^-^1]');
title('Standard deviaition of the VSF');
end
A=[angles; meanvsf];
%save vsf_maalox3.txt A -ascii
