%% spline_orientation_multi_XZ.m
% Plot tangent orientation (deg) vs arc length for multiple ParaView CSV files.
% Curve is assumed to lie in the XZ plane. Only the orientation figure is produced.

clear; clc; close all;

%% ---- User settings ----
folder   = ".";                 % folder containing the CSVs
pattern  = "*.csv";             % or e.g. "spline_*.csv"
use3DArcLength = true;          % true: s from XYZ; false: s from XZ only

smoothMethod  = "sgolay";       % "sgolay" or "spline"
sgolayWindow  = 5;             % odd; will auto-shrink if too big
splineP       = 0.999;          % used only if smoothMethod="spline"

% Resample each curve to a common s-grid so curves are comparable
resampleToCommonGrid = true;
NsCommon = 500;                 % number of samples on common grid

%% ---- Find files ----
files = dir(fullfile(folder, pattern));
if isempty(files)
    error("No files matched %s in folder %s", pattern, folder);
end

%% ---- Prepare plot ----
figure('Name','Spline tangent orientation (XZ)');
hold on; grid on;
xlabel('Arc length, s');
ylabel('\theta(s) [deg]');
title('Tangent orientation relative to +X (XZ plane)');

allResults = struct([]);

%% ---- Loop over CSVs ----
for k = 1:numel(files)
    fpath = fullfile(files(k).folder, files(k).name);

    % Read (sanitised names first)
    T = readtable(fpath);
    vn = string(T.Properties.VariableNames);

    if any(vn == "Points_0")
        x = T.Points_0; y = T.Points_1; z = T.Points_2;
    elseif any(vn == "Points0") && any(vn == "Points1") && any(vn == "Points2")
        x = T.Points0;  y = T.Points1;  z = T.Points2;
    else
        % fall back to preserved names
        Tp = readtable(fpath, "VariableNamingRule", "preserve");
        vnp = string(Tp.Properties.VariableNames);
        if any(vnp == "Points:0")
            x = Tp.("Points:0"); y = Tp.("Points:1"); z = Tp.("Points:2");
        else
            error("Could not find point columns in %s", files(k).name);
        end
    end

    x = x(:); y = y(:); z = z(:);

    % Arc length s
    if use3DArcLength
        ds = sqrt(diff(x).^2 + diff(y).^2 + diff(z).^2);
    else
        ds = sqrt(diff(x).^2 + diff(z).^2);
    end
    s = [0; cumsum(ds)];

    % Remove duplicate s
    [s, ia] = unique(s, 'stable');
    x = x(ia); z = z(ia);

    if numel(s) < 5
        warning("Skipping %s (not enough unique points).", files(k).name);
        continue;
    end

    % Smooth + tangent
    switch lower(smoothMethod)
        case "sgolay"
            w = sgolayWindow;
            if w >= numel(s)
                w = max(5, 2*floor((numel(s)-1)/2)-1); % largest odd < N
            end
            if mod(w,2) == 0, w = w+1; end

            xS = smoothdata(x, 'sgolay', w);
            zS = smoothdata(z, 'sgolay', w);

            dxds = gradient(xS, s);
            dzds = gradient(zS, s);

        case "spline"
            if ~exist('csaps','file')
                error("csaps not found. Use smoothMethod='sgolay' or install Spline Toolbox.");
            end
            ppx = csaps(s, x, splineP);
            ppz = csaps(s, z, splineP);

            dxds = fnval(fnder(ppx,1), s);
            dzds = fnval(fnder(ppz,1), s);

        otherwise
            error("Unknown smoothMethod: %s", smoothMethod);
    end

    theta_deg = rad2deg(unwrap(atan2(dzds, dxds)));

    allResults(k).file = files(k).name;
    allResults(k).s = s;
    allResults(k).theta_deg = theta_deg;

    % Plot
    if resampleToCommonGrid
        sCommon = linspace(s(1), s(end), NsCommon).';
        thetaCommon = interp1(s, theta_deg, sCommon, 'pchip');
        plot(sCommon, thetaCommon, 'LineWidth', 1.2, 'DisplayName', files(k).name);
        allResults(k).sCommon = sCommon;
        allResults(k).thetaCommon = thetaCommon;
    else
        plot(s, theta_deg, 'LineWidth', 1.2, 'DisplayName', files(k).name);
    end
end

legend('Interpreter','none','Location','best');
hold off;