% Parametric sweep script example (see section 4.3 of User Manual - Part II)
% Plot efficiency of the motor versus Tooth Width and Slot Depth
cd ../
initMXPscriptingAPI
file = getpmpath('SimFiles\example_innerrotor\example_innerrotor.mxp');   % initial mxp-project file
[motorProps,settingsMagnetostatic] = openMXP(file);    % open mxp-project
% motorProps - structure with motor model parameters (geometry, winding, mesh)
% settingsMagnetostatic - magnetostatic FEA settings

ToothWidth = [3.5 4 4.5 5 5.5 6];        % sweep through a range of Tooth Width values
SlotDepth = [15 17 19 21 23 25];         % sweep through a range of Slot Depth values

targetTorque = 100;  % N*m - all designs to be compared with the same output torque 100 N*m and speed 4500 RPM
initCurrent = 200;   % Arms - initial supply current value
settingsMagnetostatic=setParampm(settingsMagnetostatic,'Save each step solution',0);
settingsMagnetostatic=setParampm(settingsMagnetostatic,'Compute cogging torque',0);
settingsMagnetostatic = setParampm(settingsMagnetostatic,'Convergence tolerance',0.01);
settingsMagnetostatic = setParampm(settingsMagnetostatic,'Current waveform','Sinusoidal');
settingsMagnetostatic = setParampm(settingsMagnetostatic,'Advance angle',0);
settingsMagnetostatic = setParampm(settingsMagnetostatic,'Simulation time','One electrical period');
settingsMagnetostatic = setParampm(settingsMagnetostatic,'Mechanical speed',4500);           % RPM

motorProps = setParampm(motorProps,'Wire size method','Fill factor');
motorProps = setParampm(motorProps,'Coil fill factor',0.5);                                  % keep the same coil fill factor for all the designs
motorProps = setParampm(motorProps,'End winding inductance','Auto');                         % end winding inductance is automatically calculated
motorProps = setParampm(motorProps,'Phase resistance','Auto');                               % phase resistance is automatically calculated

% initialize with NaN values
Efficiency = NaN*ones(length(ToothWidth),length(SlotDepth));   
SupplyCurrent = NaN*ones(length(ToothWidth),length(SlotDepth)); 
PhaseResistance = NaN*ones(length(ToothWidth),length(SlotDepth)); 
StatorLosses = NaN*ones(length(ToothWidth),length(SlotDepth)); 
IronLosses = NaN*ones(length(ToothWidth),length(SlotDepth)); 
MagnetLosses = NaN*ones(length(ToothWidth),length(SlotDepth)); 
tic
for i=1:length(ToothWidth)
    for j=1:length(SlotDepth)
        disp(['Tooth Width: ' num2str(ToothWidth(i)) '; Slot Depth: ' num2str(SlotDepth(j))]);
        motorProps=setParampm(motorProps,'Tooth Width',ToothWidth(i));                         % Tooth Width
        motorProps=setParampm(motorProps,'Slot Depth',SlotDepth(j));                   % Slot Depth
        % assemble mxp-project
        [Geometry, Windings, Mesh, Materials, SubdomainProperty, Error, motorProps_output, Weight] = assembleMXP(file,motorProps);
        if isempty(Error) % if the mxp-project was successfully assembled
            % Determine current required to produce target torque
            iter = 1;
            Private = [];   % Private will be calculated with the first call of runMStimestepping
            supplyCurrentRMS = initCurrent;
            settingsMagnetostatic = setParampm(settingsMagnetostatic,'Number of points',1);                   % define single rotor position simulation setup
            settingsMagnetostatic = setParampm(settingsMagnetostatic,'RMS supply current',supplyCurrentRMS);  % initial current value
            while 1
                % run magnetostatic analysis (with single rotor position simulation setup - see line 48)
                [Results,Timesteppingdata,field,~,~,Private] = runMStimestepping(Private,Geometry,Mesh,Windings,settingsMagnetostatic,Materials,SubdomainProperty,0);
                if 100*abs((Results.torque-targetTorque)/targetTorque)<0.5             % 0.5% torque error tolerance
                    break
                end
                iter = iter+1;
                if iter>10
                    break         % limit maximum number of iterations
                end
                supplyCurrentRMS = (targetTorque/Results.torque)*supplyCurrentRMS;     % recalculate supply current
                settingsMagnetostatic = setParampm(settingsMagnetostatic,'RMS supply current',supplyCurrentRMS);
            end
            nPolePairs = getParampm(motorProps,'Number of pole pairs');
            Ns = getParampm(motorProps,'Number of slots');
            numberOfRotorPositions = setNumberOfPoints('Low accuracy',nPolePairs,Ns);
            settingsMagnetostatic = setParampm(settingsMagnetostatic, 'Number of points', numberOfRotorPositions);
            % run magnetostatic analysis (with multiple rotor position simulation setup - see line 66)
            [Results,Timesteppingdata] = runMStimestepping(Private,Geometry,Mesh,Windings,settingsMagnetostatic,Materials,SubdomainProperty,0);  
            if Results.SolutionConverged
                Efficiency(i,j) = Results.Efficiency;
                SupplyCurrent(i,j) = supplyCurrentRMS; 
                PhaseResistance(i,j) = Windings.Rs;
                StatorLosses(i,j) = Results.Ps;
                IronLosses(i,j) = Results.Piron; 
                MagnetLosses(i,j) = Results.MagnetLoss;
            end
        end
    end
end
toc
save('ExParSw.mat')
figure('Name','Efficiency vs. Tooth Width and Slot Depth'); surf(ToothWidth,SlotDepth,Efficiency); xlabel('Tooth Width'); ylabel('Slot Depth'); zlabel('Efficiency');

