function ResponseSpectrumGUI()
% =========================================================================
%  ResponseSpectrumGUI — Earthquake Response Spectrum Calculator
%  Linear Single-Degree-of-Freedom (SDOF) Dynamic Analysis
% =========================================================================
%
%  Author:  doc. Ing. Özgür Yurdakul, Ph.D.
%           Faculty of Transport Engineering
%           University of Pardubice
%
%  Version: 3.2 (April 2026)
%
%  License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0
%           International (CC BY-NC-SA 4.0)
%           https://creativecommons.org/licenses/by-nc-sa/4.0/
%
%  Acknowledgement:
%     This tool was developed within the framework of a project
%     co-financed from the state budget by the Technology Agency
%     of the Czech Republic under the SIGMA Programme within the
%     M.ERA-NET Cofund.
%
%  Description:
%     GUI-based tool for computing pseudo-acceleration (PSa), pseudo-
%     velocity (PSv), and spectral displacement (Sd) response spectra
%     from earthquake ground motion records. Supports single (1D) and
%     two-component (2D) analysis with geometric mean computation.
%
%  File Loading:
%     Multiple files can be loaded. Each file opens a per-file dialog:
%       - Time step dt (s)  — always required by user, never auto-filled
%       - Acceleration units (cm/s² | m/s² | g)
%       - Component label  (Single / E-W / N-S)
%       - Column selection (multi-column files only)
%       - Mode 1D / 2D     (multi-column files only)
%
%  Compute Rules:
%       - 1 file in table        → used automatically (no checkbox needed)
%       - 1 file checked (1D)    → 1D spectrum
%       - 1 file checked (2D)    → 2D from same multi-column file
%       - 2 files checked (EW+NS)→ 2D from paired single-component files
%       - >2 files checked       → error
%
%  Requirements: MATLAB R2020b or later (uifigure / App Designer)
%
% =========================================================================

close all force;

%% ── Layout ───────────────────────────────────────────────────────────────
PW  = 272;   % control panel width
PH  = 632;   % control panel height — aligns top with Sa axes top
GAP = 6;
CW  = 322;   % axes column width

XC = [PW+2*GAP,  PW+2*GAP+CW+GAP,  PW+2*GAP+2*(CW+GAP)];

AY_rec=8;   AH_rec=94;
AY_sd =110; AH_sd =99;
AY_sv =217; AH_sv =99;
AY_sa =324; AH_sa =292;

FH    = 650;
FW_1D = PW + 2*GAP + CW + GAP;
FW_2D = PW + 2*GAP + 3*CW + 3*GAP;

COL_TITLES = {'E-W','N-S','Geometric Mean'};
C_EW  = [0 0 0.7];
C_NS  = [0.8 0 0];
C_GM  = [0.5 0 0.5];
OV_COLS  = {[0 0 0.7],[0.8 0 0],[0 0.5 0]};
OV_NAMES = {'Newmark-β','Nigam-Jennings','Duhamel'};

%% ── Figure ───────────────────────────────────────────────────────────────
fig = uifigure('Name','Response Spectrum Calculator  v3.2', ...
    'Position',[40 40 FW_1D FH],'Color',[0.95 0.95 0.97]);

%% ── Control panel ────────────────────────────────────────────────────────
cp = uipanel(fig,'Title','Ground Motion Files','TitlePosition','centertop', ...
    'Position',[GAP GAP PW PH], ...
    'BackgroundColor',[0.91 0.91 0.95],'FontWeight','bold','FontSize',10);

%% ── Layout: top-down (y measured from bottom of panel) ──────────────────
%%
%%  y=688  "Loaded Files" label
%%  y=426  file table  (h=260)
%%  y=410  "Check file(s) below…" hint
%%  y=382  [Add File(s)...]  button
%%  y=354  [Remove Checked]  button
%%  ───────────────────────────────────────────────────
%%  y=328  "Spectral period range" label
%%  y=306  T_min + T_max row
%%  y=283  N points row
%%  y=259  "Damping ratio ξ" label
%%  y=236  damping field
%%  y=209  "Integration method" label
%%  y=185  method dropdown
%%  y=164  Overlay checkbox
%%  y=144  Same-lim checkbox  (hidden, 2D only)
%%  ───────────────────────────────────────────────────
%%  y=118  [COMPUTE SPECTRUM]
%%  y= 90  [Export CSV]  [Save All Figures]
%%  y=  4  Status label  (h=82)

%% ── File table section ───────────────────────────────────────────────────
uilabel(cp,'Text','Loaded Files  (check box to select)', ...
    'Position',[8 596 254 14],'FontSize',9,'FontWeight','bold', ...
    'FontColor',[0.25 0.25 0.25]);

%  Checkbox in col 1 (logical), then file info columns
tbl_files = uitable(cp,'Position',[8 414 254 180], ...
    'ColumnName',    {'','File','dt','Units','Comp.','Cols'}, ...
    'ColumnWidth',   {20, 88, 38, 44, 36, 24}, ...
    'ColumnEditable',[true false false false false false], ...
    'FontSize',9,'RowStriping','on');

uilabel(cp,'Text','1 file → used automatically.  2 files → must be E-W + N-S.', ...
    'Position',[8 392 254 14],'FontSize',8,'FontColor',[0.45 0.45 0.45], ...
    'FontAngle','italic');

uibutton(cp,'Text','➕   Add File(s)...','Position',[8 364 254 24], ...
    'BackgroundColor',[0.15 0.52 0.22],'FontColor','white','FontWeight','bold', ...
    'ButtonPushedFcn',@(~,~)cb_addfiles());

uibutton(cp,'Text','✖   Remove Checked','Position',[8 336 254 24], ...
    'BackgroundColor',[0.72 0.15 0.12],'FontColor','white','FontWeight','bold', ...
    'ButtonPushedFcn',@(~,~)cb_removefile());

%% ── Analysis parameters ──────────────────────────────────────────────────
uilabel(cp,'Text','Spectral period range  (s)', ...
    'Position',[8 308 200 14],'FontWeight','bold');
uilabel(cp,'Text','T_min','Position',[8  287 38 18]);
ef_Tmin = uieditfield(cp,'numeric','Position',[50  287 62 20], ...
    'Value',0.01,'Limits',[1e-4 100]);
uilabel(cp,'Text','T_max','Position',[126 287 38 18]);
ef_Tmax = uieditfield(cp,'numeric','Position',[168 287 82 20], ...
    'Value',4.0,'Limits',[0.01 500]);

% N points and damping on the same row
uilabel(cp,'Text','N pts','Position',[8 262 36 20]);
ef_Np = uieditfield(cp,'numeric','Position',[48 262 70 20], ...
    'Value',2000,'Limits',[10 10000],'RoundFractionalValues','on');
uilabel(cp,'Text','ξ','Position',[132 262 14 20]);
ef_xi = uieditfield(cp,'numeric','Position',[150 262 68 20], ...
    'Value',0.05,'Limits',[0 0.99]);

uilabel(cp,'Text','Integration method','Position',[8 238 200 14],'FontWeight','bold');
dd_meth = uidropdown(cp, ...
    'Items',{'Newmark-β  (Avg. Accel.)','Nigam-Jennings  (Exact)','Duhamel''s Integral'}, ...
    'Position',[8 216 254 22],'Value','Newmark-β  (Avg. Accel.)');

cb_all = uicheckbox(cp,'Text','Overlay all three methods', ...
    'Position',[8 194 254 16],'Value',false);

cb_samelim = uicheckbox(cp,'Text','Use same axis limits (2D)', ...
    'Position',[8 174 254 16],'Value',false,'Visible','off', ...
    'ValueChangedFcn',@(~,~)cb_apply_samelim());

%% ── Action buttons ───────────────────────────────────────────────────────
uibutton(cp,'Text','▶  COMPUTE SPECTRUM','Position',[8 148 254 24], ...
    'FontWeight','bold','FontSize',12, ...
    'BackgroundColor',[0.14 0.37 0.72],'FontColor','white', ...
    'ButtonPushedFcn',@(~,~)cb_compute());

uibutton(cp,'Text','Export CSV','Position',[8 120 122 24], ...
    'ButtonPushedFcn',@(~,~)cb_export());
uibutton(cp,'Text','Save All Figures','Position',[138 120 116 24], ...
    'ButtonPushedFcn',@(~,~)cb_saveallfigs());

%% ── Status label ─────────────────────────────────────────────────────────
lbl_status = uilabel(cp,'Text','Ready.  Add files to begin.', ...
    'Position',[8 4 254 82],'WordWrap','on','FontSize',8.5, ...
    'FontColor',[0.05 0.45 0.05],'VerticalAlignment','top');

%% ── Axes  (3 cols × 4 rows) ──────────────────────────────────────────────
ax_sa = gobjects(1,3);  ax_sv = gobjects(1,3);
ax_sd = gobjects(1,3);  ax_ag = gobjects(1,3);

for c = 1:3
    ax_sa(c) = uiaxes(fig,'Position',[XC(c) AY_sa CW AH_sa]);
    ax_sv(c) = uiaxes(fig,'Position',[XC(c) AY_sv CW AH_sv]);
    ax_sd(c) = uiaxes(fig,'Position',[XC(c) AY_sd CW AH_sd]);
    ax_ag(c) = uiaxes(fig,'Position',[XC(c) AY_rec CW AH_rec]);
    for ax_ = [ax_sa(c) ax_sv(c) ax_sd(c) ax_ag(c)]
        grid(ax_,'on');  box(ax_,'on');
    end
    title(ax_sa(c),COL_TITLES{c},'FontSize',10,'FontWeight','bold');
    ylabel(ax_sa(c),'PSa (g)');    xlabel(ax_sa(c),'Period T (s)');
    ylabel(ax_sv(c),'PSv (cm/s)'); xlabel(ax_sv(c),'Period T (s)');
    ylabel(ax_sd(c),'Sd (cm)');    xlabel(ax_sd(c),'Period T (s)');
    ylabel(ax_ag(c),'a_g');        xlabel(ax_ag(c),'Time t (s)');
    if c > 1, set_col_visible(c,'off'); end
end
set(ax_ag(3),'Visible','off');

%% ── Application state ────────────────────────────────────────────────────
% File struct fields:
%   .name .fullpath .rawdata .ncols
%   .dt   .units
%   .comp   : 'Single' | 'EW' | 'NS' | '2D'
%   .mode   : '1D' | '2D'   (2D = two cols from same file)
%   .col1   : column index for 1D or EW (mode=1D)
%   .colEW  : EW column (mode=2D)
%   .colNS  : NS column (mode=2D)
ST.FILES = {};
ST.T     = [];
ST.Sa    = {};  ST.Sv = {};  ST.Sd = {};

%% ═══════════════════════════════════════════════════════════════════════
%%  CALLBACKS
%% ═══════════════════════════════════════════════════════════════════════

    %% ── Add files ────────────────────────────────────────────────────────
    function cb_addfiles()
        [fnames, fpath] = uigetfile( ...
            {'*.asc;*.txt;*.csv;*.at2;*.dat','Ground Motion Files';'*.*','All'}, ...
            'Select Ground Motion File(s)','MultiSelect','on');
        if isequal(fnames,0), return; end
        if ischar(fnames), fnames = {fnames}; end

        nAdded = 0;
        for k = 1:numel(fnames)
            fullp = fullfile(fpath, fnames{k});
            try
                [rawdata, ncols] = read_raw_numeric(fullp);
            catch ME
                uialert(fig, sprintf('Cannot read  %s\n\n%s',fnames{k},ME.message),'Read Error');
                continue;
            end

            fstruct = fileAddDialog(fig, fnames{k}, rawdata, ncols);
            if isempty(fstruct), continue; end

            fstruct.name    = fnames{k};
            fstruct.fullpath = fullp;
            fstruct.rawdata = rawdata;
            fstruct.ncols   = ncols;
            ST.FILES{end+1} = fstruct;
            nAdded = nAdded + 1;
        end

        if nAdded > 0
            refreshFileTable();
            setStatus(sprintf('%d file(s) loaded.', numel(ST.FILES)),'ok');
        end
    end

    %% ── Remove checked files ─────────────────────────────────────────────
    function cb_removefile()
        if isempty(ST.FILES)
            uialert(fig,'No files loaded.','Remove');  return;
        end
        data = tbl_files.Data;
        if isempty(data)
            uialert(fig,'No files to remove.','Remove');  return;
        end

        % Find checked rows
        rows = [];
        for i = 1:size(data,1)
            if isequal(data{i,1}, true)
                rows(end+1) = i; %#ok<AGROW>
            end
        end

        if isempty(rows)
            % If nothing checked and only 1 file, remove it
            if numel(ST.FILES) == 1
                rows = 1;
            else
                uialert(fig,'Check at least one file to remove.','Remove');
                return;
            end
        end

        % Remove descending to keep indices valid
        for r = fliplr(rows)
            if r <= numel(ST.FILES)
                ST.FILES(r) = [];
            end
        end

        refreshFileTable();
        setStatus(sprintf('%d file(s) remaining.', numel(ST.FILES)),'ok');
    end

    %% ── Compute ──────────────────────────────────────────────────────────
    function cb_compute()
        if isempty(ST.FILES)
            setStatus('Add files first.','warn');  return;
        end

        n = numel(ST.FILES);

        % Determine which files to use
        if n == 1
            % Single file: use automatically regardless of checkbox
            useIdx = 1;
        else
            % Multiple files: read checkboxes
            data = tbl_files.Data;
            useIdx = [];
            for i = 1:size(data,1)
                if isequal(data{i,1}, true)
                    useIdx(end+1) = i; %#ok<AGROW>
                end
            end
            if isempty(useIdx)
                setStatus('Check at least one file in the table.','warn');  return;
            end
            if numel(useIdx) > 2
                uialert(fig, ...
                    sprintf(['Select at most 2 files.\n' ...
                    'For 1D: check 1 file.\n' ...
                    'For 2D: check 1 E-W file and 1 N-S file.']), ...
                    'Too Many Files Selected');
                return;
            end
        end

        xi   = ef_xi.Value;
        Tmin = ef_Tmin.Value;
        Tmax = ef_Tmax.Value;
        Np   = round(ef_Np.Value);
        meth = dd_meth.Value;
        doAll= cb_all.Value;
        T_vec= logspace(log10(Tmin), log10(Tmax), Np);
        ST.T = T_vec;

        setStatus('Computing...','warn');  drawnow;

        try
            t0 = tic;

            if numel(useIdx) == 1
                f = ST.FILES{useIdx};

                if strcmp(f.mode,'2D')
                    %% ── 2D from single multi-column file ─────────────────
                    agEW = convert_units(f.rawdata(:,f.colEW), f.units);
                    agNS = convert_units(f.rawdata(:,f.colNS), f.units);
                    doCompute2D(agEW, agNS, f.dt, T_vec, xi, meth, doAll, t0);

                else
                    %% ── 1D ──────────────────────────────────────────────
                    ag = convert_units(f.rawdata(:,f.col1), f.units);

                    fig.Position(3) = FW_1D;
                    set_col_visible(2,'off');  set_col_visible(3,'off');
                    cb_samelim.Visible='off';  cb_samelim.Value=false;

                    [Sd,Sv,Sa] = run_method(ag, f.dt, T_vec, xi, meth, doAll, ...
                        OV_COLS, OV_NAMES, ax_sa(1), ax_sv(1), ax_sd(1), C_EW);
                    ST.Sa={Sa};  ST.Sv={Sv};  ST.Sd={Sd};

                    tvec = (0:numel(ag)-1)' * f.dt;
                    cla(ax_ag(1));
                    plot(ax_ag(1),tvec,ag/9.80665,'Color',[0.2 0.2 0.2],'LineWidth',0.7);
                    ylabel(ax_ag(1),'a_g');  xlabel(ax_ag(1),'Time t (s)');
                    grid(ax_ag(1),'on');
                    title(ax_sa(1), sprintf('[%s]  %s', f.comp, f.name), ...
                        'FontSize',8,'FontWeight','bold','Interpreter','none');

                    elapsed = toc(t0);
                    [samax,idx] = max(Sa);
                    setStatus(sprintf('Done %.2fs  ξ=%.0f%%\nmax PSa=%.4fg @ T=%.3fs', ...
                        elapsed, xi*100, samax/9.80665, T_vec(idx)),'ok');
                end

            else
                %% ── 2 files checked: must be one EW + one NS ─────────────
                f1 = ST.FILES{useIdx(1)};
                f2 = ST.FILES{useIdx(2)};

                comps = {f1.comp, f2.comp};

                if ~(ismember('EW',comps) && ismember('NS',comps))
                    uialert(fig, sprintf([ ...
                        'For 2D analysis, select one E-W and one N-S file.\n\n' ...
                        'Current selection:\n' ...
                        '  File 1: %s  (%s)\n' ...
                        '  File 2: %s  (%s)\n\n' ...
                        'Note: U-D can only be used for 1D analysis.'], ...
                        f1.name, f1.comp, f2.name, f2.comp), ...
                        '2D Requires E-W + N-S');
                    return;
                end

                if strcmp(f1.comp,'EW'), fEW=f1; fNS=f2;
                else,                    fEW=f2; fNS=f1; end

                agEW = convert_units(fEW.rawdata(:,fEW.col1), fEW.units);
                agNS = convert_units(fNS.rawdata(:,fNS.col1), fNS.units);

                if fEW.dt ~= fNS.dt
                    uialert(fig, sprintf( ...
                        'dt mismatch!\n  E-W dt = %.5f s\n  N-S dt = %.5f s\nUsing E-W dt.', ...
                        fEW.dt, fNS.dt),'dt Mismatch');
                end

                doCompute2D(agEW, agNS, fEW.dt, T_vec, xi, meth, doAll, t0);
            end

        catch ME
            setStatus(['Error: ' ME.message],'err');
        end
    end

    function doCompute2D(agEW, agNS, dt, T_vec, xi, meth, doAll, t0)
        fig.Position(3) = FW_2D;
        set_col_visible(2,'on');  set_col_visible(3,'on');
        set(ax_ag(3),'Visible','off');
        cb_samelim.Visible = 'on';

        [SdEW,SvEW,SaEW] = run_method(agEW, dt, T_vec, xi, meth, doAll, ...
            OV_COLS, OV_NAMES, ax_sa(1), ax_sv(1), ax_sd(1), C_EW);
        ST.Sa{1}=SaEW;  ST.Sv{1}=SvEW;  ST.Sd{1}=SdEW;
        tvec = (0:numel(agEW)-1)' * dt;
        cla(ax_ag(1));
        plot(ax_ag(1),tvec,agEW/9.80665,'Color',C_EW,'LineWidth',0.7);
        ylabel(ax_ag(1),'a_g');  xlabel(ax_ag(1),'Time t (s)');  grid(ax_ag(1),'on');
        title(ax_sa(1),'E-W','FontSize',10,'FontWeight','bold');

        [SdNS,SvNS,SaNS] = run_method(agNS, dt, T_vec, xi, meth, doAll, ...
            OV_COLS, OV_NAMES, ax_sa(2), ax_sv(2), ax_sd(2), C_NS);
        ST.Sa{2}=SaNS;  ST.Sv{2}=SvNS;  ST.Sd{2}=SdNS;
        tvec2 = (0:numel(agNS)-1)' * dt;
        cla(ax_ag(2));
        plot(ax_ag(2),tvec2,agNS/9.80665,'Color',C_NS,'LineWidth',0.7);
        ylabel(ax_ag(2),'a_g');  xlabel(ax_ag(2),'Time t (s)');  grid(ax_ag(2),'on');
        title(ax_sa(2),'N-S','FontSize',10,'FontWeight','bold');

        Sa_gm = sqrt(SaEW.*SaNS);
        Sv_gm = sqrt(SvEW.*SvNS);
        Sd_gm = sqrt(SdEW.*SdNS);
        ST.Sa{3}=Sa_gm;  ST.Sv{3}=Sv_gm;  ST.Sd{3}=Sd_gm;
        plot_single(ax_sa(3),ax_sv(3),ax_sd(3),T_vec,Sa_gm,Sv_gm,Sd_gm,C_GM,xi);
        title(ax_sa(3),'Geometric Mean','FontSize',10,'FontWeight','bold');

        elapsed = toc(t0);
        [saEW,iEWp]=max(SaEW); [saNS,iNSp]=max(SaNS); [saGM,iGMp]=max(Sa_gm);
        setStatus(sprintf('Done %.2fs  ξ=%.0f%%\nE-W: %.4fg @ T=%.3fs\nN-S: %.4fg @ T=%.3fs\nGM:  %.4fg @ T=%.3fs', ...
            elapsed, xi*100, ...
            saEW/9.80665, T_vec(iEWp), ...
            saNS/9.80665, T_vec(iNSp), ...
            saGM/9.80665, T_vec(iGMp)),'ok');

        if cb_samelim.Value, cb_apply_samelim(); end
    end

    %% ── Export CSV ───────────────────────────────────────────────────────
    function cb_export()
        if isempty(ST.T)
            setStatus('Compute spectrum first.','warn');  return;
        end
        [fn,fp] = uiputfile('spectrum.csv','Save Spectrum as CSV');
        if isequal(fn,0), return; end
        if numel(ST.Sa) >= 3
            tbl2 = table(ST.T(:), ...
                ST.Sa{1}(:)/9.80665, ST.Sv{1}(:)*100, ST.Sd{1}(:)*100, ...
                ST.Sa{2}(:)/9.80665, ST.Sv{2}(:)*100, ST.Sd{2}(:)*100, ...
                ST.Sa{3}(:)/9.80665, ST.Sv{3}(:)*100, ST.Sd{3}(:)*100, ...
                'VariableNames',{'T_s', ...
                'PSa_EW_g','PSv_EW_cms','Sd_EW_cm', ...
                'PSa_NS_g','PSv_NS_cms','Sd_NS_cm', ...
                'PSa_GM_g','PSv_GM_cms','Sd_GM_cm'});
        else
            tbl2 = table(ST.T(:), ...
                ST.Sa{1}(:)/9.80665, ST.Sv{1}(:)*100, ST.Sd{1}(:)*100, ...
                'VariableNames',{'T_s','PSa_g','PSv_cms','Sd_cm'});
        end
        writetable(tbl2, fullfile(fp,fn));
        setStatus(['Saved: ' fn],'ok');
    end

    %% ── Save figures ─────────────────────────────────────────────────────
    function cb_saveallfigs()
        if isempty(ST.T)
            setStatus('Compute spectrum first.','warn');  return;
        end
        [fn,fp] = uiputfile( ...
            {'*.png','PNG (*.png)';'*.pdf','PDF (*.pdf)'},'Save All Figures');
        if isequal(fn,0), return; end
        nc    = numel(ST.Sa);
        clrs  = {C_EW, C_NS, C_GM};
        hf    = figure('Color','w','Position',[60 60 420*nc 560],'Visible','off');
        ylabs   = {'PSa (g)','PSv (cm/s)','Sd (cm)'};
        yscales = {1/9.80665, 100, 100};
        ctitles = {'E-W','N-S','Geometric Mean'};
        for r = 1:3
            for c = 1:nc
                subplot(4,nc,(r-1)*nc+c);
                y = ST.Sa{c};
                if r==2, y=ST.Sv{c}; elseif r==3, y=ST.Sd{c}; end
                col = clrs{c};  if nc==1, col=[0.2 0.2 0.8]; end
                plot(ST.T, y*yscales{r}, 'Color',col, 'LineWidth',1.3);
                ylabel(ylabs{r});  grid on;
                if r==1, title(ctitles{min(c,end)}); end
                if r==3, xlabel('Period T (s)'); end
            end
        end
        exportgraphics(hf, fullfile(fp,fn),'Resolution',300);
        close(hf);
        setStatus(['Saved: ' fn],'ok');
    end

    %% ── Same axis limits ─────────────────────────────────────────────────
    function cb_apply_samelim()
        if ~cb_samelim.Value
            for c = 1:3
                ylim(ax_sa(c),'auto'); xlim(ax_sa(c),'auto');
                ylim(ax_sv(c),'auto'); xlim(ax_sv(c),'auto');
                ylim(ax_sd(c),'auto'); xlim(ax_sd(c),'auto');
            end
            return;
        end
        if isempty(ST.T), return; end
        n = numel(ST.Sa);
        lim_sa = nice_ceil(max(cellfun(@(x)max(x)/9.80665, ST.Sa(1:n))));
        lim_sv = nice_ceil(max(cellfun(@(x)max(x)*100,     ST.Sv(1:n))));
        lim_sd = nice_ceil(max(cellfun(@(x)max(x)*100,     ST.Sd(1:n))));
        for c = 1:n
            ylim(ax_sa(c),[0 lim_sa]); xlim(ax_sa(c),[ST.T(1) ST.T(end)]);
            ylim(ax_sv(c),[0 lim_sv]); xlim(ax_sv(c),[ST.T(1) ST.T(end)]);
            ylim(ax_sd(c),[0 lim_sd]); xlim(ax_sd(c),[ST.T(1) ST.T(end)]);
        end
    end

    %% ── Helpers ──────────────────────────────────────────────────────────
    function refreshFileTable()
        n    = numel(ST.FILES);
        data = cell(n, 6);
        for i = 1:n
            f = ST.FILES{i};
            if n == 1
                data{i,1} = true;   % auto-check if only one file
            else
                % Preserve existing checkbox state if possible
                existing = tbl_files.Data;
                if ~isempty(existing) && i <= size(existing,1) && ~isempty(existing{i,1})
                    data{i,1} = existing{i,1};
                else
                    data{i,1} = false;
                end
            end
            data{i,2} = f.name;
            data{i,3} = f.dt;
            data{i,4} = f.units;
            data{i,5} = f.comp;
            data{i,6} = f.ncols;
        end
        tbl_files.Data = data;
    end

    function set_col_visible(c, vis)
        set(ax_sa(c),'Visible',vis);  set(ax_sv(c),'Visible',vis);
        set(ax_sd(c),'Visible',vis);  set(ax_ag(c),'Visible',vis);
    end

    function setStatus(msg, type)
        lbl_status.Text = msg;
        switch type
            case 'ok',   lbl_status.FontColor = [0.05 0.45 0.05];
            case 'warn', lbl_status.FontColor = [0.65 0.45 0.00];
            case 'err',  lbl_status.FontColor = [0.75 0.05 0.05];
        end
    end

end   % ← end ResponseSpectrumGUI


%% ═══════════════════════════════════════════════════════════════════════
%%  PER-FILE SETTINGS DIALOG
%%  Returns struct on OK, [] on Cancel.
%%
%%  Struct fields: dt, units, comp, mode, col1, colEW, colNS
%%
%%  comp values:  'Single' | 'EW' | 'NS' | '2D'
%%  mode values:  '1D' | '2D'
%% ═══════════════════════════════════════════════════════════════════════
function fstruct = fileAddDialog(parentFig, fname, rawdata, ncols)

fstruct  = [];
nPrev    = min(size(rawdata,1), 8);
dlgH     = 380 + 72*(ncols >= 2);    % taller for multi-column

dlg = uifigure('Name',sprintf('File Settings — %s', fname), ...
    'Position',[0 0 480 dlgH], ...
    'WindowStyle','modal','Resize','off');
movegui(dlg,'center');

%% File info
uilabel(dlg,'Text',sprintf('File:  %s', fname), ...
    'Position',[14 dlgH-28 452 18],'FontWeight','bold','FontSize',11);
uilabel(dlg,'Text',sprintf('%d column(s) of numeric data in file', ncols), ...
    'Position',[14 dlgH-48 452 16],'FontSize',10,'FontColor',[0.38 0.38 0.38]);

%% Preview table
colNames = arrayfun(@(k)sprintf('Col %d',k), 1:ncols, 'UniformOutput',false);
prevH    = 150;
uitable(dlg,'Position',[14 dlgH-48-4-prevH 452 prevH], ...
    'Data',num2cell(rawdata(1:nPrev,:)), ...
    'ColumnName',colNames,'RowName',{}, ...
    'ColumnEditable',false(1,ncols),'FontSize',9);

%% dt
baseY = dlgH - 48 - 4 - prevH - 28;
uilabel(dlg,'Text','Time step  dt  (s)  — required', ...
    'Position',[14 baseY 230 16],'FontWeight','bold','FontSize',10);
ef_dt_dlg = uieditfield(dlg,'numeric', ...
    'Position',[250 baseY-1 100 22],'Value',0,'Limits',[0 Inf]);
uilabel(dlg,'Text','(0 = not set)', ...
    'Position',[358 baseY 110 16],'FontSize',8,'FontColor',[0.55 0.55 0.55]);

%% Units
baseY = baseY - 30;
uilabel(dlg,'Text','Acceleration units', ...
    'Position',[14 baseY 180 16],'FontWeight','bold','FontSize',10);
dd_units_dlg = uidropdown(dlg,'Items',{'cm/s²','m/s²','g'}, ...
    'Position',[200 baseY-1 120 22],'Value','cm/s²');

%% ── Single-column file: just component ──────────────────────────────────
dd_mode_dlg = [];
dd_comp_dlg = [];
dd_col1     = [];  lbl_col1  = [];
dd_colEW    = [];  lbl_colEW = [];
dd_colNS    = [];  lbl_colNS = [];
lbl_comp    = [];

if ncols == 1
    % Component only — no mode or column selector needed
    baseY = baseY - 30;
    uilabel(dlg,'Text','Component', ...
        'Position',[14 baseY 90 16],'FontWeight','bold','FontSize',10);
    dd_comp_dlg = uidropdown(dlg, ...
        'Items',{'E-W','N-S','U-D'}, ...
        'Position',[110 baseY-1 120 22],'Value','E-W');
    uilabel(dlg,'Text','(single-column — 1D only)', ...
        'Position',[240 baseY 220 16],'FontSize',9, ...
        'FontColor',[0.45 0.45 0.45],'FontAngle','italic');

else
    %% ── Multi-column: Mode → Component → Column (in that order) ─────────

    % 1) Mode
    baseY = baseY - 30;
    uilabel(dlg,'Text','Mode', ...
        'Position',[14 baseY 60 16],'FontWeight','bold','FontSize',10);
    dd_mode_dlg = uidropdown(dlg,'Items',{'1D','2D (same file)'}, ...
        'Position',[80 baseY-1 130 22],'Value','1D', ...
        'ValueChangedFcn',@onModeChange);

    % 2) Component (shown for 1D; hidden for 2D — not needed when both cols selected)
    baseY = baseY - 30;
    lbl_comp = uilabel(dlg,'Text','Component', ...
        'Position',[14 baseY 90 16],'FontWeight','bold','FontSize',10);
    dd_comp_dlg = uidropdown(dlg, ...
        'Items',{'E-W','N-S','U-D'}, ...
        'Position',[110 baseY-1 120 22],'Value','E-W');

    % 3) Column selection
    baseY = baseY - 30;
    colItems = arrayfun(@(k)sprintf('Col %d',k), 1:ncols, 'UniformOutput',false);

    % 1D: single column picker
    lbl_col1 = uilabel(dlg,'Text','Column', ...
        'Position',[14 baseY 60 16],'FontSize',10);
    dd_col1  = uidropdown(dlg,'Items',colItems, ...
        'Value',colItems{min(2,ncols)}, ...
        'Position',[80 baseY-1 110 22]);

    % 2D: EW + NS column pickers (hidden initially)
    lbl_colEW = uilabel(dlg,'Text','E-W col', ...
        'Position',[14 baseY 62 16],'FontSize',10,'Visible','off');
    dd_colEW  = uidropdown(dlg,'Items',colItems, ...
        'Value',colItems{min(2,ncols)}, ...
        'Position',[82 baseY-1 108 22],'Visible','off');
    lbl_colNS = uilabel(dlg,'Text','N-S col', ...
        'Position',[210 baseY 62 16],'FontSize',10,'Visible','off');
    dd_colNS  = uidropdown(dlg,'Items',colItems, ...
        'Value',colItems{min(3,ncols)}, ...
        'Position',[278 baseY-1 108 22],'Visible','off');
end

%% OK / Cancel
uibutton(dlg,'Text','OK', ...
    'Position',[14 12 150 30], ...
    'BackgroundColor',[0.20 0.45 0.75],'FontColor','white','FontWeight','bold', ...
    'ButtonPushedFcn',@doOK);
uibutton(dlg,'Text','Cancel', ...
    'Position',[174 12 120 30], ...
    'ButtonPushedFcn',@(~,~)uiresume(dlg));

uiwait(dlg);
if isvalid(dlg)
    fstruct = dlg.UserData;
    delete(dlg);
end

    %% Mode change → show/hide component + column selectors
    function onModeChange(~,~)
        if isempty(dd_mode_dlg), return; end
        is2D = strcmp(dd_mode_dlg.Value,'2D (same file)');
        % Column pickers
        if ~isempty(lbl_col1),  lbl_col1.Visible  = onoff(~is2D); end
        if ~isempty(dd_col1),   dd_col1.Visible   = onoff(~is2D); end
        if ~isempty(lbl_colEW), lbl_colEW.Visible = onoff(is2D);  end
        if ~isempty(dd_colEW),  dd_colEW.Visible  = onoff(is2D);  end
        if ~isempty(lbl_colNS), lbl_colNS.Visible = onoff(is2D);  end
        if ~isempty(dd_colNS),  dd_colNS.Visible  = onoff(is2D);  end
        % Component: hide for 2D (not applicable — both EW+NS selected via columns)
        if ~isempty(lbl_comp),    lbl_comp.Visible    = onoff(~is2D); end
        if ~isempty(dd_comp_dlg), dd_comp_dlg.Visible = onoff(~is2D); end
    end

    function s = onoff(b)
        if b, s='on'; else, s='off'; end
    end

    %% OK: validate and store
    function doOK(~,~)
        dt_val = ef_dt_dlg.Value;
        if dt_val <= 0
            uialert(dlg,'Time step dt must be greater than 0.','dt Required');
            return;
        end

        % Map component string to code
        if ~isempty(dd_comp_dlg) && strcmp(dd_comp_dlg.Visible,'on')
            cv = dd_comp_dlg.Value;
            if     contains(cv,'E-W'), comp = 'EW';
            elseif contains(cv,'N-S'), comp = 'NS';
            else,                       comp = 'UD';
            end
        else
            comp = '2D';   % 2D mode: component not used individually
        end

        if ncols == 1
            dlg.UserData = struct('dt',dt_val,'units',dd_units_dlg.Value, ...
                'comp',comp,'mode','1D','col1',1,'colEW',1,'colNS',1);

        else
            modeVal = dd_mode_dlg.Value;
            if strcmp(modeVal,'2D (same file)')
                cEW = str2double(regexp(dd_colEW.Value,'\d+','match','once'));
                cNS = str2double(regexp(dd_colNS.Value,'\d+','match','once'));
                if cEW == cNS
                    uialert(dlg,'E-W and N-S columns must be different.','Column Error');
                    return;
                end
                dlg.UserData = struct('dt',dt_val,'units',dd_units_dlg.Value, ...
                    'comp','2D','mode','2D','col1',cEW,'colEW',cEW,'colNS',cNS);
            else
                c1 = str2double(regexp(dd_col1.Value,'\d+','match','once'));
                dlg.UserData = struct('dt',dt_val,'units',dd_units_dlg.Value, ...
                    'comp',comp,'mode','1D','col1',c1,'colEW',c1,'colNS',c1);
            end
        end
        uiresume(dlg);
    end

end   % fileAddDialog


%% ═══════════════════════════════════════════════════════════════════════
%%  RAW NUMERIC FILE READER
%%  Skips all non-numeric header lines (AFAD .asc, PEER AT2, CSV, etc.)
%%  Returns raw numeric matrix — NO time column added, NO dt detection.
%% ═══════════════════════════════════════════════════════════════════════
function [rawmat, ncols] = read_raw_numeric(filepath)
fid = fopen(filepath,'r');
if fid < 0, error('Cannot open file: %s', filepath); end
lines = {};
while ~feof(fid)
    L = fgetl(fid);
    if ischar(L), lines{end+1} = L; end %#ok<AGROW>
end
fclose(fid);

rows = {};
nc   = 0;
for i = 1:numel(lines)
    L = strtrim(lines{i});
    if isempty(L), continue; end
    if ismember(L(1),{'%','#','!',';','/'}), continue; end
    v = sscanf(L,'%f')';
    if isempty(v), continue; end
    if nc == 0, nc = numel(v); end
    if numel(v) == nc
        rows{end+1} = v; %#ok<AGROW>
    end
end

if isempty(rows)
    error('No numeric data found in file: %s', filepath);
end
rawmat = vertcat(rows{:});
ncols  = nc;
end


%%  UNIT CONVERSION  (all methods work in m/s²)
%% ═══════════════════════════════════════════════════════════════════════
function ag = convert_units(ag_raw, unit)
switch unit
    case 'cm/s²',  ag = ag_raw(:) * 1e-2;
    case 'g',      ag = ag_raw(:) * 9.80665;
    otherwise,     ag = ag_raw(:);          % m/s² — as-is
end
end


%% ═══════════════════════════════════════════════════════════════════════
%%  RUN METHOD  (single or overlay)
%% ═══════════════════════════════════════════════════════════════════════
function [Sd,Sv,Sa] = run_method(ag,dt,T_vec,xi,meth,doAll,OV_COLS,OV_NAMES,ax_sa,ax_sv,ax_sd,col)
METH_NAMES = {'Newmark-β  (Avg. Accel.)','Nigam-Jennings  (Exact)','Duhamel''s Integral'};
METH_FN    = {@nm_spectrum,@nj_spectrum,@dh_spectrum};

if doAll
    res_Sd=cell(1,3); res_Sv=cell(1,3); res_Sa=cell(1,3);
    for m = 1:3
        [res_Sd{m},res_Sv{m},res_Sa{m}] = METH_FN{m}(ag,dt,T_vec,xi);
    end
    cla(ax_sa); cla(ax_sv); cla(ax_sd);
    hold(ax_sa,'on'); hold(ax_sv,'on'); hold(ax_sd,'on');
    for m = 1:3
        plot(ax_sa,T_vec,res_Sa{m}/9.80665,'Color',OV_COLS{m},'LineWidth',1.4,'DisplayName',OV_NAMES{m});
        plot(ax_sv,T_vec,res_Sv{m}*100,    'Color',OV_COLS{m},'LineWidth',1.4);
        plot(ax_sd,T_vec,res_Sd{m}*100,    'Color',OV_COLS{m},'LineWidth',1.4);
    end
    legend(ax_sa,OV_NAMES,'Location','northeast','FontSize',7);
    hold(ax_sa,'off'); hold(ax_sv,'off'); hold(ax_sd,'off');
    Sd=res_Sd{1}; Sv=res_Sv{1}; Sa=res_Sa{1};
else
    k = find(strcmp(meth,METH_NAMES),1);
    if isempty(k), k=1; end
    [Sd,Sv,Sa] = METH_FN{k}(ag,dt,T_vec,xi);
    plot_single(ax_sa,ax_sv,ax_sd,T_vec,Sa,Sv,Sd,col,xi);
end
end


function plot_single(ax_sa,ax_sv,ax_sd,T,Sa,Sv,Sd,col,~)
lw = 1.6;
cla(ax_sa); plot(ax_sa,T,Sa/9.80665,'Color',col,'LineWidth',lw);
xlabel(ax_sa,'Period T (s)'); ylabel(ax_sa,'PSa (g)'); grid(ax_sa,'on');
cla(ax_sv); plot(ax_sv,T,Sv*100,'Color',col,'LineWidth',lw);
xlabel(ax_sv,'Period T (s)'); ylabel(ax_sv,'PSv (cm/s)'); grid(ax_sv,'on');
cla(ax_sd); plot(ax_sd,T,Sd*100,'Color',col,'LineWidth',lw);
xlabel(ax_sd,'Period T (s)'); ylabel(ax_sd,'Sd (cm)'); grid(ax_sd,'on');
end


%% ═══════════════════════════════════════════════════════════════════════
%%  METHOD 1 — NEWMARK-β  (β=1/4, γ=1/2)
%% ═══════════════════════════════════════════════════════════════════════
function [Sd,Sv,Sa] = nm_spectrum(ag,dt,T,xi)
ag=ag(:)'; beta=1/4; gamma=1/2;
a0=1/(beta*dt^2); a1=gamma/(beta*dt); a2=1/(beta*dt);
a3=1/(2*beta)-1;  a4=gamma/beta-1;
N=numel(T); wn=2*pi./T(:)'; c=2*xi*wn; Keff=wn.^2+a0+c*a1;
u=zeros(1,N); v=zeros(1,N); acc=-ag(1)*ones(1,N); umax=zeros(1,N);
for i=1:numel(ag)-1
    Peff=-ag(i+1)+(a0+c*a1).*u+(a2+c*a4).*v+a3*acc;
    u_new=Peff./Keff; acc_new=a0*(u_new-u)-a2*v-a3*acc;
    v=v+dt*((1-gamma)*acc+gamma*acc_new);
    u=u_new; acc=acc_new; umax=max(umax,abs(u));
end
Sd=umax; Sv=wn.*Sd; Sa=wn.^2.*Sd;
end


%% ═══════════════════════════════════════════════════════════════════════
%%  METHOD 2 — NIGAM-JENNINGS (1969)
%% ═══════════════════════════════════════════════════════════════════════
function [Sd,Sv,Sa] = nj_spectrum(ag,dt,T,xi)
ag=ag(:)'; N=numel(T); wn=2*pi./T(:)'; sq=sqrt(1-xi^2); wd=wn*sq;
E=exp(-xi*wn*dt); S=sin(wd*dt); C=cos(wd*dt);
A11=E.*(C+(xi/sq)*S); A12=E.*S./wd; A21=-(E.*wn/sq).*S; A22=E.*(C-(xi/sq)*S);
EC=(A11+A22)/2;
B11= A11./wn.^2+2*xi*(EC-1)./(dt*wn.^3)-A12*(1-2*xi^2)./(dt*wn.^2);
B12=-1  ./wn.^2-2*xi*(EC-1)./(dt*wn.^3)+A12*(1-2*xi^2)./(dt*wn.^2);
B21=(A21+(1-A11)/dt)./wn.^2; B22=(A11-1)./(dt*wn.^2);
u=zeros(1,N); v=zeros(1,N); umax=zeros(1,N);
for i=1:numel(ag)-1
    fn=ag(i); fn1=ag(i+1);
    un=A11.*u+A12.*v+B11*fn+B12*fn1;
    vn=A21.*u+A22.*v+B21*fn+B22*fn1;
    u=un; v=vn; umax=max(umax,abs(u));
end
Sd=umax; Sv=wn.*Sd; Sa=wn.^2.*Sd;
end


%% ═══════════════════════════════════════════════════════════════════════
%%  METHOD 3 — DUHAMEL'S INTEGRAL
%% ═══════════════════════════════════════════════════════════════════════
function [Sd,Sv,Sa] = dh_spectrum(ag,dt,T,xi)
ag=ag(:)'; N=numel(T); wn=2*pi./T(:)'; sq=sqrt(1-xi^2); wd=wn*sq;
E=exp(-xi*wn*dt); S=sin(wd*dt); C=cos(wd*dt);
fac=1./wn.^2; K1=xi*wn.*(1-E.*C)+wd.*E.*S; K2=wd.*(1-E.*C)-xi*wn.*E.*S;
I1=zeros(1,N); I2=zeros(1,N); umax=zeros(1,N);
for i=1:numel(ag)-1
    favg=(ag(i)+ag(i+1))/2;
    I1n=E.*(I1.*C-I2.*S)+fac.*K1*favg;
    I2n=E.*(I1.*S+I2.*C)+fac.*K2*favg;
    I1=I1n; I2=I2n; u=-I2./wd; umax=max(umax,abs(u));
end
Sd=umax; Sv=wn.*Sd; Sa=wn.^2.*Sd;
end


%% ═══════════════════════════════════════════════════════════════════════
%%  HELPERS
%% ═══════════════════════════════════════════════════════════════════════
function v = nice_ceil(x)
if x <= 0, v = 1; return; end
mag = 10^floor(log10(x));
v   = ceil(x / mag * 10) / 10 * mag;
end