% LCload: load a Lidar Ceilometer file into the matlab workspace
%
% USAGE:
%  [ceil_utc,ceil_date,ceil_range,ceil_b,cloud1,cloud2,cloud3]=LCload(data,ceil_filename);
%
% WHERE:
%   data=[yyyy mm dd];
%    ceil_filename is the name of the file to load
%    ceil_utc is the time in decimal hours of each ray
%    ceil_date is the date of each ray in the form yyyymmdd
%    ceil_range is the range of each gate in metres (in theory can be in feet sometimes)
%    ceil_b is the backscatter coefficient in the units 10**(-7) (srad m)**(-1)
%    cloud1, cloud2 and cloud3 are the heights of the first three cloud bases as
%         determined by the data aquisition system from the backscatter profile
%
% NOTE:
%    this program does not clean up the data


%function [ceil_utc,ceil_date,ceil_range,ceil_b,cloudcode,cloud1,cloud2,cloud3]=LCload(data,ceil_filename);
nn=4*721;

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


loop = 16;

cloud1=NaN*zeros(1,nn);
cloud2=NaN*zeros(1,nn);
cloud3=NaN*zeros(1,nn);
cloudcode=NaN*zeros(1,nn);

ceil_sc=NaN*zeros(1,nn);
ceil_mn=NaN*zeros(1,nn);
ceil_hr=NaN*zeros(1,nn);

ceil_b=NaN*zeros(256,nn);
ceil_utc=NaN*zeros(1,nn);

k=0;
fid=fopen(ceil_filename,'r');

if fid == -1
error(['Cannot read ' ceil_filename]);
end

%dmp=fgetl(fid);
%dmp1=fgetl(fid);

while ~feof(fid)
    dmp4=fscanf(fid,'%s',[1 1]);
    [data]=fscanf(fid,'%2d:%2d.%2d\n',[3 1]);
    if ~isempty(data) 
        ceil_hour=data(1);
        ceil_min=data(2);
        ceil_sec=data(3);
%k;
        temp=fgetl(fid);
        %range_flag=str2num(temp(8));
        temp=fgetl(fid); %cloud heights
        [data,count]=sscanf(temp,'%1d%*1c %d %d %d') ;
        units_flag=sscanf(temp,'%*21c %8x');

        fgetl(fid); % useless apparently


%%%%%%%%%%%%%%%%%%%% Start reader the data %%%%%%%%%%%%%%%%%%%
for a = 1:16;
    line = fgets(fid);
    kke= size(line);
    if kke(1,2)==69
        start_height = max([str2num(line(1:3))*100 NaN]); %%% in ft %%%%%%%%
        for j = 1:16;
            height((a-1)*16+j)=start_height+100*(j-1);
            ceil_b((a-1)*16+j,k+1)=hex2dec(line(3+4*(j-1)+1:3+4*j));
        end
    else
        for j = 1:16;
            height((a-1)*16+j)=NaN;
            ceil_b((a-1)*16+j,k+1)=NaN;
        end
    end

end

        %ceil_data=fscanf(fid,'%*3d%4x%4x%4x%4x%4x%4x%4x%4x%4x%4x%4x%4x%4x%4x%4x%4x',loop*16);

        k=k+1;
        %size(ceil_b)
        %ceil_b
        %k
        %ceil_hour
        %ceil_min
        %ceil_sec
        ceil_mn(k)=ceil_min;ceil_hr(k)=ceil_hour;ceil_sc(k)=ceil_sec;
        ceil_utc(k)=ceil_hour+ceil_min/60.0+ceil_sec/3600.0;
        ceil_date(k)=ceil_year*10000+ceil_month*100+ceil_day;
        %ceil_b(1:256,k)=ceil_data(1:256); % read in 1 `ray' of lidar backscatter
        
        if length(data) > 0
            cloudcode(k)=data(1);
            switch length(data)-1
            case 1, cloud1(k)=data(2);
            case 2, cloud1(k)=data(2); cloud2(k)=data(3);
            case 3, cloud1(k)=data(2); cloud2(k)=data(3); cloud3(k)=data(4);
            end
        end
        %fgetl(fid); 
        %fgetl(fid)
    end %if
    fgetl(fid); 
end % while
fclose(fid);

%clear data ceil_data

ceil_b=ceil_b(:,1:k);
index=find(ceil_b>=2^15);
if ~isempty(index), ceil_b(index)=ceil_b(index)-2^16;end
ceil_utc=ceil_utc(1:k);

% Convert metres to kilometers, assuming the units are metres not feet!
ceil_range=[1:256]*30.0/1000.0;

cloud1=cloud1(1:k)*.303;
cloud2=cloud2(1:k)*.303;
cloud3=cloud3(1:k)*.303;
cloudcode=cloudcode(1:k);
ceil_sc=ceil_sc(1:k);
ceil_mn=ceil_mn(1:k);
ceil_hr=ceil_hr(1:k);

%imagesc(log10(abs(ceil_b)));axis xy;
imagesc(ceil_utc,ceil_range,log10(abs(ceil_b)));axis xy;
height=ceil_range;
ceilo_30_sec=[cloudcode' cloud1' cloud2' cloud3' ceil_hr' ceil_mn' ceil_sc'];% this is the raw data with time
%save d:\data\ceilo\11_11_03.txt ceilo_30sec -ascii -tabs
%save d:\data\ceilo\11_11_03 ceil_utc ceil_b


