function [freq_hz, grids, EV] = parse_pch_eigenvectors(pchFile, createNamedVars)
% PARSE_PCH_EIGENVECTORS  Extract natural frequencies and eigenvectors from a Nastran .pch
%   [freq_hz, grids, EV] = parse_pch_eigenvectors(pchFile, createNamedVars)
%
% Inputs
%   pchFile         - path to .pch file (SOL103 punch with $EIGENVECTOR sections)
%   createNamedVars - (optional) logical, default true
%                     If true, creates variables in the *caller* workspace named:
%                       EV_grid<GRID>_mode<M> (1x6 row vectors)
%
% Outputs
%   freq_hz  - [nModes x 1] natural frequencies in Hz
%   grids    - [nGrids x 1] GRID IDs (consistent ordering across modes)
%   EV       - struct with fields Mode1, Mode2, ... each [nGrids x 6] eigenvector matrix
%
% Notes
% - Assumes REAL eigenvectors. Each grid’s 6-DOF components appear on two lines:
%     main line:  GRIDID  'G'  T1 T2 T3
%     -CONT-   :             R1 R2 R3
% - Natural frequency computed via f = sqrt(eigenvalue) / (2*pi)
%
% Example
%   [f, grids, EV] = parse_pch_eigenvectors('pl_fix_fix_modes.pch.3', true);

if nargin < 2 || isempty(createNamedVars)
    createNamedVars = true;
end

% Read file as lines
fid = fopen(pchFile, 'r');
if fid < 0
    error('Cannot open file: %s', pchFile);
end
C = textscan(fid, '%s', 'Delimiter', '\n', 'Whitespace', '');
fclose(fid);
lines = C{1};
nLines = numel(lines);

% Regex patterns
patEigenHdr = '^\$EIGENVALUE\s*=\s*([0-9\.\+\-Ee]+)\s*MODE\s*=\s*([0-9]+)'; % eigenvalue & mode
patGridLine = '^\s*([0-9]+)\s+G\s+([+\-]?\d\.\d+E[+\-]\d+)\s+([+\-]?\d\.\d+E[+\-]\d+)\s+([+\-]?\d\.\d+E[+\-]\d+)';
patContLine = '^\-CONT-\s+([+\-]?\d\.\d+E[+\-]\d+)\s+([+\-]?\d\.\d+E[+\-]\d+)\s+([+\-]?\d\.\d+E[+\-]\d+)';

% First pass: find all mode headers (line indices, eigenvalues, modes)
modeIdx   = [];  % line numbers where a new mode begins (header line)
eigVals   = [];  % eigenvalues (rad/s)^2
modeNums  = [];  % mode numbers (as printed)
for i = 1:nLines
    t = lines{i};
    m = regexp(t, patEigenHdr, 'tokens', 'once');
    if ~isempty(m)
        eigVals(end+1,1)  = str2double(m{1}); %#ok<AGROW>
        modeNums(end+1,1) = str2double(m{2}); %#ok<AGROW>
        modeIdx(end+1,1)  = i; %#ok<AGROW>
    end
end

if isempty(modeIdx)
    error('No $EIGENVALUE / MODE headers found. Is this a SOL103 eigenvector punch file?');
end

% Compute natural frequencies (Hz)
freq_hz = sqrt(eigVals) / (2*pi);

% Export natural frequencies to base workspace
assignin('base', 'freq_hz', freq_hz);

% Second pass: read grid blocks per mode
modesCount = numel(modeIdx);

% Keep consistent grid ordering
gridOrderMap = containers.Map('KeyType','char','ValueType','double');
nextOrdinal = 1;

% Temporary storage: per-mode grid maps
modeGridData = cell(modesCount,1);

for k = 1:modesCount
    startLine = modeIdx(k);
    if k < modesCount
        endLine = modeIdx(k+1) - 1;
    else
        endLine = nLines;
    end

    gridMap = containers.Map('KeyType','char','ValueType','any');
    i = startLine + 1;
    while i <= endLine
        L = lines{i};

        % Try main grid line
        mMain = regexp(L, patGridLine, 'tokens', 'once');
        if ~isempty(mMain)
            gridID = str2double(mMain{1});
            T1 = str2double(mMain{2});
            T2 = str2double(mMain{3});
            T3 = str2double(mMain{4});

            % Expect continuation for rotations
            if i+1 <= endLine
                Lc = lines{i+1};
                mCont = regexp(Lc, patContLine, 'tokens', 'once');
                if ~isempty(mCont)
                    R1 = str2double(mCont{1});
                    R2 = str2double(mCont{2});
                    R3 = str2double(mCont{3});
                    vec6 = [T1 T2 T3 R1 R2 R3];

                    key = sprintf('%d', gridID);
                    gridMap(key) = vec6;

                    if ~isKey(gridOrderMap, key)
                        gridOrderMap(key) = nextOrdinal;
                        nextOrdinal = nextOrdinal + 1;
                    end

                    i = i + 2; % consumed main + cont
                    continue
                end
            end

            % Missing/malformed continuation: pad rotations with NaN
            vec6 = [T1 T2 T3 NaN NaN NaN];
            key = sprintf('%d', gridID);
            gridMap(key) = vec6;
            if ~isKey(gridOrderMap, key)
                gridOrderMap(key) = nextOrdinal;
                nextOrdinal = nextOrdinal + 1;
            end
            i = i + 1;
            continue
        end

        i = i + 1;
    end

    modeGridData{k} = gridMap;
end

% Build ordered grid list
keysAll = gridOrderMap.keys;
ord  = zeros(numel(keysAll),1);
gids = zeros(numel(keysAll),1);
for j = 1:numel(keysAll)
    ord(j)  = gridOrderMap(keysAll{j});
    gids(j) = str2double(keysAll{j});
end
[~, idxSort] = sort(ord, 'ascend');
grids = gids(idxSort);

% >>> NEW: export grids to base workspace
assignin('base', 'grids', grids);

% Assemble EV struct with Mode<modeNum> fields as [nGrids x 6]
EV = struct();
for k = 1:modesCount
    gridMap = modeGridData{k};
    M = nan(numel(grids), 6);
    for j = 1:numel(grids)
        key = sprintf('%d', grids(j));
        if isKey(gridMap, key)
            M(j,:) = gridMap(key);
        end
    end
    fieldName = sprintf('Mode%d', modeNums(k));
    EV.(fieldName) = M;

    % Create named variables per grid/mode if requested
    if createNamedVars
        for j = 1:numel(grids)
            vec6 = M(j,:);
            varName = sprintf('EV_grid%d_mode%d', grids(j), modeNums(k));
            assignin('base', varName, vec6);
        end
    end
end

% (Optional) Build a 3D stack if desired:
% nGr = numel(grids);
% nMd = numel(modeNums);
% EVstack = nan(nGr, 6, nMd);
% for k = 1:nMd
%     EVstack(:,:,k) = EV.(sprintf('Mode%d', modeNums(k)));
% end

end
