%% Extract 1D arrays
salt661 = squeeze(salt(661,1:280,:));
v661 = isfinite(salt661);
salt661 = salt661(v661);
sr661 = squeeze(sr(661,1:280,:));
sr661 = sr661(v661); 
tau_remi661 = squeeze(tau_remi(661,1:280,:));
tau_remi661 = tau_remi661(v661);
tau_remi_ref661 = squeeze(tau_remi_ref(661,1:280,:));
tau_remi_ref661 = tau_remi_ref661(v661); 
tau_jmd661 = squeeze(tau_jmd(661,1:280,:));
tau_jmd661 = tau_jmd661(v661);
tau_huang661 = squeeze(tau_huang(661,1:280,:));
tau_huang661 = tau_huang661(v661);
tau_huang_ref661 = squeeze(pi_ref(661,1:280,:));
tau_huang_ref661 = tau_huang_ref661(v661);
ct661 = squeeze(ct(661,1:280,:));
ct661 = ct661(v661); 
sig1_661 = squeeze(sigma1(661,1:280,:));
sig1_661 = sig1_661(v661);
sig1_lin = 27:0.1:33;
gtlin = 22:0.1:29; 
gt_661 = squeeze(gammat(661,1:280,:));
gt_661 = gt_661(v661); 

%% Construct polynomial fit 
[p,S,mu] = polyfit(sig1_661,tau_remi661,2);
tau1_lin = polyval(p,sig1_lin,[],mu);

[pjmd,Sjmd,mujmd] = polyfit(sig1_661,tau_jmd661,2);
taujmd_lin = polyval(pjmd,sig1_lin,[],mujmd);

[phuang,Shuang,muhuang] = polyfit(sig1_661,tau_huang661,2);
tauhuang_lin = polyval(phuang,sig1_lin,[],muhuang); 

[psalt,Ssalt,musalt] = polyfit(sig1_661,salt661,2);
salt_lin = polyval(psalt,sig1_lin,[],musalt); 

[pct,Sct,muct] = polyfit(sig1_661,ct661,2);
ct_lin = polyval(pct,sig1_lin,[],muct); 

%% New way to construct polynomial fits 
[ftau,goftau] = fit(gt_661(:),tau_remi_ref661(:), 'poly2', ...
    'Normalize','on','Robust','Bisquare');
taulin = ftau(gtlin);

[fpi,gofpi] = fit(gt_661(:),tau_huang_ref661(:), 'poly2', ...
    'Normalize','on','Robust','Bisquare');
pilin = fpi(gtlin); 

[fsr,gofsr] = fit(gt_661(:),sr661(:), 'poly2', ...
    'Normalize','on','Robust','Bisquare'); 
srlin = fsr(gtlin);

[fct,gofct] = fit(gt_661(:),ct661(:), 'poly2', ...
    'Normalize','on','Robust','Bisquare'); 
ctlin = fct(gtlin); 

%% Make figures with all points for comparison
figure
subplot(2,2,1)
plot(gammat(valid),tau_remi_ref(valid),'.','Color',[0.9290 0.6940 0.1250])
hold on
plot(gt_661,tau_remi_ref661,'.','Color',[0.8500 0.3250 0.0980])
plot(gtlin,taulin,'k')
xlabel('\gamma^T_{analytic}');ylabel('\tau_{ref} (kg/m^3)')
title('(a)')
axis([22 30 -6 10])
hold off
subplot(2,2,2)
plot(gammat(valid),pi_ref(valid),'.','Color',[0.9290 0.6940 0.1250])
hold on
plot(gt_661,tau_huang_ref661,'.','Color',[0.8500 0.3250 0.0980]) 
plot(gtlin,pilin,'k')
xlabel('\gamma^T_{analytic}');ylabel('\pi_{ref} (kg/m^3)') 
title('(b)')
axis([22 30 -6 10])
subplot(2,2,3)
plot(gammat(valid),ct(valid),'.','Color',[0.9290 0.6940 0.1250])
hold on
plot(gt_661,ct661,'.','Color',[0.8500 0.3250 0.0980]) 
plot(gtlin,ctlin,'k')
xlabel('\gamma^T_{analytic}');ylabel('\Theta ({}^{\circ}C)')
title('(c)')
axis([22 30 -4 30])
hold off 
subplot(2,2,4)
plot(gammat(valid),sr(valid),'.','Color',[0.9290 0.6940 0.1250])
hold on
plot(gt_661,sr661,'.','Color',[0.8500 0.3250 0.0980])
plot(gtlin,srlin,'k')
xlabel('\gamma^T_{analytic}');ylabel('S_A (g/kg)') 
title('(d)')
axis([22 30 30 38])
print('jmr_correlation_tau_gammat','-depsc2') 
print('jmr_correlation_tau_gammat','-dpng') 
print('fig10','-dpng') 
print('fig10','-depsc2')
print('fig10','-dpdf') 
