function [cls] = rd_cls(filename)
%
%       function [cls] = rd_cls(filename)
%
%   ....This Matlab function reads cloud boundary data files generated from NASA's Cloud
%       Lidar System (CLS).  More information can be obtained on this instrument from
%       Jim Spinhirne's web page at http://virl.gsfc.nasa.gov/er2cls.html.
%
%       Input:
%         filename - filename of CLS boundary file
%
%       Output:
%         cls      - structure containing the following variables:
%                      mtime             - Matlab serial date/time
%                      seconds           - time in seconds since midnight
%                      lat               - aircraft latitude
%                      lon               - aircraft longitude
%                      alt               - aircraft altitude
%                      roll              - aircraft roll
%                      cloud_boundary    - array containing top and bottom cloud boundaries
%                                            of up to 5 layers below aircraft (NaN represent
%                                            no data.)
%                                              cls.cloud_boundary(:,1) = top of 1st cloud layer
%                                              cls.cloud_boundary(:,2) = bottom of 1st cloud layer
%                                              cls.cloud_boundary(:,3) = top of 2nd cloud layer
%                                                ...
%                      surface_elevation - surface elevation as deduced by CLS
%                      hdr               - header lines of CLS file.
%
%       'mtime' is in Matlab serial date format.  
%       To view these times, use datestr(mtime,0).
%
%       To plot using 'mtime', try 'plot(mtime,lat); datetick('x',13);'
%
%       Written by Von P. Walden
%                  28 May 1998
%

%
%   ....Opens the CLS file.
%
fp = fopen(filename);

%
%   ....Determines the number of header lines.
%
num_hdr_lines = fscanf(fp,'%d',1);
frewind(fp);

%
%   ....Reads the beginning and ending dates.
%
for k = 1:6,   % Skips 6 lines.
   fgetl(fp);
end
arr = fscanf(fp,'%d',6);
bYY = arr(1);
bMM = arr(2);
bDD = arr(3);
eYY = arr(4);
eMM = arr(5);
eDD = arr(6);
frewind(fp);

%
%   ....Reads the header lines.
%
cls.hdr = [];
for i = 1:num_hdr_lines,
   cls.hdr = [cls.hdr fgets(fp)];
end

%
%   ....Reads in all the data as integers.
%
arr = fscanf(fp,'%d');

%
%   ....Replaces missing values with Not-A-Numbers.
%
ind = find(arr == 99999);
arr(ind) = NaN;

%
%   ....Places data into specific arrays.
%
N = length(arr);
cls.seconds = arr(1:16:N);
% Converts to Matlab time.
cls.mtime = datenum(bYY,bMM,bDD) + cls.seconds/(60*60*24);
cls.lat   = arr(2:16:N)/100; % degrees
cls.lon   = arr(3:16:N)/100; % degrees
cls.alt   = arr(4:16:N);     % meters
cls.roll  = arr(5:16:N);

cls.cloud_boundary(:,1) = arr(6:16:N);
cls.cloud_boundary(:,2) = arr(7:16:N);
cls.cloud_boundary(:,3) = arr(8:16:N);
cls.cloud_boundary(:,4) = arr(9:16:N);
cls.cloud_boundary(:,5) = arr(10:16:N);
cls.cloud_boundary(:,6) = arr(11:16:N);
cls.cloud_boundary(:,7) = arr(12:16:N);
cls.cloud_boundary(:,8) = arr(13:16:N);
cls.cloud_boundary(:,9) = arr(14:16:N);
cls.cloud_boundary(:,10) = arr(15:16:N);

cls.surface_elevation = arr(16:16:N);

fclose(fp);