% LCload: load a Lidar Ceilometer file into the matlab workspace
%
% USAGE:
%  [ceil_utc,ceil_date,ceil_range,ceil_b,cloud1,cloud2,cloud3]=LCload(ceil_filename);
%
% WHERE:
%    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
%
% Robin Hogan <R.J.Hogan@reading.ac.uk>
% 16 September 1999
%
%
% NOTE:
%    this program does not clean up the data

function [ceil_utc,ceil_date,ceil_range,ceil_b,cloud1,cloud2,cloud3]=LCload_new(ceil_filename);
nn=2*720;
cloud1=NaN*zeros(1,nn);
cloud2=NaN*zeros(1,nn);
cloud3=NaN*zeros(1,nn);
ceil_b=NaN*zeros(384,nn);
ceil_utc=NaN*zeros(1,nn);

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

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

fgetl(fid);

while ~feof(fid)
[data]=fscanf(fid,'-%4d-%2d-%2d %2d:%2d:%2d\n',[6 1]);
if ~isempty(data) 
ceil_year=data(1);
ceil_month=data(2);
ceil_day=data(3);
ceil_hour=data(4);
ceil_min=data(5);
ceil_sec=data(6);

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

if range_flag==3, loop=16 ;end %0 to 7.68 Km
if range_flag==4, loop=24 ;end %0 to 11.25 Km
if range_flag==8, loop=47 ;end %0 to 22.5 Km

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;
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:loop.*16,k)=ceil_data(1:loop.*16); % read in 1 `ray' of lidar backscatter
if length(data) > 0
   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:loop.*16,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:loop.*16]*30.0/1000.0;
cloud1=cloud1(1:k)./1000.0;
cloud2=cloud2(1:k)./1000.0;
cloud3=cloud3(1:k)./1000.0;



