function varargout=mcbPMSMCharacteristics(pmsm, inverter, varargin)
% mcbPMSMCharacteristics   Draw the characteristic curves for PMSM
%   mcbPMSMCharacteristics(pmsm, inverter) accepts 2 required arguments.
%
%   mcbPMSMCharacteristics plots the current limit circle,
%    voltage limit curve, MTPA curve, constant torque curve, and MTPV
%    curve in the id-iq space. The function also plots the torque-speed,
%    power-speed, current-speed curves.
%
%   The function accepts following inputs: a pmsm object, inverter object
%   and optional name-value pairs. The possible name value pairs
%   are the following (default is mentioned first):
%       'voltageEquation'         -> 'actual' or 'approximate'
%       'speed'                   -> speed in rpm
%       'torque'                  -> load torque in N.m.
%       'driveCharacteristics'    -> 0 or 1 or 2
%       'constraintCurves'        -> 1 or 0
%       'FWCMethod'               -> 'vclmt' or 'none' or 'cvcp' or 'cccp'
%       'opacity'                 -> 1 or (1 to 0)
%       'imax'                    -> current in Amperes
%       'idqExternal'             -> [id, iq] array in Amperes, Nx2 matrix
%                                  or [id, iq, wrpm] array, Nx3 matrix.
%
%    The speed(rpm) and torque(Nm) can be provided with optional name-value
%   arguments 'speed' and 'torque' respectively. Default values for the
%   speed and torque are maximum speed and 0Nm respectively.
%    The viscous frictional load torque is added to the final load torque.
%   Static frictional load torque is not added to the computation.
%    The optional argument 'constraintCurves' followed by 0 or 1 either
%   disables or enables plotting the PMSM constraint curves in the id-iq space.
%    The optional argument 'driveCharacteristics' followed by 0 or 1 either
%   disables or enables plotting Torque-Speed, Power-Speed, Currents-Speed
%   curves for the motor-inverter combination.
%    In a case where both the driveCharacteristics and constraintCurves
%   are disabled, constraintCurves is set to 1 for plotting.
%    To use the actual voltage equations, use 'voltageEquation' and 'actual'
%   as the optional name value pair arguments. To use approximate voltage
%   equations, use 'approximate' as the value.
% 		actual equations:
% 			vd=id*r-w*Lq*iq; vq=iq*r+w*Ld*id+w*fluxPM; Vmax=Vdc/sqrt(3);
% 		approximate equations:
% 			vd=-w*Lq*iq; vq=w*Ld*id+w*fluxPM; Vmax=(Vdc/sqrt(3))-irated*r;
%    Another optional value pair argument with 'opacity' and a value
%   between 0 and 1 makes sure that multiple constraint curves can be drawn in
%   the same figure. A value of 1 ensures that the plot is drawn on a clear
%   figure window. Any value less than 1 holds the earlier plot and draws
%   again with reduced intensity of the same colors. for ex:
%   mcbPMSMCharacteristics(pmsm,inverter,'speed',4000,'opacity',0.6)
%    Another optional argument for SPMSM is 'FWCMethod' followed by
%   'cvcp' implements constant-voltage-constant-power method for FW region,
%   'cccp' implements constant-current-constant-power method for FW region,
%   'none' implements no field weakening. The (default) optional value
%   'vclmt' implements voltage-current-limited-maximum-torque method
%   for FW region. The option is set to vclmt for the IPMSM type machines.
%    The optional argument 'imax' followed by the maximum value of the
%   current in amperes plots the drive characteristics for the maximum
%   current in addition to the rated current. This value is expected
%   to be higher than the rated current to enable additional plotting.
%    The optional argument idqExternal followed by an array of Nx2 matrix
%   consisting of id, iq values provided by the user, will plot the id,iq
%   in the id-iq plane. This will be helpful in comparing test data with
%   computed data. Another possibility for idqExternal is to provide id,
%   iq, and wrpm column matrices where the wrpm corresponds to the id iq at
%   the given currents id, iq.
%    The non-linear characteristics of the PMSM can be used for plotting
%   the constraint curves. Use pmsm.PMSMLUT to provide the non-linear
%   parameter fields, idVec and iqVec in combination with LdTable, LqTable,
%   and FluxPMTable, or FluxDTable and FluxQTable. The field 'method'
%   selects the input method for motor parameters, can be assigned as
%   'lumped' or 'Ldq' or 'FluxDQ'.
%
%   Example:
%       inverter = mcb_SetInverterParameters( 'BoostXL-DRV8305');
%       pmsm     = mcb_SetPMSMMotorParameters('BLY171D');
%       mcbPMSMCharacteristics(pmsm,inverter,"voltageEquation","actual", ...
%           "opacity",1,"speed",5000,"torque",0, ...
%           "constraintCurves",1,"driveCharacteristics",1, ...
%           "FWCMethod","vclmt","imax",6, ...
%           "idqExternal",[0 0;0 0.1; 0 0.2; -0.1 0.2;-0.2 0.2])
%
%  See also mcbPMSMSpeeds

%   Copyright 2022-2024 The MathWorks, Inc.

%%



p = inputParser;
addRequired(p,'pmsm',@(x) (isstruct(x)&&min(isfield(x,{'p','Rs','Ld','Lq','FluxPM','B','I_rated'}))) || (all(cellfun(@(y) isprop(x,y),{'p','Rs','Ld','Lq','FluxPM','B','I_rated'}))));
addRequired(p,'inverter',@(x) isstruct(x)&&min(isfield(x,{'V_dc'})));
addParameter(p,'voltageEquation','actual',@(x) any(ismember(lower(x) ...
    ,{'actual','approximate'})));
addParameter(p,'speed',1i,@(x) min(isnumeric(x))&&min(isfinite(x)));
addParameter(p,'torque',0,@(x) min(isnumeric(x))&&min(isfinite(x)));
addParameter(p,'driveCharacteristics',0,@(x) any(ismember(x ,[2,1,0])));
addParameter(p,'constraintCurves',1,@(x) any(ismember(x ,[1,0])));
addParameter(p,'drawContours',0,@(x) any(ismember(x ,[1,0])));
addParameter(p,'changeDQAxes',0,@(x) any(ismember(x ,[1,0])));
addParameter(p,'opacity',1,@(x) (x>0 && x<=1));
addParameter(p,'FWCMethod','vclmt',@(x) any(ismember(lower(x) , ...
    {'vclmt','none','cccp','cvcp'})));
addParameter(p,'imax',0,@(x) isnumeric(x)&&isfinite(x));
addParameter(p,'idqExternal',[],@(x) isnumeric(x)&&min(min(isfinite(x))));

parse(p,pmsm,inverter,varargin{:});
warning('off','MATLAB:structOnObject');
pmsm=struct(p.Results.pmsm);
warning('on','MATLAB:structOnObject');
inverter=p.Results.inverter;
if ~isfield(inverter,'R_board')
    inverter.R_board=0;
end

motorSystem=pmsm;
motorSystem.Vdc=inverter.V_dc;
motorSystem.Rboard=inverter.R_board;

radps2rpm=30/pi;
elec2mech=1/pmsm.p;

if ~isfield(pmsm,'model')
    pmsm.model='pmsm';
end
if isfield(pmsm,'oldmotorType')
    pmsm.model=pmsm.oldmotorType;
end

voltageEquation = p.Results.voltageEquation;

opacity=p.Results.opacity;
if opacity<0.1
    opacity=0.1;
end

FWCMethod = lower(p.Results.FWCMethod);
%if ~isfield(pmsm,'N_max')
[pmsm.T_rated,PMSMspeeds,ID_milestone,IQ_milestone]=mcbPMSMRatedTorque(pmsm,inverter,'FWCMethod',FWCMethod, ...
    'voltageEquation',voltageEquation);
pmsm.N_max=PMSMspeeds(end);
pmsm.N_base=PMSMspeeds(1);
motorSystem.N_max=pmsm.N_max;
motorSystem.N_base=pmsm.N_base;

%end

Very_high_speed=200000;

wrpm=p.Results.speed;
T_load=p.Results.torque;

if length(T_load)~= length(wrpm)
    % if length(T_load)==1
    %     T_load=T_load*ones(1,length((wrpm)));
    if length(wrpm)==1 && wrpm==1i%#ok
        wrpm=pmsm.N_max*ones(1,length(T_load));
    end
end
if all(imag(wrpm)==0)
    wrpm=(wrpm);
else
    wrpm=pmsm.N_max;
    % wrpm=mcbPMSMMaxSpeed(pmsm,inverter,'FWCMethod',FWCMethod, ...
    %     'voltageEquation',voltageEquation);
    if wrpm>Very_high_speed % if the speeds are very high due to MTPV, set a practical limit.
        wrpm=Very_high_speed;
    end
end
drawContours=p.Results.drawContours;
changeDQAxes=p.Results.changeDQAxes;


driveCharacteristics=p.Results.driveCharacteristics;

constraintCurves=p.Results.constraintCurves;

% the imax will be used to plot the drive characteristics for the
% maximum current. When no input is provided, or when a imax lower than
% rated current is provided, the drive characteristics at this current is
% not plotted.
imax=p.Results.imax;
if imax==0
    flag_imax=0;
elseif imax<=pmsm.I_rated
    flag_imax=0;
else
    flag_imax=1;
end

idqCustom=p.Results.idqExternal;
idqsize=size(idqCustom);
if isempty(idqCustom)
    flag_idq=0;
elseif idqsize(2)==2 || idqsize(2)==3
    flag_idq=1;
else
    flag_idq=0;
end


if strcmpi(FWCMethod,'none')==1
    FWCMethodtext='no';
else
    FWCMethodtext=upper(FWCMethod);
end


% setting at least one of the plots when both options are set to zero.
if driveCharacteristics==0 && constraintCurves==0 && nargout==0
    %constraintCurves=1;
    drawContours=1;
end
motorSystem.wrpm=wrpm;

motorSystem.Tload=T_load;

motorType = 'ipmsm';
% checking if the pmsm has motorType field
if isfield(pmsm,'motorType')
    % read the motorType, if it is defined already
    motorType=pmsm.motorType;
end
saliency=max(motorSystem.Lq,motorSystem.Ld)/min(motorSystem.Lq,motorSystem.Ld);
% Setting to SPMSM when the d, q axes inductances are equal
if saliency<1.01
    motorType='spmsm';
end


motorSystem.motorType=motorType;
% the retardation torque = Bv*w (w in mech radps)
Bv=motorSystem.B;
if Bv<=0
    Bv=eps;
end
% All the following calculations are done using the Bv*w is T with w as
% electrical radps. The Bv from user input is Nm/mech.radps.
%To be able to represent the viscous torque correctly in the following
%equations, Bv is converted into equivalent units. in this case,
% new Bv becomes Bv(new)=Bv*elec2mech;
% now frictional load torque Tfriction= Bv(new) * elec.radps.
Bv=Bv*elec2mech;

%checking the stack to know if this function is called from its counterpart. If it
%is called so, do not print the LdLq warning.
stk = dbstack;
calledFromAnother=0;
if length(stk)>1
    if strcmp(stk(2).name,'mcbPMSMSpeeds') || strcmp(stk(2).name,'PMSMSpeedComputation')
        calledFromAnother=1;
    end
end
checkLdqSwap=1;
if isfield(pmsm,'PMSMLUT')
    if isfield(pmsm.PMSMLUT,'method')
        if strcmpi(pmsm.PMSMLUT.method,'Ldq') || strcmpi(pmsm.PMSMLUT.method,'FluxDQ')
            checkLdqSwap=0;
        end
    end
end

Irated=motorSystem.I_rated;
%Setting the higher inductance in Lq, and lower inductance in Ld.
if checkLdqSwap==1
    if motorSystem.Lq<motorSystem.Ld
        Lq=motorSystem.Ld;
        Ld=motorSystem.Lq;
        % if the mcbPMSMSpeeds function has already been called, this message
        % would have been displayed there
        if calledFromAnother==0
            disp(message('mcb:blocks:APILqLessThanLdSwapping').getString());
        end
    else
        Ld=motorSystem.Ld;
        Lq=motorSystem.Lq;
    end
else
    if strcmpi(pmsm.PMSMLUT.method,'Ldq') 
        Ld=mean(mean(pmsm.PMSMLUT.LdTable));
        Lq=mean(mean(pmsm.PMSMLUT.LqTable));
        pmsm.Ld=Ld;
        pmsm.Lq=Lq;
        try
            fluxPM=mean(mean(pmsm.PMSMLUT.FluxPMTable));
            pmsm.FluxPM=fluxPM;
        catch
            pmsm.PMSMLUT.FluxPMTable=0*pmsm.PMSMLUT.LqTable;
            pmsm.FluxPM=0;
        end
    else
        [Ld_derived,Lq_derived,fluxPM_derived]=mcb.internal.PMSMFluxDQ2Ldq(pmsm.PMSMLUT);
        fluxPM=mean(fluxPM_derived,'all');
        Ld=mean(Ld_derived,'all');
        Lq=mean(Lq_derived,'all');
        pmsm.Ld=Ld;
        pmsm.Lq=Lq;
        pmsm.FluxPM=fluxPM;
    end
end
R=motorSystem.Rs;

fluxPM=motorSystem.FluxPM;
Pp=motorSystem.p;

I_short=inverter.V_dc/sqrt(3)/R;

if Irated>(I_short)
    error(message('mcb:blocks:APIHighIRated',num2str(I_short)));
end

% creating a second motor variable that includes the imax as rated current.
% this will be used to generate the drive characteristics at imax.
if flag_imax==1 && driveCharacteristics>0
    pmsm1=pmsm;
    if imax>=I_short
        flag_imax=1;
        warning(message('mcb:blocks:APIHighIMax',num2str(floor(I_short))));
        imax=floor(I_short);
    end
    pmsm1.I_rated=imax;
end

goLUT=0;
if isfield(pmsm,'PMSMLUT')
    if isfield(pmsm.PMSMLUT,'method')
        if strcmpi(pmsm.PMSMLUT.method,'lumped')
            startval.goLUT=0;
        elseif strcmpi(pmsm.PMSMLUT.method,'Ldq')
            startval.goLUT=1;
        elseif strcmpi(pmsm.PMSMLUT.method,'FluxDQ')
            startval.goLUT=2;
        else
            error('Please use either ''lumped'' or ''Ldq'' or ''FluxDQ'' for pmsm.PMSMLUT.method');
        end
    else
        startval.goLUT=0;
    end
    goLUT=startval.goLUT;
end

if goLUT==2 && strcmpi(motorType,'spmsm')
    [Ld_derived,Lq_derived,fluxPM_derived]=mcb.internal.PMSMFluxDQ2Ldq(pmsm.PMSMLUT);
    pmsm.PMSMLUT.LdTable=Ld_derived;
    pmsm.PMSMLUT.LqTable=Lq_derived;
    pmsm.PMSMLUT.FluxPMTable=fluxPM_derived;
    pmsm.PMSMLUT.method='Ldq';
    startval.goLUT=1;
    goLUT=1;
end


if goLUT==1
    if ~isfield(pmsm.PMSMLUT,'TorqueTable')
        for i1=1:length(pmsm.PMSMLUT.idVec)
            for j1=1:length(pmsm.PMSMLUT.iqVec)
                TorqueTable(i1,j1)=1.5*pmsm.p*pmsm.PMSMLUT.iqVec(j1)*(pmsm.PMSMLUT.FluxPMTable(i1,j1)+(pmsm.PMSMLUT.LdTable(i1,j1)-pmsm.PMSMLUT.LqTable(i1,j1))*pmsm.PMSMLUT.idVec(i1));%#ok
            end
        end
        pmsm.PMSMLUT.TorqueTable=TorqueTable;
    end
end

if goLUT==2
    if ~isfield(pmsm.PMSMLUT,'TorqueTable')
        iq0index=find(pmsm.PMSMLUT.iqVec==0);
        if ~isempty(iq0index)
            for i1=1:length(pmsm.PMSMLUT.idVec)
                pmsm.PMSMLUT.FluxQTable(i1,:)=pmsm.PMSMLUT.FluxQTable(i1,:)-pmsm.PMSMLUT.FluxQTable(i1,iq0index);
            end
        end
        for i1=1:length(pmsm.PMSMLUT.idVec)
            for j1=1:length(pmsm.PMSMLUT.iqVec)
                TorqueTable(i1,j1)=1.5*pmsm.p*(pmsm.PMSMLUT.iqVec(j1)*pmsm.PMSMLUT.FluxDTable(i1,j1)-pmsm.PMSMLUT.idVec(i1)*pmsm.PMSMLUT.FluxQTable(i1,j1));
            end
        end
        pmsm.PMSMLUT.TorqueTable=TorqueTable;
    end
end

%
%%
% following code plots the drive characteristics of the motor
% driveCharacteristics set to 1 will plot the torque-speed, power-speed,
% id/iq-speed and ipeak-speed plots.
% when the constraintCurves is set to 1 along with driveCharacteristics to
% 1, the constraint curves will include additional curve that represents
% the control trajectory in the id-iq plane.
% driveCharacteristics set to 2 will plot the above, and also plots the
% voltage limit circle and the vd, vq trajectory in the vd-vq space.

CreateNewFigure=0;
if driveCharacteristics && opacity==1
    CreateNewFigure=1;
else
    availablefigHandles = findobj('Type', 'figure','Tag','DriveCharacteristics');
    %Picking the last created constraint curve for plotting
    if isempty(availablefigHandles)
        CreateNewFigure=1;
    else
        DrvCharFigH=availablefigHandles(1);
    end
end
if driveCharacteristics>=1 || nargout==1
    if CreateNewFigure==1 && driveCharacteristics>=1
        DrvCharFigH=figure('Name',strcat("PMSM : Drive Characteristics of ", ...
            pmsm.model), 'numberTitle','off','Tag','DriveCharacteristics');
    end
    if strcmpi(motorType,'spmsm')
        % if the motor is SPMSM type, setting the Ld and Lq equal.
        Lq=Ld;
    else
        % if the motor is IPMSM type Ld, Lq can be different.
        % if the FWCMethod is set to anything else, it is changed to VCLMT
        if strcmpi(FWCMethod,'vclmt')==0
            FWCMethod='vclmt';
            FWCMethodtext=FWCMethod;
            disp(message('mcb:blocks:APIForceSetVCLMT').getString());
        end
    end
    % this loop will be run once when imax is not provided
    % this loop will run twice when imax is provided
    for pmsmloop=1:(flag_imax+1)
        if pmsmloop==1
            motor=pmsm;% pick the original motor in the first run
        elseif pmsmloop==2
            motor=pmsm1;% pick the motor with imax in the second run
            [motor.T_rated,PMSMspeeds,ID_milestone,IQ_milestone]=mcbPMSMRatedTorque(pmsm1,inverter,'FWCMethod',FWCMethod, ...
                'voltageEquation',voltageEquation);
        end
        % [CurrentSpeedArray]=mcbPMSMSpeeds(motor,inverter,'outputAll',1, ...
        %     'voltageEquation',voltageEquation, ...
        %     'constraintCurves',0,'FWCMethod',FWCMethod);
        % base_rpm=round(CurrentSpeedArray(3,1)*(30/pi)/pmsm.p);
        base_rpm=PMSMspeeds(1);
        % mcbPMSMSpeeds with outputAll set to 1 returns an array of 3Xn size.
        % the top row is id, second row is iq, third row is w_radps. each
        % column represents the specific speed limestone. For example: corner
        % speed, speed after which the control path changes, or maximum speed
        %count=size(CurrentSpeedArray,2);
        count=length(PMSMspeeds);
        if count==1%#ok
            PMSMspeeds=[PMSMspeeds PMSMspeeds];%#ok
            ID_milestone=[ID_milestone ID_milestone];%#ok
            IQ_milestone=[IQ_milestone IQ_milestone];%#ok
            count=2;
        end
        %w_max=floor(CurrentSpeedArray(3,count)*radps2rpm*elec2mech);
        w_max=PMSMspeeds(end);
        if w_max>Very_high_speed % if the speeds are very high due to MTPV, set a practical limit.
            w_max=Very_high_speed;
        end



        if isfield(pmsm,'PMSMLUT') && isfield(pmsm.PMSMLUT,'wrpmVec')
            if size(pmsm.PMSMLUT.wrpmVec,2)==1
                pmsm.PMSMLUT.wrpmVec=pmsm.PMSMLUT.wrpmVec';
            end
            wArrayZero2Base=pmsm.PMSMLUT.wrpmVec(pmsm.PMSMLUT.wrpmVec<=PMSMspeeds(1));
            wArrayBase2MS=pmsm.PMSMLUT.wrpmVec(pmsm.PMSMLUT.wrpmVec<=PMSMspeeds(2) & pmsm.PMSMLUT.wrpmVec>PMSMspeeds(1));
            if count==3
                wArrayMS2Max=pmsm.PMSMLUT.wrpmVec(pmsm.PMSMLUT.wrpmVec>PMSMspeeds(2));
            else
                wArrayMS2Max=[];
            end
        else
            wArrayZero2Base=unique(round(linspace(1,PMSMspeeds(1),4)));
            if count==3
                wArrayBase2MS=unique(round(linspace(PMSMspeeds(1),PMSMspeeds(2),31)));
                wArrayBase2MS=wArrayBase2MS(2:end);
                wArrayMS2Max=unique(round(linspace(PMSMspeeds(2),w_max,31)));
                wArrayMS2Max=wArrayMS2Max(2:end);
            else
                wArrayBase2MS=unique(round(linspace(PMSMspeeds(1),PMSMspeeds(2),61)));
                wArrayBase2MS=wArrayBase2MS(2:end);
                wArrayMS2Max=[];
            end

        end

        wArray_Final=[wArrayZero2Base wArrayBase2MS wArrayMS2Max];


        TorqueArray=0*wArray_Final;
        PowerArray=0*wArray_Final;
        idArray=0*wArray_Final;
        iqArray=0*wArray_Final;

        %id_corner=CurrentSpeedArray(1,1);
        id_corner=ID_milestone(1);
        %iq_corner=CurrentSpeedArray(2,1);
        iq_corner=IQ_milestone(1);
        %w_corner=floor(CurrentSpeedArray(3,1)*radps2rpm*elec2mech);%mechrpm
        w_corner=PMSMspeeds(1);
        Lw1=length(wArrayZero2Base);
        Lw2=length(wArrayBase2MS);
        Lw3=length(wArrayMS2Max);
        Lw12=Lw1+Lw2;
        Lw123=Lw1+Lw2+Lw3;

        % second milestone speed. This may be max speed for few FW methods,
        % or an intermediate speed after which the control path changes
        % direction. This is termed w_c.

        %id_milestone2=CurrentSpeedArray(1,2);
        id_milestone2=ID_milestone(2);
        %iq_milestone2=CurrentSpeedArray(2,2);
        iq_milestone2=IQ_milestone(2);
        %w_c=floor(CurrentSpeedArray(3,2)*radps2rpm*elec2mech);%mechrpm
        % w_c=PMSMspeeds(2);

        if goLUT==0
            Tcorner=1.5*Pp*(fluxPM*iq_corner+(Ld-Lq)*id_corner*iq_corner);
        else
            Tcorner=local_interp2(pmsm.PMSMLUT.iqVec,pmsm.PMSMLUT.idVec,pmsm.PMSMLUT.TorqueTable,iq_corner,id_corner);
        end
        TorqueArray(1:Lw1)=Tcorner;
        PowerArray(1:Lw1)=TorqueArray(1:Lw1).*wArrayZero2Base/(radps2rpm);
        idArray(1:Lw1)=ones(1,Lw1)*id_corner;
        iqArray(1:Lw1)=ones(1,Lw1)*iq_corner;

        warray=wArrayBase2MS; % mech rpm
        warrayradps=warray/(radps2rpm*elec2mech);%electrical radps
        if strcmpi(FWCMethod,'cvcp')==1
            startval.id=0;startval.iq=0;startval.w=warrayradps;
            startval.wcorner=w_corner/(radps2rpm*elec2mech);
            %solving for Finding id, iq coordinates for CVCP curve of SPMSM
            % after corner speed, till current touches current limit circle
            % this function returns array of operating points[id;iq;wradps]
            operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(motor, ...
                inverter,startval,26,'outputAll',1, ...
                'voltageEquation',voltageEquation);
            idArray(Lw1+1:Lw12)=operatingpoint(1,:);
            iqArray(Lw1+1:Lw12)=operatingpoint(2,:);
            TorqueArray(Lw1+1:Lw12)=1.5*Pp*(fluxPM*iqArray(Lw1+1:Lw12) ...
                + (Ld-Lq)*idArray(Lw1+1:Lw12).*iqArray(Lw1+1:Lw12));
            PowerArray(Lw1+1:Lw12)=TorqueArray(Lw1 ...
                +1:Lw12).*warrayradps*elec2mech; % torque * mech radps
        elseif strcmpi(FWCMethod,'cccp')==1
            startval.id=0;startval.iq=0;startval.w=warrayradps;
            startval.wcorner=w_corner/(radps2rpm*elec2mech);
            % solving for finding id, iq coordinates for CCCP curve of SPMSM
            % after corner speed and till vd, vq touch voltage limit circle
            % in the vd-vq plane
            % this function returns array of operating points[id;iq;wradps]
            operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(motor, ...
                inverter,startval,29,'outputAll',1, ...
                'voltageEquation',voltageEquation);
            idArray(Lw1+1:Lw12)=operatingpoint(1,:);
            iqArray(Lw1+1:Lw12)=operatingpoint(2,:);
            TorqueArray(Lw1+1:Lw12)=1.5*Pp*(fluxPM*iqArray(Lw1+1:Lw12) ...
                + (Ld-Lq)*idArray(Lw1+1:Lw12).*iqArray(Lw1+1:Lw12));
            PowerArray(Lw1+1:Lw12)=TorqueArray(Lw1 ...
                +1:Lw12).*warrayradps*elec2mech; % torque * mech radps
        elseif strcmpi(FWCMethod,'none')==1
            startval.id=0;startval.iq=0;startval.w=warrayradps;
            % solving for finding id, iq coordinates for noFW line of SPMSM
            % this function returns array of operating points[id;iq;wradps]
            operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(motor, ...
                inverter,startval,23,'outputAll',1, ...
                'voltageEquation',voltageEquation);
            idArray(Lw1+1:Lw12)=operatingpoint(1,:);
            iqArray(Lw1+1:Lw12)=operatingpoint(2,:);
            TorqueArray(Lw1+1:Lw12)=1.5*Pp*fluxPM*iqArray(Lw1+1:Lw12);

            PowerArray(Lw1+1:Lw12)=TorqueArray(Lw1 ...
                +1:Lw12).*warrayradps*elec2mech; % torque * mech radps
        else
            % VCLMT method of control - a.k.a. 'optimal current vector' control
            if strcmpi(motorType,'ipmsm') || strcmpi(motorType,'synrm') || strcmpi(motorType,'pmasynrm')
                for iw=(Lw1+1):Lw12
                    w=wArray_Final(iw);
                    wradps=w/(radps2rpm*elec2mech);
                    startval.id=0;startval.iq=0;startval.w=wradps;
                    % solving for IPMSM's FW id, iq with known w
                    % this function returns array of operating points
                    % [id;iq;wradps]
                    operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(motor, ...
                        inverter,startval,10,'outputAll',1, ...
                        'voltageEquation',voltageEquation);
                    id=operatingpoint(1);
                    iq=operatingpoint(2);
                    if goLUT==0
                        torq=1.5*Pp*(fluxPM*iq+(Ld-Lq)*id*iq);
                    else
                        torq=local_interp2(pmsm.PMSMLUT.iqVec,pmsm.PMSMLUT.idVec,pmsm.PMSMLUT.TorqueTable,iq,id);
                    end
                    pow=torq*w/(radps2rpm);% torque * mech radps
                    idArray(iw)=id;
                    iqArray(iw)=iq;
                    TorqueArray(iw)=torq;
                    PowerArray(iw)=pow;
                end
            else
                for iw=(Lw1+1):Lw12
                    w=wArray_Final(iw);
                    wradps=w/(radps2rpm*elec2mech);
                    startval.id=0;startval.iq=0;startval.w=wradps;
                    operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(motor, ...
                        inverter,startval,11,'outputAll',1, ...
                        'voltageEquation',voltageEquation);
                    id=operatingpoint(1);
                    iq=operatingpoint(2);
                    if goLUT==0
                        torq=1.5*Pp*(fluxPM*iq+(Ld-Lq)*id*iq);
                    else
                        torq=local_interp2(pmsm.PMSMLUT.iqVec,pmsm.PMSMLUT.idVec,pmsm.PMSMLUT.TorqueTable,iq,id);
                    end
                    pow=torq*w/(radps2rpm);% torque * mech radps
                    idArray(iw)=id;
                    iqArray(iw)=iq;
                    TorqueArray(iw)=torq;
                    PowerArray(iw)=pow;
                end

                %{
                startval.id=0;startval.iq=0;startval.w=warrayradps;
                %solving for SPMSM FW id, iq when speed is known
                % this function returns array of operating points[id;iq;wradps]
                operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(motor, ...
                    inverter,startval,11,'outputAll',1, ...
                    'voltageEquation',voltageEquation);
                idArray(Lw1+1:Lw12)=operatingpoint(1,:);
                iqArray(Lw1+1:Lw12)=operatingpoint(2,:);
                TorqueArray(Lw1+1:Lw12)=1.5*Pp*fluxPM*iqArray(Lw1+1:Lw12);
                PowerArray(Lw1+1:Lw12)=TorqueArray(Lw1 ...
                    +1:Lw12).*warrayradps*elec2mech; % torque * mech radps
                %}
            end

        end

        % when the motor has third milestone speed, i.e. when the motor has
        % maximum speed, and the second milestone speed (w_c) after which the
        % control trajectory changes.
        if count==3
            % w_d is the maximum speed when it occurs for a motor with the
            % chosen FW control algorithm
            %w_d=floor(CurrentSpeedArray(3,3)*radps2rpm*elec2mech); % mech rpm
            % w_d=PMSMspeeds(3);
            % if w_d>Very_high_speed % if the speeds are very high due to MTPV, set a practical limit.
            %     w_d=Very_high_speed;
            % end
            %warray=(w_c+1:w_d); % mech rpm
            warray=wArrayMS2Max;
            warrayradps=warray/(radps2rpm*elec2mech); % elec radps
            if strcmpi(FWCMethod,'cvcp')==1
                startval.id=0;startval.iq=0;startval.w=warrayradps;
                startval.wcorner=w_corner/(radps2rpm*elec2mech);
                %solving for finding id, iq coordinates for CVCP curve of
                % SPMSM after current touches current limit circle, till max speed
                % this function returns array of operating points[id;iq;wradps]
                operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(motor, ...
                    inverter,startval,27,'outputAll',1, ...
                    'voltageEquation',voltageEquation);
                idArray(Lw12+1:Lw123)=operatingpoint(1,:);
                iqArray(Lw12+1:Lw123)=operatingpoint(2,:);
                TorqueArray(Lw12+1:Lw123)=1.5*Pp*(fluxPM*iqArray(Lw12+1:Lw123) ...
                    + (Ld-Lq)*idArray(Lw12+1:Lw123).*iqArray(Lw12+1:Lw123));
                PowerArray(Lw12+1:Lw123)=TorqueArray(Lw12 ...
                    +1:Lw123).*warrayradps*elec2mech; % torque * mech radps
            elseif strcmpi(FWCMethod,'cccp')==1
                startval.id=0;startval.iq=0;startval.w=warrayradps;
                startval.wcorner=w_corner/(radps2rpm*elec2mech);
                %solving for finding id, iq coordinates for CCCP curve of
                % SPMSM after  vd, vq touch voltage limit circle in vd, vq
                % plane, if id,iq are not going over the current limit
                % this function returns array of operating points[id;iq;wradps]
                operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(motor, ...
                    inverter,startval,31,'outputAll',1, ...
                    'voltageEquation',voltageEquation);
                idArray(Lw12+1:Lw123)=operatingpoint(1,:);
                iqArray(Lw12+1:Lw123)=operatingpoint(2,:);
                TorqueArray(Lw12+1:Lw123)=1.5*Pp*(fluxPM*iqArray(Lw12+1:Lw123) ...
                    + (Ld-Lq)*idArray(Lw12+1:Lw123).*iqArray(Lw12+1:Lw123));
                PowerArray(Lw12+1:Lw123)=TorqueArray(Lw12 ...
                    +1:Lw123).*warrayradps*elec2mech; % torque * mech radps
            elseif strcmpi(FWCMethod,'vclmt')==1
                % MTPV
                if strcmpi(motorType,'ipmsm') || strcmpi(motorType,'synrm') || strcmpi(motorType,'pmasynrm')
                    % since this method uses iterative solution,
                    % a sparse solution could be obtained.
                    %
                    myFcnStack=dbstack;
                    if ~(size(myFcnStack,1)>3 && strcmpi(myFcnStack(4).name,'StartPauseContinue'))
                        DisappearingMessage = fprintf( ...
                            message('mcb:blocks:APINRPleaseWait').getString());
                    end
                    id=id_milestone2;iq=iq_milestone2;
                    for iw=(Lw12+1):Lw123
                        w=wArray_Final(iw);
                        wradps=w/(radps2rpm*elec2mech);
                        startval.id=id;startval.iq=iq;startval.w=wradps;
                        % solving for IPMSM's MTPV id, iq with known w
                        % this function returns array of operating points
                        % [id;iq;wradps] and status
                        [operatingpoint, ~]=mcb.internal.PMSMSpeedCurrentsFcn( ...
                            motor, inverter,startval,9,'outputAll',1, ...
                            'voltageEquation',voltageEquation);
                        % if status_ok==0
                        %     disp(message('mcb:blocks:APINRDidnotConverge', ...
                        %         num2str(operatingpoint(1)), ...
                        %         num2str(operatingpoint(2)), ...
                        %         num2str(operatingpoint(3))).getString());
                        % end
                        id=operatingpoint(1);
                        iq=operatingpoint(2);
                        if goLUT==0
                            torq=1.5*Pp*(fluxPM*iq+(Ld-Lq)*id*iq);
                        else
                            torq=local_interp2(pmsm.PMSMLUT.iqVec,pmsm.PMSMLUT.idVec,pmsm.PMSMLUT.TorqueTable,iq,id);
                        end
                        pow=torq*w/(radps2rpm); % torque * mech radps
                        idArray(iw)=id;
                        iqArray(iw)=iq;
                        TorqueArray(iw)=torq;
                        PowerArray(iw)=pow;
                    end
                    %Retaining the message to wait when called through
                    %simulation environment. Print the message and remove
                    %it when called through scripts only.
                    if ~(size(myFcnStack,1)>3 && strcmpi(myFcnStack(4).name,'StartPauseContinue'))
                        fprintf(repmat('\b',1,DisappearingMessage));

                    end
                else
                     %{
                    startval.id=0;startval.iq=0;startval.w=warrayradps;
                    %solving for SPMSM MTPV id, iq when speed is known
                    % this function returns array of operating points
                    % [id;iq;wradps]
                    operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(motor, ...
                        inverter,startval,12,'outputAll',1, ...
                        'voltageEquation',voltageEquation);
                    idArray(Lw12+1:Lw123)=operatingpoint(1,:);
                    iqArray(Lw12+1:Lw123)=operatingpoint(2,:);
                    TorqueArray(Lw12+1:Lw123)=1.5*Pp*fluxPM*iqArray(Lw12+1:Lw123);
                    PowerArray(Lw12+1:Lw123)=TorqueArray(Lw12 ...
                        +1:Lw123).*warrayradps*elec2mech; % torque * mech radps
                    %}
                    for iw=(Lw12+1):Lw123
                        w=wArray_Final(iw);
                        wradps=w/(radps2rpm*elec2mech);
                        startval.id=0;startval.iq=0;startval.w=wradps;
                        %solving for SPMSM MTPV id, iq when speed is known
                        % this function returns array of operating points
                        % [id;iq;wradps]
                        [operatingpoint, ~]=mcb.internal.PMSMSpeedCurrentsFcn( ...
                            motor, inverter,startval,12,'outputAll',1, ...
                            'voltageEquation',voltageEquation);
                        id=operatingpoint(1);
                        iq=operatingpoint(2);
                        if goLUT==0
                            torq=1.5*Pp*fluxPM*iq;
                        else
                            torq=local_interp2(pmsm.PMSMLUT.iqVec,pmsm.PMSMLUT.idVec,pmsm.PMSMLUT.TorqueTable,iq,id);
                        end
                        pow=torq*w/(radps2rpm); % torque * mech radps
                        idArray(iw)=id;
                        iqArray(iw)=iq;
                        TorqueArray(iw)=torq;
                        PowerArray(iw)=pow;
                    end
                end
            end
        end
        %eliminating the imaginary values computed from approximate
        %equations. These are wrong, but plotting nonetheless.
        TorqueArray=TorqueArray(imag(TorqueArray)==0);
        PowerArray=PowerArray(imag(TorqueArray)==0);
        idArray=idArray(imag(TorqueArray)==0);
        iqArray=iqArray(imag(TorqueArray)==0);
        wArray_Final=wArray_Final(imag(TorqueArray)==0);

        warray=wArray_Final;warrayradps=warray/(radps2rpm*elec2mech);
        IpeakArray=sqrt(idArray.^2+iqArray.^2);
        if goLUT==1
            Ld1=local_interp2(pmsm.PMSMLUT.iqVec,pmsm.PMSMLUT.idVec,pmsm.PMSMLUT.LdTable,iqArray,idArray);
            Lq1=local_interp2(pmsm.PMSMLUT.iqVec,pmsm.PMSMLUT.idVec,pmsm.PMSMLUT.LqTable,iqArray,idArray);
            fluxPM1=local_interp2(pmsm.PMSMLUT.iqVec,pmsm.PMSMLUT.idVec,pmsm.PMSMLUT.FluxPMTable,iqArray,idArray);
            vdArray=idArray*R-Lq1.*warrayradps.*iqArray;
            vqArray=iqArray*R+Ld1.*warrayradps.*idArray+warrayradps.*fluxPM1;
        elseif goLUT==2
            fluxd1=local_interp2(pmsm.PMSMLUT.iqVec,pmsm.PMSMLUT.idVec,pmsm.PMSMLUT.FluxDTable,iqArray,idArray);
            fluxq1=local_interp2(pmsm.PMSMLUT.iqVec,pmsm.PMSMLUT.idVec,pmsm.PMSMLUT.FluxQTable,iqArray,idArray);
            vdArray=idArray*R-fluxq1.*warrayradps;
            vqArray=iqArray*R+fluxd1.*warrayradps;
        else
            vdArray=idArray*R-Lq*warrayradps.*iqArray;
            vqArray=iqArray*R+Ld*warrayradps.*idArray+warrayradps*fluxPM;
        end

        % the following script is to eliminate the data points in the plot
        % with invalid currents and voltages
        %IpeakArray=IpeakArray(IpeakArray<=Irated*1.01);
        Tfriction=Bv*warrayradps;
        indexhigh=find(TorqueArray<=Tfriction, 1 );
        if isempty(indexhigh)
            indexhigh=length(iqArray);
        end
        vm=sqrt(vdArray.^2+vqArray.^2);
        indexhigh1=find(vm>=1.001*inverter.V_dc/sqrt(3),1);
        if isempty(indexhigh1)
            indexhigh1=length(iqArray);
        end


        indexhigh2=find(IpeakArray>=motor.I_rated*1.01,1);
        if isempty(indexhigh2)
            indexhigh2=length(iqArray);
        end
        indexhigh3=find(imag(iqArray)~=0,1);
        if isempty(indexhigh3)
            indexhigh3=length(iqArray);
        end
        iwarray=1:min([indexhigh2 indexhigh3 indexhigh indexhigh1]);
        warray=warray(iwarray);
        warrayradps=warray/(radps2rpm*elec2mech);

        %
        idArray=idArray(iwarray);
        iqArray=iqArray(iwarray);
        vdArray=vdArray(iwarray);
        vqArray=vqArray(iwarray);
        if changeDQAxes==1
            tmp=idArray;
            idArray=iqArray;
            iqArray=-tmp;
        end

        if driveCharacteristics>=1
            figure(DrvCharFigH);
            if opacity==1
                sg=sgtitle(strcat("FW Control method: ",  FWCMethodtext ,", VoltageEquation: ",  voltageEquation ));
                sg.FontSize=10;
            else
                %Since the multiple figures could come from different
                %configurations, this title is not displayed
                sgtitle(" ");
                set(gcf, 'Name', 'PMSM : Drive Characteristics', 'NumberTitle', 'Off','Tag','DriveCharacteristics')
            end

            subplot(2,2,1);hold on;
            if pmsmloop==1
                thisplot=plot(warray,TorqueArray(iwarray),'b',warray,Bv*warrayradps,'r');
                thiscolor=thisplot(1).Color;
                thisplot(1).Color=[thiscolor(1:3) opacity];
                thiscolor=thisplot(2).Color;
                thisplot(2).Color=[thiscolor(1:3) opacity];
                lh=thisplot(1);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='Nm';
                lh.DataTipTemplate.DataTipRows(end+1)='Torque available';
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("irated: ",num2str(Irated));
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
                lh=thisplot(2);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='Nm';
                lh.DataTipTemplate.DataTipRows(end+1)='Torque frictional';

                title('Torque-Speed Characteristics');
                xlabel('Speed (rpm)');
                ylabel('Torque (Nm)');
                legend({'T_{available}','T_{InternalFriction}'},'location','best');
                grid on;
                if flag_imax==0
                    scale2 = axis;
                    axis([scale2(1:2) 0 1.1*scale2(4)]);
                end
            elseif pmsmloop==2
                thisplot=plot(warray,TorqueArray(iwarray),'b-.',warray,Bv*warrayradps,'r-.');
                thiscolor=thisplot(1).Color;
                thisplot(1).Color=[thiscolor(1:3) opacity];
                thiscolor=thisplot(2).Color;
                thisplot(2).Color=[thiscolor(1:3) opacity];
                lh=thisplot(1);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='Nm';
                lh.DataTipTemplate.DataTipRows(end+1)='Short-time torque available';
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("imax: ",num2str(imax));
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
                lh=thisplot(2);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='Nm';
                lh.DataTipTemplate.DataTipRows(end+1)='Torque frictional';

                legend({['T_{available}@ I_{rated}=' num2str(pmsm.I_rated) 'A'],'T_{InternalFriction}', ...
                    ['T_{available}@ I_{max}=' num2str(imax) 'A']},'location','best');
                scale2 = axis;
                axis([scale2(1:2) 0 1.1*scale2(4)]);
            end

            subplot(2,2,2);hold on;
            if pmsmloop==1
                thisplot=plot(warray,PowerArray(iwarray),'b', ...
                    warray,Bv*warrayradps.*warrayradps*elec2mech,'r');
                thiscolor=thisplot(1).Color;
                thisplot(1).Color=[thiscolor(1:3) opacity];
                thiscolor=thisplot(2).Color;
                thisplot(2).Color=[thiscolor(1:3) opacity];
                lh=thisplot(1);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='W';
                lh.DataTipTemplate.DataTipRows(end+1)='Power available';
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("irated: ",num2str(Irated));
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
                lh=thisplot(2);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='W';
                lh.DataTipTemplate.DataTipRows(end+1)='Power frictional';
                title('Power-Speed Characteristics');
                xlabel('Speed (rpm)');
                ylabel('Power (W)');
                legend({'P_{available}','P_{FrictionalLoss}'},'location','best');
                grid on;
                if flag_imax==0
                    scale2 = axis;
                    axis([scale2(1:2) 0 1.1*scale2(4)]);
                end
            elseif pmsmloop==2
                thisplot=plot(warray,PowerArray(iwarray),'b-.', ...
                    warray,Bv*warrayradps.*warrayradps*elec2mech,'r-.');
                thiscolor=thisplot(1).Color;
                thisplot(1).Color=[thiscolor(1:3) opacity];
                thiscolor=thisplot(2).Color;
                thisplot(2).Color=[thiscolor(1:3) opacity];
                lh=thisplot(1);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='W';
                lh.DataTipTemplate.DataTipRows(end+1)='Short-time power available';
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("imax: ",num2str(imax));
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
                lh=thisplot(2);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='W';
                lh.DataTipTemplate.DataTipRows(end+1)='Power frictional';
                legend({['P_{available}@ I_{rated}=' num2str(pmsm.I_rated) 'A']','P_{FrictionalLoss}', ...
                    ['P_{available}@ I_{max}=' num2str(imax) 'A']},'location','best');
                scale2 = axis;
                axis([scale2(1:2) 0 1.1*scale2(4)]);
            end

            subplot(2,2,3);hold on;
            if pmsmloop==1
                thisplot=plot(warray,idArray(iwarray),'k',warray,iqArray(iwarray),'b');
                thiscolor=thisplot(1).Color;
                thisplot(1).Color=[thiscolor(1:3) opacity];
                thiscolor=thisplot(2).Color;
                thisplot(2).Color=[thiscolor(1:3) opacity];
                lh=thisplot(1);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='id';
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("irated: ",num2str(Irated));
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
                lh=thisplot(2);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='iq';
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("irated: ",num2str(Irated));
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
                title('I_{dq} Vs Speed');
                xlabel('Speed (rpm)');
                ylabel('I_{dq} (A)');
                legend({'I_d','I_q'},'location','best');
                grid on;
                if flag_imax==0
                    scale = axis;
                    axis([scale(1:2) -Irated*1.1 Irated*1.1]);
                end
            elseif pmsmloop==2
                thisplot=plot(warray,idArray(iwarray),'k-.',warray,iqArray(iwarray),'b-.');
                thiscolor=thisplot(1).Color;
                thisplot(1).Color=[thiscolor(1:3) opacity];
                thiscolor=thisplot(2).Color;
                thisplot(2).Color=[thiscolor(1:3) opacity];
                lh=thisplot(1);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='id';
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("imax: ",num2str(imax));
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
                lh=thisplot(2);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='iq';
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("imax: ",num2str(imax));
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
                legend({['I_d@ I_{rated}=' num2str(pmsm.I_rated) 'A'], ...
                    ['I_q@ I_{rated}=' num2str(pmsm.I_rated) 'A'], ...
                    ['I_d@ I_{max}=' num2str(imax) 'A'], ...
                    ['I_q@ I_{max}=' num2str(imax) 'A']},'location','best');
                scale = axis;
                axis([scale(1:2) -imax*1.1 imax*1.1]);
            end


            subplot(2,2,4);hold on;
            if pmsmloop==1
                thisplot=plot(warray,IpeakArray(iwarray),'b', ...
                    warray,Irated*ones(1,length(iwarray)),'r');
                thiscolor=thisplot(1).Color;
                thisplot(1).Color=[thiscolor(1:3) opacity];
                thiscolor=thisplot(2).Color;
                thisplot(2).Color=[thiscolor(1:3) opacity];
                lh=thisplot(1);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='ipeak';
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
                lh=thisplot(2);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='irated';
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
                title('I_{peak} Vs Speed');
                xlabel('Speed (rpm)');
                ylabel('I_{peak} (A)');
                legend({'I_{peak}','I_{rated}'},'location','best');
                grid on;
                if flag_imax==0
                    scale = axis;
                    axis([scale(1:2) 0 Irated*1.1]);
                end
            elseif pmsmloop==2
                thisplot=plot(warray,IpeakArray(iwarray),'b-.', ...
                    warray,imax*ones(1,length(warray)),'r-.');
                thiscolor=thisplot(1).Color;
                thisplot(1).Color=[thiscolor(1:3) opacity];
                thiscolor=thisplot(2).Color;
                thisplot(2).Color=[thiscolor(1:3) opacity];
                lh=thisplot(1);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='ipeak with imax';
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
                lh=thisplot(2);
                lh.DataTipTemplate.DataTipRows(1).Label='rpm';
                lh.DataTipTemplate.DataTipRows(2).Label='imax';
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
                lh.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
                legend({['I_{peak}@ I_{rated}=' num2str(pmsm.I_rated) 'A'], ...
                    ['I_{rated}(' num2str(pmsm.I_rated) 'A)'], ...
                    ['I_{peak}@ I_{max}=' num2str(imax) 'A'], ...
                    ['I_{max}(' num2str(imax) 'A)']},'location','best');
                scale = axis;
                axis([scale(1:2) 0 imax*1.1]);
            end
        end

        if pmsmloop==1
            % store the voltage and current arrays for the original motor
            vdArray0=vdArray;
            vqArray0=vqArray;
            idArray0=idArray;
            iqArray0=iqArray;
            wArray_Final0=wArray_Final;
        end
    end
    % restoring the currents and voltages for the original motor.
    % these will be used when plotting along with the constraint curves
    vdArray=vdArray0;
    vqArray=vqArray0;
    idArray=idArray0;
    iqArray=iqArray0;
    wArray_Final=wArray_Final0;
end
%%

% The critical current point from the approximate equations, used primarily
% for finding the scale for plotting
Icritical=-fluxPM/Ld;
% define the second quadrant only
xmax=max(Irated,abs(Icritical));
w=motorSystem.wrpm/(radps2rpm*elec2mech);
if length(motorSystem.Tload)>length(w)
    T=motorSystem.Tload(1:length(w))+Bv*w;
    T1=motorSystem.Tload((length(w)+1):end)+Bv*w(end);
    T=[T(:);T1(:)];
else
    T=motorSystem.Tload+Bv*w(1:length(motorSystem.Tload));
end

% The following code plots the constraint curves for the PMSM
CreateNewFigure=0;
if constraintCurves==1 && opacity==1
    CreateNewFigure=1;
else
    availablefigHandles = findobj('Type', 'figure','Tag','ConstraintCurve');
    %Picking the last created constraint curve for plotting
    if isempty(availablefigHandles)
        CreateNewFigure=1;
    else
        CnstCrvFigH=availablefigHandles(1);
    end
end
if constraintCurves==1

    if  goLUT==1 || goLUT==2
        xmax=min(xmax,max(max(abs(pmsm.PMSMLUT.idVec)),max(abs(pmsm.PMSMLUT.iqVec))));
        % end
        % if isfield(pmsm,'PMSMLUT')
        startval.ParamTableData=pmsm.PMSMLUT;
    end
    if CreateNewFigure==1
        CnstCrvFigH=figure('Name',strcat("Constraint Curves of ", ...
            pmsm.model), 'numberTitle','off','Tag','ConstraintCurve');
    end
    figure(CnstCrvFigH);
    if opacity==1
        sg=sgtitle(strcat("FW Control method: ",  FWCMethodtext ,", VoltageEquation: ",  voltageEquation ));
        sg.FontSize=10;
    else
        %Since the multiple figures could come from different
        %configurations, this title is not displayed
        sgtitle(" ");
        set(gcf, 'Name', 'PMSM : Constraint Curves', 'NumberTitle', 'Off','Tag','ConstraintCurve');
    end
    axis equal;
    hold on;
    curvesplotted=0;

    if drawContours==1
        linesplotted=drawPMSMContours(pmsm,inverter,changeDQAxes);
        for lgd=1:linesplotted
            curvesplotted=curvesplotted+1;
            legendtxt(curvesplotted)="";%#ok
        end
    end

    wrpm_unique=unique(wrpm);
    w_unique=unique(w);

    if length(T_load)*length(w_unique)>1
        sgtitle(" ");
        set(gcf, 'Name', 'PMSM : Constraint Curves', 'NumberTitle', 'Off','Tag','ConstraintCurve');
    end

    % current limiting circle
    %Creating the x-axis points
    % large number of points result in smooth curve
    x=linspace(-Irated*2,Irated*2,10000);
    %passing the id values to the function
    startval.id=x;startval.iq=0;startval.w=0;
    %solving for finding id, iq coordinates for Current limit circle
    % this function returns array of operating points[id;iq;wradps]
    operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(pmsm, ...
        inverter,startval,13,'outputAll',1, ...
        'voltageEquation',voltageEquation);
    id=operatingpoint(1,:);
    iq=operatingpoint(2,:);
    if goLUT==0
        id(end+1)=id(1);
        iq(end+1)=iq(1);
    else
        id=id(iq>=0);%returns only second quadrant results
        iq=iq(iq>=0);
        iq=iq(id<=0);
        id=id(id<=0);
    end
    if changeDQAxes==1
        tmp=id;
        id=iq;
        iq=-tmp;
    end
    lh=plot(id,iq,'k');
    thiscolor=lh.Color;
    lh.Color=[thiscolor(1:3) opacity];
    lh.DataTipTemplate.DataTipRows(1).Label='id';
    lh.DataTipTemplate.DataTipRows(2).Label='iq';
    lh.DataTipTemplate.DataTipRows(end+1)='Current limit curve';
    lh.DataTipTemplate.DataTipRows(end+1)=strcat("I_{rated} ",string(Irated),"A");
    curvesplotted=curvesplotted+1;
    if opacity<1
        legendtxt(curvesplotted)="Current Limit";
    else
        legendtxt(curvesplotted)=strcat("Current Limit(I_{rated}=",string(Irated),"A)");
    end
    CnstCurve.current.id=id;
    CnstCurve.current.iq=iq;
    %drawnow;

    L_w=length(unique(wrpm));
    L_T=length(T_load);
    if L_w < L_T
        Opacityratio=T_load/(max(T_load));
        motorSystem_wrpm(1:L_w)=motorSystem.wrpm;
        motorSystem_wrpm((L_w+1):L_T)=motorSystem_wrpm(L_w);
        motorSystem_Tload=motorSystem.Tload;
    elseif L_w > L_T
        Opacityratio=wrpm/(max(wrpm));
        motorSystem_Tload(1:L_T)=motorSystem.Tload;
        motorSystem_Tload((L_T+1):L_w)=motorSystem_Tload(L_T);
        motorSystem_wrpm=motorSystem.wrpm;
    else
        Opacityratio=wrpm/(max(wrpm));
        motorSystem_wrpm=motorSystem.wrpm;
        motorSystem_Tload=motorSystem.Tload;
    end
    if isnan(Opacityratio)
        Opacityratio=1;
    end
    Opacityratio(Opacityratio<0.1)=0.1;

    for iw=1:length(T_load)
        % constant torque curve
        % creating the x-axis values
        % large number of points result in smooth curve
        if goLUT==0
            x=linspace(-xmax,0,10000);
        else
            x=linspace(-xmax,0,1000);
        end
        %passing the id values and torque to the function;
        startval.id=x;startval.iq=0;startval.t=T(iw);
        if strcmpi(motorSystem.motorType,'spmsm')
            %solving for finding id, iq coordinates for
            % constant torque hyperbola for SPMSM
            % this function returns array of operating points[id;iq;wradps]
            operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(pmsm, ...
                inverter,startval,14,'outputAll',1, ...
                'voltageEquation',voltageEquation);
        else
            %solving for finding id, iq coordinates for
            % constant torque hyperbola for IPMSM
            % this function returns array of operating points[id;iq;wradps]
            operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(pmsm, ...
                inverter,startval,15,'outputAll',1, ...
                'voltageEquation',voltageEquation);
        end
        if ~isempty(operatingpoint)
            id=operatingpoint(1,:);
            iq=operatingpoint(2,:);
            if changeDQAxes==1
                tmp=id;
                id=iq;
                iq=-tmp;
            end
            lh=plot(id,iq);
            
            lh.Color=[1 0 0 opacity*Opacityratio(iw)];
            lh.DataTipTemplate.DataTipRows(1).Label='id';
            lh.DataTipTemplate.DataTipRows(2).Label='iq';
            lh.DataTipTemplate.DataTipRows(end+1)='Constant torque curve';
            lh.DataTipTemplate.DataTipRows(end+1)=strcat(string(round(motorSystem_wrpm(iw)))," rpm");
            lh.DataTipTemplate.DataTipRows(end+1)=strcat(string(motorSystem_Tload(iw))," N.m. load");
            curvesplotted=curvesplotted+1;
            if opacity<1
                legendtxt(curvesplotted)="Constant Torque";
            else
                legendtxt(curvesplotted)=strcat("Constant Torque(T_{load}=" ...
                    ,string(motorSystem.Tload(iw)),"N.m, T_{friction}=" ...
                    ,string(T(iw)-motorSystem.Tload(iw)),"N.m)");
            end
            CnstCurve.torque.id=id;
            CnstCurve.torque.iq=iq;
        else
            CnstCurve.torque.id=[];
            CnstCurve.torque.iq=[];
        end
    end
    %drawnow;


    % MTPA curve
    if goLUT==0
        x=linspace(-xmax,0,10000);
    else
        x=linspace(-xmax,0,100);
    end
    % large number of points result in smooth curve
    startval.id=x;startval.iq=0;startval.w=0;
    if strcmpi(motorSystem.motorType,'spmsm')
        %solving for finding id, iq coordinates for MTPA line of SPMSM
        % this function returns array of operating points[id;iq;wradps]
        operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(pmsm, ...
            inverter,startval,17,'outputAll',1, ...
            'voltageEquation',voltageEquation);
    else
        %solving for finding id, iq coordinates for MTPA curve of IPMSM
        % this function returns array of operating points[id;iq;wradps]
        operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(pmsm, ...
            inverter,startval,18,'outputAll',1, ...
            'voltageEquation',voltageEquation);
    end
    id=operatingpoint(1,:);
    iq=operatingpoint(2,:);
    if changeDQAxes==1
        tmp=id;
        id=iq;
        iq=-tmp;
    end
    lh=plot(id,iq,'b--');
    thiscolor=lh.Color;
    lh.Color=[thiscolor(1:3) opacity];
    lh.DataTipTemplate.DataTipRows(1).Label='id';
    lh.DataTipTemplate.DataTipRows(2).Label='iq';
    lh.DataTipTemplate.DataTipRows(end+1)='MTPA curve';
    curvesplotted=curvesplotted+1;
    legendtxt(curvesplotted)="MTPA";
    CnstCurve.mtpa.id=id;
    CnstCurve.mtpa.iq=iq;
    %drawnow;

    Opacityratio=wrpm_unique/(max(wrpm_unique));
    if isnan(Opacityratio)
        Opacityratio=1;
    end
    Opacityratio(Opacityratio<0.1)=0.1;
    for iw=1:length(w_unique)
        % voltage limit curve
        if goLUT==0
            x=linspace(-xmax*4,xmax*4,20000);
        else
            x=linspace(-xmax*4,xmax*4,4000);
        end
        % large number of points result in smooth curve
        startval.id=x;startval.iq=0;startval.w=w_unique(iw);
        %solving for finding id, iq coordinates for voltage limit curve
        % for IPMSM, SPMSM
        % this function returns array of operating points[id;iq;wradps]
        operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(pmsm, ...
            inverter,startval,16,'outputAll',1, ...
            'voltageEquation',voltageEquation);
        if ~isempty(operatingpoint)
            id=operatingpoint(1,:);%id(end+1)=id(1);
            iq=operatingpoint(2,:);%iq(end+1)=iq(1);
            if changeDQAxes==1
                tmp=id;
                id=iq;
                iq=-tmp;
            end
            lh=plot(id,iq,'b:','LineWidth',2);
            thiscolor=lh.Color;

            lh.Color=[thiscolor(1:3) opacity*Opacityratio(iw)];
            lh.DataTipTemplate.DataTipRows(1).Label='id';
            lh.DataTipTemplate.DataTipRows(2).Label='iq';
            lh.DataTipTemplate.DataTipRows(end+1)='Voltage limit curve';
            lh.DataTipTemplate.DataTipRows(end+1)=strcat(string(round(wrpm_unique(iw)))," rpm");
            lh.DataTipTemplate.DataTipRows(end+1)=strcat(string(motorSystem.Vdc)," V");
            lh.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
            lh.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
            critical.y=(max(iq)+min(iq))/2;
            critical.x=(max(id)+min(id))/2;
            curvesplotted=curvesplotted+1;
            if opacity<1
                legendtxt(curvesplotted)='Voltage Limit';
            else
                legendtxt(curvesplotted)=strcat("Voltage Limit(V_{DC}=", ...
                    string(motorSystem.Vdc) ,"V, w=", ...
                    string(round(wrpm_unique(iw))),"rpm)");
            end
            CnstCurve.voltage.id=id;
            CnstCurve.voltage.iq=iq;
        else
            CnstCurve.voltage.id=[];
            CnstCurve.voltage.iq=[];
        end
        %when the FW control mode is CVCP, also draw the cvcp fw trajectory
        if strcmpi(FWCMethod,'cvcp')
            startval.id=x;startval.iq=0;startval.w=w_unique(iw);
            %solving for finding id, iq coordinates for CVCP trajectory for
            %plotting
            % this function returns array of operating points[id;iq;wradps]
            operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(pmsm, ...
                inverter,startval,32,'outputAll',1, ...
                'voltageEquation',voltageEquation);
            id=operatingpoint(1,:);
            iq=operatingpoint(2,:);
            if changeDQAxes==1
                tmp=id;
                id=iq;
                iq=-tmp;
            end
            lh=plot(id,iq,'g');
            lh.DataTipTemplate.DataTipRows(1).Label='id';
            lh.DataTipTemplate.DataTipRows(2).Label='iq';
            lh.DataTipTemplate.DataTipRows(end+1)=strcat(string(round(wrpm_unique(iw)))," rpm");
            lh.DataTipTemplate.DataTipRows(end+1)='CVCP line';
            curvesplotted=curvesplotted+1;
            legendtxt(curvesplotted)="CVCP Line";
            CnstCurve.cvcp.id=id;
            CnstCurve.cvcp.iq=iq;
        end
    end
    %drawnow;


    % MTPV curve
    % large number of points result in smooth curve
    if goLUT==0
        x1=linspace(-xmax,0,1001);
    else
        x1=linspace(-xmax,0,1000);
    end
    if strcmpi(FWCMethod,'vclmt')
        startval.id=x1;startval.iq=0;startval.w=w(1);
        if strcmpi(motorSystem.motorType,'spmsm')
            %solving for finding id, iq coordinates for MTPV line of SPMSM
            % this function returns array of operating points[id;iq;wradps]
            operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(pmsm, ...
                inverter,startval,19,'outputAll',1, ...
                'voltageEquation',voltageEquation);
        else
            %solving for finding id, iq coordinates for MTPV curve of IPMSM
            % this function returns array of operating points[id;iq;wradps]
            operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(pmsm, ...
                inverter,startval,20,'outputAll',1, ...
                'voltageEquation',voltageEquation);
        end
        if ~isempty(operatingpoint)
            id=operatingpoint(1,:);
            iq=operatingpoint(2,:);
            if changeDQAxes==1
                tmp=id;
                id=iq;
                iq=-tmp;
            end
            lh=plot(id,iq,'m-.');
            thiscolor=lh.Color;
            lh.Color=[thiscolor(1:3) opacity];
            curvesplotted=curvesplotted+1;
            legendtxt(curvesplotted)="MTPV";
            lh.DataTipTemplate.DataTipRows(1).Label='id';
            lh.DataTipTemplate.DataTipRows(2).Label='iq';
            lh.DataTipTemplate.DataTipRows(end+1)='MTPV curve';
            %lh.DataTipTemplate.DataTipRows(end+1)=strcat(string(motorSystem.wrpm)," rpm");
            lh.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
            lh.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
            legendtxt(curvesplotted)="MTPV";
            CnstCurve.mtpv.id=id;
            CnstCurve.mtpv.iq=iq;
        else
            CnstCurve.mtpv.id=[];
            CnstCurve.mtpv.iq=[];
        end

    end
    %drawnow;

    if goLUT==0
        % center of voltage limit curve
        thisarray=[1-opacity 1-opacity 1-opacity];
        if changeDQAxes==1
            tmp=critical.x;
            critical.x=critical.y;
            critical.y=-tmp;
        end
        lh=plot(critical.x,critical.y,'*','color',thisarray);
        curvesplotted=curvesplotted+1;
        legendtxt(curvesplotted)="";
        lh.DataTipTemplate.DataTipRows(1).Label='id';
        lh.DataTipTemplate.DataTipRows(2).Label='iq';
        lh.DataTipTemplate.DataTipRows(end+1)=strcat(string(round(motorSystem.wrpm(iw)))," rpm");
    end

    if flag_idq==1
        thisplot=plot(idqCustom(:,1),idqCustom(:,2),'m-*');
        curvesplotted=curvesplotted+1;
        legendtxt(curvesplotted)="user-provided id,iq";
        thisplot.DataTipTemplate.DataTipRows(1).Label='id';
        thisplot.DataTipTemplate.DataTipRows(2).Label='iq';
        thisplot.DataTipTemplate.DataTipRows(end+1)='User-provided id,iq';
        if idqsize(2)==3
            thisplot.DataTipTemplate.DataTipRows(end+1)=dataTipTextRow('rpm', string(idqCustom(:,3)));
        end
    end

    %

    if changeDQAxes==1
        xlim([0 xmax]);
        ylim([0 xmax]);
        set(gca,'YAxisLocation','left');
    else
        xlim([-xmax 0]);
        ylim([0 xmax]);
        set(gca,'YAxisLocation','right');
    end


    if driveCharacteristics>=1
        wrpm_currentplot=1:length(wArray_Final);
        wArray_Final1=wArray_Final;
        wArray_Final1(1:Lw1)=base_rpm;
        if opacity==1
            thistxt="  \leftarrow Possible Control Trajectory";
        else
            thistxt=" ";
        end
        % if changeDQAxes==1
        %     tmp=idArray;
        %     idArray=iqArray;
        %     iqArray=-tmp;
        % end
        thisplot=plot(idArray,iqArray,'c-','LineWidth',2);
        % thisplot.MarkerIndices = floor(linspace(1,length(iqArray),20));
        thisplot.DataTipTemplate.DataTipRows(1).Label='id';
        thisplot.DataTipTemplate.DataTipRows(2).Label='iq';
        thisplot.DataTipTemplate.DataTipRows(end+1)='id,iq trajectory';
        thisplot.DataTipTemplate.DataTipRows(end+1)=dataTipTextRow('rpm', round(wArray_Final1(wrpm_currentplot)));
        thisplot.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
        thisplot.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
        text(idArray(round(median(iwarray))),iqArray(round(median(iwarray))),thistxt,"Color",'c');
        curvesplotted=curvesplotted+1;
        legendtxt(curvesplotted)="";
    end
    legend(legendtxt ,"location","best");

    grid on
    xlabel('d-axis current (A)');
    ylabel('q-axis current (A)');
    drawnow;


    % the following code will plot the voltage limit curve in the
    % vd-vq plane and the corresponding control trajectory for the
    % chosen FW control algorithm
    CreateNewFigure=0;
    if driveCharacteristics>1 && opacity==1
        CreateNewFigure=1;
    else
        availablefigHandles = findobj('Type', 'figure','Tag','vdvqCurve');
        %Picking the last created constraint curve for plotting
        if isempty(availablefigHandles)
            CreateNewFigure=1;
        else
            vdvqCrvFigH=availablefigHandles(1);
        end
    end

    if driveCharacteristics>1
        if CreateNewFigure==1
            vdvqCrvFigH=figure('Name',strcat("PMSM: voltage trajectory in the" + ...
                " vd vq space for ", ...
                pmsm.model), 'numberTitle','off','Tag','vdvqCurve');
        end
        figure(vdvqCrvFigH);
        if opacity==1
            sg=sgtitle(strcat("FW Control method: ",  FWCMethodtext ,", VoltageEquation: ",  voltageEquation ));
            sg.FontSize=10;
        else
            %Since the multiple figures could come from different
            %configurations, this title is not displayed
            sgtitle(" ");
            set(gcf, 'Name', 'PMSM : voltage trajectory in the vd vq space', 'NumberTitle', 'Off','Tag','vdvqCurve');
        end
        axis equal;
        hold on;
        % voltage limiting circle
        %Creating the x-axis points
        x=linspace(-(motorSystem.Vdc/sqrt(3)),(motorSystem.Vdc/sqrt(3)),1000);
        %passing the id values to the function
        startval.vd=x;startval.vq=0;startval.w=0;
        % finding the vd, vw coordinates for circle diagram in vd, vq space
        % this function returns array of operating points[id;iq;wradps]
        operatingpoint=mcb.internal.PMSMSpeedCurrentsFcn(pmsm, ...
            inverter,startval,21,'outputAll',1);
        % using the actual voltage limit circle in the above computation
        vd=operatingpoint(1,:);
        vq=operatingpoint(2,:);
        %pg=polyshape(vd,vq);

        if changeDQAxes==1
            tmp=vd;
            vd=vq;
            vq=-tmp;
        end
        plot(vd,vq,'Color','black');

        wrpm_voltageplot=1:length(wArray_Final);
        thisplot=plot(vdArray,vqArray,'m*-','LineWidth',2);
        thisplot.MarkerIndices = floor(linspace(1,length(vqArray),20));
        thiscolor=thisplot.Color;
        thisplot.Color=[thiscolor(1:3) opacity];
        thisplot.DataTipTemplate.DataTipRows(1).Label='vd';
        thisplot.DataTipTemplate.DataTipRows(2).Label='vq';
        thisplot.DataTipTemplate.DataTipRows(end+1)='Vd,Vq trajectory';
        thisplot.DataTipTemplate.DataTipRows(end+1)=dataTipTextRow('rpm', round(wArray_Final(wrpm_voltageplot)));
        thisplot.DataTipTemplate.DataTipRows(end+1)=strcat("FW Control method: ",FWCMethodtext);
        thisplot.DataTipTemplate.DataTipRows(end+1)=strcat("voltageEquation: ",voltageEquation);
        txt1(1)="Voltage limit circle";
        txt1(2)="Possible Control Trajectory";
        legend(txt1 ,"location","best");

        grid on
        xlabel('d-axis voltage (V)');
        ylabel('q-axis voltage (V)');
    end

elseif drawContours==1
    thisMotorType=pmsm.model;
    if isfield(pmsm,'oldmotorType')
        thisMotorType=pmsm.oldmotorType;
    end
    figure('Name',strcat("Torque and voltage limit contours for ", ...
        thisMotorType), 'numberTitle','off')
    axis equal;
    hold on;

    drawPMSMContours(pmsm,inverter,changeDQAxes);

end
if nargout==1
    if constraintCurves==1
        outStructure=CnstCurve;
    end
    outStructure.idArray=idArray;
    outStructure.iqArray=iqArray;
    outStructure.vdArray=vdArray;
    outStructure.vqArray=vqArray;
    outStructure.wArray=warray;
    outStructure.TArray=TorqueArray(iwarray);
    outStructure.PArray=PowerArray(iwarray);
    outStructure.pmsm=pmsm;
    outStructure.inverter=inverter;
    outStructure.FWCMethod=FWCMethod;
    outStructure.voltageEquation=voltageEquation;
    outStructure.speed_milestone=PMSMspeeds;
    outStructure.id_milestone=ID_milestone;
    outStructure.iq_milestone=IQ_milestone;
    varargout{1}=outStructure;
end

end


%% a local function to draw contours
function linesplotted=drawPMSMContours(pmsm,inverter,changeDQAxes,varargin)
id=linspace(-pmsm.I_rated,0,100);
iq=linspace(0,pmsm.I_rated,100);



goLUT=0;
Tsteps=20;
wsteps=20;
verbose=0;
Very_high_speed=200000;
MeshGridSize=100;
if isfield(pmsm,'PMSMLUT') && isfield(pmsm.PMSMLUT,'method')
    if strcmpi(pmsm.PMSMLUT.method,'Ldq')
        if ~(isfield(pmsm.PMSMLUT,'idVec') && isfield(pmsm.PMSMLUT,'iqVec') ...
                && isfield(pmsm.PMSMLUT,'LdTable') ...
                && isfield(pmsm.PMSMLUT,'LqTable') ...
                && isfield(pmsm.PMSMLUT,'FluxPMTable'))
            goLUT=0;
            if verbose>0
                disp('#### One or more of idvec, iqvec, LqTable, LqTable, and FluxPMTable are missing. Proceeding with lumped parameters.');
            end
        else
            goLUT=1;
        end
    elseif strcmpi(pmsm.PMSMLUT.method,'FluxDQ')
        if ~(isfield(pmsm.PMSMLUT,'idVec') && isfield(pmsm.PMSMLUT,'iqVec') ...
                && isfield(pmsm.PMSMLUT,'FluxDTable') ...
                && isfield(pmsm.PMSMLUT,'FluxQTable'))
            goLUT=0;
            if verbose>0
                disp('#### FluxD or FluxQ LUT is missing. Proceeding with lumped parameters.');
            end
        else
            goLUT=2;
        end
    end
end
if goLUT>0
    [id_m1,iq_m1]=meshgrid(linspace(max(-pmsm.I_rated,min(pmsm.PMSMLUT.idVec)),min(0,max(pmsm.PMSMLUT.idVec)),MeshGridSize),linspace(max(0,min(pmsm.PMSMLUT.iqVec)),min(pmsm.I_rated,max(pmsm.PMSMLUT.iqVec)),MeshGridSize));

    id_mesh=id_m1;
    iq_mesh=iq_m1;
    id=pmsm.PMSMLUT.idVec;
    iq=pmsm.PMSMLUT.iqVec;
    if goLUT==1
        pmsm.PMSMLUT.LdTable=interp2(iq,id,pmsm.PMSMLUT.LdTable,iq_mesh,id_mesh)';
        pmsm.PMSMLUT.LqTable=interp2(iq,id,pmsm.PMSMLUT.LqTable,iq_mesh,id_mesh)';
        pmsm.PMSMLUT.FluxPMTable=interp2(iq,id,pmsm.PMSMLUT.FluxPMTable,iq_mesh,id_mesh)';
    elseif goLUT==2
        iq0index=find(iq==0);
        if ~isempty(iq0index)
            for i1=1:length(id)
                pmsm.PMSMLUT.FluxQTable(i1,:)=pmsm.PMSMLUT.FluxQTable(i1,:)-pmsm.PMSMLUT.FluxQTable(i1,iq0index);
            end
        end
        pmsm.PMSMLUT.FluxDTable=interp2(iq,id,pmsm.PMSMLUT.FluxDTable,iq_mesh,id_mesh)';
        pmsm.PMSMLUT.FluxQTable=interp2(iq,id,pmsm.PMSMLUT.FluxQTable,iq_mesh,id_mesh)';
    end
    if isfield(pmsm.PMSMLUT,'TorqueTable')
        pmsm.PMSMLUT.TorqueTable=interp2(iq,id,pmsm.PMSMLUT.TorqueTable,iq_mesh,id_mesh)';
    end

    id=id_mesh(1,:);
    iq=iq_mesh(:,1);
end
if goLUT==1
    % id=pmsm.PMSMLUT.idVec;
    % iq=pmsm.PMSMLUT.iqVec;

    if isfield(pmsm.PMSMLUT,'TorqueTable')
        T_computed=pmsm.PMSMLUT.TorqueTable;
    else
        % id=id(id<=pmsm.I_rated);
        % id=id(id>=-pmsm.I_rated);
        % iq=iq(iq<=pmsm.I_rated);
        % iq=iq(iq>=-pmsm.I_rated);
        % id=id(id<=0);
        % iq=iq(iq>=0);
        for ix=1:length(id)
            for iy=1:length(iq)
                T_computed(ix,iy)=1.5*pmsm.p*iq(iy)*(pmsm.PMSMLUT.FluxPMTable(ix,iy)+(pmsm.PMSMLUT.LdTable(ix,iy)-pmsm.PMSMLUT.LqTable(ix,iy))*id(ix));%#ok
            end
        end
    end
elseif goLUT==2
    % id=pmsm.PMSMLUT.idVec;
    % iq=pmsm.PMSMLUT.iqVec;

    if isfield(pmsm.PMSMLUT,'TorqueTable')
        T_computed=pmsm.PMSMLUT.TorqueTable;
    else
        iq0index=find(pmsm.PMSMLUT.iqVec==0);
        if ~isempty(iq0index)
            for i1=1:length(pmsm.PMSMLUT.idVec)
                pmsm.PMSMLUT.FluxQTable(i1,:)=pmsm.PMSMLUT.FluxQTable(i1,:)-pmsm.PMSMLUT.FluxQTable(i1,iq0index);
            end
        end
        for ix=1:length(id)
            for iy=1:length(iq)
                T_computed(ix,iy)=1.5*pmsm.p*(iq(iy)*pmsm.PMSMLUT.FluxDTable(ix,iy)-id(ix)*pmsm.PMSMLUT.FluxQTable(ix,iy));%#ok
            end
        end
    end
else

    for ix=1:length(id)
        for iy=1:length(iq)
            T_computed(ix,iy)=1.5*pmsm.p*iq(iy)*(pmsm.FluxPM+(pmsm.Ld-pmsm.Lq)*id(ix));%#ok
        end
    end
end


hold on
linesplotted=0;
PMSM_NMAX=pmsm.N_max;

if goLUT==1 || goLUT==2
    if isfield(pmsm.PMSMLUT,'wrpmVec')
        wgrid=pmsm.PMSMLUT.wrpmVec;
        %wgrid=wgrid(wgrid<=PMSM_NMAX);
        wgrid=wgrid(wgrid>0);
    else
        N_max=min(PMSM_NMAX,Very_high_speed);
        wgrid=linspace(100,N_max,wsteps);
    end
else
    if isfield(pmsm,'PMSMLUT')
        if isfield(pmsm.PMSMLUT,'wrpmVec')
            wgrid=pmsm.PMSMLUT.wrpmVec;
            %wgrid=wgrid(wgrid<=PMSM_NMAX);
            wgrid=wgrid(wgrid>0);
        else
            N_max=min(PMSM_NMAX,Very_high_speed);
            wgrid=linspace(100,N_max,wsteps);
        end
    else
        N_max=min(PMSM_NMAX,Very_high_speed);
        wgrid=linspace(100,N_max,wsteps);
    end
end

for iw=1:length(wgrid)
    w=wgrid(iw)*pmsm.p*pi/30;
    for ix=1:length(id)
        for iy=1:length(iq)
            if goLUT==0
                vd=id(ix)*pmsm.Rs-w*iq(iy)*pmsm.Lq;
                vq=iq(iy)*pmsm.Rs+w*(id(ix)*pmsm.Ld+pmsm.FluxPM);
            elseif goLUT==1
                vd=id(ix)*pmsm.Rs-w*iq(iy)*pmsm.PMSMLUT.LqTable(ix,iy);
                vq=iq(iy)*pmsm.Rs+w*(id(ix)*pmsm.PMSMLUT.LdTable(ix,iy)+pmsm.PMSMLUT.FluxPMTable(ix,iy));
            elseif goLUT==2
                vd=id(ix)*pmsm.Rs-w*pmsm.PMSMLUT.FluxQTable(ix,iy);
                vq=iq(iy)*pmsm.Rs+w*pmsm.PMSMLUT.FluxDTable(ix,iy);
            end
            vs_cpmputed(ix,iy)=sqrt(vd^2+vq^2);%#ok
        end
    end

    % [M,c]=contour(X,Y,vs_cpmputed',[inverter.V_dc/sqrt(3) inverter.V_dc/sqrt(3)],'b:','ShowText','off','DisplayName',string(wgrid(iw)));
    % c.EdgeAlpha=0.3;
    M=contourc(id,iq,vs_cpmputed',[inverter.V_dc/sqrt(3) inverter.V_dc/sqrt(3)]);
    if ~isempty(M)
        M=M';
        if 2*(M(1,2)+1)~=numel(M)
            idtmp=[];
            iqtmp=[];i=1;
            while i<size(M,1)
                tmplen=M(i,2);
                idtmp1=M((i+1):(i+tmplen),1);
                iqtmp1=M((i+1):(i+tmplen),2);
                if ~isempty(idtmp1)
                    while(idtmp1(1)~=min(idtmp1))
                        idtmp1=[idtmp1(2:end); idtmp1(1)];
                        iqtmp1=[iqtmp1(2:end); iqtmp1(1)];
                    end
                end
                idtmp=[idtmp ; idtmp1];%#ok
                iqtmp=[iqtmp ; iqtmp1];%#ok
                i=i+tmplen+1;
            end
        else
            M=M(2:end,:);
            idtmp=M(:,1);
            iqtmp=M(:,2);
        end
        if changeDQAxes==1
            tmp=idtmp;
            idtmp=iqtmp;
            iqtmp=-tmp;
        end
        lh=plot(idtmp,iqtmp,'b:');
        thiscolor=lh.Color;
        lh.Color=[thiscolor(1:3) 0.3];
        linesplotted=linesplotted+1;
    end
end


PMSMRatedTorque=pmsm.T_rated;

if goLUT==1 || goLUT==2
    if isfield(pmsm.PMSMLUT,'trefVec')
        Tgrid=pmsm.PMSMLUT.trefVec;
    else
        TMax=PMSMRatedTorque;
        Tgrid=linspace(-TMax,TMax,Tsteps);
    end
else
    if isfield(pmsm,'PMSMLUT')
        if isfield(pmsm.PMSMLUT,'trefVec')
            Tgrid=pmsm.PMSMLUT.trefVec;
        else
            TMax=PMSMRatedTorque;
            Tgrid=linspace(-TMax,TMax,Tsteps);
        end
    else
        TMax=PMSMRatedTorque;
        Tgrid=linspace(-TMax,TMax,Tsteps);
    end
end

% linesplotted=linesplotted+2;
% [~,c]=contour(X,Y,T_computed',Tgrid,"r:",'ShowText','off','DisplayName',"");
% c.EdgeAlpha=0.3;
for it=1:length(Tgrid)
    thisT=Tgrid(it);
    M=contourc(id,iq,T_computed',[thisT thisT]);
    if ~isempty(M)
        M=M';
        if 2*(M(1,2)+1)~=numel(M)
            idtmp=[];
            iqtmp=[];i=1;
            while i<size(M,1)
                tmplen=M(i,2);
                idtmp1=M((i+1):(i+tmplen),1);
                iqtmp1=M((i+1):(i+tmplen),2);
                if ~isempty(idtmp1)
                    while(idtmp1(1)~=min(idtmp1))
                        idtmp1=[idtmp1(2:end); idtmp1(1)];
                        iqtmp1=[iqtmp1(2:end); iqtmp1(1)];
                    end
                end
                idtmp=[idtmp ; idtmp1];%#ok
                iqtmp=[iqtmp ; iqtmp1];%#ok
                i=i+tmplen+1;
            end
        else
            M=M(2:end,:);
            idtmp=M(:,1);
            iqtmp=M(:,2);
        end
        if changeDQAxes==1
            tmp=idtmp;
            idtmp=iqtmp;
            iqtmp=-tmp;
        end
        lh=plot(idtmp,iqtmp,'r:');
        thiscolor=lh.Color;
        lh.Color=[thiscolor(1:3) 0.3];
        linesplotted=linesplotted+1;
    end

end


if nargin==5
    x1=varargin{1};
    y1=varargin{2};
    plot(x1,y1,'r*')
    r = sqrt(x1^2+y1^2);
    c = [0 0];
    pos = [c-r 2*r 2*r];
    rectangle('Position',pos,'Curvature',[1 1])
end
r = pmsm.I_rated;
c = [0 0];
pos = [c-r 2*r 2*r];
rectangle('Position',pos,'Curvature',[1 1])
axis equal
if changeDQAxes==0
    xlim([-pmsm.I_rated 0]);
    ylim([0 pmsm.I_rated]);
    set(gca,'YAxisLocation','right');
else
    xlim([0 pmsm.I_rated]);
    ylim([0 pmsm.I_rated]);
    set(gca,'YAxisLocation','left');
end

xlabel('d-axis current (A)');
ylabel('q-axis current (A)');

end

%% a local function to interpolate
function valout=local_interp2(vec2,vec1,table1,val2,val1)
% a local function to limit the inputs to the bounds of the breakpoint vectors
% val1=max(min(vec1),min(max(vec1),val1));
if val1<vec1(1)
    val1=vec1(1);
elseif val1>vec1(end)
    val1=vec1(end);
end
% val2=max(min(vec2),min(max(vec2),val2));
if val2<vec2(1)
    val2=vec2(1);
elseif val2>vec2(end)
    val2=vec2(end);
end
valout=interp2(vec2,vec1,table1,val2,val1);
end
