%reads in ceilometer files

% program written by C. Fairall.  Modified by D. E. Lane
% on 5 March 2001.  Added 10-min statistics.

disp('Reading ceilo files')
clear all;
fclose('all');

%enter path of the processed data:
processed_data_path='C:\Data\ICEALOT\data\ceilometer\';
%enter path where to write the datafile with all of the pertinent info in it:
final_data_path='C:\Data\ICEALOT\data\ceilometer\';

cruise_name='ICEALOT';
year='2008';
tit = cruise_name;
lc = length(cruise_name) + 1;
tit(lc:lc+) = year;
s=[processed_data_path '*08.txt'];
dr=dir(s);			%lists all files in directory
[n m]=size(dr);
nx=n;%				%# of files in data directory (includes . and ..)     


%******************  setup read data loop  **********
ss=processed_data_path;

dd=1;
jamx=1;
for ibg = 1:n;%nx 			%major read loop
    dr(ibg).name
   e=[ss dr(ibg).name]; %gets names of files in directory
   stux=load(e);   		%loads all data from files into stux
   
%******   decode date  ************
   if isempty(stux)		% if no data, end the program
      break;
   end;
    [nl ml]=size(stux);  		%nl=file length - m1 is number of cols
    ff=dr(ibg).name;		 		%names of files
    jdx=str2num(ff(1:3));
    month=str2num(ff(5:6));	%get month
    day=str2num(ff(8:9));		%get day
    year=str2num(ff(11:12));		%get year
    jd0=datenum(2000+year,month,day)-datenum(2000+year-1,12,31); % calculates julian date
    jd=jd0+((stux(:,5)+stux(:,6)/60+stux(:,7)/3600)/24); % calculate decimal jdate
   
    happ(jamx:jamx+nl-1,1:5)=[jd stux(:,1:4)];  % make data array
    jamx=jamx+nl  
end; %  data line loop
happ(:,3)=happ(:,3);%convert ft to meters
% plot 15-sec data
ik=find(happ(:,3)==0);happ(ik,2)=5;
ik=find(happ(:,2)==4 & happ(:,3)<0);happ(ik,2)=5;
ill=find(happ(:,2)<5 & happ(:,2)>0);
figure(1);
plot(happ(:,1),happ(:,2),'o'); xlabel('Year Day'); ylabel('Number of Cloud Layers');
title(cruise_name, year)
print_buffer = [ss , cruise_name, '_nolayers.png'];
print('-dpng ', print_buffer);

figure(2);
plot(happ(ill,1),happ(ill,3)/1000,'.');xlabel('Year Day (2007)');ylabel('Cloud Base (km)');
title('GOMECC 2007')
print_buffer = [ss , cruise_name, '_cball.png'];
print('-dpng ', print_buffer);

% save 15-sec data

ceilo_30sec=happ;					% this is the raw data with date
% save d:\data\ceilo\neaqs04_ceilo_30s_a.txt ceilo_30sec -ascii -tabs
save([final_data_path, cruise_name 'ceilo_30s_a.txt'], 'ceilo_30sec', '-ascii', '-tabs')

%******************  setup read data loop  ********** 

ny=jamx-1;
lam=0:24;  			  % hours
jdx=lam/24;			  % fraction of day
jdd=happ(:,1);		  % date vector
jd0=floor(jdd(1));  % first Year Day of data series
jxm=floor(jdd(ny)); % last Year Day of data series
jx=1;					  % counts number of hours of data
kx=1;					  % counts number of 10-min blocks of data
min=0:6;				  % 10-min parts of an hour
jdmx=min/6/24;			  % 10-min fraction of hour

while jd0<=jxm+1  % the first day is l.t.e. last day
   ii=find(jdd>jd0 & jdd<(jd0+1));  % find all data
   if isempty(ii)  % if no data
      jd0=jd0+1;   % increment the start date
   else
   	jj=1;			 % if not empty start first hour calc.
      while jj<25  % not the end of the day
         pp=jdd(ii)-floor(jdd(ii)); % fractional part of day
      	ij=find(pp>jdx(jj) & pp<jdx(jj+1)); % select hour
      	if isempty(ij)
            frac(jx)=NaN;
            hgt(jx,1:2)=NaN;
      	else
            kk=1;			 % if hourly data exists start 10-min calc.
            while kk<7
               kj=find(pp(ij)>(pp(ij(1))+jdmx(kk)) & pp(ij)<(pp(ij(1))+jdmx(kk+1)));  % select 10-min section
               if isempty(kj)
                  fracm(kx)=NaN;
                  hgtm(kx,1:2)=NaN;
                  fracmin(kx)=NaN;
                  fracmino(kx)=NaN;
                  mcloudy(kx)=NaN;
                  hgtmin(kx)=NaN;
                  hgtsm1(kx)=NaN;
                  hgtsm2(kx)=NaN;
                  mindata(kx)=0;
                  clearm(kx)=NaN;
                  onecldm(kx)=NaN;
                  mltcldm(kx)=NaN;
                  obscurem(kx)=NaN;
                  pobsm(kx)=NaN;
                  fracclr(kx)=NaN;
                  
               else   
                  [mindata(kx) c]=size(pp(kj));  % how much data in 10-min section
                  mincloud=find(happ(ii(ij(kj)),2)>0 & happ(ii(ij(kj)),2)<4);  % if cloud
                  minobscure=find(happ(ii(ij(kj)),2)==4);  % if vertical vis
                  if isempty(mincloud);
                     mcloud(kx)=0;
                  else
                     [mcloud(kx),dx]=size(mincloud);
                  end;
                  if isempty(minobscure)
                     obscurem(kx)=0;
                  else
                     [obscurem(kx),dy]=size(minobscure);
                  end;
                  mcloudy(kx)=mcloud(kx)+obscurem(kx);	% total possible cloud
                  fracmin(kx)=mcloud(kx)/mindata(kx);		% cloud fraction w/o obscure
                  fracmino(kx)=mcloudy(kx)/mindata(kx);	% cloud fraction 2/ obscure
                  cleardata=find(happ(ii(ij(kj)),2)==0);	% find clear samples
                  [clearm(kx),dz]=size(cleardata);			% number of clear samples
                  fracclr(kx)=clearm(kx)/mindata(kx);		% clear fraction
                  monelay=find(happ(ii(ij(kj)),2)==1);	% find times with one layer
                  [onecldm(kx),ddx]=size(monelay);			% number times one layer
                  mmltcld=find(happ(ii(ij(kj)),2)==2 | happ(ii(ij(kj)),2)==3);
                  [mltcldm(kx),ddy]=size(mmltcld);			% number time multiple layers
                  ispobsm=find(happ(ii(ij(kj)),2)==5);
                  [pobsm(kx),ddz]=size(ispobsm);			% number of partially obscured times
                  if mcloudy(kx)==0
                     hgtmin(kx)=NaN;
                     hgtsm1(kx)=NaN;
                     hgtsm2(kx)=NaN;
                  else
                     if mcloudy(kx)>0
                        allmincover=[mincloud' minobscure'];  	%all possible cloud samples  
                 			hgtmin(kx)=median(happ(ii(ij(kj(allmincover))),3));  % take median hgt
                        chgtsm=happ(ii(ij(kj(allmincover))),3);	%get hgts all poss. clouds
                        yy=sort(chgtsm);									%sort
                        [nny mmy]=size(yy);
                        jj1=max(floor(nny*.16),1);		%calculates 1-sigma for height
                        jj2=max(floor(nny*.84),1);
                        hgtsm1(kx)=yy(jj1);
                        hgtsm2(kx)=yy(jj2);
                     end;                       
            		end;  % end if not cloudy
               end; 
               jdmin(kx)=jd0+jdx(jj)+jdmx(kk); % integrate timestep
               kk=kk+1;
               kx=kx+1;
            end  %while kk      
            [alldata(jx) b]=size(pp(ij));		% all data in selected hour            
            iscloud=find(happ(ii(ij),2)>0 & happ(ii(ij),2)<4);%if cloud 
            isobscure=find(happ(ii(ij),2)==4);%if vertical vis
            if isempty(iscloud)
               cloud(jx)=0;
            else
               [cloud(jx),dx]=size(iscloud);
            end;
            if isempty(isobscure)
              	obscure(jx)=0;
           	else
              	[obscure(jx),dy]=size(isobscure);
              end;
            cloudy(jx)=cloud(jx)+obscure(jx);%add cloud and 'totally obscured' cases
            ceilofrac(jx)=cloud(jx)/alldata(jx);  % cloud fraction cloud only
            ceilofraco(jx)=cloudy(jx)/alldata(jx);% cloud fraction cloud+obscure
            isclr=find(happ(ii(ij),2)==0); % if clear
            [clr(jx),c]=size(isclr);       % number of clear points
            clrfrac(jx)=clr(jx)/alldata(jx); % fraction of time clear
            isonecld=find(happ(ii(ij),2)==1); % one cloud layer
            [onecld(jx) d]=size(isonecld);    % number of point one layer
            ismltcld=find(happ(ii(ij),2)==2 | happ(ii(ij),2)==3);
            [mltcld(jx) e]=size(ismltcld);    % number of points multiple layer
            ispobs=find(happ(ii(ij),2)==5);   % is part obscured
            [pobs(jx) f]=size(ispobs);        % amt part obscured
         	if cloudy(jx)==0
               cloudhgt(jx)=NaN;
               hgts1(jx)=NaN;
               hgts2(jx)=NaN;
            else
               if cloudy(jx)>0
                  allcover=[iscloud' isobscure'];	% all possible cloud
                  cloudhgt(jx)=median(happ(ii(ij(allcover)),3));	%median hgt 
                  chgts=happ(ii(ij(allcover)),3);	% all hgts poss. cloud
      				y=sort(chgts);
                  [ny my]=size(y);
                  j1=max(floor(ny*.16),1);	% calculate +/- 1-sigma hgts
                  j2=max(floor(ny*.84),1);
                  hgts1(jx)=y(j1);
                  hgts2(jx)=y(j2);
               end;
            end;  % end if not cloudy
         end;%end if isempty ij
         jdh(jx)=jd0+jdx(jj);
     	   jj=jj+1;
      	jx=jx+1;
      end;%end while jj
      jd0=jd0+1
   end;%end if isempty ii
end;%end while jx

% plot rest of data
ii=find(cloudhgt==0);cloudhgt(ii)=NaN;ceilofrac(ii)=NaN;hgts1(ii)=NaN;hgts2(ii)=NaN;
jdh=jdh(1:length(ceilofrac));

figure(3);
plot(jdh,ceilofrac,'o');   %hold on;plot(jdmin,fracmin,'r.')
xlabel(['Year Day (' year ')']);ylabel('Ceilo Cloud Fraction');
title('GOMECC 2007')
print_buffer = [ss , cruise_name, '_fraction.png'];
print('-dpng ', print_buffer);

figure(4);plot(jdh,cloudhgt,'o',jdmin,hgtmin,'.');
xlabel(['Year Day (' year ')']);ylabel('Ceilo Cloudbase Hts');legend('Hourly','10-min');
title('GOMECC 2007')
print_buffer = [ss , cruise_name, '_base.png'];
print('-dpng ', print_buffer);

figure(5);
subplot(2,1,1);
plot(jdh,cloudhgt/1000,'.',jdh,hgts1/1000,'o',jdh,hgts2/1000,'x');
hold on;
errorbar_dana(jdh,cloudhgt/1000,hgts1/1000,hgts2/1000,'.');
ylabel('Hourly Hts (km:+/- sigma)');
title('GOMECC 2007 Ceilo Cloud Base Heights');
xlim([jdh(1) jdh(length(jdh))])
subplot(2,1,2);
plot(jdmin,hgtmin/1000,'.',jdmin,hgtsm1/1000,'o',jdmin,hgtsm2/1000,'x');
xlim([jdmin(1) jdmin(length(jdmin))])
hold on;
errorbar_dana(jdmin,hgtmin/1000,hgtsm1/1000,hgtsm2/1000,'.');ylabel('10-min Hts (km:+/- sigma)');
xlabel(['Year Day (2007)']);
print_buffer = [ss , cruise_name, '_base_1h_10min.png'];
print('-dpng ', print_buffer);

% save the data in hourly and 10-min files
ceilo_hourly=[jdh' alldata' clr' onecld' mltcld' obscure' pobs' clrfrac' ceilofrac' ceilofraco'  ...  
              cloudhgt'  hgts1' hgts2'];
save([final_data_path, cruise_name 'ceilo_h_a.txt'], 'ceilo_hourly', '-ascii', '-tabs')

ceilo_10min=[jdmin' mindata' clearm' onecldm' mltcldm' obscurem' pobsm' fracclr' fracmin' ...
             fracmino' hgtmin' hgtsm1' hgtsm2'];
save([final_data_path, cruise_name 'ceilo_10min_a.txt'], 'ceilo_10min', '-ascii', '-tabs')