close all
clear all


way_rawdata='C:\Data';


fluxes_jd=[];

for ddd=216;
    jd=num2str(ddd);
    if  ddd<100;
        jd=['0' num2str(ddd)];	
        end;
    for hhh=0:13;
        hr=num2str(hhh);
        if hhh<10;
           hr=['0' num2str(hhh)];	
        end;
        
        
        Licor_imp=importdata([way_rawdata '\GOMECC\June08\O3WLic_q\Licor_H2O',jd,hr,'.txt']);
        r=Licor_imp./1000;
        
        Wcat_imp=importdata([way_rawdata '\GOMECC\June08\O3WLic_q\Wcat_',jd,hr,'.txt']);
        Wcat=Wcat_imp(3,:);
        
        O3ppb_imp=importdata([way_rawdata '\GOMECC\June08\O3WLic_q\O3ppb_real_interp',jd,hr,'.txt']);
        O3ppb_real=O3ppb_imp;
   
        disp([jd,hr]);
%fluxes uncorrected for quenching

[c1, lags1]=xcov(Wcat(1:6000),O3ppb_real(1:6000),100,'unbiased');
    p = polyfit(lags1,c1,26);
    covEpol = polyval(p,lags1);
    ll=find(-100<lags1 & lags1<-35);
    jj=find(covEpol(ll)==min(covEpol(ll)));
    lagtime_1=-lags1(jj);  %lagtime in points
    flux_1_o=min(covEpol(jj));

[c2, lags2]=xcov(Wcat(6001:12000),O3ppb_real(6001:12000),100,'unbiased');
    p = polyfit(lags2,c2,26);
    covEpol = polyval(p,lags2);
    ll=find(-100<lags2 & lags2<-35);
    jj=find(covEpol(ll)==min(covEpol(ll)));
    lagtime_2=-lags2(jj);  %lagtime in points
    flux_2_o=min(covEpol(jj));

[c3, lags3]=xcov(Wcat(12001:18000),O3ppb_real(12001:18000),100,'unbiased');
    p = polyfit(lags3,c3,26);
    covEpol = polyval(p,lags3);
    ll=find(-100<lags3 & lags3<-35);
    jj=find(covEpol(ll)==min(covEpol(ll)));
    lagtime_3=-lags3(jj);  %lagtime in points
    flux_3_o=min(covEpol(jj));

[c4, lags4]=xcov(Wcat(18001:24000),O3ppb_real(18001:24000),100,'unbiased');
    p = polyfit(lags4,c4,26);
    covEpol = polyval(p,lags4);
    ll=find(-100<lags4 & lags4<-35);
    jj=find(covEpol(ll)==min(covEpol(ll)));
    lagtime_4=-lags4(jj);  %lagtime in points
    flux_4_o=min(covEpol(jj));
    
[c5, lags5]=xcov(Wcat(24001:30000),O3ppb_real(24001:30000),100,'unbiased');
    p = polyfit(lags5,c5,26);
    covEpol = polyval(p,lags5);
    ll=find(-100<lags5 & lags5<-35);
    jj=find(covEpol(ll)==min(covEpol(ll)));
    lagtime_5=-lags5(jj);  %lagtime in points
    flux_5_o=min(covEpol(jj));

[c6, lags6]=xcov(Wcat(30001:36000),O3ppb_real(30001:36000),100,'unbiased');
    p = polyfit(lags6,c6,26);
    covEpol = polyval(p,lags6);
    ll=find(-100<lags6 & lags6<-35);
    jj=find(covEpol(ll)==min(covEpol(ll)));
    lagtime_6=-lags6(jj);  %lagtime in points
    flux_6_o=min(covEpol(jj));

O3mean_1_o=mean(O3ppb_real(1:6000));
O3mean_2_o=mean(O3ppb_real(6001:12000));
O3mean_3_o=mean(O3ppb_real(12001:18000));
O3mean_4_o=mean(O3ppb_real(18001:24000));
O3mean_5_o=mean(O3ppb_real(24001:30000));
O3mean_6_o=mean(O3ppb_real(30001:36000));

O3stdev_1_o=std(O3ppb_real(1:6000));
O3stdev_2_o=std(O3ppb_real(6001:12000));
O3stdev_3_o=std(O3ppb_real(12001:18000));
O3stdev_4_o=std(O3ppb_real(18001:24000));
O3stdev_5_o=std(O3ppb_real(24001:30000));
O3stdev_6_o=std(O3ppb_real(30001:36000));



%Ozone new method
O3ppb_real=O3ppb_imp(52:35991); 

O3mean_1=mean(O3ppb_real(1:5990));
O3mean_2=mean(O3ppb_real(5991:11980));
O3mean_3=mean(O3ppb_real(11981:17970));
O3mean_4=mean(O3ppb_real(17971:23960));
O3mean_5=mean(O3ppb_real(23961:29950));
O3mean_6=mean(O3ppb_real(29951:35940));

O3stdev_1=std(O3ppb_real(1:5990));
O3stdev_2=std(O3ppb_real(5991:11980));
O3stdev_3=std(O3ppb_real(11981:17970));
O3stdev_4=std(O3ppb_real(17971:23960));
O3stdev_5=std(O3ppb_real(23961:29950));
O3stdev_6=std(O3ppb_real(29951:35940));

flux_1=cov(Wcat(1:5990),O3ppb_real(1:5990));
flux_1=flux_1(1,2);

flux_2=cov(Wcat(5991:11980),O3ppb_real(5991:11980));
flux_2=flux_2(1,2);

flux_3=cov(Wcat(11981:17970),O3ppb_real(11981:17970));
flux_3=flux_3(1,2);

flux_4=cov(Wcat(17971:23960),O3ppb_real(17971:23960));
flux_4=flux_4(1,2);

flux_5=cov(Wcat(23961:29950),O3ppb_real(23961:29950));
flux_5=flux_5(1,2);

flux_6=cov(Wcat(29951:35940),O3ppb_real(29951:35940));
flux_6=flux_6(1,2);


%            [Sw,F]=psd2(detrend(O3ppb_real),length(O3ppb_real),10.);
%            [Sws,Fsw]=specsmoo(Sw,10.);
%            
%            figure; loglog(Fsw(length(Fsw)-31:length(Fsw)),((Fsw(length(Fsw)-31:length(Fsw))).^(-2/3))./(Fsw(length(Fsw)-18).^(-2/3)).*Sws(length(Fsw)-18).*Fsw(length(Fsw)-18),'r');
%            hold on;loglog(Fsw,Sws.*Fsw,'g');
%            
%            title(['O3 Spectra:   Day ',jd,hr]);
%            xlabel('Frequency in Hz');
%            ylabel('Spectrum * Frequency');

%Fluxes corrected for quenching
O3ppb_real=O3ppb_imp;
a=5;
O3ppb_real=O3ppb_real.*(1+a*r');


[c1, lags1]=xcov(Wcat(1:6000),O3ppb_real(1:6000),100,'unbiased');
    p = polyfit(lags1,c1,26);
    covEpol = polyval(p,lags1);
    ll=find(-100<lags1 & lags1<-35);
    jj=find(covEpol(ll)==min(covEpol(ll)));
    lagtime_1_q=-lags1(jj);  %lagtime in points
    flux_1_o_q=min(covEpol(jj));

[c2, lags2]=xcov(Wcat(6001:12000),O3ppb_real(6001:12000),100,'unbiased');
    p = polyfit(lags2,c2,26);
    covEpol = polyval(p,lags2);
    ll=find(-100<lags2 & lags2<-35);
    jj=find(covEpol(ll)==min(covEpol(ll)));
    lagtime_2_q=-lags2(jj);  %lagtime in points
    flux_2_o_q=min(covEpol(jj));

[c3, lags3]=xcov(Wcat(12001:18000),O3ppb_real(12001:18000),100,'unbiased');
    p = polyfit(lags3,c3,26);
    covEpol = polyval(p,lags3);
    ll=find(-100<lags3 & lags3<-35);
    jj=find(covEpol(ll)==min(covEpol(ll)));
    lagtime_3_q=-lags3(jj);  %lagtime in points
    flux_3_o_q=min(covEpol(jj));

[c4, lags4]=xcov(Wcat(18001:24000),O3ppb_real(18001:24000),100,'unbiased');
    p = polyfit(lags4,c4,26);
    covEpol = polyval(p,lags4);
    ll=find(-100<lags4 & lags4<-35);
    jj=find(covEpol(ll)==min(covEpol(ll)));
    lagtime_4_q=-lags4(jj);  %lagtime in points
    flux_4_o_q=min(covEpol(jj));
    
[c5, lags5]=xcov(Wcat(24001:30000),O3ppb_real(24001:30000),100,'unbiased');
    p = polyfit(lags5,c5,26);
    covEpol = polyval(p,lags5);
    ll=find(-100<lags5 & lags5<-35);
    jj=find(covEpol(ll)==min(covEpol(ll)));
    lagtime_5_q=-lags5(jj);  %lagtime in points
    flux_5_o_q=min(covEpol(jj));

[c6, lags6]=xcov(Wcat(30001:36000),O3ppb_real(30001:36000),100,'unbiased');
    p = polyfit(lags6,c6,26);
    covEpol = polyval(p,lags6);
    ll=find(-100<lags6 & lags6<-35);
    jj=find(covEpol(ll)==min(covEpol(ll)));
    lagtime_6_q=-lags6(jj);  %lagtime in points
    flux_6_o_q=min(covEpol(jj));

O3mean_1_o_q=mean(O3ppb_real(1:6000));
O3mean_2_o_q=mean(O3ppb_real(6001:12000));
O3mean_3_o_q=mean(O3ppb_real(12001:18000));
O3mean_4_o_q=mean(O3ppb_real(18001:24000));
O3mean_5_o_q=mean(O3ppb_real(24001:30000));
O3mean_6_o_q=mean(O3ppb_real(30001:36000));

O3stdev_1_o_q=std(O3ppb_real(1:6000));
O3stdev_2_o_q=std(O3ppb_real(6001:12000));
O3stdev_3_o_q=std(O3ppb_real(12001:18000));
O3stdev_4_o_q=std(O3ppb_real(18001:24000));
O3stdev_5_o_q=std(O3ppb_real(24001:30000));
O3stdev_6_o_q=std(O3ppb_real(30001:36000));


%Ozone new methodO3pp
O3ppb_real=O3ppb_imp;
a=5;
O3ppb_real=O3ppb_real.*(1+a*r');
O3ppb_real=O3ppb_real(52:35991); 

O3mean_1_q=mean(O3ppb_real(1:5990));
O3mean_2_q=mean(O3ppb_real(5991:11980));
O3mean_3_q=mean(O3ppb_real(11981:17970));
O3mean_4_q=mean(O3ppb_real(17971:23960));
O3mean_5_q=mean(O3ppb_real(23961:29950));
O3mean_6_q=mean(O3ppb_real(29951:35940));

O3stdev_1_q=std(O3ppb_real(1:5990));
O3stdev_2_q=std(O3ppb_real(5991:11980));
O3stdev_3_q=std(O3ppb_real(11981:17970));
O3stdev_4_q=std(O3ppb_real(17971:23960));
O3stdev_5_q=std(O3ppb_real(23961:29950));
O3stdev_6_q=std(O3ppb_real(29951:35940));

flux_1_q=cov(Wcat(1:5990),O3ppb_real(1:5990));
flux_1_q=flux_1_q(1,2);

flux_2_q=cov(Wcat(5991:11980),O3ppb_real(5991:11980));
flux_2_q=flux_2_q(1,2);

flux_3_q=cov(Wcat(11981:17970),O3ppb_real(11981:17970));
flux_3_q=flux_3_q(1,2);

flux_4_q=cov(Wcat(17971:23960),O3ppb_real(17971:23960));
flux_4_q=flux_4_q(1,2);

flux_5_q=cov(Wcat(23961:29950),O3ppb_real(23961:29950));
flux_5_q=flux_5_q(1,2);

flux_6_q=cov(Wcat(29951:35940),O3ppb_real(29951:35940));
flux_6_q=flux_6_q(1,2);



           % Power Spectra corrected for quenching    
%            [Sw,F]=psd2(detrend(O3ppb_real),length(O3ppb_real),10.);
%            [Sws,Fsw]=specsmoo(Sw,10.);
%            
%            figure; loglog(Fsw(length(Fsw)-31:length(Fsw)),((Fsw(length(Fsw)-31:length(Fsw))).^(-2/3))./(Fsw(length(Fsw)-18).^(-2/3)).*Sws(length(Fsw)-18).*Fsw(length(Fsw)-18),'r');
%            hold on;loglog(Fsw,Sws.*Fsw,'g');
%            
%            title(['O3 Spectra (quenching corrected:   Day ',jd,hr]);
%            xlabel('Frequency in Hz');
%            ylabel('Spectrum * Frequency');
%            print([way_rawdata,'\GOMECC\Analysis\spectraO3_' jd, hr,  '.jpg'],'-djpeg');
            
%            disp(['Calculating cross spectra for ',num2str(hhh)])
%            [C_wO3,Fson_wO3]=csd2(detrend(Wcat(1:35940)),detrend(O3ppb_real),length(O3ppb_real),10,hamming(length(O3ppb_real)));
%            [C_wO3s,Fson_wO3s]=specsmoo(C_wO3,10);
%            figure;semilogx(Fson_wO3s,C_wO3s.*Fson_wO3s,'g');
%            
%            title(['wO3 Cross Spectra:   Day ',jd,hr]);
%            xlabel('Frequency in Hz');
%            ylabel('Spectrum * Frequency');
           %print([way_rawdata,'\GOMECC\Analysis\cross_spectraWO3_' jd, hr,  '.jpg'],'-djpeg');
           
           
%humidity flux
wq_1=cov(Wcat(1:6000),r(1:6000));
format long g
wq_1=wq_1(1,2);

wq_2=cov(Wcat(6001:12000),r(6001:12000));
wq_2=wq_2(1,2);

wq_3=cov(Wcat(12001:18000),r(12001:18000));
wq_3=wq_3(1,2);

wq_4=cov(Wcat(18001:24000),r(18001:24000));
wq_4=wq_4(1,2);

wq_5=cov(Wcat(24001:30000),r(24001:30000));
wq_5=wq_5(1,2);

wq_6=cov(Wcat(30001:36000),r(30001:36000));
wq_6=wq_6(1,2);
           

wO3_1=(1+a*nanmean(r(1:6000))).*flux_1+a*O3mean_1*wq_1;
wO3_2=(1+a*nanmean(r(6001:12000))).*flux_2+a*O3mean_2*wq_2;
wO3_3=(1+a*nanmean(r(12001:18000))).*flux_3+a*O3mean_3*wq_3;
wO3_4=(1+a*nanmean(r(18001:24000))).*flux_4+a*O3mean_4*wq_4;
wO3_5=(1+a*nanmean(r(24001:30000))).*flux_5+a*O3mean_5*wq_5;
wO3_6=(1+a*nanmean(r(30001:36000))).*flux_6+a*O3mean_6*wq_6;

fluxes=[ddd+hhh/24+00/24/60 lagtime_1 O3mean_1_o O3stdev_1_o flux_1_o O3mean_1 O3stdev_1 flux_1 lagtime_1_q O3mean_1_o_q O3stdev_1_o_q flux_1_o_q O3mean_1_q O3stdev_1_q flux_1_q mean(r(1:6000)) wq_1 wO3_1
         ddd+hhh/24+10/24/60 lagtime_2 O3mean_2_o O3stdev_2_o flux_2_o O3mean_2 O3stdev_2 flux_2 lagtime_2_q O3mean_2_o_q O3stdev_2_o_q flux_2_o_q O3mean_2_q O3stdev_2_q flux_2_q mean(r(6001:12000)) wq_2 wO3_2
         ddd+hhh/24+20/24/60 lagtime_3 O3mean_3_o O3stdev_3_o flux_3_o O3mean_3 O3stdev_3 flux_3 lagtime_3_q O3mean_3_o_q O3stdev_3_o_q flux_3_o_q O3mean_3_q O3stdev_3_q flux_3_q mean(r(12001:18000)) wq_3 wO3_3
         ddd+hhh/24+30/24/60 lagtime_4 O3mean_4_o O3stdev_4_o flux_4_o O3mean_4 O3stdev_4 flux_4 lagtime_4_q O3mean_4_o_q O3stdev_4_o_q flux_4_o_q O3mean_4_q O3stdev_4_q flux_4_q mean(r(18001:24000)) wq_4 wO3_4
         ddd+hhh/24+40/24/60 lagtime_5 O3mean_5_o O3stdev_5_o flux_5_o O3mean_5 O3stdev_5 flux_5 lagtime_5_q O3mean_5_o_q O3stdev_5_o_q flux_5_o_q O3mean_5_q O3stdev_5_q flux_5_q mean(r(24001:30000)) wq_5 wO3_5
         ddd+hhh/24+50/24/60 lagtime_6 O3mean_6_o O3stdev_6_o flux_6_o O3mean_6 O3stdev_6 flux_6 lagtime_6_q O3mean_6_o_q O3stdev_6_o_q flux_6_o_q O3mean_6_q O3stdev_6_q flux_6_q mean(r(30001:36000)) wq_6 wO3_6];

% save([way_rawdata, '\GOMECC\fluxes_q_' , jd, hr,'.txt'], 'fluxes', '-ascii')% 
fluxes_jd=[fluxes_jd;fluxes];
     
    end;
end;

 save([way_rawdata, '\GOMECC\fluxes_q_' , jd, '.txt'], 'fluxes_jd', '-ascii')%   

