%  7th Edition of NLIA_SUB  (03/23/01)
%
%	Algorithm to invert a selection of scans from a file,
%
%	usage:
%
%      [vd,tau,dias32,scat,mzsc, press,temp,aux,bvolt,vdate]=nlia_sub('datafile','zscfile','dcalfile',is,ie,navg,cal_const,instrument_type);
%
%  e.g. 
%      [vd,tau,dias32,scat,mzsc, press,temp,aux,bvolt,vdate]=nlia_sub('c1041_12.dat','z1041_04.dat','dcal_c.mat',23,95,10,1,2);
%
%	where analysis begins at scan # is, ending at ie.[ in the example above is=23, ie=95]
%  If the ending scan no. is not known, use ie=-1, program will find the final scan automatically
%	vd is volume distribution; tau is optical transmission;	scat is fully calibrated scattering
%	instrument type is 1 for the 5-500 micron version, 2 for 1.2-250 micron version,3 for 2.5-500
%	navg is number of scans to use in the average
%	beam_c can be computed by c=-20*log(transmission) in m^-1
%	tau is attenuation time series,
%	dias32 are the 32 sizes (logspaced)
%	scat is the fully corrected scattering signature
%	mzsc is the mean zscat profile
%
%  press is the pressure in rough psi
%  temp is the temperature in deg C
%  Aux is the Aux analog channel
%  bvolt is the battery voltage
%  date is a two column variable holding date is whatever fashion was specified at deployment
%
%  *** Needs files m32x32_X.asc, dias32_X.asc, and H32_X.asc in \matlab directory
%  *** dcalfile is detector calibration, also called ringarea.asc
%  *** dcalfile may be a .mat file containing variable dcal, or an ascii file with 32 elements.
%  Modified for Version 5 by HCP on 7/13/98
%  Indexing fixed, change to ASCII input, and z made row vector by LKS on 6/26/00
%  String extraction fixed by LKS 9/18/00 to allow paths to be included
%  Added five new output variables that pass the aux values out of the array, by lks 9/29/00
%
%  Added option to use new 2k matricies for type "b" and "c" instruments 

	function [vd,tau,dias32,scat,mzsc, press,temp,aux,bvolt,vdate]=nlia_sub(datafile,zscfile,dcalfile,is,ie,navg,cal_const,inst_type)
   
   %  constants used for the aux params
   %  these are values for a specific instrument and are "typical"
   pscale = 0.08368;
   poffset = -0.5;
   raux = 5/4096;
   
   %  increase indexes of frames when the LISST software starts with Data Frame #0
   %    Note:  intension is to eventually modify the LISST program to index from 1 at
   %           which time this coding should be removed
   if(ie>=is) 
      ie = ie + 1;
   else
      ie = -1;
   end
   is = is + 1;
   %
   [is ie]
%%   pause
   
   %  load the appropriate files; use ASCII versions to be in perfect alignment with LISST software
   	if inst_type==1          
      	load M32x32_a.asc;
      	load Dias32_a.asc;
         m2=M32x32_a;dias32=Dias32_a';
         load H32_a.asc;
 		   H=H32_a;
        elseif inst_type==2 
      	load m2kb_v.asc;
      	load H2kb_v.asc;
      	load Dias32_b.asc;
         m2=m2kb_v;dias32=Dias32_b';
 		 	H=H2kb_v;    
		elseif inst_type==3 
      	load m2kc_v.asc;
      	load H2kc_v.asc;
      	load Dias32_c.asc;
         m2=m2kc_v;dias32=Dias32_c';
 		 	H=H2kc_v;    
      end
        
% ____I____Open ringarea file, also called dcalfile
      
   eval(['load ',dcalfile]);
	s=num2str(dcalfile);
   k=find(s=='.');
   k2 = 0;
   for j=1:k-1,
      if (s(j:j)=='/' | s(j:j)=='\')
         k2 = j;
      end;
   end;
   ftype=s(k+1:k+3);
   if (ftype=='asc')
      eval (['dcal=',s(k2+1:k-1);]);
   end
   
   
   %  make sure decal is a row vector
   [r,c]=size(dcal);
   if (r>c)
      dcal=dcal';
   end
   
figure(1),hold off,plot(dcal),title('The Detector Calibration Array'),pause(1)
% ____II____ Open zsc file
   
   s=num2str(zscfile);
   k=find(s=='.');
   k2 = 0;
   for j=1:k-1,
      if (s(j:j)=='/' | s(j:j)=='\')
         k2 = j;
      end;
   end;
   ftype=s(k+1:k+3);
   if (ftype=='asc')
      eval(['load ',zscfile;]);
      eval (['z=',s(k2+1:k-1);]);
   else
   %   load zscfile in binary format
     zsc=tt2mat(zscfile,40);
     z=mean(zsc);
   end

%% make sure z is a row vector (must match dim of data(i,1:32) below)
[rz, cz] = size(z);
if (cz<rz)
   z = z';
end
%
figure(2),hold off,plot(z(1:32)),title('The Current ZSCAT'),pause(0.5)
% determine the ratio, r, of the water-transmitted and reference laser power from zscat.
		r=z(33)/z(36);

% ____III____ Open datafile, binary format

   s=num2str(datafile);
   k=find(s=='.');
   k2 = 0;
   for j=1:k-1,
      if (s(j:j)=='/' | s(j:j)=='\')
         k2 = j;
      end;
   end;
   ftype=s(k+1:k+3);

   if (ftype=='log' | ftype=='asc')
      dm=load(datafile);
   else
      dm=tt2mat(datafile,40);
   end
   [row,col]=size(dm);
   
% ____IV____ Get data corrected for zsc and with dcal applied; then invert
%  first, if the ending scan number is not known, then set the ending scan no.
	if (ie<=0)
     ie=row;
	end
   
   if (ie==is)
      data(1,:)=dm(is,:);
      tau(1)=data(1,33)/r./data(1,36);
      dc(1,:)=(z(36)/data(1,36)*data(1,1:32)./tau(1))-z(1:32);
      dc_cal(1,:)=dc(1,:).*dcal;
      
      %  calculate aux params for this group
      press(1) = poffset + pscale*data(1,37);
      temp(1) = 0.01 * data(1,38);
      aux(1) = raux * data(1,35);
      bvolt(1) = 0.01 * data(1,34);
      vdate(1,:)=[dm(is,39) dm(1,40)];	% store date of first scan in averaged group
   else   
      for i=1:(ie-is+1)/navg
         i1=is + navg*(i-1);
         i2=is + i*navg - 1;
         if (i1==i2)
            data(i,:)=dm(i1,:);
         else
            data(i,:)=mean(dm(i1:i2,:));
         end
         %  calculate aux params for this group
         press(i) = poffset + pscale*data(i,37);
         temp(i) = 0.01 * data(i,38);
         aux(i) = raux * data(i,35);
         bvolt(i) = 0.01 * data(i,34);
         vdate(i,:)=[dm(i1,39) dm(i1,40)];	% store date of first scan in averaged group
         
         %  continue with calculations of PDF
         tau(i)=data(i,33)/r./data(i,36);
         dc(i,:)=(z(36)/data(i,36)*data(i,1:32)./tau(i))-z(1:32);
         dc_cal(i,:)=dc(i,:).*dcal;
      end
   end
   
%   
scat=dc_cal;
%%scat(:,18)=0.5*(scat(:,17)+scat(:,19)); % remove this at some time
mzsc=z;
%
% ____V____________ starting inverse____________

	[r,c]=size(scat);

	   for ik=1:r
         ad(ik,:)=nlia(scat(ik,:),m2,H)/cal_const;
         if inst_type > 1 
            vd(ik,:)=ad(ik,:);
         else
            vd(ik,:)=ad(ik,:).*dias32;
         end
	   end
      
figure(1)
clf,hold off
%	subplot(222),semilogx(dias32,ad'),title('Areas')
subplot(221),plot(scat'),axis([0 33 min(min(scat)) 1.2*max(max(scat))])
title(['Scattering ',datafile] )
subplot(222),semilogx(dias32,vd'),axis([min(dias32) max(dias32) 0 1.2*max(max(vd))])
if (is==ie)
   title(strcat('Volume Distribution for Frame #',num2str(is-1)))
else
   title(strcat('Volume Distribution for Frames #',num2str(is-1),'-',num2str(ie-1)))
end
	subplot(223),plot(z(1:32)),axis([0 33 0 1000]),title('ZSCAT'),text(10,800,['laser =',num2str(z(33))])
%________

	
