% This MATLAB script was used to extract the PMSM flux linkage data and
% other parameters from Motor-CAD that are required by example model
% ee_import_fem_motorcad.slx. Lines that call Motor-CAD are commented
% out and return data loaded from a previously-created mat file.

% Copyright 2018 The MathWorks, Inc. Portions Copyright Motor Design Ltd 2018.

%% Start MotorCAD ActiveX server and load design data
%mcad = actxserver('motorcad.appautomation');
%invoke (mcad,'LoadFromFile','C:\PMSM.mot' );
clear;

% Clear workspace
close all
clear
clc

% Load model parameters
Ts                  = 5e-7;                         % [s] Model sampling time (200 KHz)
Ts_ctrl             = 2e-5;                     % [s] Controller sampling time (50us)5e-5
f_speed_ctrl        = 500;                          % [Hz] Speed/torque Controller frequency = (500 Hz)
Ts_Spd_ctl          = 1/f_speed_ctrl;
PWM_Cnt             = 6000;
PWM_Half_Cnt        = PWM_Cnt/2;
DeatTime = 0;
DT_PWM_CNT = PWM_Cnt/(Ts_ctrl)*DeatTime;
Torque_Mode = 1;
Speed_Mode = 2;
Ctrl_Mode = Torque_Mode;
%Simulink provider Motor parameters
n_polePairs  = 4;        % [-] Number of motor pole pairs
PM           = 0.01688;   % Permanent magnet flux linkage, 
Ld           = 0.7e-4/1.4;%0.07e-3;% d-axis inductance, 
Lq           = 0.178e-3/1.4; % q-axis inductance,
Rs           = 0.011;   % Stator resistance,
J            = 0.03945; % Moment of inertia,
bandwith     = 1600*2*pi;
i_Udc        = 96;      % DCbus max voltage
Rload = 1e6;

% A1 6匝电机数据
%Current=xlsread('E:\works\项目\MC100\电机仿真数据\PhiD-I230405.xlsx','Sheet1','B2:N2');
%FluxD = xlsread('E:\works\项目\MC100\电机仿真数据\PhiD-I230405.xlsx','Sheet1','B4:N40');
%FluxQ = xlsread('E:\works\项目\MC100\电机仿真数据\PhiQ-I230405.xlsx','Sheet1','B4:N40');
%Torque = xlsread('E:\works\项目\MC100\电机仿真数据\T-I230405.xlsx','Sheet1','B4:N40');

% A1 5匝电机数据
% FluxD = xlsread('E:\works\项目\MC100\电机仿真数据\5N-phiD-Current.xlsx','Sheet1','B4:N40');
% FluxQ = xlsread('E:\works\项目\MC100\电机仿真数据\5N-phiQ-Current.xlsx','Sheet1','B4:N40');
% Torque = xlsread('E:\works\项目\MC100\电机仿真数据\5N-Torque-Current.xlsx','Sheet1','B4:N40');

% 五菱电机数据
% FluxD = xlsread('E:\works\项目\MC100\电机仿真数据\96-9000\M155-phiD-IS.xlsx','Sheet1','B4:N40');
% FluxQ = xlsread('E:\works\项目\MC100\电机仿真数据\96-9000\M155-phiQ-IS.xlsx','Sheet1','B4:N40');
% Torque = xlsread('E:\works\项目\MC100\电机仿真数据\96-9000\M155-Torque-IS.xlsx','Sheet1','B4:N40');
%% Define operating point and flux linkage table breakpoints
PointsPerCycle = 36;    % Number of angle data points
Shaft_Speed = 3500;     % Nominal shaft speed (RPM)
DCBusVoltage = 96;     % DC bus voltage
magVec   = linspace(0,300,13);                  % Current magnitude vector
gammaVec  = linspace(0,360,PointsPerCycle+1); % Current advance angle
angleVec = linspace(0,360/n_polePairs,PointsPerCycle+1);  % Rotor angle vector

nMag = length(magVec);
nGamma = length(gammaVec);
nAngle = length(angleVec);
fluxDmat = zeros(nMag,nGamma,nAngle);
fluxQmat = zeros(nMag,nGamma,nAngle);
torqueMat = zeros(nMag,nGamma,nAngle);
% for i = 1:nMag
%    for j = 1:nGamma
%        for m = 1:nAngle
%            fluxDmat(i,j,m)    = FluxD(j, i);
%            fluxQmat(i,j,m)    = FluxQ(j, i);
%            torqueMat(i,j,m)   = Torque(j, i);
%        end
%    end
% end
load A1_N5_motor_flux_data
%load A1_motor_flux_data
%load M155_motor_flux_data
fluxDmat=fluxDmat.*0.707;
fluxQmat=fluxQmat.*0.707;
torqueMat=torqueMat.*0.707;
%% Extract flux linkdage data from Motor-CAD
% nMag = length(magVec);
% nGamma = length(gammaVec);
% nX = length(angleVec);
% fluxDmat = zeros(nMag,nGamma,nX);
% fluxQmat = zeros(nMag,nGamma,nX);
% torqueMat = zeros(nMag,nGamma,nX);
% fluxAmat = zeros(nMag,nGamma,nX);
% for i = 1:nMag
%     for j = 1:nGamma
%         invoke(mcad,'SetVariable','PeakCurrent',magVec(i));
%         invoke(mcad,'SetVariable','DCBusVoltage',DCBusVoltage);
%         invoke(mcad,'SetVariable','Shaft_Speed_[RPM]',Shaft_Speed);
%         invoke(mcad,'SetVariable','TorqueCalculation',1);
%         invoke(mcad,'SetVariable','TorquePointsPerCycle',PointsPerCycle);
%         invoke(mcad,'SetVariable','PhaseAdvance',gammaVec(j));
%         invoke(mcad,'DoMagneticCalculation');
%         for m = 1:nX % Loop over rotor angles
%             xe = N*angleVec(m)*pi/180; % Electrical angle in radians
%             [res1, x1, torque] = invoke(mcad, 'GetMagneticGraphPoint', 'TorqueVW',m-1);
%             [res2, x2, fluxD]  = invoke(mcad, 'GetMagneticGraphPoint', 'FluxLinkageLoadTotalD',m-1); 
%             [res3, x3, fluxQ]  = invoke(mcad, 'GetMagneticGraphPoint', 'FluxLinkageLoadTotalQ',m-1); 
%             fluxDmat(i,j,m)    = fluxD;
%             fluxQmat(i,j,m)    = fluxQ;
%             torqueMat(i,j,m)   = torque;
%             % Optional calculation of A-phase flux linkage. Note Park's
%             % transform for MotorCAD is Q leads D, rotor angle measured
%             % from A-phase flux axis to Q-axis.
%             invP = [sin(xe+[0;-2*pi/3;2*pi/3]),cos(xe+[0;-2*pi/3;2*pi/3])];
%             fluxes = invP*[fluxD;fluxQ];
%             fluxAmat(i,j,m) = fluxes(1);
%         end
%     end
% end