function [jeoldata]=mjs_jeolimport()
%   [jeoldata]=jeolimport()
%   Imports PFG-NMR data in Jeol delta format
%   Useage: Point to the *.jdf file that contains the raw data. The
%           imported data will be returned in the structure jeoldata
%           containing the following members:
%               procpar: structure containing the information in the *.hdr
%                        file (process and acqusition parameters)
%               ngrad: number of gradient levels
%               sw: spectral width (in ppm)
%               sp: start of spectrum (in ppm)
%               filename: original file name and path
%               np: number of complex data points per gradient level
%               sfrq: spectrometer frequency (im MHz)
%               at: acquisition time (in seconds)
%               gamma: magnetogyric ratio of the nucleus
%               Gzlvl: gradient strengths
%               DELTA: diffusion time
%               delta: diffusion encoding time
%               dosyconstant: gamma.^2*delts^2*DELTAprime
%               FID: Free induction decays

%   Example:
%   See also: GNAT, dosy_mn, score_mn, decra_mn, mcr_mn, varianimport,
%             brukerimport, jeolimport, peakpick_mn, dosyplot_mn,
%             dosyresidual, dosyplot_gui, scoreplot_mn, decraplot_mn,
%             mcrplot_mn
%
%   This is a part of the GNAT
%   Copyright 2017  <Mathias Nilsson>

%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License along
%   with this program; if not, write to the Free Software Foundation, Inc.,
%   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
%   Dr. Mathias Nilsson
%   School of Chemistry, University of Manchester,
%   Oxford Road, Manchester M13 9PL, UK
%   Telephone: +44 (0) 161 306 4465
%   Fax: +44 (0)161 275 4598
%   mathias.nilsson@manchester.ac.uk

[file, path]=uigetfile('*.jdf','Choose the JEOL file (*.jdf)');
clc
byte_format='b'; %always Big Endian for Header, History and List
fileid=fopen([path file],'r',byte_format);
jeoldata.filename=[path file];

%% Read in Header
File_Identifier=fread(fileid,8,'char=>char')'; %#ok<*NASGU> % JEOL.NMR otherwise incorrect data
Header.Endian=fread(fileid,1,'uint8');  % 0 for Big, 1 for Little
Header.Major_Version=fread(fileid,1,'uint8');   % must be 1
Header.Minor_Version=fread(fileid,2,'uint8');   % must be 1
Header.Data_Dimension_Number=fread(fileid,1,'uint8');   % 1-8
Header.Data_Dimension_Exist=fread(fileid,1,'uint8');
Header.Data_Type=fread(fileid,1,'ubit2');
Header.Data_Format=fread(fileid,1,'bit6');
Header.Instrument=fread(fileid,1,'uint8');
Header.Translate=fread(fileid,8,'uint8');
Header.Data_Axis_Type=fread(fileid,8,'uint8'); % 1=real, 3=complex
for k=1:8
    Header.Data_Units(k).Power=fread(fileid,1,'ubit4');
    Header.Data_Units(k).SIprefix=fread(fileid,1,'ubit4');
    Header.Data_Units(k).Base=fread(fileid,1,'int8');
end
Header.Title=fread(fileid,124,'char=>char')';
Header.Data_Axis_Ranged=fread(fileid,8,'ubit4')';
Header.Data_Points=fread(fileid,8,'int32')';
Header.Data_Offset_Start=fread(fileid,8,'int32')';
Header.Data_Offset_Stop=fread(fileid,8,'int32')';
Header.Data_Axis_Start=fread(fileid,8,'double')';
Header.Data_Axis_Stop=fread(fileid,8,'double')';
Header.Creation_Time=fread(fileid,8,'ubit4');
Header.Revision_Time=fread(fileid,8,'ubit4');
Header.Node_Name=fread(fileid,16,'char=>char')';
Header.Site=fread(fileid,128,'char=>char')';
Header.Author=fread(fileid,128,'char=>char')';
Header.Comment=fread(fileid,128,'char=>char')';
Header.Data_Axis_Titles=fread(fileid,256,'char=>char')';
Header.Base_Freq=fread(fileid,8,'double')';
Header.Zero_Point=fread(fileid,8,'float64')';
Header.Reversed=fread(fileid,8,'ubit1')';
fread(fileid,11);
Header.History_Used=fread(fileid,1,'int32')';
Header.History_Length=fread(fileid,1,'int32')';
Header.Param_Start=fread(fileid,1,'int32')';
Header.Param_Length=fread(fileid,1,'int32')';
Header.List_Start=fread(fileid,8,'int32')';
Header.List_Length=fread(fileid,8,'int32')';
Header.Data_Start=fread(fileid,1,'int32')';
Header.Data_Length=fread(fileid,1,'int64')';
Header.Context_Start=fread(fileid,1,'int64')';
Header.Context_Length=fread(fileid,1,'int32')';
Header.Annote_Start=fread(fileid,1,'int64')';
Header.Annote_Length=fread(fileid,1,'int32')';
Header.Total_Size=fread(fileid,1,'int64')';
Header.Unit_Location=fread(fileid,8,'int8')';
Header.Compound_Units=fread(fileid,2,'bit12')';

%% Read in History
fclose(fileid);
byte_format='b';
fileid=fopen([path file],'r',byte_format);
fread(fileid,1360);
History=fread(fileid,Header.History_Length,'char=>char')'; %#ok<NASGU>

%% Read in Parameters

if Header.Endian==0
    byte_format='b'; % Big Endian
else
    byte_format='l'; % Little Endian
end

fclose(fileid);
fileid=fopen([path file],'r',byte_format);
fread(fileid,Header.Param_Start);

% Read parameter section header
Parameter_Size=fread(fileid,1,'int32')'; %#ok<NASGU>
Low_Index=fread(fileid,1,'int32')';      %#ok<NASGU>
High_Index=fread(fileid,1,'int32')';
Total_Size=fread(fileid,1,'int32')';     %#ok<NASGU>

nparams=High_Index+1;

for k=1:nparams
    fread(fileid,4); % skip undocumented Class Structure
    Unit_Scaler=fread(fileid,1,'int16')';
    for m=1:5
        Units(m).Power=fread(fileid,1,'ubit4');   %#ok<*AGROW>
        Units(m).SIprefix=fread(fileid,1,'ubit4');
        Units(m).Base=fread(fileid,1,'int8');
    end
    fread(fileid,16);
    Value_Type=fread(fileid,1,'int32')';
    fseek(fileid,-20,'cof');

    switch Value_Type
        case 0 % string
            Value=fread(fileid,16,'char=>char')';
            fread(fileid,4);
        case 1 % integer
            Value=fread(fileid,1,'int32')';
            fread(fileid,16);
        case 2 % float
            Value=fread(fileid,1,'float64')';
            fread(fileid,12);
        case 3 % complex
            rl=fread(fileid,1,'float64')';
            im=fread(fileid,1,'float64')';
            Value=complex(rl,im);
            fread(fileid,4);
        case 4 % Infinity
            disp('Infinity - not sure how to handle this')
            Value=fread(fileid,1,'int32')';
            fread(fileid,16);
        otherwise
            disp('Unknown Value_Type')
            fread(fileid,20);
    end
    Name=fread(fileid,28,'char=>char')';
    Name(isspace(Name))=[];
    Name=strrep(Name,'.','_'); % Replace '.' with '_'
    procpar.(Name).Value=Value;
    procpar.(Name).Units=Units;
    procpar.(Name).Unit_Scaler=Unit_Scaler;
    procpar.(Name).Value_Type=Value_Type;
end


%% Read in Data

switch Header.Data_Format

    case 1 % 1D
        disp('1D')
        % Per the JDF spec, for One_D format:
        %   submatrix edge = 8
        %   Data_Points[1] must be a multiple of 8
        %   For Complex data there are 2 sections: Real then Imaginary
        %   Each section is Data_Points[1] doubles laid out sequentially
        %   (the 8-point submatrix blocking is contiguous so order is preserved)
        %
        % Data_Offset_Start/Stop define the VALID point range within Data_Points.

        stored_points = Header.Data_Points(1);  % total stored (multiple of 8)
        readpoints    = Header.Data_Offset_Stop(1) - Header.Data_Offset_Start(1) + 1;
        offset_start  = Header.Data_Offset_Start(1);

        if Header.Data_Type==0 % 64-bit float
            jeoldata.arraydim = 1;
            jeoldata.ngrad    = 1;

            if Header.Data_Axis_Type(1) == 3 % Complex
                % Section 1: Real, Section 2: Imaginary
                % Each section is stored_points doubles
                fseek(fileid, Header.Data_Start, 'bof');
                REAL_all = fread(fileid, stored_points, 'float64');
                IMAG_all = fread(fileid, stored_points, 'float64');
                % Extract only the valid range
                REAL = REAL_all(offset_start+1 : offset_start+readpoints);
                IMAG = IMAG_all(offset_start+1 : offset_start+readpoints);
                DATA = complex(REAL, IMAG);

            else % Real only (Data_Axis_Type == 1)
                fseek(fileid, Header.Data_Start, 'bof');
                REAL_all = fread(fileid, stored_points, 'float64');
                DATA = REAL_all(offset_start+1 : offset_start+readpoints);
            end

        elseif Header.Data_Type==1 % 32-bit float
            disp('32bit data - untested import')
            fseek(fileid, Header.Data_Start, 'bof');
            REAL_all = fread(fileid, stored_points, 'float32');
            IMAG_all = fread(fileid, stored_points, 'float32');
            REAL = REAL_all(offset_start+1 : offset_start+readpoints);
            IMAG = IMAG_all(offset_start+1 : offset_start+readpoints);
            DATA = complex(REAL, IMAG);
        else
            error('Unknown Data_Type')
        end

    case {2, 12} % 2D (Two_D submatrix edge=32) or Small_Two_D (edge=4)
        if Header.Data_Format==2
            disp('2D')
            submatrix_edge = 32;
        else
            disp('small 2D')
            submatrix_edge = 4;
        end

        jeoldata.arraydim = procpar.Y_POINTS.Value;
        jeoldata.ngrad    = procpar.Y_POINTS.Value;

        stored_x   = Header.Data_Points(1);
        stored_y   = Header.Data_Points(2);
        readpoints = Header.Data_Offset_Stop(1) - Header.Data_Offset_Start(1) + 1;
        twoDsize   = Header.Data_Offset_Stop(2) - Header.Data_Offset_Start(2) + 1;

        % Total stored points per section (Real or Imag)
        section_points = stored_x * stored_y;

        % Per spec: data is stored as submatrices of size submatrix_edge x submatrix_edge
        % laid out in row-major order. Read the whole section then extract valid data.
        fseek(fileid, Header.Data_Start, 'bof');
        REAL_section = fread(fileid, section_points, 'float64');
        IMAG_section = fread(fileid, section_points, 'float64');

        % Reshape into [stored_x, stored_y] accounting for submatrix layout
        % For Small_Two_D the submatrix ordering means we reshape as:
        % [submatrix_edge, submatrix_edge, nx_sub, ny_sub] then permute
        nx_sub = stored_x / submatrix_edge;
        ny_sub = stored_y / submatrix_edge;
        REAL_mat = reshape(REAL_section, submatrix_edge, submatrix_edge, nx_sub, ny_sub);
        REAL_mat = permute(REAL_mat, [1,3,2,4]);
        REAL_mat = reshape(REAL_mat, stored_x, stored_y);
        IMAG_mat = reshape(IMAG_section, submatrix_edge, submatrix_edge, nx_sub, ny_sub);
        IMAG_mat = permute(IMAG_mat, [1,3,2,4]);
        IMAG_mat = reshape(IMAG_mat, stored_x, stored_y);

        % Extract valid range
        x0 = Header.Data_Offset_Start(1) + 1;
        y0 = Header.Data_Offset_Start(2) + 1;
        DATA = complex(REAL_mat(x0:x0+readpoints-1, y0:y0+twoDsize-1), ...
                       IMAG_mat(x0:x0+readpoints-1, y0:y0+twoDsize-1));

    otherwise
        % General N-D reader using the submatrix unscrambling approach from the JDF spec.
        % Submatrix edge lengths per Data_Format:
        %   Two_D(2)=32, Three_D(3)=8, Four_D(4)=8
        %   Small_Two_D(12)=4, Small_Three_D(13)=4, Small_Four_D(14)=4
        switch Header.Data_Format
            case {3},        submatrix_edge = 8;
            case {4},        submatrix_edge = 8;
            case {13, 14},   submatrix_edge = 4;
            otherwise,       submatrix_edge = 4;
        end

        ndim         = Header.Data_Dimension_Number;
        stored_pts   = Header.Data_Points(1:ndim);   % row vector of stored sizes
        total_stored = prod(stored_pts);
        fprintf('%dD data (Data_Format=%d), submatrix_edge=%d, stored_pts=%s\n', ...
            ndim, Header.Data_Format, submatrix_edge, mat2str(stored_pts))

        % Read the two flat sections (Real, then Imaginary)
        fseek(fileid, Header.Data_Start, 'bof');
        REAL_flat = fread(fileid, total_stored, 'float64');
        IMAG_flat = fread(fileid, total_stored, 'float64');

        % Unscramble submatrix tiling for each section.
        % The spec stores data as row-major submatrices of size edge^ndim.
        % reshape into [edge, edge, ...(ndim times)..., nsub1, nsub2, ...]
        % then permute to interleave the within-submatrix and between-submatrix dims
        % and reshape back to [stored_pts].
        n_sub   = stored_pts / submatrix_edge;   % number of submatrices per dim
        % Build reshape size: [edge edge ... edge  nsub1 nsub2 ... nsubN]
        reshape_sz = [repmat(submatrix_edge, 1, ndim), n_sub];
        % Build permutation to interleave: dim1_within, dim1_between, dim2_within, ...
        % i.e. [1, ndim+1,  2, ndim+2,  3, ndim+3, ...]
        perm = reshape([1:ndim; ndim+1:2*ndim], 1, []);

        REAL_nd = reshape(REAL_flat, reshape_sz);
        REAL_nd = permute(REAL_nd, perm);
        REAL_nd = reshape(REAL_nd, stored_pts);

        IMAG_nd = reshape(IMAG_flat, reshape_sz);
        IMAG_nd = permute(IMAG_nd, perm);
        IMAG_nd = reshape(IMAG_nd, stored_pts);

        % Extract valid range using Data_Offset_Start/Stop for each dim
        idx = cell(1, ndim);
        for d = 1:ndim
            idx{d} = Header.Data_Offset_Start(d)+1 : Header.Data_Offset_Stop(d)+1;
        end
        DATA = complex(REAL_nd(idx{:}), IMAG_nd(idx{:}));
        fprintf('DATA size after extraction: %s\n', mat2str(size(DATA)))

end

fclose(fileid);

%% Set up the data for the DOSY Toolbox

jeoldata.procpar = procpar;
jeoldata.lp = 0;
jeoldata.rp = 0;
jeoldata.np = procpar.X_POINTS.Value;

% Convert all frequency/time parameters to base SI first, then derive ppm.
% jeol_unitstring multiplier converts raw value -> base SI (Hz for freq, s for time).

% Spectrometer frequency
[sfrq_unit, sfrq_mult] = jeol_unitstring(procpar.X_FREQ);
sfrq_hz       = procpar.X_FREQ.Value * sfrq_mult;   % Hz
jeoldata.sfrq = sfrq_hz * 1e-6;                      % MHz
fprintf('X_FREQ: %.6g %s -> %.6g MHz\n', procpar.X_FREQ.Value, sfrq_unit, jeoldata.sfrq)

% Spectral width: sw(ppm) = sw(Hz) / sfrq(MHz)  [since 1 ppm = sfrq_MHz Hz]
[sw_unit, sw_mult] = jeol_unitstring(procpar.X_SWEEP);
sw_hz       = procpar.X_SWEEP.Value * sw_mult;       % Hz
jeoldata.sw = sw_hz / jeoldata.sfrq;                  % ppm  (Hz / MHz = ppm)
fprintf('X_SWEEP: %.6g %s -> %.6g Hz -> %.6g ppm\n', procpar.X_SWEEP.Value, sw_unit, sw_hz, jeoldata.sw)

% Acquisition time
[at_unit, at_mult] = jeol_unitstring(procpar.X_ACQ_DURATION);
jeoldata.at = procpar.X_ACQ_DURATION.Value * at_mult; % s
fprintf('X_ACQ_DURATION: %.6g %s -> %.6g s\n', procpar.X_ACQ_DURATION.Value, at_unit, jeoldata.at)

jeoldata.tr = procpar.transition_ratio.Value;

% Magnetogyric ratio based on observed nucleus (rad/s/T)
if isfield(procpar, 'X_DOMAIN')
    nucleus = strtrim(procpar.X_DOMAIN.Value);
    switch nucleus
        case 'Proton',       jeoldata.gamma =  267524618.573;
        case 'Deuterium',    jeoldata.gamma =   41066279.1;
        case 'Carbon13',     jeoldata.gamma =   67282840.0;
        case 'Nitrogen15',   jeoldata.gamma =  -27116000.0;
        case 'Fluorine19',   jeoldata.gamma =  251814800.0;
        case 'Phosphorus31', jeoldata.gamma =  108394000.0;
        case 'Lithium7',     jeoldata.gamma =  103962000.0;
        case 'Boron11',      jeoldata.gamma =   85847000.0;
        case 'Sodium23',     jeoldata.gamma =   70761000.0;
        case 'Silicon29',    jeoldata.gamma =  -53190000.0;
        case 'Copper63',     jeoldata.gamma =   71117930.0;
        case 'Oxygen17',     jeoldata.gamma =  -36264000.0;
        case 'Nitrogen14',   jeoldata.gamma =   19331000.0;
        case 'Chlorine35',   jeoldata.gamma =   26241000.0;
        case 'Potassium39',  jeoldata.gamma =   12500000.0;
        case 'Cesium133',    jeoldata.gamma =   35271000.0;
        otherwise
            warning('Unknown nucleus ''%s'' - defaulting to 1H gamma', nucleus)
            jeoldata.gamma = 267524618.573;
    end
else
    warning('X_DOMAIN parameter not found - defaulting to 1H gamma')
    jeoldata.gamma = 267524618.573;
end

%% Defaults for diffusion/gradient fields (overwritten below if DOSY data)
jeoldata.Gzlvl        = [];
jeoldata.ngrad        = 1;
jeoldata.dosyconstant = 'non existing';
jeoldata.DELTA        = 'non existing';
jeoldata.delta        = 'non existing';
jeoldata.DAC_to_G     = 'non existing';

%% ===================== READ JEOL LISTS (ALL ARRAYED AXES) ==================
% Collect every non-trivial listed axis into jeoldata.axes, a struct array.
% Each element has: .index .name .values .unit_str .multiplier
% Axis 1 (X / FID) is always the direct dimension and is skipped.
% This handles 1D, 2D, 3D, ... uniformly.

jeoldata.axes = struct('index',{},'name',{},'values',{},'unit_str',{},'multiplier',{});
raw_titles = reshape(Header.Data_Axis_Titles, 32, [])';

try
    for ax = 2:Header.Data_Dimension_Number

        if Header.Data_Axis_Ranged(ax) ~= 0 && Header.List_Start(ax) > 0

            % Lists are ALWAYS big endian regardless of file endianness
            fid = fopen([path file], 'r', 'b');
            fseek(fid, Header.List_Start(ax), 'bof');
            full_list = fread(fid, Header.Data_Points(ax), 'float64');
            fclose(fid);

            % Trim to valid range
            list = full_list( ...
                Header.Data_Offset_Start(ax)+1 : ...
                Header.Data_Offset_Stop(ax)+1 );

            % Skip trivial (constant) axes
            if numel(unique(list)) <= 1
                continue
            end

            % Decode axis name from header title string
            name = strtrim(raw_titles(ax, :));
            name = regexprep(name, '[^a-zA-Z0-9_]', '');
            name = lower(name);
            if isempty(name)
                name = sprintf('axis_%d', ax);
            end

            % Decode axis units from Header.Data_Units(ax).
            % Construct a minimal procpar-compatible struct so we can reuse
            % jeol_unitstring without duplicating the lookup tables.
            axis_unit_struct.Units = repmat( ...
                struct('SIprefix',0,'Base',0,'Power',0), 1, 5);
            axis_unit_struct.Units(1)    = Header.Data_Units(ax);
            axis_unit_struct.Unit_Scaler = 0;
            [unit_str, multiplier] = jeol_unitstring(axis_unit_struct);

            % Append to axes array
            n = numel(jeoldata.axes) + 1;
            jeoldata.axes(n).index      = ax;
            jeoldata.axes(n).name       = name;
            jeoldata.axes(n).values     = list(:)';
            jeoldata.axes(n).unit_str   = unit_str;
            jeoldata.axes(n).multiplier = multiplier;

            fprintf('Axis %d (%s): %d points, units=%s, SI multiplier=%.6g\n', ...
                ax, name, numel(list), unit_str, multiplier)
        end
    end

catch ME
    warning('jeolimport: error reading axis lists: %s', ME.message)
end

% ---- Convenience aliases for single-array-dimension files (backward compat) -
if numel(jeoldata.axes) >= 1
    jeoldata.array_param_name       = jeoldata.axes(1).name;
    jeoldata.array_param_axis       = jeoldata.axes(1).index;
    jeoldata.array_param_values     = jeoldata.axes(1).values;
    jeoldata.array_param_unit_str   = jeoldata.axes(1).unit_str;
    jeoldata.array_param_multiplier = jeoldata.axes(1).multiplier;
else
    jeoldata.array_param_name       = 'y';
    jeoldata.array_param_axis       = 2;
    jeoldata.array_param_values     = [];
    jeoldata.array_param_unit_str   = 'dimensionless';
    jeoldata.array_param_multiplier = 1;
end

%% ===================== INTERPRET ARRAYED PARAMETERS ======================
% Loop over every detected axis and classify it by name.
% Each axis's values are converted to SI using the multiplier decoded from
% the header unit structure, so the logic is identical for 2D, 3D, 4D...

for ai = 1:numel(jeoldata.axes)
    ax_name  = jeoldata.axes(ai).name;
    ax_vals  = jeoldata.axes(ai).values;
    ax_mult  = jeoldata.axes(ai).multiplier;
    ax_unit  = jeoldata.axes(ai).unit_str;

    % Convert raw list values to SI units using the header-decoded multiplier
    ax_vals_si = ax_vals * ax_mult;

    if strcmp(ax_name, 'g') || contains(ax_name, 'grad')
        % ---- Diffusion gradient axis ------------------------------------
        disp(['Gradient axis detected: ' ax_name ...
              '  unit=' ax_unit ...
              sprintf('  max=%.4g T/m', max(abs(ax_vals_si)))])
        jeoldata.Gzlvl = ax_vals_si;
        if isfield(procpar, 'grad_shape')
        jeoldata.gpnam=procpar.grad_shape.Value;
        disp(['Grad_shape is: ' jeoldata.gpnam])
        end
        if strcmp(jeoldata.gpnam,'SQUARE')
            jeoldata.gpnam="RECT";
        end 
        if isfield(procpar, 'alpha')
            jeoldata.alpha=procpar.alpha.Value;
        end 
        % Diffusion timing from procpar
        if isfield(procpar, 'diffusion_time') && isfield(procpar, 'delta')
            [~, delta_mult] = jeol_unitstring(procpar.delta);
            [~, DELTA_mult] = jeol_unitstring(procpar.diffusion_time);
            jeoldata.delta        = procpar.delta.Value          * delta_mult;
            jeoldata.DELTA        = procpar.diffusion_time.Value * DELTA_mult;
            jeoldata.dosyconstant = jeoldata.gamma^2 * jeoldata.delta^2 * ...
                                    (jeoldata.DELTA - jeoldata.delta/3);
            fprintf('delta=%.4g s  DELTA=%.4g s  dosyconstant=%.4g\n', ...
                jeoldata.delta, jeoldata.DELTA, jeoldata.dosyconstant)
            if isfield(procpar, 'tau')
                [~, tau_mult] = jeol_unitstring(procpar.tau);
                jeoldata.tau = procpar.tau.Value * tau_mult;
            end
            % if isfield(procpar, 'st_coef')
            %     jeoldata.DELTA = procpar.st_coef.Value*DELTA_mult;
            % end
        else
            disp('Could not find diffusion parameters (diffusion_time / delta)')
            jeoldata.DELTA    = 'non existing';
            jeoldata.delta    = 'non existing';
            jeoldata.DAC_to_G = 'non existing';
        end

    elseif strcmp(ax_name, 'delay_array_pre') || strcmp(ax_name, 'tau_interval') || strcmp(ax_name, 'delay_array')
        % ---- Variable delay axis ----------------------------------------
        disp(['Variable delay axis detected: ' ax_name '  unit=' ax_unit])
        jeoldata.vdlist = ax_vals_si;

    elseif strcmp(ax_name, 'loop_number')
        % ---- Variable count axis ----------------------------------------
        disp(['Variable count axis detected: ' ax_name])
        jeoldata.vclist     = ax_vals_si;
        jeoldata.vcconstant = 5*procpar.x_pulse.Value + 4*procpar.tau_interval.Value;
        jeoldata.vdlist     = jeoldata.vclist * jeoldata.vcconstant;

    else
        % ---- Unrecognised axis - store raw SI values --------------------
        fprintf('Unrecognised axis: %s  unit=%s  (stored as axes(%d).values_si)\n', ...
            ax_name, ax_unit, ai)
        jeoldata.axes(ai).values_si = ax_vals_si;
    end
end

% ngrad = number of points along the first detected indirect axis (Y)
if numel(jeoldata.axes) >= 1
    jeoldata.ngrad = numel(jeoldata.axes(1).values);
elseif isfield(procpar, 'Y_POINTS')
    jeoldata.ngrad = procpar.Y_POINTS.Value;
else
    jeoldata.ngrad = 1;
end

% ngrad = number of points along the first detected indirect axis (Y)
if numel(jeoldata.axes) >= 2
    jeoldata.ngrad2 = numel(jeoldata.axes(2).values);
elseif isfield(procpar, 'Z_POINTS')
    jeoldata.ngrad2 = procpar.Z_POINTS.Value;
else
    jeoldata.ngrad2 = 1;
end
%jeoldata.arraydim = jeoldata.ngrad * jeoldata.ngrad2;
data_sz = size(DATA);
jeoldata.np      = data_sz(1);
jeoldata.arraydim = prod(data_sz(2:end));   % collapse all indirect dims
jeoldata.FID     = reshape(DATA, jeoldata.np, jeoldata.arraydim);

% Spectral offset (ppm) - must come after jeoldata.gamma is assigned.
% sp = X_OFFSET(ppm) - sw(ppm)/2
% X_OFFSET is the centre of the spectrum in ppm, sw is already in ppm above.
jeoldata.sp = procpar.X_OFFSET.Value - jeoldata.sw / 2;
fprintf('sp = %.4g ppm  (X_OFFSET=%.4g ppm, sw=%.4g ppm)\n', ...
    jeoldata.sp, procpar.X_OFFSET.Value, jeoldata.sw)

%jeoldata.FID = DATA;

%% Digital filter group-delay correction (digshift)

total   = 0;
factors = str2num(procpar.factors.Value); %#ok<ST2NM>
orders  = str2num(procpar.orders.Value);  %#ok<ST2NM>
for i = 1:orders(1)
    total = total + (orders(i+1) - 1) / prod(factors(i:end));
end
total = 0.5 * total;
% Use sw_hz (already in Hz) rather than raw X_SWEEP.Value which may be in kHz etc.
delay = total / sw_hz;
jeoldata.digshift = ceil(delay / jeoldata.at * (procpar.X_POINTS.Value - 1));

end

%% =========================================================================
function [unit_str, multiplier] = jeol_unitstring(param)
%  [unit_str, multiplier] = jeol_unitstring(param)
%
%  Decodes the unit structure of a JEOL JDF parameter into a human-readable
%  string and an SI multiplier to convert the raw value to base SI units.
%
%  Input:
%    param      - one parameter struct as read by jeolimport, e.g.
%                 procpar.X_FREQ  or  procpar.X_SWEEP
%                 Must contain fields: .Value, .Units, .Unit_Scaler
%                 Units is an array of up to 5 unit structs, each with:
%                   .Base      (integer, see JDF spec)
%                   .SIprefix  (integer, see JDF spec)
%                   .Power     (integer, exponent on the unit)
%
%  Outputs:
%    unit_str   - string such as 'MHz', 'ms', 'mT/m', 'Hz', 's', 'ppm'
%    multiplier - multiply param.Value by this to get the value in base
%                 SI units (e.g. Hz, s, T/m).
%                 Includes both the SI prefix and Unit_Scaler.
%
%  Example:
%    [u, m] = jeol_unitstring(procpar.X_FREQ)
%    % u = 'MHz',  m = 1e6  ->  freq_hz = procpar.X_FREQ.Value * m
%
%    [u, m] = jeol_unitstring(procpar.X_ACQ_DURATION)
%    % u = 'ms',   m = 1e-3  ->  at_s = procpar.X_ACQ_DURATION.Value * m

% SI prefix lookup  (SIprefix is a 4-bit signed integer stored in the file)
% Negative codes = super-unit (k, M, G...), positive = sub-unit (m, u, n...)
prefix_codes = [-8,-7,-6,-5,-4,-3,-2,-1, 0, 1, 2, 3, 4, 5, 6, 7];
prefix_names = {'Y','Z','E','P','T','G','M','k','','m','u','n','p','f','a','z'};
prefix_exp   = [ 24, 21, 18, 15, 12,  9,  6,  3,  0,-3,-6,-9,-12,-15,-18,-21];

% SI base unit lookup  (Base field, 0-indexed in spec, 1-indexed here)
base_names = { ...
    '',          ...  %  0  None
    'abundance', ...  %  1
    'A',         ...  %  2  Ampere
    'cd',        ...  %  3  Candela
    'C',         ...  %  4  Celsius
    'C',         ...  %  5  Coulomb
    'deg',       ...  %  6  Degree
    'eV',        ...  %  7  Electronvolt
    'F',         ...  %  8  Farad
    'Sv',        ...  %  9  Sievert
    'g',         ...  % 10  Gram
    'Gy',        ...  % 11  Gray
    'H',         ...  % 12  Henry
    'Hz',        ...  % 13  Hertz
    'K',         ...  % 14  Kelvin
    'J',         ...  % 15  Joule
    'L',         ...  % 16  Liter
    'lm',        ...  % 17  Lumen
    'lx',        ...  % 18  Lux
    'm',         ...  % 19  Meter
    'mol',       ...  % 20  Mole
    'N',         ...  % 21  Newton
    'Ohm',       ...  % 22  Ohm
    'Pa',        ...  % 23  Pascal
    '%',         ...  % 24  Percent
    'pt',        ...  % 25  Point
    'ppm',       ...  % 26  Ppm
    'rad',       ...  % 27  Radian
    's',         ...  % 28  Second
    'S',         ...  % 29  Siemens
    'sr',        ...  % 30  Steradian
    'T',         ...  % 31  Tesla
    'V',         ...  % 32  Volt
    'W',         ...  % 33  Watt
    'Wb',        ...  % 34  Weber
    'dB',        ...  % 35  Decibel
    'Da',        ...  % 36  Dalton
    'Th',        ...  % 37  Thompson
    ''           ...  % 38  Ugeneric (treated as None)
};

unit_str   = '';
multiplier = 1;

for k = 1:5
    u = param.Units(k);

    if u.Base == 0 || u.Base == 38 || u.Power == 0
        continue   % unused slot or zero exponent
    end

    % Look up prefix
    idx = find(prefix_codes == u.SIprefix, 1);
    if isempty(idx)
        warning('jeol_unitstring: unknown SIprefix %d', u.SIprefix)
        pfx_name = '?';  pfx_exp = 0;
    else
        pfx_name = prefix_names{idx};
        pfx_exp  = prefix_exp(idx);
    end

    % Look up base name
    if u.Base + 1 > numel(base_names)
        warning('jeol_unitstring: unknown Base %d', u.Base)
        base_name = '?';
    else
        base_name = base_names{u.Base + 1};
    end

    % Build string component
    if u.Power == 1
        component = [pfx_name base_name];
    elseif u.Power == -1
        component = ['/' pfx_name base_name];
    else
        component = [pfx_name base_name '^' num2str(u.Power)];
    end
    unit_str = [unit_str component]; %#ok<AGROW>

    % Accumulate SI multiplier
    % Special case: Gram - base SI is kg, so effective prefix is pfx-3
    % Dimensionless units (ppm, %, point) contribute factor 1
    switch u.Base
        case 10  % Gram -> kg
            base_factor = 10^((pfx_exp - 3) * u.Power);
        case {24, 25, 26, 38}  % %, point, ppm, Ugeneric
            base_factor = 1;
        otherwise
            base_factor = 10^(pfx_exp * u.Power);
    end
    multiplier = multiplier * base_factor;
end

% Apply Unit_Scaler: raw value is scaled by 10^Unit_Scaler
if param.Unit_Scaler ~= 0
    multiplier = multiplier * 10^double(param.Unit_Scaler);
end

if isempty(unit_str)
    unit_str = 'dimensionless';
end

end