% LCload: load a Lidar Ceilometer file into the matlab workspace
%d
% 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);
%%% For the CL31
%resolution=10; %10m resolution
%gates_number=770;  %number of gates
%%% For the CT25k
resolution=30; %100ft resolution or 30m
gates_number=250;  %number of gates

%nn=resolution*gates_number;   %10m x 770 samples, range 7700m (ms1_10x770)  CL31
nn=resolution*gates_number;   %100ft x 250 samples, range 25000ft  CT25k 30m x 250 = 7500m


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

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');
%%%%%%%%%  CT25k ICEALOT ceilometer test data
%-Ceilometer Logfile
%-File created: 2/14/2008 19:39:30
%-2008-02-14 19:40:03
%
CT22023
%10 00490 ///// ///// 00000F00
%100 N 101 +24 ///  197  +0  202 LF7HN1 132
%00000480022001A00190017001A001A0017001600150019001A001A007F0740061D
%0160185003E00130009000100080002FFFDFFF80009000D00000002FFFF000F0016
%0320000001800120021001C0036001A001F00000008FFFFFFF6000900090041FFFD
%048001F000A000F0015FFF2002000390014004E0044006F00750023FFEC00230051
%0640061007F002D00840073008F006600300058004F0065007400480030FFFF007E
%080006D002D000E009C0020001F00AB00AC00AD00ADFFB40036008B00F7FFC4FF99
%09600FBFFD9FFD000F900950023FFC6FF29FFC3005FFF9700C8FFEB00AA00EEFFF1
%112003F004000ADFFF8FF7D0084005900100023004EFF8AFF13FF0EFFAE00300025
%1280043FF20FFC3FF48FE8100F7FFEEFEC5FF2C0112015A0025001600BFFE64FF61
%144004D00FD00D50066FE8200C700E3FF7B012FFF90008BFDE7FFDEFEC9FE6800CB
%16001500008FF7FFEA700820217FF1C0108020A00C200EB00F0FF52FE5BFDE9FE6E
%176FD3AFF8FFF6E01120106FF07FEB901C2014301190042029F01A3FE9700DEFC88
%19202AA01ADFFBCFD3FFFA3007D0080FFD2028FFED900B1FC6CFFF6FF0BFBCD02D1
%208FF3DFF60014D000700EA0279FCB2020EFC08028EFCF40380FF9101E1FECA01C8
%224006DFF8FFA00FF3AFB39FF7A0094FC81038C07B8036BFE5F01BE02060423FE2F
%240FD40FD74FE23FE47FED305E9FCC10512FB49FDB700D600000000000000000000
%
%
%-2008-02-14 19:40:18
%
CT22023
%1W 00490 ///// ///// 00040F00
%100 N 101 +24 ///  205  +0  202 LF7HN1 143
%0000009001F001B001B001B001E001C001D001800150020001F001D013C07210699
%01601BD0041000A0008FFFE00040001FFF1FFF800000000FFF1000C0001FFFA0000
%%%%%%%%%%%%%%%%%%%%%%%
if fid == -1
    error(['Cannot read ' ceil_filename]);
end

ict = 0;
while ~feof(fid)
    %for pp = 1:8256 % short loop for testing
    ict = ict + 1;
    if ict == 1
        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
        ceil_hour=data(1);
        ceil_min=data(2);
        ceil_sec=data(3);
        temp1=fgetl(fid);  %get rid of ''
        temp2=fgetl(fid);  %get rid of _CT22023_
        range_flag=str2num(temp2(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);
        %temp=fgetl(fid);  %cloud heights 1W 23570 ///// ///// 000020000000
        if ict == 4047
        tmp = fgetl(fid) % dump 00100 10 0770 100 +28 097 02 0005 L0016HN15 238
        end
        else
        dmp4=fscanf(fid,'%s',[1 1]);  %dump -2007-02-24
        if ict >404
            temp
            dmp4 % test line
        end
        [data]=fscanf(fid,'%2d:%2d:%2d \n',[3 1]); %read 00:00:12
        if length(data) == 0;
            break
        end
        data;
        ceil_hour=data(1);
        ceil_min=data(2);
        ceil_sec=data(3);
        %temp=fgetl(fid)  %get rid of ''
        temp2=fgetl(fid);  %get rid of 
CL013011
        range_flag=str2num(temp2(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); % 00 ///// ///// ///// 000000000080
        tmp = fgetl(fid);  %dump 00100 10 0770 100 +28 097 02 0005 L0016HN15 238
    end

    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;
    mm = 1;
    %Next line commented out 5/3/2007 DEW
    %fgetl(fid) % dump 00100 10 0770 100 +28 097 02 0005 L0016HN15 238

    %%%%%%%%%%%%%%%%%%%% Start reading the data %%%%%%%%%%%%%%%%%%%
    start_height = 0; %%% in m %%%%%%%%
    %%%%%%%%%%%
    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

    %  convert ht to meters if necessary
    height = height * mm;
    %%%%%%%%%%%%%%

    if 0 % skip this section DEW 2182008  for CL31 ceilometer
        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;
    end  % skip section
    %
    %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) > 1
        %This is a ht line from CT25 K08032000.DAT
        %50 ///// ///// ///// 00000900
        tempht;
        data(2);
        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 data not empty
    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;
%figure;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