% Reads 1-minute raw data b.c. it is high resolution before substitution,
% adjustment, or QC.
% Tsea(jd_pc) = snake
% tsgm(jd_scs) = ship thermosalinograph
% Compute the bias of the sea snake relative to the ship thermosalinograph
% during the night time hours, when there is no solar warm layer.

% cd E:\VOCALS_2008\RHB\Scientific_analysis\programs\flux\
read_parameters
prtit=false;
saveit=false;
plotit=false;

stday=288; % yeardays
enday=302;
for ddd=stday:enday
    if ~exist([way_proc_data_flux 'trh_pc_1min' sprintf('%03i',ddd) '.txt'],'file')
        read_pc_day;
        % save 1-minute data
        fid=fopen([way_proc_data_flux 'trh_pc_1min' sprintf('%03i',ddd) '.txt'],'wt');
        fprintf(fid,'%f\t%f\t%f\n',[jd_pc' Tvais' Rhvais']');
        fclose(fid);
    end
    if ~exist([way_proc_data_flux 'trh_scs' sprintf('%03i',ddd) '.txt'],'file')
        read_scs_day;
        fid=fopen([way_proc_data_flux 'trh_scs' sprintf('%03i',ddd) '.txt'],'wt');
        fprintf(fid,'%f\t%f\t%f\n',[ddd+gpstim'/24 ta' rh']'); %Fing inconsistent time stamps
        fclose(fid);
    end
end

if true
    % concatenate files
    filecat([way_proc_data_flux 'trh_scs*.txt'],[way_proc_data_flux 'cat/trh_scs.txt'])
    filecat([way_proc_data_flux 'trh_pc_1min*.txt'],[way_proc_data_flux 'cat/trh_pc_1min.txt'])

    % load concatenated data
    pcdata=load([way_proc_data_flux 'cat/trh_pc_1min.txt']);
    jdpc=pcdata(:,1);
    Tapc_0=pcdata(:,2);
    RHpc_0=pcdata(:,3);
    scsdata=load([way_proc_data_flux 'cat/trh_scs.txt']); % slow
    jdscs=scsdata(:,1);
    Tascs_0=scsdata(:,2); % at 2s time intervals
    RHscs_0=scsdata(:,3);

    jd=min(jdpc):(1/24/60):max(jdpc);
    Tascs=intavg(jdscs,Tascs_0,jd); % SLOW
    RHscs=intavg(jdscs,RHscs_0,jd);
    Tapc=intavg(jdpc,Tapc_0,jd);
    RHpc=intavg(jdpc,RHpc_0,jd);
    save([way_proc_data_flux 'trh_1min.mat'], 'jd', 'Tapc', 'RHpc', 'Tascs', 'RHscs');
else
    load([way_proc_data_flux 'trh_1min.mat']);
end

% temperature difference and lowpass filter of difference
dT=Tascs-Tapc;
width=60;
offset=1:width/2:length(Tascs);
weight=hanning(width);
for i=1:length(offset)-2
    dTlo(i)=nansum(weight.*dT(offset(i):min(length(dT),offset(i)+width-1)))/sum(weight.*isfinite(dT(offset(i):min(length(dT),offset(i)+width-1))));
end

% repeat Nan's plot
figure
subplot(3,1,1,'align')
plot(jd,RHscs,'b')
hold on
plot(jd,RHpc,'k')
axis([stday enday 55 110])
subplot(3,1,2,'align')
plot(jd,Tascs,'b')
hold on
plot(jd,Tapc,'k')
axis([stday enday 16 30])
% plot temperature difference
subplot(3,1,3,'align')
plot(jd,dT)
hold on
plot(jd(offset(1:end-2)),dTlo,'r','linewidth',1.5)
set(gca,'ylim',[-0.25 1.0])
ylabel('T ship-PSD (\circ C)')
xlabel('yearday');
orient tall
print('-depsc2',[way_images_flux 'ta_rh_compare_ts.eps'])


[Tascslc,lags,nobs]=lagcov(Tascs,Tascs,1000); %lagged autocovariance
Tapclc=lagcov(Tapc,Tapc,1000);
RHscslc=lagcov(RHscs,RHscs,1000);
RHpclc=lagcov(RHpc,RHpc,1000);
Tlxc=lagcov(Tascs,Tapc,1000); %lagged cross covariance - positive lag means a lags b
RHlxc=lagcov(RHscs,RHpc,1000);
TRHscslxc=lagcov(Tascs,RHscs,1000);
TRHpclxc=lagcov(Tapc,RHpc,1000);

% plot lag auto and cross correlations
figure
subplot(2,1,1)
plot(lags/60,RHscslc/nancov(RHscs),'-')
hold on
plot(lags/60,RHpclc/nancov(RHpc),'--')
plot(lags/60,RHlxc/nanstd(RHpc)/nanstd(RHscs),'r','linewidth',1.5) % lags by 1.5 hour
axis tight
title('RH'); ylabel('correlation');
legend('ship','PSD','cross')
subplot(2,1,2)
plot(lags/60,Tascslc/nancov(Tascs),'-')
hold on
plot(lags/60,Tapclc/nancov(Tapc),'--')
plot(lags/60,Tlxc/max(Tlxc),'r','linewidth',1.5) % lags by 1.5 hour
axis tight
title('T_{air}'); ylabel('correlation');
xlabel('ship leads <           lag (hours)             > PSD leads')
print('-depsc2',[way_images_flux 'ta_rh_compare_lag.eps'])