%% MAKE_ITRAIN_FILE
%
% Create compact Preisach training-current histories for a magnet operated
% from -3 A to +3 A.
%
% Every history:
%   1. Begins at -3 A.
%   2. Moves monotonically to a local operating current.
%   3. Makes an approximately +/-20 percent minor loop around that current.
%   4. Includes smaller nested reversals.
%
% The output file is formatted similarly to:
%   Minor Hysteresis Currents.txt
%
% Output:
%   Preisach Training Currents.txt
%
% Each training data set contains 24 current points.

clear;
clc;

%% User-adjustable settings

outputFileName = 'Preisach Training Currents.txt';

% Magnet current limits.
Imin = -3.0;
Imax =  3.0;

% Centers of the expected local/random minor-loop operating regions.
% There are 9 training sets.
operatingCenters = [ ...
    -2.4, ...
    -1.8, ...
    -1.2, ...
    -0.6, ...
     0.0, ...
     0.6, ...
     1.2, ...
     1.8, ...
     2.4 ];

% For nonzero centers, outer loop amplitude is:
%   outerAmplitude = 20% * abs(operatingCenter)
%
% At zero current, use a fixed loop amplitude because 20% of zero is zero.
percentExcursion = 0.20;
zeroCurrentExcursion = 0.30;     % A

% Number of current points used for the initial path:
%   -3 A -> operating center
%
% Including the initial -3 A and the center current.
numberOfRampPoints = 12;

% Nested-loop excursion expressed as a fraction of the outer excursion.
% 0.50 means the nested reversal is half of the +/-20% outer excursion.
nestedFraction = 0.50;

%% Create current-history cell array

numberOfSets = numel(operatingCenters);

Itrain = cell(numberOfSets, 1);

for k = 1:numberOfSets

    Icenter = operatingCenters(k);

    % Determine the outer local-loop amplitude.
    if abs(Icenter) < eps
        outerAmplitude = zeroCurrentExcursion;
    else
        outerAmplitude = percentExcursion * abs(Icenter);
    end

    % Nested-loop amplitude.
    nestedAmplitude = nestedFraction * outerAmplitude;

    % Outer current limits around the operating center.
    Iupper = Icenter + outerAmplitude;
    Ilower = Icenter - outerAmplitude;

    % Nested current limits around the operating center.
    InestedUpper = Icenter + nestedAmplitude;
    InestedLower = Icenter - nestedAmplitude;

    % Protect against requested values outside the magnet current range.
    Iupper = min(Iupper, Imax);
    Ilower = max(Ilower, Imin);

    InestedUpper = min(InestedUpper, Imax);
    InestedLower = max(InestedLower, Imin);

    % ---------------------------------------------------------------
    % Initial monotonic current ramp:
    %
    % -3 A --> operating-center current
    % ---------------------------------------------------------------
    rampToCenter = linspace(Imin, Icenter, numberOfRampPoints);

    % ---------------------------------------------------------------
    % Local outer and nested minor-loop sequence.
    %
    % The center current is already the final point of rampToCenter,
    % so it is not repeated as the first entry below.
    %
    % Path:
    %
    % center
    %   -> midpoint toward upper limit
    %   -> upper limit
    %   -> midpoint toward upper limit
    %   -> center
    %   -> midpoint toward lower limit
    %   -> lower limit
    %   -> midpoint toward lower limit
    %   -> center
    %   -> nested upper reversal
    %   -> center
    %   -> nested lower reversal
    %   -> center
    % ---------------------------------------------------------------

    upperMidpoint = 0.5 * (Icenter + Iupper);
    lowerMidpoint = 0.5 * (Icenter + Ilower);

    localLoop = [ ...
        upperMidpoint, ...
        Iupper, ...
        upperMidpoint, ...
        Icenter, ...
        lowerMidpoint, ...
        Ilower, ...
        lowerMidpoint, ...
        Icenter, ...
        InestedUpper, ...
        Icenter, ...
        InestedLower, ...
        Icenter ];

    % Complete current history.
    Itrain{k} = [rampToCenter, localLoop];

    % Safety check.
    if numel(Itrain{k}) > 25
        error('Itrain{%d} has more than 25 points.', k);
    end

    if any(Itrain{k} < Imin) || any(Itrain{k} > Imax)
        error('Itrain{%d} contains a current outside [%g, %g] A.', ...
            k, Imin, Imax);
    end
end

%% Display current sets in the MATLAB Command Window

fprintf('\n');
fprintf('Preisach Training Current Sets\n');
fprintf('===============================================\n');

for k = 1:numberOfSets

    fprintf('I_train%d = ', k);

    for n = 1:numel(Itrain{k})

        if n < numel(Itrain{k})
            fprintf('%.5g,', Itrain{k}(n));
        else
            fprintf('%.5g', Itrain{k}(n));
        end
    end

    fprintf(', Length = %d\n', numel(Itrain{k}));
end

fprintf('\n');

%% Write current sets to text file

fileID = fopen(outputFileName, 'w');

if fileID == -1
    error('Could not open output file: %s', outputFileName);
end

fprintf(fileID, 'Preisach Training Current Loops\n');
fprintf(fileID, 'Current range: %.3f A to %.3f A\n', Imin, Imax);
fprintf(fileID, 'All loops begin at %.3f A\n', Imin);
fprintf(fileID, 'Outer local excursion: +/- %.1f %% of operating current\n', ...
    100 * percentExcursion);
fprintf(fileID, 'Zero-current outer excursion: +/- %.3f A\n\n', ...
    zeroCurrentExcursion);

for k = 1:numberOfSets

    fprintf(fileID, 'I_train%d = ', k);

    for n = 1:numel(Itrain{k})

        if n < numel(Itrain{k})
            fprintf(fileID, '%.5f,', Itrain{k}(n));
        else
            fprintf(fileID, '%.5f', Itrain{k}(n));
        end
    end

    fprintf(fileID, ', Length = %d\n', numel(Itrain{k}));
end

fclose(fileID);

fprintf('Training currents written to:\n');
fprintf('   %s\n\n', outputFileName);

%% Optional plot of all training current histories

figure('Name', 'Preisach Training Current Histories', 'Color', 'w');
hold on;
grid on;
box on;

colors = lines(numberOfSets);

for k = 1:numberOfSets

    pointNumber = 1:numel(Itrain{k});

    plot(pointNumber, Itrain{k}, ...
        '-o', ...
        'Color', colors(k,:), ...
        'LineWidth', 1.2, ...
        'MarkerSize', 4, ...
        'DisplayName', sprintf('I\\_train%d, center = %.1f A', ...
        k, operatingCenters(k)));
end

xlabel('Measurement Point Number');
ylabel('Current [A]');
title('Recommended Preisach Training Current Histories');
legend('Location', 'eastoutside');
ylim([Imin - 0.1, Imax + 0.1]);

%% Save data to MAT file also

save('Preisach Training Currents.mat', ...
    'Itrain', ...
    'operatingCenters', ...
    'Imin', ...
    'Imax', ...
    'percentExcursion', ...
    'zeroCurrentExcursion', ...
    'nestedFraction');

fprintf('MATLAB data also saved to:\n');
fprintf('   Preisach Training Currents.mat\n');