%% Example of building Simscape FEM-parameterized PMSM block tables (version 10-June-2024)
% https://www.mathworks.com/help/sps/ref/femparameterizedpmsm.html
% Pleace this script to \CustomScripts

% Tabulate DQ-axes flux linkage data (resolved into direct (D) and
% quadrature (Q) axes) or A-phase flux linkage data (assuming the B and C
% phase data is the same except shifted in phase)
% For DQ-axes flux linkage, the rotor angles in the range 0 to 360/nPolePairs/3 degrees are required. 
% For A-phase flux linkage, the rotor angle range must be 0 to 360/nPolePairs degrees. 
% Four parameterization approaches (see below paramApproach):
% 1) D and Q axes flux linkages as a function of D-axis current (iD), Q-axis current (iQ), and rotor angle (theta)
% 2) D and Q axes flux linkages as a function of peak current magnitude (I), current advance angle (B), and rotor angle (theta)
% 3) A-phase flux linkage as a function of D-axis current (iD), Q-axis current (iQ), and rotor angle (theta)
% 4) A-phase flux linkage as a function of peak current magnitude (I), current advance angle (B), and rotor angle (theta)
    
% Note that MotorXP uses Option 1 for the Park’s convention: Q leads D, rotor angle measured from A-phase to D-axis, 
% as described in the Simscape FEM-parameterized PMSM block documentation https://www.mathworks.com/help/sps/ref/femparameterizedpmsm.html

%% Step 1: Initialization
cd ../
initMXPscriptingAPI;
maxNumberOfProcesses = feature('numcores');                      % maximum number of processes working in parallel (equal to the number of physical cpu cores)
% maxNumberOfProcesses = 7;                                      % or defined directly
save2filename = 'example_innerrotor_simscape_multipleCPU.mat';   % save flux and losses tables to this mat-file
accuracyRotorAngle = 'Low';                                      % 'Low' | 'Medium' | 'High'
Irmssupply_max = 200;                                            % maximum RMS supply current in Amperes
Irmssupply_step = Irmssupply_max/10;                             % current step in Amperes
gamma_step = 10;                                                 % advance angle step in electrical degrees
Speed_max = 10000;                                               % maximum speed in RPM
Speed_step = 100;                                                % speed step in RPM
paramApproach = 4;                                               % see above four parameterization approaches
includeIronLoss = 1;                                             % calculate iron and magnet losses

%% Step 2: Specify design mxp-file and setup initial settings
designfile = getpmpath('SimFiles\example_innerrotor\example_innerrotor.mxp'); % mxp-project file
[motorProps,settingsMagnetostatic] = openMXP(designfile);
[Geometry, Windings, Mesh, Materials, SubdomainProperty, Error, motorProps_test, Weight] = assembleMXP(designfile,motorProps);
settingsMagnetostatic=setParampm(settingsMagnetostatic,'Convergence tolerance',0.01);           % convergance tolerance of FEA solution
settingsMagnetostatic=setParampm(settingsMagnetostatic,'Current waveform','Sinusoidal');
settingsMagnetostatic=setParampm(settingsMagnetostatic,'Save each step solution',0);
settingsMagnetostatic=setParampm(settingsMagnetostatic,'Compute cogging torque',0);

%% Step 3: Rotor angle accuracy
nPolePairs = getParampm(motorProps,'Number of pole pairs');
Ns = getParampm(motorProps,'Number of slots');
% determine number of rotor positions per one electrical period
if strcmp(accuracyRotorAngle,'Low')                 % Low accuracy: 2 points per cogging torque period
    numberOfRotorPositions = setNumberOfPoints('Low accuracy',nPolePairs,Ns);
elseif strcmp(accuracyRotorAngle,'Medium')          % Medium accuracy: 4 points per cogging torque period
    numberOfRotorPositions = setNumberOfPoints('Medium accuracy',nPolePairs,Ns);
elseif strcmp(accuracyRotorAngle,'High')            % High accuracy: 8 points per cogging torque period
    numberOfRotorPositions = setNumberOfPoints('High accuracy',nPolePairs,Ns);
else
    error('Undefined rotor angle accuracy');
end           

f1 = 1;     % define supply frequency 1 Hz
settingsMagnetostatic = setParampm(settingsMagnetostatic, 'Mechanical speed', (f1/nPolePairs)*60);
if paramApproach==1 || paramApproach==2
    if ~includeIronLoss
        % For DQ-axes flux linkage, the rotor angles in the range 0 to 360/nPolePairs/3 degrees are required (but iron losses require full electrical period). 
        settingsMagnetostatic = setParampm(settingsMagnetostatic,'Simulation time','Define time interval');
        settingsMagnetostatic = setParampm(settingsMagnetostatic,'Time from',0);
        settingsMagnetostatic = setParampm(settingsMagnetostatic,'Time to',1/3);    % 1/3 of electrical period
        numberOfRotorPositions = round(numberOfRotorPositions/3);
    else
        % iron losses require full electrical period
        settingsMagnetostatic=setParampm(settingsMagnetostatic,'Simulation time','One electrical period');
    end
else  % paramApproach==3 || paramApproach==4
    % For A-phase flux linkage, the rotor angle range must be 0 to 360/nPolePairs degrees. 
    settingsMagnetostatic=setParampm(settingsMagnetostatic,'Simulation time','One electrical period');
end
settingsMagnetostatic = setParampm(settingsMagnetostatic, 'Number of points', numberOfRotorPositions);

%% Step 4: Calculate flux tables and initial iron and magnet losses tables
% Supply current RMS values (Amperes)
Irmssupply = 0:Irmssupply_step:Irmssupply_max;
if Irmssupply(end)~=Irmssupply_max
    Irmssupply = [Irmssupply Irmssupply_max];
end
% Supply current peak values (Amperes)
Ipeak = Irmssupply*sqrt(2);
% Advance angle values (electrical degrees)
Gamma = 0:gamma_step:360;
if Gamma(end)==360
    Gamma(end) = [];
end
% Speed values (RPM)
Speed = 0:Speed_step:Speed_max;
% Rotor angle values (degrees)
if strcmp(getParampm(settingsMagnetostatic,'Simulation time'),'One electrical period')
    RotorRotationAngle = 360/nPolePairs;
    RotorAngle =linspace(-180/nPolePairs,RotorRotationAngle-180/nPolePairs,numberOfRotorPositions+1);
else  % 1/3 of electrical period
    RotorRotationAngle = 360/nPolePairs/3;
    RotorAngle =linspace(-180/nPolePairs,RotorRotationAngle-180/nPolePairs,numberOfRotorPositions);
end

% initialize with NaN values
FluxLinkage3D_a = NaN*ones(length(Irmssupply),length(Gamma),numberOfRotorPositions);
FluxLinkage3D_q = NaN*ones(length(Irmssupply),length(Gamma),numberOfRotorPositions);
FluxLinkage3D_d = NaN*ones(length(Irmssupply),length(Gamma),numberOfRotorPositions);
Torque3D = NaN*ones(length(Irmssupply),length(Gamma),numberOfRotorPositions);
Piron_eddy_stator0 = NaN*ones(length(Irmssupply),length(Gamma));
Piron_eddy_rotor0 = NaN*ones(length(Irmssupply),length(Gamma));
Piron_hyst_stator0 = NaN*ones(length(Irmssupply),length(Gamma));
Piron_hyst_rotor0 = NaN*ones(length(Irmssupply),length(Gamma));
MagnetLoss0 = NaN*ones(length(Irmssupply),length(Gamma));

initMStimestepping_par(); % init running magnetostatic FE simulations in parallel
numberOfProcesses = 0;       % current number of processes
Private = [];
nresultij =[];
% [i j] pair for each result number
for i = 2:length(Irmssupply)
    for j = 1:length(Gamma)
        nresultij =[nresultij;
                    [i j]];
    end
end
for i = 1:length(Irmssupply)
    settingsMagnetostatic = setParampm(settingsMagnetostatic,'RMS supply current',Irmssupply(i));  % supply current value
    if i==1    % for zero current flux and losses do not depend on anvance angle
        settingsMagnetostatic = setParampm(settingsMagnetostatic,'Advance angle',0);
        [Results,Timesteppingdata,~,~,~,Private] = runMStimestepping([],Geometry,Mesh,Windings,settingsMagnetostatic,Materials,SubdomainProperty,[]);
        for j = 1:length(Gamma)
            FluxLinkage3D_a(i,j,:) = Timesteppingdata.Fluxlinkage_a;
            FluxLinkage3D_q(i,j,:) = Timesteppingdata.Fluxlinkage_q;
            FluxLinkage3D_d(i,j,:) = Timesteppingdata.Fluxlinkage_d;
            Torque3D(i,j,:) = Timesteppingdata.Torque_maxwell;
            Piron_eddy_stator0(i,j) = Results.Piron_tooltip.Piron_eddy_stator;
            Piron_eddy_rotor0(i,j) = Results.Piron_tooltip.Piron_eddy_rotor;
            Piron_hyst_stator0(i,j) = Results.Piron_tooltip.Piron_hyst_stator;
            Piron_hyst_rotor0(i,j) = Results.Piron_tooltip.Piron_hyst_rotor;
            MagnetLoss0(i,j) = Results.MagnetLoss;
        end
    else
        for j = 1:length(Gamma)
            settingsMagnetostatic = setParampm(settingsMagnetostatic,'Advance angle',Gamma(j));
            % run magnetostatic FE simulation in parallel
            if numberOfProcesses<maxNumberOfProcesses
                fprintf('run simulation #[%i %i]', i, j);
                % add one process
                runMStimestepping_par(Private,Geometry,Mesh,Windings,settingsMagnetostatic,Materials,SubdomainProperty,[]);
                numberOfProcesses = numberOfProcesses+1;
            end
            if numberOfProcesses == maxNumberOfProcesses
                % wait when one of processes finish calculation
                while 1
                    % check if any results are ready
                    [cell_Results, cell_Timesteppingdata, ~, ~, numberOfProcesses] = getMStimestepping_par();   % get magnetostatic analysis results
                    if numberOfProcesses<maxNumberOfProcesses
                        break;
                    end
                    pause(1); % wait for 1 second
                end
                % put results
                for nresult = 1:length(cell_Results)
                    if ~isempty(cell_Results{nresult})
                        Results = cell_Results{nresult};
                        Timesteppingdata = cell_Timesteppingdata{nresult};
                        ii = nresultij(nresult,1);
                        jj = nresultij(nresult,2);
                        FluxLinkage3D_a(ii,jj,:) = Timesteppingdata.Fluxlinkage_a;
                        FluxLinkage3D_q(ii,jj,:) = Timesteppingdata.Fluxlinkage_q;
                        FluxLinkage3D_d(ii,jj,:) = Timesteppingdata.Fluxlinkage_d;
                        Torque3D(ii,jj,:) = Timesteppingdata.Torque_maxwell;
                        Piron_eddy_stator0(ii,jj) = Results.Piron_tooltip.Piron_eddy_stator;
                        Piron_eddy_rotor0(ii,jj) = Results.Piron_tooltip.Piron_eddy_rotor;
                        Piron_hyst_stator0(ii,jj) = Results.Piron_tooltip.Piron_hyst_stator;
                        Piron_hyst_rotor0(ii,jj) = Results.Piron_tooltip.Piron_hyst_rotor;
                        MagnetLoss0(ii,jj) = Results.MagnetLoss;
                    end
                end
            end
        end
    end
end

% get all (the rest) results
numberOfProcesses0 = numberOfProcesses;
while numberOfProcesses>0
    % check if any results are ready
    [cell_Results, cell_Timesteppingdata, ~, ~, numberOfProcesses] = getMStimestepping_par();
    if numberOfProcesses ~= numberOfProcesses0  % if new results are ready (number of processes changed)
        % put results
        for nresult = 1:length(cell_Results)
            if ~isempty(cell_Results{nresult})
                Results = cell_Results{nresult};
                Timesteppingdata = cell_Timesteppingdata{nresult};
                ii = nresultij(nresult,1);
                jj = nresultij(nresult,2);
                FluxLinkage3D_a(ii,jj,:) = Timesteppingdata.Fluxlinkage_a;
                FluxLinkage3D_q(ii,jj,:) = Timesteppingdata.Fluxlinkage_q;
                FluxLinkage3D_d(ii,jj,:) = Timesteppingdata.Fluxlinkage_d;
                Torque3D(ii,jj,:) = Timesteppingdata.Torque_maxwell;
                Piron_eddy_stator0(ii,jj) = Results.Piron_tooltip.Piron_eddy_stator;
                Piron_eddy_rotor0(ii,jj) = Results.Piron_tooltip.Piron_eddy_rotor;
                Piron_hyst_stator0(ii,jj) = Results.Piron_tooltip.Piron_hyst_stator;
                Piron_hyst_rotor0(ii,jj) = Results.Piron_tooltip.Piron_hyst_rotor;
                MagnetLoss0(ii,jj) = Results.MagnetLoss;
            end
        end
        numberOfProcesses0 = numberOfProcesses;
    end
    pause(1); % wait for 1 second
end

%% check if any simulation failed and try again
for i = 2:length(Irmssupply)
    for j = 1:length(Gamma)
        if isnan(MagnetLoss0(i,j))
            settingsMagnetostatic = setParampm(settingsMagnetostatic,'RMS supply current',Irmssupply(i));  % supply current value
            settingsMagnetostatic = setParampm(settingsMagnetostatic,'Advance angle',Gamma(j));
            [Results,Timesteppingdata] = runMStimestepping(Private,Geometry,Mesh,Windings,settingsMagnetostatic,Materials,SubdomainProperty,[]);
            FluxLinkage3D_a(i,j,:) = Timesteppingdata.Fluxlinkage_a;
            FluxLinkage3D_q(i,j,:) = Timesteppingdata.Fluxlinkage_q;
            FluxLinkage3D_d(i,j,:) = Timesteppingdata.Fluxlinkage_d;
            Torque3D(i,j,:) = Timesteppingdata.Torque_maxwell;
            Piron_eddy_stator0(i,j) = Results.Piron_tooltip.Piron_eddy_stator;
            Piron_eddy_rotor0(i,j) = Results.Piron_tooltip.Piron_eddy_rotor;
            Piron_hyst_stator0(i,j) = Results.Piron_tooltip.Piron_hyst_stator;
            Piron_hyst_rotor0(i,j) = Results.Piron_tooltip.Piron_hyst_rotor;
            MagnetLoss0(i,j) = Results.MagnetLoss;
        end
    end
end
            
%% Step 5: Calculate full iron and magnet losses tables
alfa_stator = Private.MaterialProperties(4).ironloss.alfa;     % hysteresis iron loss Steinmetz coefficient (frequency power) for stator core material
alfa_rotor = Private.MaterialProperties(5).ironloss.alfa;      % hysteresis iron loss Steinmetz coefficient (frequency power) for rotor core material
Piron_eddy_stator = zeros(length(Irmssupply),length(Gamma),length(Speed));
Piron_eddy_rotor = zeros(length(Irmssupply),length(Gamma),length(Speed));
Piron_hyst_stator = zeros(length(Irmssupply),length(Gamma),length(Speed));
Piron_hyst_rotor = zeros(length(Irmssupply),length(Gamma),length(Speed));
MagnetLoss = zeros(length(Irmssupply),length(Gamma),length(Speed));
for n = 1:length(Speed)
    f1 = Speed(n)*nPolePairs/60;
    for i = 1:length(Irmssupply)
        for j = 1:length(Gamma)
            Piron_eddy_stator(i,j,n) = Piron_eddy_stator0(i,j)*(f1^2);
            Piron_eddy_rotor(i,j,n) = Piron_eddy_rotor0(i,j)*(f1^2);
            Piron_hyst_stator(i,j,n) = Piron_hyst_stator0(i,j)*(f1^alfa_stator);
            Piron_hyst_rotor(i,j,n) = Piron_hyst_rotor0(i,j)*(f1^alfa_rotor);
            MagnetLoss(i,j,n) = MagnetLoss0(i,j)*(f1^2);
        end
    end
end
Piron_stator = Piron_eddy_stator+Piron_hyst_stator;
Piron_rotor = Piron_eddy_rotor+Piron_hyst_rotor;
Pcore_rotor = Piron_rotor+MagnetLoss;
%% Step 6: Stator resistance per phase measured on the motor’s terminals 
PhaseResistance = Windings.Rs;

%% Step 7: Prepare and save Simscape model data
N = nPolePairs; % Number of pole pairs
magVec = Ipeak; % Peak current magnitude vector [A]
gammaVec = Gamma; % Current advance angle vector [deg]
angleVec = RotorAngle; % Rotor position vector [deg]
fluxAmat = FluxLinkage3D_a; % A-phase flux linkage matrix [Wb]
fluxDmat = FluxLinkage3D_d; % D-axis flux linkage matrix [Wb]
fluxQmat = FluxLinkage3D_q; % Q-axis flux linkage matrix [Wb]
torqueMat = Torque3D; % Torque matrix [N*m]
Rs = PhaseResistance; % Phase resistance [Ohm]
rpmVec = Speed; % Rotor speed vector [rpm]
rotorCoreLossMat = Pcore_rotor; % Rotor core loss matrix [W]
statorIronLossMat = Piron_stator; % Stator iron loss matrix [W]

if strcmp(getParampm(settingsMagnetostatic,'Simulation time'),'One electrical period')
    % Rotor position vector must cover a range of 360/nPolePairs deg and be cyclic.
    fluxAmat(:,:,end+1) = fluxAmat(:,:,1);
    fluxDmat(:,:,end+1) = fluxDmat(:,:,1);
    fluxQmat(:,:,end+1) = fluxQmat(:,:,1);
    torqueMat(:,:,end+1) = torqueMat(:,:,1);
end

% Current advance angle vector must cover a range of 360 deg and be cyclic.
gammaVec(end+1) = 360;
fluxAmat(:,end+1,:) = fluxAmat(:,1,:);
fluxDmat(:,end+1,:) = fluxDmat(:,1,:);
fluxQmat(:,end+1,:) = fluxQmat(:,1,:);
torqueMat(:,end+1,:) = torqueMat(:,1,:);
rotorCoreLossMat(:,end+1,:) = rotorCoreLossMat(:,1,:);
statorIronLossMat(:,end+1,:) = statorIronLossMat(:,1,:);

save(save2filename,'N','magVec','gammaVec','angleVec','fluxAmat','fluxDmat','fluxQmat','torqueMat','Rs','rpmVec','rotorCoreLossMat','statorIronLossMat');





