function OUT = pch_transmissibility_T123_ref(pchFile, inputGrid, varargin)
    % PCH_TRANSMISSIBILITY_T123_REF
    % (docstring omitted for brevity; use your full doc again)

    % ---------- parse options ----------
    p = inputParser;
    addParameter(p,'Result','ACCELERATION',@(s)ischar(s)||isstring(s));
    addParameter(p,'Subcase',[],@(x) isempty(x) || (isscalar(x)&&isnumeric(x)));
    addParameter(p,'ForceBaseAxis',[],@(x) isempty(x) || (isscalar(x)&&ismember(x,1:3)));
    addParameter(p,'RelDenTol',1e-14,@(x)isscalar(x)&&x>=0);
    addParameter(p,'MaskDen',true,@(x)islogical(x)||ismember(x,[0 1]));
    addParameter(p,'AssignToBase',true,@(x)islogical(x)||ismember(x,[0 1]));
    addParameter(p,'VarPrefix','',@(s)ischar(s)||isstring(s));
    parse(p, varargin{:});
    wantType      = upper(strtrim(string(p.Results.Result)));
    wantSC        = p.Results.Subcase;
    forceBaseAxis = p.Results.ForceBaseAxis;
    relTol        = p.Results.RelDenTol;
    doMask        = logical(p.Results.MaskDen);
    doAssign      = logical(p.Results.AssignToBase);
    varPrefix     = char(string(p.Results.VarPrefix));

    % ---------- open file ----------
    fid = fopen(pchFile, 'r');
    if fid < 0
        error('PCH:FileOpenFailed', 'Cannot open PCH file: %s', pchFile);
    end
    cleanup = onCleanup(@() fclose(fid));

    % ---------- helpers & init ----------
    num_tokens = @(s) local_num_tokens(s);
    Rows = struct('f', {}, 'sc', {}, 'grid', {}, 'vals', {});
    formsBySC = containers.Map('KeyType','double','ValueType','char');
    gridsSet = [];
    gridMap = containers.Map('KeyType','int32','ValueType','int32');
    inBlock = false;
    curSC = NaN;
    curF = NaN;
    curForm = 'RECT';

    % ---------- parse lines ----------
    while true
        ln = fgetl(fid);
        if ~ischar(ln)
            break;
        end
        u = upper(strtrim(ln));
        if startsWith(u, '$')
            % control block logic (same as before) ...
            if startsWith(u, '$DISPLACEMENT') || startsWith(u, '$DISPLACEMENTS')
                inBlock = any(strcmp(wantType, ["DISPLACEMENT","DISPLACEMENTS"]));
            elseif startsWith(u, '$VELOCITY') || startsWith(u, '$VELOCITIES')
                inBlock = any(strcmp(wantType, ["VELOCITY","VELOCITIES"]));
            elseif startsWith(u, '$ACCELERATION') || startsWith(u, '$ACCELERATIONS')
                inBlock = any(strcmp(wantType, ["ACCELERATION","ACCELERATIONS"]));
            elseif contains(u,'$SPCFORCE') || contains(u,'$ELFORCE') ...
                    || contains(u,'$FORCE') || contains(u,'$STRESS') || contains(u,'$STRAIN')
                inBlock = false;
            end
            if contains(u, 'MAGNITUDE-PHASE OUTPUT') || contains(u, '$PHASE')
                curForm = 'POLAR';
            elseif contains(u, 'REAL OUTPUT') || contains(u, 'RECTANGULAR') ...
                    || contains(u, 'REAL/IMAG')
                curForm = 'RECT';
            end
            if contains(u, '$SUBCASE ID')
                v = regexp(u, '=\s*([0-9D+Ee\-\.+]+)', 'tokens', 'once');
                if ~isempty(v)
                    curSC = local_str2num(v{1});
                end
            end
            if contains(u, '$FREQUENCY')
                v = regexp(u, '=\s*([0-9D+Ee\-\.+]+)', 'tokens', 'once');
                if ~isempty(v)
                    curF = local_str2num(v{1});
                end
            end
            continue;
        end
        if ~inBlock
            continue;
        end
        if ~isempty(wantSC) && ~isnan(curSC) && curSC ~= wantSC
            continue;
        end

        m = regexp(ln, '^\s*(?:GRID\s+)?(\d+)\s+([GS])\s+(.+)$', 'tokens', 'once');
        if isempty(m)
            continue;
        end
        if m{2} == 'S'
            continue;
        end

        grid = int32(str2double(m{1}));
        nums = num_tokens(m{3});

        while numel(nums) < 12
            pos = ftell(fid);
            ln2 = fgetl(fid);
            if ~ischar(ln2)
                break;
            end
            if contains(ln2, '-CONT-')
                k = strfind(ln2, '-CONT-');
                rest = ln2(k + 6 : end);
                nums = [nums, num_tokens(rest)]; %#ok<AGROW>
            else
                fseek(fid, pos, 'bof');
                break;
            end
        end
        if numel(nums) < 6
            continue;
        end

        if strcmp(curForm, 'POLAR')
            if numel(nums) >= 12
                mag = nums(1:6);
                phdeg = nums(7:12);
            else
                mag = nums(1:6);
                phdeg = zeros(1,6);
            end
            vals = mag .* exp(1j * (pi/180) * phdeg);
        else
            if numel(nums) >= 12
                vals = nums(1:6) + 1i * nums(7:12);
            else
                vals = nums(1:6) + 1i * zeros(1,6);
            end
        end

        if ~isKey(gridMap, grid)
            gridMap(grid) = numel(gridsSet) + 1;
            gridsSet(end+1,1) = grid; %#ok<AGROW>
        end
        Rows(end+1) = struct('f', curF, 'sc', curSC, 'grid', grid, 'vals', vals); %#ok<AGROW>
        if ~isnan(curSC) && ~isKey(formsBySC, curSC)
            formsBySC(curSC) = curForm;
        end
    end  % end while

    if isempty(Rows)
        error('PCH:NoDataFound', ...
            'No %s data found in PCH file %s. Check case control or PUNCH format.', ...
            wantType, pchFile);
    end

    % ---------- assemble data ----------
    SCs = unique([Rows.sc]);
    if isempty(SCs) || any(isnan(SCs))
        SCs = 0;
        for i = 1:numel(Rows)
            if isnan(Rows(i).sc)
                Rows(i).sc = 0;
            end
        end
    end
    if ~isempty(wantSC)
        SCs = SCs(SCs == wantSC);
    end
    sc = SCs(1);

    r = Rows([Rows.sc] == sc);
    freqs = unique([r.f]);
    freqs = freqs(~isnan(freqs));
    if isempty(freqs)
        for i = 1:numel(r)
            r(i).f = i;
        end
        freqs = (1 : numel(r)).';
    else
        freqs = sort(freqs(:));
    end

    nF = numel(freqs);
    nG = numel(gridsSet);
    data = complex(zeros(nF, nG, 6));
    for i = 1:numel(r)
        fi = find(freqs == r(i).f, 1, 'first');
        gi = gridMap(r(i).grid);
        data(fi, gi, :) = reshape(r(i).vals, [1 1 6]);
    end

    % ---------- convert to acceleration domain ----------
    f = freqs(:);
    omega = 2 * pi * f;
    switch upper(wantType)
        case {'ACCELERATION', 'ACCELERATIONS'}
            Aacc = data;
        case {'VELOCITY', 'VELOCITIES'}
            Aacc = (1i * omega) .* data;
        case {'DISPLACEMENT', 'DISPLACEMENTS'}
            Aacc = -(omega .^ 2) .* data;
        otherwise
            error('PCH:UnsupportedResult', 'Unsupported Result type: %s', wantType);
    end

    % ---------- locate input grid (denominator) ----------
    g_in = find(gridsSet == inputGrid, 1);
    if isempty(g_in)
        if isempty(gridsSet)
            gridListStr = '(none parsed)';
        else
            countShow = min(numel(gridsSet), 10);
            listed = sprintf('%d ', double(gridsSet(1:countShow)));
            if numel(gridsSet) > countShow
                listed = [listed, '...'];
            end
            gridListStr = listed;
        end
        error('PCH:GridNotFound', ...
            'Input grid %d not found in PCH file %s.\nAvailable grids: %s', ...
            inputGrid, pchFile, gridListStr);
    end

    % ---------- detect base axis ----------
    if ~isempty(forceBaseAxis)
        baseAxis = forceBaseAxis;
    else
        rmsVals = zeros(1,3);
        for d = 1:3
            v = abs(squeeze(Aacc(:, g_in, d)));
            tol = max(1e-16, 1e-8 * max(v));
            v(v < tol) = 0;
            rmsVals(d) = sqrt(mean(v .^ 2));
        end
        [~, baseAxis] = max(rmsVals);
    end
    baseAxisName = sprintf('T%d', baseAxis);

    % ---------- denominator / masking ----------
    den = squeeze(Aacc(:, g_in, baseAxis));
    denMag = abs(den);
    finiteMask = isfinite(denMag);
    if ~any(finiteMask) || max(denMag(finiteMask)) == 0
        warning('PCH:ZeroDenominator', ...
            'Input grid %d, %s has zero accel across the band. Ratios will be NaN.', ...
            inputGrid, baseAxisName);
    else
        if doMask
            dmax = max(denMag(finiteMask));
            denTol = max(1e-16, relTol * dmax);
            badDen = ~isfinite(den) | (denMag <= denTol);
            if nnz(badDen) / numel(den) > 0.95
                badDen = ~isfinite(den) | (denMag == 0);
                denTol = 0;
                warning('PCH:DenMaskRelaxed', ...
                    'Denominator nearly all below threshold; relaxing mask to exact zeros only.');
            end
            den(badDen) = NaN;
        else
            badDen = ~isfinite(den) | (denMag == 0);
            den(badDen) = NaN;
        end
    end

    % ---------- compute transmissibility for other grids ----------
    idxOther = setdiff(1:nG, g_in);
    gridsOther = gridsSet(idxOther);

    OUT = struct();
    OUT.info = struct( ...
        'pch', pchFile, ...
        'subcase', sc, ...
        'resultIn', char(wantType), ...
        'inputGrid', inputGrid, ...
        'grids', gridsSet, ...
        'baseAxis', baseAxis, ...
        'baseAxisName', baseAxisName ...
    );

    for g = reshape(gridsOther, 1, [])
        go = find(gridsSet == g, 1);
        for d = 1:3
            nume = squeeze(Aacc(:, go, d));
            H = nume ./ den;
            OUT.(sprintf('mag_T%d_%d', d, g)) = [f, abs(H)];
            OUT.(sprintf('cpx_T%d_%d', d, g)) = [f, H];
        end
    end

    % ---------- debug info ----------
    OUT.debug.n_freq = nF;
    OUT.debug.n_grids = numel(gridsSet);
    OUT.debug.sample_max = max(abs(Aacc(:)));
    OUT.debug.detected_base_axis = baseAxis;

    % ---------- optional assign to base workspace ----------
    if doAssign
        try
            assignin('base', [varPrefix 'OUT'], OUT);
            assignin('base', [varPrefix 'pch_grids'], double(gridsSet));
            assignin('base', [varPrefix 'pch_freq_hz'], f);
            assignin('base', [varPrefix 'pch_input_grid'], double(inputGrid));
            assignin('base', [varPrefix 'pch_base_axis'], double(baseAxis));
            fn = fieldnames(OUT);
            for k = 1:numel(fn)
                name = fn{k};
                if startsWith(name, 'mag_') || startsWith(name, 'cpx_')
                    assignin('base', [varPrefix name], OUT.(name));
                end
            end
        catch ME
            warning('PCH:AssignBaseFailed', ...
                'assignin to base workspace failed: %s', ME.message);
        end
    end
end

% ===== helper subfunctions =====

function arr = local_num_tokens(s)
    pat = '[-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[DEde][+\-]?\d+)?';
    tok = regexp(s, pat, 'match');
    if isempty(tok)
        arr = [];
        return;
    end
    tok = regexprep(tok, '[dD]', 'E');
    arr = str2double(tok);
end

function x = local_str2num(s)
    x = str2double(regexprep(s, '[dD]', 'E'));
end
