% 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.
read_parameters
prtit=false;
saveit=false;
plotit=false;

stday=294; % yeardays the snake was in the water tracking TSG all night
enday=299;
for ddd=stday:enday
    read_pc_day;
    read_scs_day;
    % append and save 1-minute data
    fid=fopen([way_proc_data_flux 'snake_1min' sprintf('%03i',ddd) '.txt'],'wt');
    fprintf(fid,'%f\t%f\n',[jd_pc' Tsea']');
    fclose(fid);
    fid=fopen([way_proc_data_flux 'tsg_1min' sprintf('%03i',ddd) '.txt'],'wt');
    fprintf(fid,'%f\t%f\n',[jd_scs tsgm]');
    fclose(fid);
end
% concatenate data outside matlab

% custom axis limits
xlim=[...
   24   28
   18   25
   19   21
   17   22
   17   20
   18   19 ];
ylim=[...
   26   30
   20   27
   21   23
   19   24
   19   22
   20   21 ];

% load concatenated data
snakedata=load([way_proc_data_flux 'snake_1min.txt']);
jd_pc=snakedata(:,1);
Tsea=snakedata(:,2);
tsgdata=load([way_proc_data_flux 'tsg_1min.txt']);
jd_scs=tsgdata(:,1);
tsgm=tsgdata(:,2);
Ttsg=interp1(jd_scs,tsgm,jd_pc);

for ddd=stday:enday
    ind=ddd-stday+1;
     & ... % night
       floor(jd_pc)==ddd;                            % select day
    Ttsg_mean(ind)=mean(Ttsg(ii));
    Tsea_mean(ind)=mean(Tsea(ii));
    P=polyfit(Ttsg(ii)-Ttsg_mean(ind),Tsea(ii)-Tsea_mean(ind),1);
    Ttsg_std(ind)=std(Ttsg(ii));
    Tsea_std(ind)=std(Tsea(ii));
    covtmp=cov(Ttsg(ii),Tsea(ii));
    Cov(ind)=covtmp(2,1);
    bias(ind)=Tsea_mean(ind)-Ttsg_mean(ind);   % bias
    gain(ind)=P(1); % gain
    
    if true % true to replot
        ax(ind)=subplot(3,2,ind,'align');
        xwidth=[-6 2];
        
        plot(Ttsg,Tsea,'.')
        hold on
        plot(Ttsg(ii),Tsea(ii),'r.')
        plot(Ttsg_mean(ind),Tsea_mean(ind),'k+')
        plot(Ttsg_mean(ind)+xwidth,Ttsg_mean(ind)+bias(ind)+xwidth*gain(ind),'k')
        %axis([18 28 20 30])
        axis square
        axis([xlim(ind,:) ylim(ind,:)])
        %axis tight
        xlabel('SCS thermosalinograph');
        ylabel('PSD snake');
        %     text(18.5, 19.5,sprintf('bias at mean = %+5.2f\\circC',bias))
        %     text(18.5, 19.2,sprintf('gain = %3.3f', gain))
        %     text(19.2, 21.5,'night','color','r')
        title(['yearday ' sprintf('%3i',ddd)])
    end % plot/don't plot
end
orient tall
print('-depsc2',[way_raw_images_flux 'snake_tsg_comparison' sprintf('%03i-%03i',stday,enday) '.eps'])

Cor=Cov./(Tsea_std.*Ttsg_std);
unexpl_Tsnake_dev=Tsea_std.*sqrt(1-Cor.^2); % degree C

save([way_proc_data_flux 'snake_cal.mat'], 'bias','gain','Cov','Cor','unexpl_Tsnake_dev');

% plot concatenated timeseries of tsg and snake
figure
ii=rem(jd_pc,1)<0.45 & isfinite(Ttsg+Tsea);
plot(jd_pc,Tsea-Ttsg,'b.','markersize',3)
hold on
plot(jd_pc(ii),Tsea(ii)-Ttsg(ii),'r.','markersize',3)
set(gca,'ylim',[1.5 2])
orient landscape
print('-depsc2',[way_raw_images_flux 'snake_tsg_timeseries294-299.eps'])

figure
plot(Ttsg(ii),Tsea(ii)-Ttsg(ii),'b.','markersize',2)
axis square; axis tight
xlabel('TSG (C)')
ylabel('SeaSnake-TSG (C)')