conn = connect_penguin();
query_k0_seawater_lab = "SELECT * FROM k0_seawater_lab WHERE gliderfet_id IN ('10', '16') AND datetime > '2026-03-11' ORDER BY datetime ASC";
query_FT_k0_spec_ph = "SELECT * FROM FT_k0_spec_ph WHERE LOWER(comments) != 'junk'";
data = fetch(conn, query_k0_seawater_lab);
spec = fetch(conn, query_FT_k0_spec_ph);

data.datetime = datetime(data.datetime);
data.datetime.TimeZone = 'UTC';

for i = 1:height(spec)
    datetimestr = strsplit(string(spec.Comments(i)),'_');
    spec.datetime_sample(i) = datetime(strjoin(datetimestr(2:3)),InputFormat='uuuuMMdd HHmmSS');
end
spec.datetime_sample.TimeZone = 'UTC';
spec.datetime_sample.Format = "dd/MM/uuuu HH:mm:ss";
spec.datetime_sample = spec.datetime_sample - minutes(2);
spec.pHin = NaN(height(spec), 1);  % pre-initialize pHin column

gfets = ["10", "16"];
k2s = [0.000371965, 0.000481629];

k0_output = table('Size',[0,3], ...
    'VariableTypes',{'string','double','datetime'}, ...
    'VariableNames',{'gliderfet_id','k0_seawater','datetime_sample'});
k0_output.datetime_sample.TimeZone = 'UTC';

for g = 1:length(gfets)
    d = data(ismember(data.gliderfet_id, gfets(g)),:);
    tolerance = seconds(60);
    time_a = spec.datetime_sample;
    time_b = d.datetime;
    idx_a = zeros(length(time_a), 1);
    idx_b = zeros(length(time_a), 1);

    for i = 1:length(time_a)
        diffs = abs(time_b - time_a(i));
        [minDiff, idx] = min(diffs);
        if minDiff <= tolerance
            idx_a(i) = i;
            idx_b(i) = idx;
        end
    end

    valid = idx_a > 0;
    idx_a = idx_a(valid);
    idx_b = idx_b(valid);

    for j = 1:length(idx_b)
        trex = CO2SYS(2270, spec.pH_tot(idx_a(j)), 1, 3, d.salinity(idx_b(j)), spec.Temp_Samp(idx_a(j)), d.temperature(idx_b(j)), 0, 0, 0, 0, 1, 10, 1);
        pHin = trex(18);
        spec.pHin(idx_a(j)) = pHin;
        k0_val = calc_k0_ext(d.vrse(idx_b(j)), d.temperature(idx_b(j)), d.salinity(idx_b(j)), pHin, k2s(g));
        new_row = {gfets(g), k0_val, spec.datetime_sample(idx_a(j))};
        k0_output = [k0_output; new_row];
    end

    T = k0_output(ismember(k0_output.gliderfet_id, gfets(g)),:);
    k0 = mean(T.k0_seawater);
    disp(k0)
    d.pHin = calc_pHext(d.vrse, d.temperature, d.salinity, k0, k2s(g));
    plot(d.datetime,d.pHin)
    hold on;
    scatter(spec.datetime_sample,spec.pHin,'+')
end

disp(k0_output)

% d10 = data(ismember(data.gliderfet_id,"10"),:);
% d10.pHin = calc_pHext(d10.vrse, d10.temperature, d10.salinity, k0, k2s(g));
% d16 = data(ismember(data.gliderfet_id,"16"),:);

% writetable(k0_output, 'k0_output.csv')  % save to CSV