close all; clear; % Not needed in production, but useful for debugging.

%%
datapth = "../../../../../Data";

lbsync = {load("models/lbsync.mat")};
lbsync = [lbsync, load(fullfile(datapth, "2024-07-01 XFEL Characterization", "data/lbsync.mat"))];
lbsync = [lbsync, load(fullfile(datapth, "2024-10-10 XFEL Characterization", "data/lbsync.mat"))];
lbsync = [lbsync, load(fullfile(datapth, "2024-12-11 XFEL Characterization", "data/lbsync.mat"))];
lbsync = [lbsync, load(fullfile(datapth, "2024-12-12 XFEL Characterization", "data/lbsync.mat"))];

lbsync = merge_structs(lbsync{:});

%%
export_mdl(lbsync.xfel.mo.Fd, 1e6, "mo_pn");
export_mdl(lbsync.xfel.laser.mlo1.Fd, 1e4, "mlo_pn");
export_mdl(lbsync.xfel.laser.slo1.Fd, 1e4, "slo_pn");
export_mdl(lbsync.xfel.laser.sase3sxp.Fd, 1e4, "sxp_pn");

%%
export_mdl(lbsync.xfel.link.slo1.Fd, 1e5, "lonk_tn");
export_mdl(lbsync.xfel.link.relonk.Fd, 1e5, "relonk_tn");
export_mdl(lbsync.xfel.link.sase3sxp.Fd, 1e5, "sxp_tn");

function export_mdl(sys, floor_start, name)
    [t, t_smooth] = sys.UserData.raw.data{1:2};
    assert(isequal(t.Frequency, t_smooth.Frequency));
    f = t.Frequency;
    mask = f < floor_start;

    t.ReverseIntegratedTimingJitterNoFloor = nan(height(t), 1);
    t.ReverseIntegratedTimingJitterNoFloor(mask) = flip(sqrt(-cumtrapz(flip(f(mask)), flip(db2pow(t.TimingJitter(mask))))));

    t_smooth.ReverseIntegratedTimingJitterNoFloor = nan(height(t_smooth), 1);
    t_smooth.ReverseIntegratedTimingJitterNoFloor(mask) = flip(sqrt(-cumtrapz(flip(f(mask)), flip(db2pow(t_smooth.TimingJitter(mask))))));

    fresp = abs(squeeze(freqresp(sys, f, "Hz"))).^2;
    t_fit = t(:,"Frequency");
    t_fit.TimingJitter = pow2db(fresp);
    t_fit.IntegratedTimingJitter = sqrt(cumtrapz(f, fresp));
    t_fit.ReverseIntegratedTimingJitter = flip(sqrt(-cumtrapz(flip(f), flip(fresp))));
    t_fit.ReverseIntegratedTimingJitterNoFloor = nan(height(t_fit), 1);
    t_fit.ReverseIntegratedTimingJitterNoFloor(mask) = flip(sqrt(-cumtrapz(flip(f(mask)), flip(fresp(mask)))));

    t = reducepoly_table(t, 1e-3);
    t_smooth = reducepoly_table(t_smooth, 1e-3);
    t_fit = reducepoly_table(t_fit, 1e-3);

    writetable(t, "xfel_" + name + ".dat", Delimiter='\t');
    writetable(t_smooth, "xfel_" + name + "_smooth.dat", Delimiter='\t');
    writetable(t_fit, "xfel_" + name + "_fit.dat", Delimiter='\t');
end
