function [ceil_utc,ceil_date,ceil_hr,ceil_mn,ceil_sc,cloudcode,cloud1,cloud2,cloud3,ceil_range,ceil_b]=load_CT25_2009(date,ceil_filename)


% USAGE:
%  [ceil_utc,ceil_date,ceil_hr,ceil_mn,ceil_sc,cloudcode,cloud1,cloud2,cloud3,ceil_range,ceil_b]=load_CT25_2009(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 
%    cloud1, cloud2 and cloud3 are the heights of the first three cloud bases 
%    ceil_range is the range of each gate in metres
%    ceil_b is the backscatter coefficient in the units (srad km)^(-1)
% 
% Revision history, original version 02/02/2009
% put bacscatter unit to (srad km)^(-1) instead of (srad m)^(-1) --> to be
% consistent with CL31 program

% EXAMPLES:
%  clear; close all;
%  date=[2008 11 17];
%  ceil_filename='C:\A8111700DAT.txt';

nn=7000; %should be exactly 5760 samples (86400s/15s), but put more space just in case
gates_number=256; %number of gates
resolution=30; %30m resolution
nbits=16; % number of bits for reflectivity
start_height = 0; % in m

% Declare Arrays
tempcloud1=ones(nn,1)*NaN;
tempcloud2=ones(nn,1)*NaN;
tempcloud3=ones(nn,1)*NaN;
tempcloudcode=ones(nn,1)*NaN;
tempceil_sec=ones(nn,1)*NaN;
tempceil_min=ones(nn,1)*NaN;
tempceil_hour=ones(nn,1)*NaN;
tempceil_b=NaN+zeros(gates_number,nn);

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

% dump header
dmp=fgetl(fid);   %dump line -Ceilometer Logfile
dmp1=fgetl(fid);  %dump line -File created: 17.11.2008 0:00:00

%start read loop
k=0; % count time samples
while ~feof(fid)
    k=k+1;
    time=fscanf(fid,'%*11c %2d:%2d:%2d',[3 1]); %dump date, read hh:mm:ss
    if ~isempty(time);
        tempceil_hour(k)=time(1);
        tempceil_min(k)=time(2);
        tempceil_sec(k)=time(3);

        fgetl(fid);       %go to end of line
        fgetl(fid);       %get rif of CT02023 
        
        tempht=fgetl(fid); %cloud base heights  10 00860 ///// ///// 00000F00
        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+8-1); % hexadecimal digits
            units_bit=sscanf(Detection_status(6),'%1x',1)>=1;
            % 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
        
        line3=fgetl(fid); % 100 N 100 +21  74  219 +13    0 LF7HN1  38
        if k<2 % update scale and resolution once per file
            scaleres=sscanf(line3,'%5d');
            scale=scaleres(1)/100;
%             resolution=scaleres(2);
%             gates_number=scaleres(3);
        end

              
        % read the 16 lines of reflectivity data...
        for a = 1:16;
            line = fgets(fid);
            if size(line,2)==69              
                tempceil_b((a-1)*16+1:a*16,k)=sscanf(line(4:69),'%4x');
            else
                tempceil_b((a-1)*16+1:a*16,k)=NaN;
            end;    
        end;        
        
    end %if ~isempty(time)
    fgetl(fid); % dump 
    fgetl(fid);  %get rid of blank line
end % while

fclose(fid);


ceil_b=tempceil_b(:,1:k);
% wrap 2s complement for nbits=16
ceil_b2=ceil_b;
index=find(ceil_b>=2^(nbits-1));
if ~isempty(index), ceil_b2(index)=ceil_b(index)-2^nbits;end

if scale==1
    ceil_b=ceil_b2*1e-4; % scaling factor of 1e4/(km srad), see p37 of CL25 manual, data message no.2        
%     ceil_b=ceil_b2*1e-7; % scaling factor of 1e7/(m srad), see p37 of CL25 manual, data message no.2     
else
    ceil_b=ceil_b2*scale*1e-4;
%     ceil_b=ceil_b2*scale*1e-7;
end
% ceil_b is now VOLUME BACKSCATTER COEFFICIENT [units: 1/(km*steradians)].


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

% 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+(0:gates_number-1)*resolution); % m (profile range is 0 to 7500 m)

%do some tests to see if last points are not from the next hour...
tt=find(diff(ceil_hr)<0);
cloudcode(tt+1)=[]; 
cloud1(tt+1)=[]; 
cloud2(tt+1)=[];  
cloud3(tt+1)=[];  
ceil_hr(tt+1)=[];  
ceil_mn(tt+1)=[];  
ceil_sc(tt+1)=[]; 
ceil_b(:,tt+1)=[]; 
ceil_utc(tt+1)=[]; 
ceil_date(tt+1)=[]; 

%does the same for extra NaNs...
tt=find(~isfinite(ceil_hr));
cloudcode(tt)=[]; 
cloud1(tt)=[]; 
cloud2(tt)=[];  
cloud3(tt)=[];  
ceil_hr(tt)=[];  
ceil_mn(tt)=[];  
ceil_sc(tt)=[]; 
ceil_b(:,tt)=[]; 
ceil_utc(tt)=[]; 
ceil_date(tt)=[]; 
