%evaluate 7200lag time for HiWinGS 2013
clear all
cruise='HiWinGS_2013';     % Cruise name
if ~exist('yyyy','var'); yyyy=input('Input year to evaluate as yyyy : ','s'); end;
if ~exist('ddd','var'); ddd=input('Input yearday to evaluate as ddd : '); end;
% jd=sprintf('%03i',ddd);
% [m,d]=yd2md(str2num(yyyy), ddd);
% Vdate=[str2num(yyyy) m d];
way_raw_data_flux='C:\Data\CRUISES\HiWinGS_2013\Knorr\flux\Raw\';
prtit=true;  % for printing plots--ineffective unless plotit is also true going to a temp directory
graphformat='.png';  %select graphics format files
graphdevice='-dpng'; %select graphic device

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Look at lag time from Puff on 7200
% Supposedely on top of each hour
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for ddd=ddd
    jd=sprintf('%03i',ddd);
    [m,d]=yd2md(str2num(yyyy), ddd);
    Vdate=[str2num(yyyy) m d];
    for hhh=00:23
        hr=sprintf('%02i',hhh);
    %     read_son_DYNAMO2011;     %
        licor_prefix='li0';
        read_7200_HiWinGS_2013;
    % end

    % figure(1);clf;
    % plot(t7200,LI7200_CO2mfd);xlim([0 20]);grid
    % xlabel('Time (s)');
    % CO2base=nanmedian(LI7200_CO2mfd(t7200<5));

    ti=t7200((t7200<20));
    CO2i=LI7200_CO2mfd((t7200<20));

    % ti=0:.001:20;
    % CO2i=interp1(t7200((t7200<20)),LI7200_CO2mfd((t7200<20)),ti);

    CO2min=min(CO2i);
    LI7200_CO2mf2=CO2i-CO2min;
    CO2base=nanmedian(LI7200_CO2mf2(t7200<5));
    LI7200_CO2mf2=LI7200_CO2mf2/CO2base;
    figure(1);clf;
    plot(ti,LI7200_CO2mf2);grid
    xlabel('Time (s)');ylabel('Normalized CO2')
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i, hour%02i).  Puff and CO2',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hhh),'FontWeight','Bold','Interpreter','none')

    hold on;plot([0 20],[1-1/exp(1) 1-1/exp(1)],'k--')


    %finds lag
    detpic=[0 diff(CO2i)];
    % ii=find(detpic<min(detpic(t7200<5)));
    % ii=find(detpic<-.5);
    ii=find(detpic<(min(detpic(t7200<5)))+(min(detpic(t7200<5)))/2);
    lag=ti(ii(1)-1);  %assumes DSA clock synchronized with UH sustem... need to check that
    hold on;plot(ti(1:ii(1)-1),LI7200_CO2mf2(1:ii(1)-1),'r.')

    %finds time constant for filter function
    tref=[0:.001:20];
    LI7200_CO2mf2i=interp1(ti,LI7200_CO2mf2,tref);
    % hold on;plot(tref,LI7200_CO2mf2i,'k-')
    ti1=tref(tref<8);
    T1=ti1(abs(LI7200_CO2mf2i(tref<8)-(1-1/exp(1)))==min(abs(LI7200_CO2mf2i(tref<8)-(1-1/exp(1))))); %first time it crosses 1/e on its way down
    hold on;plot(T1,(1-1/exp(1)),'ro')
    pp=find(CO2i==CO2min);
    T2=ti(pp(end)); hold on;plot(T2,0,'rs'); %time when it starts to return up
    ti1=tref(tref>T2);
    T3=ti1(abs(LI7200_CO2mf2i(tref>T2)-(1-1/exp(1)))==min(abs(LI7200_CO2mf2i(tref>T2)-(1-1/exp(1))))); hold on;plot(T3,(1-1/exp(1)),'ro')

    % T1=ti(abs(LI7200_CO2mf2-(1-1/exp(1)))==min(abs(LI7200_CO2mf2-(1-1/exp(1))))); %first time it crosses 1/e on its way down
    % hold on;plot(T1,(1-1/exp(1)),'ro')
    % T2=ti(CO2i==CO2min); hold on;plot(T2,0,'ro'); %time when it starts to return up
    % ti1=ti((ti>T2));
    % T3=ti1(abs(LI7200_CO2mf2(ti>T2)-(1-1/exp(1)))==min(abs(LI7200_CO2mf2(ti>T2)-(1-1/exp(1))))); hold on;plot(T3,(1-1/exp(1)),'ro')

    % ii=find(detpic>max(abs(detpic(t7200<5))));T2=ti(ii(1)-1);hold on;plot(T2,0,'ro'); 
    ii=find(detpic>0.5);iii=find(ti(ii)>10);
    % T2=ti(ii(1)-1);hold on;plot(T2,0,'ro'); 
    T2=ti(ii(iii(1))-1);hold on;plot(T2,0,'ro'); 
    % figure;plot(ti(ti>=T2 & ti<=T3),[0 diff(CO2i(ti>=T2 & ti<=T3))],'+')
    % <<--to explore if the above does not work 

    % [lag T1 T3-T2]
    tauc=T3-T2; % time constant estimation.

    oo=find(ti==lag);
    tau=ti(oo);
    st=ones(1,length(ti));
    st(1,oo:length(ti))=exp((tau-ti(oo:length(ti)))./tauc); %
    % figure;plot(ti,st,'r--');
    oo2=find(ti==T2);
    tau2=ti(oo2+2);
    st2=st;
    st2(1,oo2+2:length(ti))=1-exp((tau2-ti(oo2+2:length(ti)))./tauc); %
    % st2(st2<0)==0
    hold on;plot(ti,st2,'r--');


    oo2=find(ti>=T2);
    % figure;plot(ti(oo2(4:end)),LI7200_CO2mf2(oo2(4:end)),ti(oo2(4:end)),1-exp((tau2-ti(oo2(4:end)))./tauc),'r')
    % figure;plot(LI7200_CO2mf2(oo2(4:end)),1-exp((tau2-ti(oo2(4:end)))./tauc),'.',[0 1.2],[0 1.2],'k--')
    % r=corrcoef(LI7200_CO2mf2(oo2(4:end)),1-exp((tau2-ti(oo2(4:end)))./tauc));r(1,2)

    taul=0.1:0.01:0.7;
    r2=[];
    for ll=1:length(taul)
        r=corrcoef(LI7200_CO2mf2(oo2(4:end)),1-exp((tau2-ti(oo2(4:end)))./taul(ll)));
        r2(ll)=r(1,2);
    end    
    tauc2=taul(r2==max(r2));%figure;plot(taul,r2)

    oo=find(ti==lag);
    tau=ti(oo);
    st=ones(1,length(ti));
    st(1,oo:length(ti))=exp((tau-ti(oo:length(ti)))./tauc2); %
    % figure;plot(ti,st,'r--');
    oo2=find(ti==T2);
    tau2=ti(oo2+2);
    st2=st;
    st2(1,oo2+2:length(ti))=1-exp((tau2-ti(oo2+2:length(ti)))./tauc2); %
    % st2(st2<0)==0
    hold on;plot(ti,st2,'g--');

    if prtit; 
        print(graphdevice,['C:\Data\CRUISES\HiWinGS_2013\7200\puff\' cruise '_puff_' sprintf('%04i_%02i_%02i_%03i_hour%02i',Vdate(1),Vdate(2),Vdate(3),ddd,hhh) graphformat]);
    %     saveas(gcf,['C:\Data\CRUISES\DYNAMO_2011\CO2\puff_v2', jd,hr '.fig'])
    end


    %%Now look at synch between UH and PSD systems
    %start by using the cross-correlation approach on the vertical
    %accelerations

    %**********************************************
    %gets PSD motion pack
    read_motion_HiWinGS_2013;
    %**********************************************
    %gets UH motion
    read_dms_HiWinGS_2013
    time_dms = dmsdata(1,:)-dmsdata(1,1); % time in s
    UH_accz=9.8*dmsdata(11,:)./nanmedian(dmsdata(11,:));   
    UH_Lag_Pulse=dmsdata(12,:);
    clear dmsdata

    % figure;plot(tref,accz,time_dms,UH_accz,'r');grid
    ClockDrift = findlag(UH_accz,accz,20,0);  %positive means PSD delayed WRT UH
    % figure;plot(tref+ClockDrift/10,accz,time_dms,UH_accz,'r');grid
    % figure;plot(time_dms,UH_Lag_Pulse,'r+--');grid;xlim([-1 10])
    uhi=find(UH_Lag_Pulse>4);[time_dms(uhi(1)) time_dms(uhi(end))]


    puffmat=[ddd hhh lag T1 T2 T3 tauc tauc2 ClockDrift time_dms(uhi(1)) time_dms(uhi(end))]
    puffmat=double(puffmat);
    sdg=['save ', 'C:\Data\CRUISES\HiWinGS_2013\7200\puff\HiWinGS2013_pufflag_',jd hr,'.prn puffmat -ascii'];
    eval(sdg)

    end
    filecat(['C:\Data\CRUISES\HiWinGS_2013\7200\puff\HiWinGS2013_pufflag_',jd, '*.prn'],['HiWinGS2013pufflag_',jd, '.prn'])
end
% 
% % filecat(['C:\Data\CRUISES\HiWinGS_2013\7200\puff\HiWinGS2013pufflag_*.prn'],['HiWinGS2013pufflag_all.prn'])
% 
CO2lagall=load(['C:\Data\CRUISES\HiWinGS_2013\7200\puff\HiWinGS2013pufflag_all.prn']);

%sometimes the outliers are du to the fact that the clockdrift is for the
%entire hour 0 wheras the lag should be corrected from previous hour
%drifts...
as=nanmedian(CO2lagall(:,3)-CO2lagall(:,9)/10-CO2lagall(:,10))+2*nanstd(CO2lagall(:,3)-CO2lagall(:,9)/10-CO2lagall(:,10));
am=nanmedian(CO2lagall(:,3)-CO2lagall(:,9)/10-CO2lagall(:,10))-2*nanstd(CO2lagall(:,3)-CO2lagall(:,9)/10-CO2lagall(:,10));
aai=find(CO2lagall(:,3)-CO2lagall(:,9)/10-CO2lagall(:,10) > as | CO2lagall(:,3)-CO2lagall(:,9)/10-CO2lagall(:,10) < am);

truelag=(CO2lagall(:,3)-CO2lagall(:,9)/10-CO2lagall(:,10));
truelag(aai)=CO2lagall(aai,3)-CO2lagall(aai-1,9)/10-CO2lagall(aai,10);

figure(2);clf;
set(gca,'fontsize',14);
subplot(2,2,1);
plot(CO2lagall(:,1)+CO2lagall(:,2)/24,CO2lagall(:,3),'b.')
% hold on;plot(CO2lagall(:,1)+CO2lagall(:,2)/24,CO2lagall(:,3)-CO2lagall(:,9)/10-CO2lagall(:,10),'g.')
hold on;plot(CO2lagall(:,1)+CO2lagall(:,2)/24,truelag,'g.')
hold on;plot(CO2lagall(:,1)+CO2lagall(:,2)/24,CO2lagall(:,5)-CO2lagall(:,9)/10-CO2lagall(:,11),'m.')
grid;xlabel('DOY (UTC)');ylabel('Lag (s)')
legend(sprintf('%.1fs',nanmedian(CO2lagall(:,3))),sprintf('%.1f +/- %.1fs',nanmedian(CO2lagall(:,3)-CO2lagall(:,9)/10-CO2lagall(:,10)),nanstd(CO2lagall(:,3)-CO2lagall(:,9)/10-CO2lagall(:,10))),sprintf('%.1f +/- %.1fs',nanmedian(CO2lagall(:,5)-CO2lagall(:,9)/10-CO2lagall(:,11)),nanstd(CO2lagall(:,5)-CO2lagall(:,9)/10-CO2lagall(:,11))),'Location','NorthWest')
xlim([283 CO2lagall(end,1)+1])
title(sprintf('%s .  Puff results',cruise),'FontWeight','Bold','Interpreter','none')

subplot(2,2,2);plot(CO2lagall(:,1)+CO2lagall(:,2)/24,CO2lagall(:,9)/10,'b.')
grid;xlabel('DOY (UTC)');ylabel('ClockDrift (s)')
% legend('Positive means PSD delayed WRT UH','Location','SouthWest')
title('Positive Drift = PSD delayed WRT UH')
xlim([283 CO2lagall(end,1)+1])

subplot(2,2,3);plot(CO2lagall(:,1)+CO2lagall(:,2)/24,CO2lagall(:,7),'b.',CO2lagall(:,1)+CO2lagall(:,2)/24,CO2lagall(:,8),'g.')
grid;xlabel('DOY (UTC)');ylabel('Time constant (s)')
xlim([283 CO2lagall(end,1)+1])

subplot(2,2,4);plot(CO2lagall(:,1)+CO2lagall(:,2)/24,1./(2*pi*CO2lagall(:,7)),'b.',CO2lagall(:,1)+CO2lagall(:,2)/24,1./(2*pi*CO2lagall(:,8)),'g.')
grid;xlabel('DOY (UTC)');ylabel('Cutoff frequency (Hz)')
xlim([283 CO2lagall(end,1)+1])

if prtit; print(graphdevice,['C:\Data\CRUISES\HiWinGS_2013\7200\puff\HiWinGS2013_puffresults' graphformat]); end
