% Reads in raw ceilometer .DAT files as written by the Vaisala logging program
%Writes out the old style cloud base height file and the backscatter file
% This version is looking for data on  E:\
%It ignores files smaller than 8k to avoid the log files
%Uses LCload_1 adapted from Pavlos Kollias on Stratus 03 Revelle cruise
%Uses LCload_CL31_2007.m adapted from Ludovic Bariteau for the new
%ceilometer CL31.

% program written by C. Fairall.
% on 8 Aug 2004.

dv=datevec(dn);
year = num2str(dv(1));
filename=sprintf('%sa%c%02d%02d*.dat', way_raw_data_ceilo,year(4), dv(2), dv(3) );
fname=sprintf('a%c%02d%02d*.dat',year(4), dv(2), dv(3) );

dr=dir(filename);			%lists all files in directory
[nfiles m]=size(dr);

if nfiles == 0
    return
end

%******************  setup read data loop  **********
wrtit=1;
nbits=20; % number of bits for reflectivity
nn=5760; %largest possible number of samples --should be exactly 5760 samples, but who knows?
start_height = 0; % in m % SPdeS ?should be height of ceilometer on deck?
gates_number=770;  %number of gates, to be updated from file
resolution=10; %10m resolution, to be updated from file

tempcloud1=NaN+zeros(nn,1);
tempcloud2=NaN+zeros(nn,1);
tempcloud3=NaN+zeros(nn,1);
tempcloudcode=NaN+zeros(nn,1);

tempceil_sec=NaN+zeros(nn,1);
tempceil_min=NaN+zeros(nn,1);
tempceil_hour=NaN+zeros(nn,1);
tempceil_b=NaN+zeros(gates_number,nn);

for ibg = 1:nfiles 			%major read loop (no appending!)
    if dr(ibg).bytes>8000
        %ff=dr(ibg).name;		 		%names of files
        %ceil_filename=[way_raw_data_ceilo ff];
        fprintf(1,'Reading ceilo file %s\n', dr(ibg).name);
        ceil_filename=sprintf('%s%s', way_raw_data_ceilo,dr(ibg).name);
        mon=dr(ibg).name(3:4);	%get month
        %month=str2double(mon);	%get month
        dy=dr(ibg).name(5:6);		%get day
        %day=str2double(dy);		%get day
        yr=year(3:4);
        %fyear=str2double(year);       
        %date=[dv(1) dv(2) dv(3)];%get year %SPdeS changed data to date
        jd=datenum(dv(1),dv(2),dv(3))-datenum(dv(1)-1,12,31);
        jdn=sprintf('%03d',jd);
        
        % and now a very slow read of a horrendously formatted Vaisala
        % ceilometer file...
        LCload_CL31_BAO_2010a; % SPdeS re-implement as function
    end;     % if file is big enough
end;

% truncate data  x=tempx(1:k)
ceil_b=tempceil_b(:,1:k);
% wrap 2s complement for nbits=20
%ceil_b2=ceil_b;
index=find(ceil_b>=2^(nbits-1));
if ~isempty(index)
    ceil_b(index)=ceil_b(index)-2^nbits;
end
% % for comparison puposes, compute backscatter with wrong 16-bit 2s complement
% ceil_bwrong=ceil_b;
% ceil_bwrong(ceil_b>=2^15)=ceil_b(ceil_b>=2^15)-2^16;
% ceil_bwrong=ceil_bwrong*1e-5;

%if scale==1
%    ceil_b=ceil_b2*1e-5; % scaled by default in units 1e5/(km srad)
%else
    ceil_b=ceil_b*scale*1e-5;
%end
% ceil_b is now VOLUME BACKSCATTER COEFFICIENT [units 1/(km*steradians)].
% It is the extinction-normalized power backward-scattered per
% atmospheric scattering volume drdOmega [m*srad]. This is
% related to the power transmitted and recieved at the instrument
% by the lidar equation (CL31 manual p. 75).
% In practice it is sometimes negative, presumably because of noise.

% convert units to meters if necessary
%if ~units_bit
    cloud1=tempcloud1(1:k)*mm;    %feet to meters
    cloud2=tempcloud2(1:k)*mm;
    cloud3=tempcloud3(1:k)*mm;
%else
%    cloud1=tempcloud1(1:k);
%    cloud2=tempcloud2(1:k);
%    cloud3=tempcloud3(1:k);
%end
cloudcode=tempcloudcode(1:k);
ceil_sc=tempceil_sec(1:k);
ceil_mn=tempceil_min(1:k);
ceil_hr=tempceil_hour(1:k);
ceil_utc=tempceil_hour(1:k)+tempceil_min(1:k)/60.0+tempceil_sec(1:k)/3600.0;
ceil_date(1:k,1)=dv(1)*10000+dv(2)*100+dv(3);

ceil_range=(start_height+(1:gates_number)*resolution)/1000.0; % km

% % plot
% %subplot(2,1,1)
% imagesc(ceil_utc,ceil_range,ceil_b); axis xy; colorbar;
% xlabel('hour UTC'); ylabel('range (km)'); title({num2str([date(1) date(2) date(3)]),'CL31 log_{10}backscatter (m^{-1} srad^{-1})'});
% % plot the wrong 2s complement ceil_b:
% % subplot(2,1,2)
% % imagesc(ceil_utc,ceil_range,log10(max(0,ceil_bwrong))); axis xy; colorbar;
% % xlabel('hour UTC'); ylabel('range (km)'); title('CL31 log_{10}backscatter (m^{-1} srad^{-1})');

% ceilo_30_sec is data in the form required by the calling program
ceilo_30_sec=[cloudcode cloud1 cloud2 cloud3 ceil_hr ceil_mn ceil_sc];% this is the raw 1-D time series data

if wrtit
    %fn=[mon '_' dy '_' yr];
    %fn=sprintf('%02d_%02d_%02d', mon, dy, yr);
    %sta=[way_proc_data_ceilo jdn '_' fn '.txt'];
    sta=sprintf('%scloudbase_%s_%s_%s_%s.txt',way_proc_data_ceilo,jdn, mon, dy, yr);
    save( sta, 'ceilo_30_sec', '-ascii','-tabs')
    
    %stb=[way_proc_data_ceilo 'backscatter_' jdn '_' fn '.mat'];
    stb=sprintf('%sbackscatter_%s_%s_%s_%s.mat',way_proc_data_ceilo,jdn,mon,dy,yr);
    %save(stb,'ceil_utc','ceil_b','-ascii','-tabs');
    save(stb,'ceil_utc','ceil_range','ceil_b');  % height --> ceil_range SPdeS
end; % write
