function [ceil_utc,ceil_date,ceil_range,ceil_b,cloudcode,cloud1,cloud2,cloud3,ceilo_30_sec]=LCload_CL31(date,ceil_filename)
% LCload_CL31: load a Lidar Ceilometer CL31 file and return time series of
% cloud heights and a time-height series of backscatter.
%
% USAGE:
%  [ceil_utc,ceil_date,ceil_range,ceil_b,cloud1,cloud2,cloud3,ceilo_30_sec]=LCload_CL31(date,ceil_filename);
%
% INPUTS: 
%    date=[yyyy mm dd];
%    ceil_filename is the full-path name of the file to load
% OUTPUTS:
%    ceil_utc is the time in decimal hours of each ray
%    ceil_date is the date of each ray as a float that looks like yyyymmdd 
%    ceil_range is the range of each gate in metres
%    ceil_b is the backscatter coefficient in the units (srad km)**(-1)
%    cloud1, cloud2 and cloud3 are the heights of the first three cloud
%         bases as determined by the CL31 data aquisition system from
%         the backscatter profile.
%
% AUTHORS:
% Dana, Chris, Dan, Ludovic, and Simon
% 
% 2008 September 17
% -eliminated hex2dec calls, sped code up ~30x
% -declared temporary arrays in fast orientation, truncated at end of read
% -fixed feet/meters unit conversion
% -cleaned up loop and logic
% -fixed 2's complement unwrapping for 20-bit integers
% -reads range gates from file
% -returns ceil_b in VOLUME BACKSCATTER COEFFICIENT (p.75 CL31 manual)
% -alarms flagged with cloudcode=-1

% %for testing:
%  clear; close all;
%  date=[2007 10 17];
%  ceil_filename='C:\Data\STRATUS_2007\RHB\ceilometer\Raw\A7101700.DAT';

nbits=20; % number of bits for reflectivity
resolution=10; %10m resolution, to be updated from file
gates_number=770;  %number of gates, to be updated from file
%resolution*gates_number;   %10m x 770 samples, range 7700m (ms1_10x770)
nn=7000; %largest possible number of samples --should be exactly 5760 samples, but who knows?
start_height = 0; % in m % SPdeS ?should be height of ceilometer on deck?

% fix single digit years in date (e.g. 7-->2007):
if date(1)<70
    date(1)=date(1)+2000;
elseif date(1)>=70 && date(1)<1000
    date(1)=date(1)+1900;
end

ceil_year=date(1);
ceil_month=date(2);
ceil_day=date(3);

% Declare Arrays
% transposed for faster indexing, SPdeS
% Since we don't know the final size, 
% use temp arrays and truncate copy to final arrays at end.
tempcloud1=NaN+zeros(nn,1);
tempcloud2=NaN+zeros(nn,1);
tempcloud3=NaN+zeros(nn,1);
tempcloudcode=NaN+zeros(nn,1);

tempceil_sec=NaN+zeros(nn,1);
tempceil_min=NaN+zeros(nn,1);
tempceil_hour=NaN+zeros(nn,1);
tempceil_b=NaN+zeros(gates_number,nn);

% open file
fid=fopen(ceil_filename,'r');
if fid == -1
    error(['Cannot read ' ceil_filename]);
end

% frewind(fid) % start reading at top of file here
% dump header
dmp=fgetl(fid);   %dump line -Ceilometer Logfile
dmp1=fgetl(fid);  %dump line -File created: 2/24/2007 12:00:00 AM

% Start gargantuan read loop
dmp
k=0; % count time samples
while ~feof(fid)
    %for pp = 1:20  % short loop for testing
    if strcmp(fscanf(fid,'%s',[1 1]),'-File'); break; end; %dump -2007-02-24 and break on last '-File Closed:' line
    k=k+1;
    time=fscanf(fid,'%2d:%2d:%2d',[3 1]); %dump date, read hh:mm:ss
    if ~isempty(time);% catch other errors and empty at last line of file
        tempceil_hour(k)=time(1);
        tempceil_min(k)=time(2);
        tempceil_sec(k)=time(3);

        fgetl(fid);       %go to end of line
        temp=fgetl(fid);  %get rid of 
CL013011
        %range_flag=str2num(temp(8));
        tempht=fgetl(fid); %cloud base heights
        if tempht(1)=='/' || tempht(2)=='A' % Alarms present!
            tempcloudcode(k)=-1; % flag error
        else % no alarms present
            tempcloudcode(k)=sscanf(tempht(1:2),'%1d%*1c'); % data is number of clouds detected: 0,1,2, or 3
            if tempcloudcode(k) < 4 % codes 4:fully obscured, no base detected; 5:some obscuration, but determined transparent
                if tempcloudcode(k) > 0
                    tempcloud1(k)=sscanf(tempht(4:20),'%5d',1);
                    if tempcloudcode(k) > 1
                        tempcloud2(k)=sscanf(tempht(10:20),'%5d',1);
                        if tempcloudcode(k) == 3
                            tempcloud3(k)=sscanf(tempht(16:20),'%5d',1);
                        end
                    end
                end
            end
        end % no alarms present
        % position 21 is a space before the detection status bits in 12
        if k<2 % only check detection status once per file
            Detection_status=tempht(22:22+12-1); % hexadecimal digits
            units_bit=sscanf(Detection_status(11),'%1x',1)>7;
            % units_bit true if meters; see table of hexadecimal digits and detection bits in ceilo manual
            if units_bit  %units in cloud heights line are meters
                mm=1;     %no conversion coefficient
            else          %units in cloud heights line are feet
                mm=0.303; %conversion coefficient from feet to meters
            end;
        end
        %line2p5=fgetl(fid); % line2.5: -1 ///  0 ///  0 ///  0 ///  0 ///  This is the extra line put in for CL31
        line3=fgetl(fid);   % line3: 00100 10 0770 100 +28 097 02 0005 L0016HN15 238
        if k<2 % update scale and resolution once per file
            scaleres=sscanf(line3,'%5d %2d %4d');
            scale=scaleres(1)/100;
            resolution=scaleres(2);
            gates_number=scaleres(3);
        end
        % read the long line of ranged data...
        line = fgetl(fid);
        if size(line,2)==5*gates_number
            tempceil_b(:,k)=sscanf(line,'%5x');
        else
            tempceil_b(:,k)=NaN;
        end;
    end %if ~isempty(time)
    fgetl(fid); % dump 6 character line
    fgetl(fid);  %get rid of blank line (if k>1?)
end % while

fclose(fid);

% truncate data  x=tempx(1:k)
ceil_b=tempceil_b(:,1:k);
% wrap 2s complement for nbits=20
ceil_b2=ceil_b;
index=find(ceil_b>=2^(nbits-1));
if ~isempty(index), ceil_b2(index)=ceil_b(index)-2^nbits;end
% % for comparison puposes, compute backscatter with wrong 16-bit 2s complement
% ceil_bwrong=ceil_b;
% ceil_bwrong(ceil_b>=2^15)=ceil_b(ceil_b>=2^15)-2^16;
% ceil_bwrong=ceil_bwrong*1e-5;

if scale==1
    ceil_b=ceil_b2*1e-5; % scaled by default in units 1e5/(km srad)
else
    ceil_b=ceil_b2*scale*1e-5;
end
% ceil_b is now VOLUME BACKSCATTER COEFFICIENT [units 1/(km*steradians)].
% It is the extinction-normalized power backward-scattered per
% atmospheric scattering volume drdOmega [m*srad]. This is
% related to the power transmitted and recieved at the instrument
% by the lidar equation (CL31 manual p. 75).
% In practice it is sometimes negative, presumably because of noise.

% convert units to meters if necessary
if ~units_bit
    cloud1=tempcloud1(1:k)*mm;    %feet to meters
    cloud2=tempcloud2(1:k)*mm;
    cloud3=tempcloud3(1:k)*mm;
else
    cloud1=tempcloud1(1:k);
    cloud2=tempcloud2(1:k);
    cloud3=tempcloud3(1:k);
end
cloudcode=tempcloudcode(1:k);
ceil_sc=tempceil_sec(1:k);
ceil_mn=tempceil_min(1:k);
ceil_hr=tempceil_hour(1:k);
ceil_utc=tempceil_hour(1:k)+tempceil_min(1:k)/60.0+tempceil_sec(1:k)/3600.0;
ceil_date(1:k,1)=ceil_year*10000+ceil_month*100+ceil_day;

ceil_range=(start_height+(1:gates_number)*resolution)/1000.0; % km

% % plot
% %subplot(2,1,1)
% imagesc(ceil_utc,ceil_range,ceil_b); axis xy; colorbar;
% xlabel('hour UTC'); ylabel('range (km)'); title({num2str([date(1) date(2) date(3)]),'CL31 log_{10}backscatter (m^{-1} srad^{-1})'});
% % plot the wrong 2s complement ceil_b:
% % subplot(2,1,2)
% % imagesc(ceil_utc,ceil_range,log10(max(0,ceil_bwrong))); axis xy; colorbar;
% % xlabel('hour UTC'); ylabel('range (km)'); title('CL31 log_{10}backscatter (m^{-1} srad^{-1})');

% ceilo_30_sec is data in the form required by the calling program
ceilo_30_sec=[cloudcode cloud1 cloud2 cloud3 ceil_hr ceil_mn ceil_sc];% this is the raw 1-D time series data

return