

clear
close all


% import table

g = readtable('GDF6_initial_cal.csv', 'headerlines', 24);

g.Properties.VariableNames = {'DT', 'TC', 'tc_bath', 'bath_sp', 'pres', 'maxflow', 'flowrate', 'Vpump', ...
    'crc', 'id', 'msc_dt', 'utc_sec', 'ctd_depth', 'ctd_tc', 'ctd_salinity', 'scount', 'reset_count', 'error_count', ...
    'tc_msc', 'humid', 'Vinput', 'Iinput', 'foobar', 'Vbackup', 'Vrs', 'Vrs_std', 'Vk', 'Vk_std', ...
    'Ik', 'Ib'};

% calculate k2, f(P)
g.meangamma = -0.000113.*g.TC + 0.9067;
g.delVo = 16.38 + 0.0896*g.TC - 0.001264.*g.TC.^2;
g.TK = g.TC + 273.15;

g.Vct = -R().*g.TK./F() .* log((g.meangamma.*0.01).^2);
g.Vcp = -g.delVo./23061.*0.0242.*(g.pres-10)/10*1/1.01;

g.Ktp = g.Vrs - g.Vct - g.Vcp;

% do k2
ik2 = g.pres < 40;

c = polyfit(g.TC(ik2), g.Ktp(ik2), 1);
k2 = c(1); k0 = c(2);

% calculate fP
g.fP = (g.Ktp - k2.*g.TC - k0).*1000; % do it in mV

% indicies for fP
iPincrease = g.flowrate > 0.01;
iPdecrease = g.flowrate < -0.01;

i20C = inrange(g.TC, [19 21]);
i15C = inrange(g.TC, [14 16]);
i12C = inrange(g.TC, [11 13]);

% bin data
Xpres = [30:5:2000];

[Xup, iXup] = sort(g.pres(iPincrease));
    Yup = g.fP(iPincrease); Yup = Yup(iXup);
fPup = bin_data_to_X(Xup, Yup, Xpres);

[Xdn, iXdn] = sort(g.pres(iPdecrease));
    Ydn = g.fP(iPdecrease); Ydn = Ydn(iXdn);
fPdn = bin_data_to_X(Xdn, Ydn, Xpres);

dfP = fPup - fPdn; 

%% plot for k2

f1 = figure; 
hold on; box on;
plot(g.TC, g.Ktp, 'k.');
plot(g.TC(ik2), g.Ktp(ik2), 'r.');
plot(g.TC, polyval(c, g.TC), 'b');
xlabel('Temperature');
ylabel('Vrs - Vct - Vcp');
legend('all data', '1atm only');
title(['k2 = ',num2str(k2*1000,'%1.3f'),' mV/C, k0 = ',num2str(k0,'%1.6f'),' V']);
saveas(f1, 'k2_k0_1atm.jpg');

%% fP related plots


f2 = figure;
hold on; box on;
plot(g.pres(iPincrease), g.fP(iPincrease), 'k.');
plot(g.pres(iPdecrease), g.fP(iPdecrease), 'r.');
xlabel('pressure (dBar)');
ylabel('f(P) [mV]');
xlim([0 2100]);
legend('Increasing P', 'Decreasing P');
title('All temperatures');
saveas(f2, 'fP_hysteresis_alltemp.jpg');

%%

f3 = figure('unit', 'inch', 'position', [1 1 16 9]);

tiledlayout(2,2,'TileSpacing','compact');

nexttile;
    hold on; box on;
    plot(g.pres(iPincrease), g.fP(iPincrease), 'k.');
    plot(g.pres(iPdecrease), g.fP(iPdecrease), 'r.');
    xlabel('pressure (dBar)');
    ylabel('f(P) [mV]');
    xlim([0 2100]);
    legend('Increasing P', 'Decreasing P');
    title('All temperatures');  
    
nexttile;
    hold on; box on;
    scatter(g.pres(iPincrease), g.fP(iPincrease), 10, g.TC(iPincrease), 'filled');
    colormap(cbrewer('div', 'RdYlBu', 1000, 'PCHIP'));   
    colorbar;
    xlabel('Pressure (dBar)');
    ylabel('f(P) [mV]');
    title('Increasing P');
    
nexttile;
    hold on; box on;
    scatter(g.pres(iPdecrease), g.fP(iPdecrease), 10, g.TC(iPdecrease), 'filled');
    colormap(cbrewer('div', 'RdYlBu', 1000, 'PCHIP'));   
    colorbar;
    xlabel('Pressure (dBar)');
    ylabel('f(P) [mV]');
    title('Decreasing P');

nexttile;
    hold on; box on;
    plot(Xpres, dfP, 'k.');
    ylabel('\DeltafP (Pinc - Pdec)');
    xlabel('Pres');
    title('all data binned');
    
saveas(f3, 'fP_hysteresis_summary.jpg');




%% Environment and diagnostic plots

f3 = figure('unit', 'inch', 'position', [1 1 16 9]);

tiledlayout(2,1,'TileSpacing','compact');

% pres and temp
nexttile;
hold on; box on;

yyaxis left
    plot(g.DT, g.pres, 'b');
    ylabel('Pressure (dBar)');

yyaxis right
    plot(g.DT, g.TC, 'r');
    ylabel('Temperature');
%     ylim([14.75 14.85]); % comment out and file name below to get zoom in
    datetick;
    
nexttile;
hold on; box on;

yyaxis left
    plot(g.DT, g.Vrs, 'b');
    ylabel('Vrs');

yyaxis right
    plot(g.DT, g.Vk, 'r');
    ylabel('Vk');
    datetick;

saveas(f3, 'PT_cycle_and_Vrs_Vk.jpg');
% saveas(f3, 'PT_cycle_and_Vrs_Vk_TCzoomin.jpg');

% Std and Ik/Ib
f4 = figure('unit', 'inch', 'position', [1 1 16 9]);

tiledlayout(2,1,'TileSpacing','compact');

% pres and temp
nexttile;
hold on; box on;

yyaxis left
    plot(g.DT, g.Ik.*1E9, 'b');
    ylabel('Ik (nA)');

yyaxis right
    plot(g.DT, g.Ib.*1E9, 'r');
    ylabel('Ib (nA)');
    datetick;
    
nexttile;
hold on; box on;

yyaxis left
    plot(g.DT, g.Vrs_std.*1E6, 'b');
    ylabel('Vrs std (uV)');

yyaxis right
    plot(g.DT, g.Vk_std.*1E6, 'r');
    ylabel('Vk std (uV)');
    datetick;

saveas(f4, 'Diagnostics_stacked.jpg');

    
%%


% Std and Ik/Ib
f4 = figure('unit', 'inch', 'position', [1 1 16 9]);

tiledlayout(2,1,'TileSpacing','compact');

% pres and temp
nexttile;
hold on; box on;
    
yyaxis left
    plot(g.DT, g.Vk, 'b');
    ylabel('Vk');
yyaxis right
    plot(g.DT, g.flowrate, 'r');
    ylabel('flow rate (mL/min)');
    datetick;

nexttile;
yyaxis left
    plot(g.DT, g.Vk, 'b');
    ylabel('Vk');
yyaxis right
    plot(g.DT, g.Ik.*1E9, 'r');
    ylabel('Ik (nA)');
    datetick;
    
saveas(f4, 'Vkspikes_with_others.jpg');
    
%% thermal response

f5 = figure;
hold on; box on;
plot(g.DT, g.TC, 'k');
plot(g.DT, ones(size(g.DT)).*14.83, 'r');
plot(g.DT, ones(size(g.DT)).*14.98, 'r:');
ylabel('temperature');
xlim([datetime('4/21/2021 14:00') datetime('4/21/2021 16:00')]);
datetickzoom;

saveas(f5, '20C_to_15C_zoomin.jpg');

f5 = figure;
hold on; box on;
plot(g.DT, g.TC, 'k');
plot(g.DT, ones(size(g.DT)).*11.88, 'r');
plot(g.DT, ones(size(g.DT)).*11.98, 'r:');


ylabel('temperature');
xlim([datetime('4/22/2021 7:00') datetime('4/22/2021 8:30')]);
datetickzoom('x', 'keeplimits');

saveas(f5, '15C_to_12C_zoomin.jpg');

% make basic figures
% 
% figure
% hold on; box on;



