% 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);
resolution=10; %10m resolution
gates_number=770;  %number of gates
nn=resolution*gates_number;   %10m x 770 samples, range 7700m (ms1_10x770)

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_date=NaN*zeros(1,nn);

ceil_b=NaN*zeros(gates_number,nn);
ceil_utc=NaN*zeros(1,nn);
height=[];

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

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

%dmp=fgetl(fid);   %dump line -Ceilometer Logfile 
%dmp1=fgetl(fid);  %dump line -File created: 2/24/2007 12:00:00 AM
%%DEW
%dmp=fgetl(fid)  %-File Closed: 5/1/2007 12:22:15 PM  
dmp=fgetl(fid)   %dump line -Ceilometer Logfile 
dmp1=fgetl(fid)  %dump line -File created: 2/24/2007 12:00:00 AM
%dmp4=fscanf(fid,'%s',[1 1])  %dump -2007-02-24
%[data]=fscanf(fid,'%2d:%2d:%2d \n',[3 1]); %read 00:00:12 temp=fgetl(fid);  %get rid of ''
%temp=fgetl(fid);  %get rid of _CL013011_
%range_flag=str2num(temp(8));
%temp=fgetl(fid);  %cloud heights 1W 23570 ///// ///// 000020000000 
%fgetl(fid); % dump 00100 10 0770 100 +28 097 02 0005 L0016HN15 238
%%DEW
ict = 0;
while ~feof(fid)
%for pp = 1:20  % short loop for testing
    ict = ict + 1;
    dmp4=fscanf(fid,'%s',[1 1])  ;%dump -2007-02-24
    [data]=fscanf(fid,'%2d:%2d:%2d \n',[3 1])  %read hh:mm:ss
    if ~isempty(data) 
        ceil_hour=data(1);
        ceil_min=data(2);
        ceil_sec=data(3);
%k;
        temp=fgetl(fid);  %get rid of 
CL013011
        %range_flag=str2num(temp(8));
        tempht=fgetl(fid); %cloud heights
%         [data,count]=sscanf(temp,'%1d%*1c %d %d %d') ;    error if status bit contains no characters...(like 00001000)
        [data,count]=sscanf(tempht(1:(length(tempht)-12)),'%1d%*1c %d %d %d'); 
        data;
%         units_flag=sscanf(temp,'%*21c %12x');
        Detection_status=tempht((length(tempht)-12):(length(tempht))-1);
        units_flag=dec2bin(hex2dec(Detection_status(11)),4);
        fgetl(fid); % dump 00100 10 0770 100 +28 097 02 0005 L0016HN15 238

        %%%***********************  This does not appear to correctly
        %%%identify the resolution based on comparison of backscatter and
        %%%cloud 
        if strcmp(units_flag,'1000')==1  %units in cloud heights line are meters
            mm=1;   %no convertion coefficient
        elseif strcmp(units_flag,'1000')==0  %units in cloud heights line are feet
             mm=0.303; %convertion coefficient for feet to meters
        end;
        
        %Next line commented out 5/3/2007 DEW
        %fgetl(fid) % dump 00100 10 0770 100 +28 097 02 0005 L0016HN15 238

%%%%%%%%%%%%%%%%%%%% Start reader the data %%%%%%%%%%%%%%%%%%%
start_height = 0; %%% in m %%%%%%%%
line = fgets(fid);
kke= size(line); 
if kke(1,2)==3852
    kke(1,2);  %Print out line length
    for i=1:770
        i;
        height(i)=start_height+resolution*(i);
        ceil_b(i,k+1)=hex2dec(line((5*(i-1)+1):(5*(i-1)+1+4)));
    end
else
    for i=1:770
        height(i)=NaN;
        ceil_b(i,k+1)=NaN;
    end;
end;
% 
% 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
            data(1);
            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); 
    if ict > 1
temp=fgetl(fid);  %get rid of ''
end
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:gates_number)*resolution/1000.0;

cloud1=cloud1(1:k)*mm;    %output is meters
cloud2=cloud2(1:k)*mm;
cloud3=cloud3(1:k)*mm;
% cloud1=cloud1(1:k);
% cloud2=cloud2(1:k);
% cloud3=cloud3(1:k);
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