addpath C:\NEAQS_2004\matlab_code\
pathdata = 'c:\NEAQS_2004\ETLceilometer\';
outpath  = 'c:\NEAQS_2004\ETLceilometer\results\';
eval(['cd ', pathdata]);

aa = dir('b4*00.dat');

k = length(aa);
[num] = num2str(k);
disp('Number of ceilometer daily *.raw files : '); disp( [num]); 



%%%%%%%%%%%%%%%%%% last file == k %%%%%%%%%%%%%%%

for w = 1:18;

[numf] = num2str(w);

disp('# of *.day file in process :'); disp ( [numf]); 

%%%%%%% k is the # of *.cdf files need to be processed. %%%%%%%%

ceil_filename = [aa(w).name]


%%%%%%%%%%%%%%%%%% INPUT DATA FILE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ceil_year     =  2004;
ceil_month    =  str2num(ceil_filename(3:4));
ceil_day      =  str2num(ceil_filename(5:6));
%%%%% READ HEX CEILOMETER FILE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
LCload_1

%%%% SAVE DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


jd0=(datenum(2000+ceil_year,ceil_month,ceil_day)-datenum(2000+ceil_year-1,12,31));
jdd= num2str(jd0)

ceilo_30_sec=[cloudcode' cloud1' cloud2' cloud3' ceil_hr' ceil_mn' ceil_sc'];% this is the raw data with time
year = '04';
if ceil_month < 10
    month = ['0' num2str(ceil_month)];
else
    month = num2str(ceil_month);
end
if ceil_day < 10
    day = ['0' num2str(ceil_day)];
else
    day = num2str(ceil_day);
end
out_file = [ outpath jdd  '_' month '_' day '_' year '.txt ceilo_30_sec -ascii -tabs'];
eval(['save ',out_file]);
out_file2 = [outpath 'backscatter_' jdd  '_' month '_' day '_' year ' ceil_b ceil_hr ceil_mn ceil_sc height'];
eval(['save ',out_file2]);
clear cloudcode cloud1 cloud2 cloud3 ceil_hr ceil_mn ceil_sc ceilo_30_sec ceil*;

%%%%% 2 LEVEL OF DATA ANALYSIS -- STATISTICS %%%%%%%%%


end

clear all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% PLOT CEILOMETER DAILY AND COMPOSITE TIME SERIES AND BACKSCATTER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% INTENSITY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pathdata = 'C:\NEAQS_2004\ETLceilometer\results\';

eval(['cd ', pathdata]);

a = dir('*.txt');

k = length(a);


[num] = num2str(k);
disp('Number of ceilometer *.txt files : '); disp( [num]); 


for w = 15:k;

[numf] = num2str(w);

disp('# of *.txt file in process :'); disp ( [numf]); 

%%%%%%% k is the # of *.cdf files need to be processed. %%%%%%%%

data = [a(w).name]
scatter_data = ['backscatter_' a(w).name(1:12)];
ceil_series = load(data);
load(scatter_data);
ceil_time = ceil_series(:,5)+ceil_series(:,6)/60+ceil_series(:,7)/3600;

figure(5)
set(gcf,'DefaultAxesFonts',12);
plot(ceil_time,ceil_series(:,2),'b.');
xlim([0 24]); set(gca,'xtick',[0:4:24]);
%ylim([0 4000]);
xlabel('Time (UTC)','FontName','times','FontSize',16);
ylabel('Height (m)','FontName','times','FontSize',12);
title(['NEAQS 2004 ' data(5:6) '/' data(8:9) ', Ceilometer Cloud Base'],'FontName','times','FontSize',12);
print_buffer = ['C:\NEAQS_2004\ETLceilometer\results\NEAQS2004_ceilometer_' data(5:12) '_base.jpg'];
print('-djpeg90 ', print_buffer);

pause(5);
figure(6)
set(gcf,'DefaultAxesFonts',12);
imagesc(ceil_time,height*.303,log10(abs(ceil_b)));
axis xy; colorbar;
xlim([0 24]); set(gca,'xtick',[0:4:24]);
%ylim([0 4000]);
xlabel('Time (UTC)','FontName','times','FontSize',12);
ylabel('Height (m)','FontName','times','FontSize',12);
title(['NEAQS 2004 ' data(5:6) '/' data(8:9) ', Ceilometer Backscatter'],'FontName','times','FontSize',12);
print_buffer = ['C:\NEAQS_2004\ETLceilometer\results\NEAQS2004_ceilometer_' data(5:12) '_backscatter.jpg'];
print('-djpeg90 ', print_buffer);
%%% START WITH THE RADIOMETER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clear cei_b ceil* 
pause(5);
end


clear all


%%%%%% COMPOSITE OF CEILOMETER BACKSCATTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



pathdata = 'C:\NEAQS_2004\ETLceilometer\results\';

eval(['cd ', pathdata]);

a = dir('*.txt');

k = length(a);


[num] = num2str(k);
disp('Number of ceilometer *.txt files : '); disp( [num]); 

st=1%%%%%%%%%%%%%%%%%%%%%%%by changing st you can plot they dates separetely

for w = st:k;

[numf] = num2str(w);

disp('# of *.txt file in process :'); disp ( [numf]); 

%%%%%%% k is the # of *.cdf files need to be processed. %%%%%%%%

data = [a(w).name]
scatter_data = ['backscatter_' a(w).name(1:12)];
load(scatter_data);
ceil_series = load(data);

ceil_time = str2num(data(1:3))+(ceil_hr+ceil_mn/60+ceil_sc/3600)/24;

if w == st 

sf = length(ceil_time);
time_all = ceil_time;
height_all = height;
backscatter = ceil_b;
base        = ceil_series(:,2);
else
    
    gf = length(time_all);
    sf = length(ceil_time);
    time_all(gf+1:gf+sf)=ceil_time;
    backscatter(1:256,gf+1:gf+sf)=ceil_b;
    base(gf+1:gf+sf)=ceil_series(:,2);
end


end

figure(6)
set(gcf,'DefaultAxesFonts',12);
imagesc(time_all(1:3:length(time_all)),height(1:100)*.303,log10(abs(backscatter(1:100,1:3:length(time_all)))),[1 3.5]);
axis xy; colorbar;
%ylim([0 3000]);
xlabel('Day','FontName','times','FontSize',12);
ylabel('Height (m)','FontName','times','FontSize',12);
title('NEAQS 2004 Ceilometer Backscatter','FontName','times','FontSize',12);
print_buffer = ['C:\NEAQS_2004\ETLceilometer\results\NEAQS2004_ceilometer_alldays_backscatter_b.jpg'];
print('-djpeg90 ', print_buffer);

pause(5);
figure(5)
set(gcf,'DefaultAxesFonts',12);
plot(time_all(1:3:length(time_all)),base(1:3:length(time_all)),'b.');
xlim([min(time_all) max(time_all)]);
%ylim([0 4000]);
xlabel('Day','FontName','times','FontSize',16);
ylabel('Height (m)','FontName','times','FontSize',12);
title('NEAQS 2004, Ceilometer Cloud Base','FontName','times','FontSize',12);
print_buffer = ['C:\NEAQS_2004\ETLceilometer\results\NEAQS2004_ceilometer_alldays_base_b.jpg'];
print('-djpeg90 ', print_buffer);
pause(5);
