 clear all;
 
 % calibration procedure with Balloon Sonde
 
 file_ball=uipickfiles('FilterSpec','*.txt','Prompt','Balloon Sonde File Selector');
 sonde=read_sonde(file_ball);
 nsondes=length(sonde);
 
 
 % select Raman Data from Big Lidar
 
 file_raman=uipickfiles('FilterSpec','*.txt','Prompt','HU Big Raman File Selector');
 niles=length(file_raman);
 ts=[sonde.time];
 
 
 % determine time stamps
 for j=1:niles
     [datadir,s,sext]=fileparts(char(file_raman(j)));
     yy=str2num(s(2:3));
     mm=str2num(s(4:4));
     dd=str2num(s(5:6));
     hh=str2num(s(7:15));
     timestamp(j)=datenum([yy,mm,dd,hh+4,0,0]);
 end
 
 [C,is]=min(abs([sonde.time]-timestamp(1)));
 kb=1.38064852e-23;
 calib_s=sonde(is);
 idx=find(calib_s.PRESS>0);
 calib_TEMP=calib_s.TEMP(idx)+276; % convert temperature from Celsius to Kelvin
 calib_PRESS=calib_s.PRESS(idx)*100; % convert pressure from millibars to Pa
 calib_Height=calib_s.GPH(idx);
 calib_mol_dens=calib_PRESS./(kb*calib_TEMP); % units are number of molecules/m^3
 
 clear sonde calib_s
 ch=zeros(10000,9);
 ch_filt=ch;

for k = 1:niles
              
    lidar=read_HUBIG(file_raman(k));
    duration=(lidar.time_e-lidar.time_s)*24*60*60; %in sec
    val=2.998e8/2/7.5/(duration*20);
    % convert photon counting vals to rates
    lidar.ch(2).data=lidar.ch(2).data*val;
    lidar.ch(4).data=lidar.ch(4).data*val;
    lidar.ch(6).data=lidar.ch(6).data*val;
    lidar.ch(8).data=lidar.ch(8).data*val;
    

    % pull signals for raman channels and average;
    pp=[lidar.ch.data];
    ch=ch+pp/niles;
end

M=40; %set range binning to 40*7.5=300m


    
% remove background noise from signals and filter using median filter

bkr=mean(ch(6666:7900,:),1);
for kk=1:9
    ch_filt(:,kk)=ch(:,kk)-bkr(kk);
    ch_filt(:,kk)=medfilt2(ch_filt(:,kk),[40 1]);
end  
 
    Tc=(cl6(1:4000)./tm3(1:4000).*exp(cumsum(((355/353.5)^4+1)*tm3(1:4000))*dz));
    Tc1=(cl5(1:4000)./tm3(1:4000).*exp(cumsum(((355/386)^4+1)*tm3(1:4000))*dz));
    Tc2=Tc1.^(2/(1+355/386));
    Tc3=Tc2;
    Tc3(1001:length(Tc3))=Tc3(1001);
   
 
   ntem=469./(log(abs(cl5)./abs(cl6))+2.35);
   mtem(kk,:)=ntem;
  
   gff=temw(1:600)./mtem(1,1:600);
   ntem1=ntem;
   ntem1(1:600)=ntem(1:600).*gff(1:600);
   cmtem(kk,:)=smooth(ntem1,30);
   
   time(kk)=sscanf(Files(k).name(5:6), '%f')+(sscanf(Files(k).name(7:8), '%f')+ sscanf(Files(k).name(10:11), '%f')/60)/24;
  nff=char(Files(k).name(2:13));
 name1(kk,:)=nff;
 kk=kk+1;





 
mt1=cmtem;
for k = 3:kk-3
    mt1(k,:)=mean(cmtem(k-2:k+2,:));
end



alt_new = 0.0075:0.0075:length(mt1)*0.0075;
figure
imagesc(time,alt_new,abs(mt1'));
set(gca,'YDir','normal');
caxis([200 310]);
colormap([0 0 0;jet(100)]);
title(' The relative ozone density ','FontSize',12);
xlabel(' UTC Time (July.31 2017)','FontSize',12);
ylabel('Altitude (km)','FontSize',12);
colorbar;   
