MATLAB Examples

Contents

function OUTPUT = mfeval(PARAMETER_SOURCE,INPUTS,USE_MODE) %#codegen
%MFEVAL evaluates Magic Formula 6.1.2 in steady state for series of input
%variables.
%
% OUTPUT = mfeval(PARAMETER_SOURCE,INPUTS,USE_MODE)
%
% The formulation includes combined force/moment and turn slip calculation
% ISO coordinate system is used in all calculations.
%
% PARAMETER_SOURCE refers to a MF-Tyre tyre property file (.TIR) containing
% all the Magic Formula coefficients or to a structure with all the
% parameters.
%
% INPUTS = [FZ KAPPA ALPHA GAMMA PHIT VX P* OMEGA*], where
% FZ     = normal load on the tyre  [N]
% KAPPA  = longitudinal slip        [dimensionless, -1: locked wheel]
% ALPHA  = side slip angle          [rad]
% GAMMA  = inclination angle        [rad]
% PHIT   = turn slip                [1/m]
% VX     = forward velocity         [m/s]
% P*     = pressure                 [Pa]
% OMEGA* = rotational speed         [rad/s]
%
% P* and OMEGA* are optional inputs. If they are not included pressure is
% constant and equal to the nominal pressure on the TIR file and the
% rotational speed is approximated.
%
% USE_MODE specifies the type of calculation performed:
%      1: combined force/moment calculation
%      2: combined force/moment calculation + turn slip
%    +10: revoke alpha_star definition
%    +20: include alpha_star definition
%   +100: include limit checks
%   +200: ignore limit checks
%
% For example: USE_MODE = 121 implies:
%   -combined force/moment
%   -include alpha_star
%   -include limit checks
%
% For normal situations turn slip may be neglected, meaning that the radius
% of the path is close to infinity (straight line).
% Alpha_star improves the accuracy of the model for very large slip angles
% and possibly backward running of the wheel.
% The limit checks verify that the inputs are inside the stable range of
% the model and activates the low speed model.
%
% OUTPUT consists of 20 columns:
% 1 - Fx: longitudinal force         11 - Vx: longitudinal velocity
% 2 - Fy: lateral force              12 - P: pressure
% 3 - Fz: normal force               13 - Re: effective rolling radius
% 4 - Mx: overturning moment         14 - rho: tyre deflection
% 5 - My: rolling resistance moment  15 - a: contact patch length
% 6 - Mz: self aligning moment       16 - t: pneumatic trail
% 7 - kappa: longitudinal slip       17 - longitudinal friction coefficient
% 8 - alpha: side slip angle         18 - lateral friction coefficient
% 9 - gamma: inclination angle       19 - omega: rot. speed
% 10 - phit: turn slip               20 - Rl : loaded radius
%
% The equations here presented are published in the book:
% Title:    	Tire and Vehicle Dynamics
% Author:       Hans Pacejka
% Edition:      3, revised
% Publisher:	Elsevier, 2012
% ISBN:         0080970176, 9780080970172
% Length:       672 pages
%
% And in the following paper:
% Besselink, I. J. M. , Schmeitz, A. J. C. and Pacejka, H. B.(2010) 'An
% improved Magic Formula/Swift tyre model that can handle inflation
% pressure changes', Vehicle System Dynamics, 48: 1, 337 — 352
% Link to this Article: DOI: 10.1080/00423111003748088
% URL: http://dx.doi.org/10.1080/00423111003748088
%
% Script version:        1.4
% Date:                  03/07/2017
% Code developed by:     Marco Furlan
% Email address 1:       mfurlant@jaguarlandrover.com
% Email address 2:       marcofurlan92@gmail.com
%
% See the <a href="matlab:web(fullfile(mfevalRoot, 'html','index.html'))">documentation</a> for more details and examples

Versions change log

% Version	|   Date       |    Description
% -------       ----            -----------
% V1.0          25/01/2017      -First implementation
%
% V1.1          31/01/2017      -Eqn (4.E40) Bt modified.
%                               -Allow a structure of parameters to be
%                               passed in instead of a TIR file name.
%                               -Replace for loops in the coeff checks with
%                               vector expressions.
%                               -Rho calculation reviewed avoiding (NaNs)
%                               in external effects.
%                               -Coefficient checks can be disabled.
%                               -Coefficients of friction are checked
%                               before being exported.
%
% V1.2          16/02/2017      -Typo in the help text corrected.
%                               -Mux and Muy warnings can be deactivated
%                               with the coefficient checks.
%                               -Rolling resistance (4.E70), pneumatic
%                               trail (4.E73), Bt (4.E40), Dt (4.E43) and s
%                               (4.E76) equations have been modified to
%                               match TNO solver.
%
% V1.3          06/03/2017      -Code modified for code generation
%                               -Simulink blocks with transients have been
%                               included
%                               -Use mode implementation modified
%                               -Turn Slip equations modified
%                               -Pneumatic trail equation fixed
%
% V1.4          03/07/2017      -First implementation of the low speed mode
%                               -Publication of mfeval in MathWorks File
%                               Exchange

Future work

- Low speed model is not perfect, review the effect on pneumatic trail (t) and aligning torque (Mz) - Turn slip equations needs to be reviewed, specially for Mz

Declare extrinsic functions

If mfeval is used in a Simulink block, the code generation software will generate code for the call to an extrinsic function. Functions 'warning' and 'loadTIR' are not supported for code generation. By adding coder.extrinsic('warning') we can to bypass code generation. Type : coder.screener('mfeval') into the command window to check

coder.extrinsic('warning')
coder.extrinsic('loadTIR')

Check the arguments of the function

If the USE_MODE (3rd argument) is not specified, set it to 111

if nargin < 3
    USE_MODE = 111;
    warning('USE_MODE not defined, default mode 111 selected')
end

Apply Use Mode

Read the USE_MODE and apply the changes

if isempty(USE_MODE)
    USE_MODE = 111;
    warning('USE_MODE is empty, default mode 111 selected')
end

% Check the number of digits of USE_MODE
numdigits_USE_MODE = max(ceil(log10(abs(USE_MODE))),1);

if numdigits_USE_MODE~=3
    error('USE_MODE should contain exactly 3 digits')
end

% Save the elements of the string into separate variables
hundreds = floor(USE_MODE/100);
tens = floor((USE_MODE-hundreds*100)/10);
units = floor((USE_MODE-hundreds*100-tens*10));

switch hundreds
    case 1
        % 1: include limits checks
        useLimitsCheck = true;
    case 2
        % 2: revoke limits checks
        useLimitsCheck = false;
    otherwise
        error('USE_MODE not contemplated')
end

switch tens
    case 1
        % 1: revoke alpha_star definition
        useAlphaStar = false;
    case 2
        % 2: include alpha_star definition
        useAlphaStar = true;
    otherwise
        error('USE_MODE not contemplated')
end

switch units
    case 1
        % 1: combined force/moment calculation
        useTurnSlip = false;
    case 2
        % 2: combined force/moment calculation + turn slip
        useTurnSlip = true;
    otherwise
        error('USE_MODE not contemplated')
end

Inputs (Parameters)

Load the TIR file parameters

if(ischar(PARAMETER_SOURCE)) % File location
    TIRparam = loadTIR(PARAMETER_SOURCE); % Call an external function to
    % read all the parameters of the TIR file.

elseif(isstruct(PARAMETER_SOURCE))
    TIRparam = PARAMETER_SOURCE; % Structure of parameter
else
    error('PARAMETER_SOURCE not recognised')
end

% Low speed singularity
epsilon  = 1e-6; % [Eqn (4.E6a) Page 178 - Book]
epsilonv = epsilon;
epsilonx = epsilon;
epsilonk = epsilon;
epsilony = epsilon;
epsilonr = epsilon;

Inputs (Variables)

Unpack the input matrix into separate input variables

Fz       = INPUTS(:,1); % vertical load            (N)
kappa    = INPUTS(:,2); % longitudinal slip        (-) (-1 = locked wheel)
alpha    = INPUTS(:,3); % side slip angle          (radians)
gamma    = INPUTS(:,4); % inclination angle        (radians)
phit     = INPUTS(:,5); % turn slip                (1/m)
Vcx      = INPUTS(:,6); % forward velocity         (m/s)

% Check the size of the inputs
[mrows,ncolumns] = size(INPUTS);

% If the pressure is specified by the user, grab it.
% If is not, use the Inflation pressure of the TIR file
if ncolumns > 6
    p = INPUTS(:,7);                        % pressure     (Pa)
else
    p = ones(mrows,1).*TIRparam.INFLPRES;   % pressure     (Pa)
end

% IMPORTANT NOTE: Vx = Vcx [Eqn (7.4) Page 331 - Book]
% It is assumed that the difference between the wheel centre longitudinal
% velocity Vx and the longitudinal velocity Vcx of the contact centre is
% negligible

Limits

This code applies the appropriate limits to the input values to the model based on the MF6.1 limits specified in the TIR File

% Save unlimited variables for other calculations
Fz_unlimited       = Fz;
kappa_unlimited    = kappa;

if useLimitsCheck
    % Minimum Speed
    VXLOW = TIRparam.VXLOW;

    % Inflation Pressure Range
    PRESMIN = TIRparam.PRESMIN;
    PRESMAX = TIRparam.PRESMAX;

    % Vertical Force Range
    FZMIN = TIRparam.FZMIN;
    FZMAX = TIRparam.FZMAX;

    % Slip Angle Range
    ALPMIN = TIRparam.ALPMIN;
    ALPMAX = TIRparam.ALPMAX;

    % Inclination Angle Range
    CAMMIN = TIRparam.CAMMIN;
    CAMMAX = TIRparam.CAMMAX;

    % Long Slip Range
    KPUMIN = TIRparam.KPUMIN;
    KPUMAX = TIRparam.KPUMAX;

    % Speed limits (apply low slip mode)
    isLowSpeed = abs(Vcx) < VXLOW;

    % If Vcx is close to zero then SR is also 0 (linear digression)
    kappa(isLowSpeed) = kappa(isLowSpeed).*(Vcx(isLowSpeed)./VXLOW).*...
        sign(Vcx(isLowSpeed));

    % Create a reduction factor for forces and moments to get zero output
    % at zero speed
    fraction = Vcx(isLowSpeed)./TIRparam.VXLOW;
    reduction_factor = sin(fraction.*1.5708); % 90 deg into rad

    isLowAndNegative = Vcx(isLowSpeed)<0;
    reduction_factor(isLowAndNegative) = -reduction_factor(isLowAndNegative);

    % Turn slip modifications
    phit = phit.*cos(alpha); % Empirically discovered

    % If the speed is negative, the turn slip is also negative
    isNegativeSpeed = Vcx<0;
    phit(isNegativeSpeed) = -phit(isNegativeSpeed);

    % If Vcx is close to zero, use linear digression
    phit(isLowSpeed) = phit(isLowSpeed).*abs(fraction);

    % Check Slip Angle limits
    isLowSlip = alpha < ALPMIN;
    alpha(isLowSlip) = ALPMIN;

    isHighSlip = alpha > ALPMAX;
    alpha(isHighSlip) = ALPMAX;

    % Check camber limits
    isLowCamber = gamma < CAMMIN;
    gamma(isLowCamber) = CAMMIN;

    isHighCamber = gamma > CAMMAX;
    gamma(isHighCamber) = CAMMAX;

    % Check Fz limits
    Fz(Fz < 0) = 0;

    isLowLoad = Fz < FZMIN & Fz > 0;
    Fz(isLowLoad) = FZMIN;

    isHighLoad = Fz > FZMAX & Fz > 0;
    Fz_unlimited = Fz; % Some formulas use the unsaturated Fz in TNO
    Fz(isHighLoad) = FZMAX;

    % Check pressure limits
    isLowPressure = p < PRESMIN;
    p(isLowPressure) = PRESMIN;

    isHighPressure = p > PRESMAX;
    p(isHighPressure) = PRESMAX;

    % Check slip ratio limits
    isLowSlipR = kappa < KPUMIN;
    kappa(isLowSlipR) = KPUMIN;

    isHighSlipR = kappa > KPUMAX;
    kappa(isHighSlipR) = KPUMAX;

    % Flag if anything is out of range.
    chk =   isLowSlip  | isHighSlip | isLowCamber   | ...
        isHighCamber   | isLowLoad  | isHighLoad | isLowPressure | ...
        isHighPressure | isLowSlipR | isHighSlipR;

    if(any(chk))
        warning (['Limits of the input values to the model exceeded.',...
            'Values have been saturated.']);
    end

    if any(isLowSpeed)
        warning(['Speed input VX below the limit. ',...
            'Low speed mode activated.']);
    end
else % Not using limits checks
    % isLowSpeed = logical(zeros(length(Fz),1));
    isLowSpeed = false(length(Fz),1);
    reduction_factor = 1;
end

First calculations

% Basic calculations
[Fzo_prime, dfz, dpi, alpha_star, gamma_star, alpha_prime, LMUX_star, LMUY_star, LMUX_prime, LMUY_prime] = calculateBasic(TIRparam,useAlphaStar, Fz, kappa, alpha, gamma, Vcx, p, epsilonv);

% Effective rolling radius
[Re, Romega, omega] = calculateRe(TIRparam, ncolumns, INPUTS, dpi, Fz_unlimited, kappa_unlimited); % Unlimited inputs used to match TNO results

% Contact patch lengths
[a, ~] = calculateContactPatch(TIRparam, Fz_unlimited, dpi); % Unlimited Fz used to match TNO results

Turn slip calculations

Extension of the Model for Turn Slip, Page 183

[zetaStr] = calculateZeta(TIRparam, useTurnSlip, useLimitsCheck, Fz, alpha, gamma, kappa, Vcx, omega, phit, alpha_star, gamma_star, Fzo_prime, dfz, dpi, LMUY_star, LMUY_prime, epsilonr, epsilonk);

Fxo - Pure Longitudinal Force

(alpha = 0) Page 179 of the book

[Fxo, mux, Kxk] = calculateFxo(TIRparam, useLimitsCheck, Fz, kappa, gamma, Vcx, dfz, dpi, epsilonx, LMUX_star, LMUX_prime, zetaStr, isLowSpeed, reduction_factor);

Fyo - Pure Lateral Force

(kappa = 0) Page 179 of the book

[Fyo, muy, Kya, Kygo, SHy, SVy, By, Cy] = calculateFyo(TIRparam, useAlphaStar, useLimitsCheck, Fz, alpha_star, gamma_star, Fzo_prime, dfz, dpi, Vcx, epsilony, epsilonk, LMUY_star, LMUY_prime, zetaStr);

Mzo - Pure Aligning Torque

(kappa = 0) Page 180 of the book

[~, alphar, alphat, Dr, Cr, Br, Dt, Ct, Bt, Et, Kya_prime] = calculateMzo(TIRparam, useAlphaStar, useLimitsCheck, Fz, gamma, alpha_star, alpha_prime, gamma_star, Fzo_prime, dfz, dpi, Vcx, epsilonk, epsilony, Kya, Kygo, SHy, SVy, By, Cy, LMUY_star, LMUY_prime, zetaStr);

Fx - Combined Longitudinal Force

(Combined Slip), Page 181 of the book

[Fx] = calculateFx(TIRparam, useLimitsCheck, kappa, alpha_star, gamma_star, dfz, Fxo);

Fy - Combined Lateral Force

(Combined Slip), Page 181 of the book

[Fy, ~] = calculateFy(TIRparam, useLimitsCheck, Fz, alpha_star, kappa, gamma_star, dfz, Fyo, muy, zetaStr, isLowSpeed, reduction_factor);

Mx - Overturning Couple

(Also see Section 4.3.5), Page 182 of the book and Page 12 of the paper

[Mx] = calculateMx(TIRparam, Fz_unlimited, Fy, gamma, dpi);

My - Rolling Resistance Moment

(See Eqns (9.236, 9.230, 9.231)), Page 182 of the book

[My] = calculateMy(TIRparam, Fz_unlimited, Fx, Vcx, kappa, gamma, p, isLowSpeed, reduction_factor);

Mz - Aligning Torque

(Combined slip), Page 182 of the book

[Mz, t] = calculateMz(TIRparam, useAlphaStar, useLimitsCheck, alphar, alphat, Kxk, Kya_prime, kappa, Fz, Fy, Fx, Fzo_prime, dfz, gamma, Dr, Cr, Br, Dt, Ct, Bt, Et, alpha_star, Fyo, muy, alpha_prime, Vcx, dpi, epsilony, epsilonk, LMUY_star, LMUY_prime, isLowSpeed, reduction_factor);

Tyre Deflection and Loaded Radius

Methodology explained by Henning Olsson (henning.olsson@calspan.com)

[rho, Rl] = calculateRhoRl(TIRparam, omega, Fz_unlimited, Fx, Fy, dpi, Romega); % Unlimited Fz used to match TNO results

Generate the output matrix

% Check the sign of the coefficient of friction
% The calculation of Fx and Fy is not affected by the sign of mux and muy
if(useLimitsCheck)
    if any(mux < 0)
        mux = abs(mux);
        warning('Negative longitudinal coefficient of friction forced to be positive')
    end

    if any(muy < 0)
        muy = abs(muy);
        warning('Negative lateral coefficient of friction forced to be positive')
    end
end

% Preallocate OUTPUT variable
OUTPUT = zeros(mrows,20);

% Pack all the outputs
OUTPUT(:,1) = Fx;
OUTPUT(:,2) = Fy;
OUTPUT(:,3) = Fz_unlimited;
OUTPUT(:,4) = Mx;
OUTPUT(:,5) = My;
OUTPUT(:,6) = Mz;
OUTPUT(:,7) = kappa;
OUTPUT(:,8) = alpha;
OUTPUT(:,9) = gamma;
OUTPUT(:,10) = phit;
OUTPUT(:,11) = Vcx;
OUTPUT(:,12) = p;
OUTPUT(:,13) = Re;
OUTPUT(:,14) = rho;
OUTPUT(:,15) = 2*a;
OUTPUT(:,16) = t;
OUTPUT(:,17) = mux;
OUTPUT(:,18) = muy;
OUTPUT(:,19) = omega;
OUTPUT(:,20) = Rl;
end

Supplementary functions

Functions to calculate forces and moments

function [Fzo_prime, dfz, dpi, alpha_star, gamma_star, alpha_prime, LMUX_star, LMUY_star, LMUX_prime, LMUY_prime] = calculateBasic(TIRparam,useAlphaStar, Fz, kappa, alpha, gamma, Vcx, p, epsilonv)%#codegen

%[MODEL]
LONGVL 	= TIRparam.LONGVL 	; %Nominal speed

%[OPERATING_CONDITIONS]
NOMPRES	= TIRparam.NOMPRES	; %Nominal tyre inflation pressure

%[VERTICAL]
FNOMIN	= TIRparam.FNOMIN 	; %Nominal wheel load

%[SCALING_COEFFICIENTS]
LFZO	= TIRparam.LFZO     ; % Scale factor of nominal (rated) load
LMUX	= TIRparam.LMUX   	; % Scale factor of Fx peak friction coefficient
LMUY	= TIRparam.LMUY    	; % Scale factor of Fy peak friction coefficient

% New scaling factor in Pacejka 2012 with it's default value. This
% scaling factor is not present in the standard MF6.1 TIR files.
LMUV    = 0 ; % Scale factor with slip speed Vs decaying friction

% Rename the TIR file variables in the Pacejka style
Fzo = FNOMIN;               % Nominal (rated) wheel load
pio = NOMPRES;              % Reference pressure
Vo = LONGVL;                % Nominal speed (LONGVL)

% Some basic calculations are done before starting to calculate forces and
% moments.

% (1) Velocities in point S (slip)
Vsx = -kappa.*abs(Vcx); % [Eqn (4.E5) Page 181 - Book]
Vsy = -tan(alpha).*abs(Vcx); % [Eqn (2.12) Page 67 - Book] and [(4.E3) Page 177 - Book]
Vs = sqrt(Vsx.^2+Vsy.^2); % [Eqn (3.39) Page 102 - Book] -> Slip velocity of the slip point S

% (2) Velocities in point C (contact)
Vcy = Vsy; % Assumption from page 67 of the book, paragraph above Eqn (2.11)
Vc = sqrt(Vcx.^2+Vcy.^2); % Velocity of the wheel contact centre C, Not described in the book but is the same as [Eqn (3.39) Page 102 - Book]

% (3) Effect of having a tire with a different nominal load
Fzo_prime =  LFZO.*Fzo; % [Eqn (4.E1) Page 177 - Book]

% (4) Normalized change in vertical load
dfz = (Fz - Fzo_prime)./Fzo_prime; % [Eqn (4.E2a) Page 177 - Book]

% (5) Normalized change in inflation pressure
dpi = (p - pio)./pio; % [Eqn (4.E2b) Page 177 - Book]

% (6) Use of alpha_star definition (Large slip angles correction)
if useAlphaStar
    alpha_star = tan(alpha).*sign(Vcx); % [Eqn (4.E3) Page 177 - Book]
else
    alpha_star = alpha; % To match TNO results
end

% (7) Spin due to the camber angle
gamma_star = sin(gamma); % [Eqn (4.E4) Page 177 - Book]

% (8) For the aligning torque at high slip angles
Vc_prime = Vc + epsilonv.*sign(Vc); % [Eqn (4.E6a) Page 178 - Book] [sign(Vc) term explained on page 177]
alpha_prime = acos(Vcx./Vc_prime); % [Eqn (4.E6) Page 177 - Book]

% (9) Slippery surface with friction decaying with increasing (slip) speed
LMUX_star = LMUX./(1 + LMUV.*Vs./Vo); % [Eqn (4.E7) Page 179 - Book]
LMUY_star = LMUY./(1 + LMUV.*Vs./Vo); % [Eqn (4.E7) Page 179 - Book]

% (10) Digressive friction factor
% On Page 179 of the book is suggested Amu=10, but after comparing the use
% of the scaling factors against TNO, Amu = 1 was giving perfect match
Amu = 1;
LMUX_prime = Amu.*LMUX_star./(1+(Amu-1).*LMUX_star); % [Eqn (4.E8) Page 179 - Book]
LMUY_prime = Amu.*LMUY_star./(1+(Amu-1).*LMUY_star); % [Eqn (4.E8) Page 179 - Book]
end

function [zetaStr] = calculateZeta(TIRparam, useTurnSlip, useLimitsCheck, Fz, alpha, gamma, kappa, Vcx, omega, phit, alpha_star, gamma_star, Fzo_prime, dfz, dpi, LMUY_star, LMUY_prime, epsilonr, epsilonk)%#codegen
% Declare extrinsic functions
coder.extrinsic('warning')

if useTurnSlip

    %[SCALING_COEFFICIENTS]
    LKY                     = TIRparam.LKY             		; % Scale factor of Fy cornering stiffness
    LYKA                    = TIRparam.LYKA             	; % Scale factor of alpha influence on Fx
    LKYC                    = TIRparam.LKYC             	; % Scale factor of camber force stiffness
    LKZC                    = TIRparam.LKZC             	; % %Scale factor of camber torque stiffness
    LMP                     = TIRparam.LMP             		; % Scale factor of Parking Moment

    %[TURNSLIP_COEFFICIENTS]
    PDXP1                   = TIRparam.PDXP1           		; %Peak Fx reduction due to spin parameter
    PDXP2                   = TIRparam.PDXP2            	; %Peak Fx reduction due to spin with varying load parameter
    PDXP3                   = TIRparam.PDXP3             	; %Peak Fx reduction due to spin with kappa parameter
    PKYP1                   = TIRparam.PKYP1             	; %Cornering stiffness reduction due to spin
    PDYP1                   = TIRparam.PDYP1           		; %Peak Fy reduction due to spin parameter
    PDYP2                   = TIRparam.PDYP2             	; %Peak Fy reduction due to spin with varying load parameter
    PDYP3                   = TIRparam.PDYP3             	; %Peak Fy reduction due to spin with alpha parameter
    PDYP4                   = TIRparam.PDYP4            	; %Peak Fy reduction due to square root of spin parameter
    PHYP1                   = TIRparam.PHYP1             	; %Fy-alpha curve lateral shift limitation
    PHYP2                   = TIRparam.PHYP2          		; %Fy-alpha curve maximum lateral shift parameter
    PHYP3                   = TIRparam.PHYP3            	; %Fy-alpha curve maximum lateral shift varying with load parameter
    PHYP4                   = TIRparam.PHYP4            	; %Fy-alpha curve maximum lateral shift parameter
    PECP1                   = TIRparam.PECP1           		; %Camber w.r.t. spin reduction factor parameter in camber stiffness
    PECP2                   = TIRparam.PECP2            	; %Camber w.r.t. spin reduction factor varying with load parameter in camber stiffness
    QDTP1                   = TIRparam.QDTP1            	; %Pneumatic trail reduction factor due to turn slip parameter
    QCRP1                   = TIRparam.QCRP1           		; %Turning moment at constant turning and zero forward speed parameter
    QCRP2                   = TIRparam.QCRP2           		; %Turn slip moment (at alpha=90deg) parameter for increase with spin
    QBRP1                   = TIRparam.QBRP1           		; %Residual (spin) torque reduction factor parameter due to side slip
    QDRP1                   = TIRparam.QDRP1             	; %Turn slip moment peak magnitude parameter

    if isfield(TIRparam,'QDRP2')==0 %Coefficient not found in the TIR file
        QDRP2               = 0                             ; %Turn slip moment peak position parameter
    else
        QDRP2               = TIRparam.QDRP2             	; %Turn slip moment peak position parameter
    end

    %[LATERAL_COEFFICIENTS]
    PDY1                    =  TIRparam.PDY1    			; %Lateral friction Muy
    PDY2                    =  TIRparam.PDY2  				; %Variation of friction Muy with load
    PDY3                    =  TIRparam.PDY3     			; %Variation of friction Muy with squared camber
    PKY1                    =  TIRparam.PKY1    			; %Maximum value of stiffness Kfy./Fznom
    PKY2                    =  TIRparam.PKY2     			; %Load at which Kfy reaches maximum value
    PKY3                    =  TIRparam.PKY3    			; %Variation of Kfy./Fznom with camber
    PKY4                    =  TIRparam.PKY4     			; %Curvature of stiffness Kfy
    PKY5                    =  TIRparam.PKY5   				; %Peak stiffness variation with camber squared
    PKY6                    =  TIRparam.PKY6   				; %Fy camber stiffness factor
    PKY7                    =  TIRparam.PKY7   				; %Vertical load dependency of camber stiffness
    PVY3                    =  TIRparam.PVY3   				; %Variation of shift Svy./Fz with camber
    PVY4                    =  TIRparam.PVY4  				; %Variation of shift Svy./Fz with camber and load
    PPY1                    =  TIRparam.PPY1    			; %influence of inflation pressure on cornering stiffness
    PPY2                    =  TIRparam.PPY2     			; %influence of inflation pressure on dependency of nominal tyre load on cornering stiffness
    PPY3                    =  TIRparam.PPY3   				; %linear influence of inflation pressure on lateral peak friction
    PPY4                    =  TIRparam.PPY4    			; %quadratic influence of inflation pressure on lateral peak friction
    PPY5                    =  TIRparam.PPY5   				; %Influence of inflation pressure on camber stiffness
    RBY1                    =  TIRparam.RBY1     			; %Slope factor for combined Fy reduction
    RBY2                    =  TIRparam.RBY2    			; %Variation of slope Fy reduction with alpha
    RBY3                    =  TIRparam.RBY3  				; %Shift term for alpha in slope Fy reduction
    RBY4                    =  TIRparam.RBY4     			; %Influence of camber on stiffness of Fy combined
    RCY1                    =  TIRparam.RCY1    			; %Shape factor for combined Fy reduction
    REY1                    =  TIRparam.REY1    			; %Curvature factor of combined Fy
    REY2                    =  TIRparam.REY2    			; %Curvature factor of combined Fy with load
    RHY1                    =  TIRparam.RHY1  				; %Shift factor for combined Fy reduction
    RHY2                    =  TIRparam.RHY2   				; %Shift factor for combined Fy reduction with load

    %[ALIGNING_COEFFICIENTS]
    QDZ8                    =  TIRparam.QDZ8   				; %Variation of peak factor Dmr" with camber
    QDZ9                    =  TIRparam.QDZ9   				; %Variation of peak factor Dmr" with camber and load
    QDZ10                   =  TIRparam.QDZ10     			; %Variation of peak factor Dmr with camber squared
    QDZ11                   =  TIRparam.QDZ11     			; %Variation of Dmr with camber squared and load

    %[DIMENSION]
    Ro                      = TIRparam.UNLOADED_RADIUS      ; %Free tyre radius

    % Speed limits (avoid zero speed)
    Vc_prime = Vcx;
    isLowSpeed = abs(Vcx) < TIRparam.VXLOW;
    signVcx = sign(Vcx(isLowSpeed));
    signVcx(signVcx == 0) = 1;
    Vc_prime(isLowSpeed) = TIRparam.VXLOW * signVcx; % Singularity protected velocity, text on Page 184

    psi_dot = phit.*Vc_prime; % Rearrange [Eqn (4.75) Page 183 - Book]

    % IMPORTANT NOTE: Eqn (4.75) has been modified
    % In chapter 2.2 "Definition of tire input quantities" in the Pacejka
    % book, it is assumed that the z-axis of the road coordinate system
    % "points downward normal to the road plane" (p. 62). Due to this
    % definition, Pacejka introduces the minus sign for the spin slip so
    % that a positive spin leads to a positive torque Mz (p. 68).But in
    % CarMaker (and other MBS software), the road coordinate system is
    % orientated differently. The z-axis points upward to the
    % road plane. Thus switching the signs is not needed here.

    epsilong = PECP1.*(1 + PECP2.*dfz); % [Eqn (4.90) Page 186 - Book] Camber reduction factor
    phi = (-1./Vc_prime).*(psi_dot-(1-epsilong).*omega.*sin(gamma)); % [Eqn (4.76) Page 184 - Book]

    zeta3 = cos(atan(PKYP1.*Ro.^2.*phi.^2)); % [Eqn (4.79) Page 185 - Book]

    Kygo = Fz.*(PKY6 + PKY7 .*dfz).*(1 + PPY5.*dpi).*LKYC; % (=dFyo./dgamma at alpha = gamma = 0) (= CFgamma) (4.E30)
    Kya = PKY1.*Fzo_prime.*(1 + PPY1.*dpi).*(1 - PKY3.*abs(gamma_star)).*sin(PKY4.*atan((Fz./Fzo_prime)./((PKY2+PKY5.*gamma_star.^2).*(1+PPY2.*dpi)))).*zeta3.*LKY; % (= ByCyDy = dFyo./dalphay at alphay = 0) (if gamma =0: =Kyao = CFa) (PKY4=2)(4.E15)
    Kyao = PKY1.*Fzo_prime.*(1 + PPY1.*dpi).*(1 - PKY3.*abs(0)).*sin(PKY4.*atan((Fz./Fzo_prime)./((PKY2+PKY5.*0.^2).*(1+PPY2.*dpi)))).*zeta3.*LKY;

    % IMPORTANT NOTE: Explanation of the above equation, Kyao
    % Kya0 is the cornering stiffness when the camber angle is zero
    % (gamma=0) which is again the product of the coefficients By, Cy and
    % Dy at zero camber angle. Information from Kaustub Ragunathan, email:
    % carmaker-service-uk@ipg-automotive.com

    signKya = sign(Kya);
    signKya(signKya == 0) = 1; % If [Kya = 0] then [sign(0) = 0]. This is done to avoid [num / 0 = NaN] in Eqn 4.E39

    signKyao = sign(Kyao);
    signKyao(signKyao == 0) = 1; % If [Kyao = 0] then [sign(0) = 0]. This is done to avoid [num / 0 = NaN]

    Kya_prime = Kya + epsilonk.*signKya; % (4.E39) [sign(Kya) term explained on page 177]
    Kyao_prime = Kyao + epsilonk.*signKyao; % epsilonk is a small factor added to avoid the singularity condition during zero velocity (equation 308, CarMaker reference Manual).

    muy = (PDY1 + PDY2 .* dfz).*(1 + PPY3.*dpi + PPY4 .*dpi.^2).*(1 - PDY3.*gamma_star.^2).*LMUY_star; % (4.E23)
    Cyk = RCY1; % (4.E63)
    Byk = (RBY1 + RBY4.*gamma_star.^2).*cos(atan(RBY2.*(alpha_star - RBY3))).*LYKA; % (> 0) (4.E62)
    Eyk = REY1 + REY2.*dfz; % (<=1) (4.E64)
    SHyk = RHY1 + RHY2.*dfz; % (4.E65)
    Gyko = cos(Cyk.*atan(Byk.*SHyk - Eyk.*(Byk.*SHyk - atan(Byk.*SHyk)))); % (4.E60)
    Gyk_k = cos(Cyk.*atan(Byk.*kappa - Eyk.*(Byk.*kappa-atan(Byk.*kappa))))./Gyko; % (> 0)(4.E59)

    Byp = PDYP1.*(1 + PDYP2.*dfz).*cos(atan(PDYP3.*tan(alpha))); % [Eqn (4.78) Page 185 - Book]
    Bxp = PDXP1.*(1 + PDXP2.*dfz).*cos(atan(PDXP3.*kappa)); % [Eqn (4.106) Page 188 - Book]
    CHyp = PHYP1; % (>0) [Eqn (4.85) Page 186 - Book]
    DHyp = (PHYP2 + PHYP3.*dfz).*sign(Vcx); % [Eqn (4.86) Page 186 - Book]
    EHyp = PHYP4; % (<=1) [Eqn (4.87) Page 186 - Book]

    % Limits check
    if(useLimitsCheck)
        if EHyp > 1
            EHyp=1;
            warning('Ex over limit (>1), Eqn(4.E14)')
        end
    end

    KyRpo = Kygo./(1-epsilong); % Eqn (4.89)
    BHyp = KyRpo./(CHyp.*DHyp.*Kyao_prime); %[Eqn (4.88) Page 186 - Book]
    SHyp = DHyp.*sin(CHyp .* atan(BHyp.*Ro.*phi - EHyp.*(BHyp.*Ro.*phi - atan(BHyp.*Ro.*phi)))).*sign(Vcx); % [Eqn (4.80) Page 185 - Book]

    Mzp_inf = QCRP1.*muy.*Ro.*Fz.*sqrt(Fz./Fzo_prime).*LMP; % [Eqn (4.95) Page 187 - Book]

    % Mzp_inf should be always > 0
    negative_Mzp_inf = Mzp_inf < 0;
    Mzp_inf(negative_Mzp_inf) = 1e-6;

    Mzp90 = Mzp_inf.*(2./pi).*atan(QCRP2.*Ro.*abs(phit)).*Gyk_k; % [Eqn (4.103) Page 188 - Book]

    CDrp = QDRP1; % (>0) [Eqn (4.96) Page 187 - Book]
    EDrp = QDRP2; % (<=1) [Eqn (4.97) Page 187 - Book]
    DDrp = Mzp_inf./sin(0.5.*pi.*CDrp); % [Eqn (4.94) Page 187 - Book]
    Kzgro = Fz.*Ro.*(QDZ8.*QDZ9.*dfz + (QDZ10 + QDZ11.*dfz.*abs(gamma))).*LKZC; %[Eqn (4.99) Page 187 - Book]
    BDrp = Kzgro./(CDrp.*DDrp.*(1 - epsilong) + epsilonr.*sign(1 - epsilong)); % [Eqn (4.98) Page 187 - Book] [sign(1 - epsilong) term is explained on page 177]
    Drp = DDrp.*sin(CDrp.*atan(BDrp.*Ro.*(-phi) - EDrp.*(BDrp.*Ro.*(-phi) - atan(BDrp.*Ro.*(-phi))))); % Eqn (4.93) Page 187 - Book]

    % IMPORTANT NOTE:
    % Drp equation modified to operate with -phi, so that a positive spin
    % leads to a positive torque Mz (p. 68). Original equation below:
    % Drp = DDrp.*sin(CDrp.*atan(BDrp.*Ro.*(phi) - EDrp.*(BDrp.*Ro.*(phi) - atan(BDrp.*Ro.*(phi))))); % Eqn (4.93) Page 187 - Book]

    zeta0 = 0; % [Eqn (4.83) Page 186 - Book]
    zeta1 = cos(atan(Bxp.*Ro.*phi)); % [Eqn (4.105) Page 188 - Book]
    zeta2 = cos(atan(Byp.*(Ro.*abs(phi)+ PDYP4.*sqrt(Ro.*abs(phi))))); % [Eqn (4.77) Page 184 - Book]
    SVyg = Fz.*(PVY3 + PVY4.*dfz).*gamma_star.* LKYC .* LMUY_prime .* zeta2; % (4.E28)
    zeta4 = 1 + SHyp - SVyg./Kya_prime; % [Eqn (4.84) Page 186 - Book]
    zeta5 = cos(atan(QDTP1.*Ro.*phi)); %[Eqn (4.91) Page 186 - Book]
    zeta6 = cos(atan(QBRP1.*Ro.*phi)); % [Eqn (4.102) Page 188 - Book]

    signDrp = sign(Drp);
    signDrp(signDrp == 0) = 1; % If [Drp = 0] then [sign(0) = 0]. This is done to avoid [num / 0 = NaN] in Eqn 4.104

    x = Mzp90./(-abs(Drp)+epsilonr.*signDrp); % Argument of Eqn (4.104) [sign(Drp) term is explained on page 177]

    % Text from Pacejka book page 188: (In any case, the argument of
    % (4.104) should remain <1.) This is because the inputs of "acos"
    % shpuld be between +1 to -1 to aboid complex solutions.

    if any(x > 1) || any(x < -1)
        isAbove = x>1;
        isBelow = x<-1;
        x(isAbove) = 1;
        x(isBelow) = -1;
    end

    zeta7 = (2./pi).*acos(x); % [Eqn (4.104) Page 188 - Book]

    % IMPORTANT NOTE:
    % zeta7 equation modified to operate with -(abs(Drp)) so that a
    % positive spin leads to a positive torque Mz (p. 68). Original
    % equation below:
    % zeta7 = (2./pi).*acos(Mzp90./(abs(Drp)+epsilonr.*signDrp)); % [Eqn (4.104) Page 188 - Book]

    zeta8 = 1 + Drp;

else
    % No turn slip and small camber angles
    % First paragraph on page 178 of the book
    zeta0 = 1;
    zeta1 = 1;
    zeta2 = 1;
    zeta3 = 1;
    zeta4 = 1;
    zeta5 = 1;
    zeta6 = 1;
    zeta7 = 1;
    zeta8 = 1;
end

% Pack zeta in a structure
zetaStr.zeta0 = zeta0;
zetaStr.zeta1 = zeta1;
zetaStr.zeta2 = zeta2;
zetaStr.zeta3 = zeta3;
zetaStr.zeta4 = zeta4;
zetaStr.zeta5 = zeta5;
zetaStr.zeta6 = zeta6;
zetaStr.zeta7 = zeta7;
zetaStr.zeta8 = zeta8;

end

function [Fxo, mux, Kxk] = calculateFxo(TIRparam, useLimitsCheck, Fz, kappa, gamma, Vcx, dfz, dpi, epsilonx, LMUX_star, LMUX_prime, zetaStr, isLowSpeed, reduction_factor)%#codegen
% Declare extrinsic functions
coder.extrinsic('warning')
% Unpack zeta structure
zeta1 = zetaStr.zeta1;

%[SCALING_COEFFICIENTS]
LCX                     = TIRparam.LCX             		; % Scale factor of Fx shape factor
LEX                     = TIRparam.LEX             		; % Scale factor of Fx curvature factor
LKX                     = TIRparam.LKX             		; % Scale factor of Fx slip stiffness
LHX                     = TIRparam.LHX             		; % Scale factor of Fx horizontal shift
LVX                     = TIRparam.LVX            		; % Scale factor of Fx vertical shift

%[LONGITUDINAL_COEFFICIENTS]
PCX1                    =  TIRparam.PCX1     			; %Shape factor Cfx for longitudinal force
PDX1                    =  TIRparam.PDX1     			; %Longitudinal friction Mux at Fznom
PDX2                    =  TIRparam.PDX2   				; %Variation of friction Mux with load
PDX3                    =  TIRparam.PDX3   				; %Variation of friction Mux with camber squared
PEX1                    =  TIRparam.PEX1  				; %Longitudinal curvature Efx at Fznom
PEX2                    =  TIRparam.PEX2  				; %Variation of curvature Efx with load
PEX3                    =  TIRparam.PEX3  				; %Variation of curvature Efx with load squared
PEX4                    =  TIRparam.PEX4    			; %Factor in curvature Efx while driving
PKX1                    =  TIRparam.PKX1     			; %Longitudinal slip stiffness Kfx./Fz at Fznom
PKX2                    =  TIRparam.PKX2     			; %Variation of slip stiffness Kfx./Fz with load
PKX3                    =  TIRparam.PKX3   				; %Exponent in slip stiffness Kfx./Fz with load
PHX1                    =  TIRparam.PHX1 				; %Horizontal shift Shx at Fznom
PHX2                    =  TIRparam.PHX2 				; %Variation of shift Shx with load
PVX1                    =  TIRparam.PVX1   				; %Vertical shift Svx./Fz at Fznom
PVX2                    =  TIRparam.PVX2   				; %Variation of shift Svx./Fz with load
PPX1                    =  TIRparam.PPX1 				; %linear influence of inflation pressure on longitudinal slip stiffness
PPX2                    =  TIRparam.PPX2				; %quadratic influence of inflation pressure on longitudinal slip stiffness
PPX3                    =  TIRparam.PPX3				; %linear influence of inflation pressure on peak longitudinal friction
PPX4                    =  TIRparam.PPX4				; %quadratic influence of inflation pressure on peak longitudinal friction

Cx = PCX1.*LCX; % (> 0) (4.E11)
mux = (PDX1 + PDX2.*dfz).*(1 + PPX3.*dpi + PPX4.*dpi.^2).*(1 - PDX3.*gamma.^2).*LMUX_star; % (4.E13)
Dx = mux.*Fz.*zeta1; % (> 0) (4.E12)
Kxk = Fz.*(PKX1 + PKX2.*dfz).*exp(PKX3.*dfz).*(1 + PPX1.*dpi + PPX2.*dpi.^2).*LKX;  % (= BxCxDx = dFxo./dkx at kappax = 0) (= Cfk) (4.E15)

signDx = sign(Dx);
signDx(signDx == 0) = 1; % If [Dx = 0] then [sign(0) = 0]. This is done to avoid [Kxk / 0 = NaN] in Eqn 4.E16

Bx = Kxk./(Cx.*Dx + epsilonx.*signDx); % (4.E16) [sign(Dx) term explained on page 177]
SHx = (PHX1 + PHX2.*dfz).*LHX; % (4.E17)
SVx = Fz.*(PVX1 + PVX2.*dfz).*LVX.*LMUX_prime.*zeta1; % (4.E18)
kappax = kappa + SHx; % (4.E10)
Ex = (PEX1 + PEX2.*dfz + PEX3.*dfz.^2).*(1 - PEX4.*sign(kappax)).*LEX; % (<=1) (4.E14)

% Limits check
if(useLimitsCheck)
    if(any(Ex > 1))
        warning('Ex over limit (>1), Eqn(4.E14)')
        Ex(Ex > 1) = 1;
    end
end

Fxo = Dx.*sin(Cx.*atan(Bx.*kappax-Ex.*(Bx.*kappax-atan(Bx.*kappax))))+SVx; % (4.E9)

% Backward speed check
Fxo(Vcx < 0) = -Fxo(Vcx < 0);

% Low speed model
if any(isLowSpeed)
    Fxo(isLowSpeed) = Fxo(isLowSpeed).* reduction_factor ;
end
end

function [Fyo, muy, Kya, Kygo, SHy, SVy, By, Cy] = calculateFyo(TIRparam, useAlphaStar, useLimitsCheck, Fz, alpha_star, gamma_star, Fzo_prime, dfz, dpi, Vcx, epsilony, epsilonk, LMUY_star, LMUY_prime, zetaStr)%#codegen
% Declare extrinsic functions
coder.extrinsic('warning')
% Unpack zeta structure
zeta0= zetaStr.zeta0;
zeta2= zetaStr.zeta2;
zeta3= zetaStr.zeta3;
zeta4= zetaStr.zeta4;

%[SCALING_COEFFICIENTS]
LCY                     = TIRparam.LCY             		; % Scale factor of Fy shape factor
LEY                     = TIRparam.LEY             		; % Scale factor of Fy curvature factor
LKY                     = TIRparam.LKY             		; % Scale factor of Fy cornering stiffness
LHY                     = TIRparam.LHY             		; % Scale factor of Fy horizontal shift
LVY                     = TIRparam.LVY             		; % Scale factor of Fy vertical shift
LKYC                    = TIRparam.LKYC             	; % Scale factor of camber force stiffness

%[LATERAL_COEFFICIENTS]
PCY1                    =  TIRparam.PCY1     			; %Shape factor Cfy for lateral forces
PDY1                    =  TIRparam.PDY1    			; %Lateral friction Muy
PDY2                    =  TIRparam.PDY2  				; %Variation of friction Muy with load
PDY3                    =  TIRparam.PDY3     			; %Variation of friction Muy with squared camber
PEY1                    =  TIRparam.PEY1   				; %Lateral curvature Efy at Fznom
PEY2                    =  TIRparam.PEY2    			; %Variation of curvature Efy with load
PEY3                    =  TIRparam.PEY3   				; %Zero order camber dependency of curvature Efy
PEY4                    =  TIRparam.PEY4    			; %Variation of curvature Efy with camber
PEY5                    =  TIRparam.PEY5				; %Variation of curvature Efy with camber squared
PKY1                    =  TIRparam.PKY1    			; %Maximum value of stiffness Kfy./Fznom
PKY2                    =  TIRparam.PKY2     			; %Load at which Kfy reaches maximum value
PKY3                    =  TIRparam.PKY3    			; %Variation of Kfy./Fznom with camber
PKY4                    =  TIRparam.PKY4     			; %Curvature of stiffness Kfy
PKY5                    =  TIRparam.PKY5   				; %Peak stiffness variation with camber squared
PKY6                    =  TIRparam.PKY6   				; %Fy camber stiffness factor
PKY7                    =  TIRparam.PKY7   				; %Vertical load dependency of camber stiffness
PHY1                    =  TIRparam.PHY1  				; %Horizontal shift Shy at Fznom
PHY2                    =  TIRparam.PHY2  				; %Variation of shift Shy with load
PVY1                    =  TIRparam.PVY1   				; %Vertical shift in Svy./Fz at Fznom
PVY2                    =  TIRparam.PVY2    			; %Variation of shift Svy./Fz with load
PVY3                    =  TIRparam.PVY3   				; %Variation of shift Svy./Fz with camber
PVY4                    =  TIRparam.PVY4  				; %Variation of shift Svy./Fz with camber and load
PPY1                    =  TIRparam.PPY1    			; %influence of inflation pressure on cornering stiffness
PPY2                    =  TIRparam.PPY2     			; %influence of inflation pressure on dependency of nominal tyre load on cornering stiffness
PPY3                    =  TIRparam.PPY3   				; %linear influence of inflation pressure on lateral peak friction
PPY4                    =  TIRparam.PPY4    			; %quadratic influence of inflation pressure on lateral peak friction
PPY5                    =  TIRparam.PPY5   				; %Influence of inflation pressure on camber stiffness

Kya = PKY1.*Fzo_prime.*(1 + PPY1.*dpi).*(1 - PKY3.*abs(gamma_star)).*sin(PKY4.*atan((Fz./Fzo_prime)./((PKY2+PKY5.*gamma_star.^2).*(1+PPY2.*dpi)))).*zeta3.*LKY; % (= ByCyDy = dFyo./dalphay at alphay = 0) (if gamma =0: =Kyao = CFa) (PKY4=2)(4.E25)
SVyg = Fz.*(PVY3 + PVY4.*dfz).*gamma_star.* LKYC .* LMUY_prime .* zeta2; % (4.E28)
Kygo = Fz.*(PKY6 + PKY7 .*dfz).*(1 + PPY5.*dpi).*LKYC; % (=dFyo./dgamma at alpha = gamma = 0) (= CFgamma) (4.E30)

signKya = sign(Kya);
signKya(signKya == 0) = 1; % If [Kya = 0] then [sign(0) = 0]. This is done to avoid [num / 0 = NaN] in Eqn 4.E27

SHy = (PHY1 + PHY2.*dfz).* LHY + ((Kygo .*gamma_star - SVyg)./(Kya + epsilonk.*signKya)).*zeta0 +zeta4 -1; % (4.E27) [sign(Kya) term explained on page 177]
SVy = Fz.*(PVY1 + PVY2.*dfz).*LVY.*LMUY_prime.*zeta2 + SVyg; % (4.E29)
alphay = alpha_star + SHy; % (4.E20)
Cy = PCY1.*LCY; % (> 0) (4.E21)
muy = (PDY1 + PDY2 .* dfz).*(1 + PPY3.*dpi + PPY4 .*dpi.^2).*(1 - PDY3.*gamma_star.^2).*LMUY_star; % (4.E23)
Dy = muy.*Fz.*zeta2; % (4.E22)
signAlphaY = sign(alphay);
signAlphaY(signAlphaY == 0) = 1;
Ey = (PEY1 + PEY2.*dfz).*(1 + PEY5.*gamma_star.^2 - (PEY3 + PEY4.*gamma_star).*signAlphaY).*LEY; % (<=1)(4.E24)

% Limits check
if(useLimitsCheck)
    if(any(Ey > 1))
        warning('Ey over limit (>1), Eqn(4.E24)')
        Ey(Ey > 1) = 1;
    end
end

signDy = sign(Dy);
signDy(signDy == 0) = 1; % If [Dy = 0] then [sign(0) = 0]. This is done to avoid [Kya / 0 = NaN] in Eqn 4.E26

By = Kya./(Cy.*Dy + epsilony.*signDy); % (4.E26) [sign(Dy) term explained on page 177]
Fyo = Dy .* sin(Cy.*atan(By.*alphay-Ey.*(By.*alphay - atan(By.*alphay))))+ SVy; % (4.E19)

% Backward speed check
% Only for alpha_star
if(useAlphaStar)
    Fyo(Vcx < 0) = -Fyo(Vcx < 0);
end

end

function [Mzo, alphar, alphat, Dr, Cr, Br, Dt, Ct, Bt, Et, Kya_prime] = calculateMzo(TIRparam, useAlphaStar, useLimitsCheck, Fz, gamma, alpha_star, alpha_prime, gamma_star, Fzo_prime, dfz, dpi, Vcx, epsilonk, epsilony, Kya, Kygo, SHy, SVy, By, Cy, LMUY_star, LMUY_prime, zetaStr)%#codegen
% Declare extrinsic functions
coder.extrinsic('warning')

% Unpack zeta structure
zeta0= zetaStr.zeta0;
zeta2= zetaStr.zeta2;
zeta5= zetaStr.zeta5;
zeta6= zetaStr.zeta6;
zeta7= zetaStr.zeta7;
zeta8= zetaStr.zeta8;

%[DIMENSION]
Ro                      = TIRparam.UNLOADED_RADIUS      ; %Free tyre radius

%[SCALING_COEFFICIENTS]
LKY                     = TIRparam.LKY             		; % Scale factor of Fy cornering stiffness
LTR                     = TIRparam.LTR             		; % Scale factor of peak of pneumatic trail
LRES                    = TIRparam.LRES             	; % Scale factor for offset of residual torque
LKZC                    = TIRparam.LKZC             	; % %Scale factor of camber torque stiffness

%[ALIGNING_COEFFICIENTS]
QBZ1                    =  TIRparam.QBZ1    			; %Trail slope factor for trail Bpt at Fznom
QBZ2                    =  TIRparam.QBZ2     			; %Variation of slope Bpt with load
QBZ3                    =  TIRparam.QBZ3     			; %Variation of slope Bpt with load squared
QBZ4                    =  TIRparam.QBZ4				; %Variation of slope Bpt with camber
QBZ5                    =  TIRparam.QBZ5     			; %Variation of slope Bpt with absolute camber
QBZ9                    =  TIRparam.QBZ9     			; %Factor for scaling factors of slope factor Br of Mzr
QBZ10                   =  TIRparam.QBZ10    			; %Factor for dimensionless cornering stiffness of Br of Mzr
QCZ1                    =  TIRparam.QCZ1   				; %Shape factor Cpt for pneumatic trail
QDZ1                    =  TIRparam.QDZ1    			; %Peak trail Dpt = Dpt.*(Fz./Fznom.*R0)
QDZ2                    =  TIRparam.QDZ2   				; %Variation of peak Dpt" with load
QDZ3                    =  TIRparam.QDZ3   				; %Variation of peak Dpt" with camber
QDZ4                    =  TIRparam.QDZ4     			; %Variation of peak Dpt" with camber squared
QDZ6                    =  TIRparam.QDZ6  				; %Peak residual torque Dmr" = Dmr./(Fz.*R0)
QDZ7                    =  TIRparam.QDZ7  				; %Variation of peak factor Dmr" with load
QDZ8                    =  TIRparam.QDZ8   				; %Variation of peak factor Dmr" with camber
QDZ9                    =  TIRparam.QDZ9   				; %Variation of peak factor Dmr" with camber and load
QDZ10                   =  TIRparam.QDZ10     			; %Variation of peak factor Dmr with camber squared
QDZ11                   =  TIRparam.QDZ11     			; %Variation of Dmr with camber squared and load
QEZ1                    =  TIRparam.QEZ1    			; %Trail curvature Ept at Fznom
QEZ2                    =  TIRparam.QEZ2     			; %Variation of curvature Ept with load
QEZ3                    =  TIRparam.QEZ3    			; %Variation of curvature Ept with load squared
QEZ4                    =  TIRparam.QEZ4   				; %Variation of curvature Ept with sign of Alpha-t
QEZ5                    =  TIRparam.QEZ5    			; %Variation of Ept with camber and sign Alpha-t
QHZ1                    =  TIRparam.QHZ1  				; %Trail horizontal shift Sht at Fznom
QHZ2                    =  TIRparam.QHZ2 				; %Variation of shift Sht with load
QHZ3                    =  TIRparam.QHZ3    			; %Variation of shift Sht with camber
QHZ4                    =  TIRparam.QHZ4   				; %Variation of shift Sht with camber and load
PPZ1                    =  TIRparam.PPZ1				; %effect of inflation pressure on length of pneumatic trail
PPZ2                    =  TIRparam.PPZ2				; %Influence of inflation pressure on residual aligning torque

SHt = QHZ1 + QHZ2.*dfz + (QHZ3 + QHZ4.*dfz).*gamma_star; % (4.E35)

signKya = sign(Kya);
signKya(signKya == 0) = 1; % If [Kya = 0] then [sign(0) = 0]. This is done to avoid [num / 0 = NaN] in Eqn 4.E38

Kya_prime = Kya + epsilonk.*signKya; % (4.E39) [sign(Kya) term explained on page 177]
SHf = SHy + SVy./Kya_prime; % (4.E38)
alphar = alpha_star + SHf; % = alphaf (4.E37)
alphat = alpha_star + SHt; % (4.E34)
Dto = Fz.*(Ro./Fzo_prime).*(QDZ1 + QDZ2.*dfz).*(1 - PPZ1.*dpi).* LTR.*sign(Vcx); % (4.E42)
% Dt = Dto.*(1 + QDZ3.*abs(gamma_star) + QDZ4.*gamma_star.^2).*zeta5; % (4.E43)
%
% IMPORTANT NOTE: The above original equation (4.E43) was not matching the
% TNO solver. The coefficient Dt affects the pneumatic trail (t) and the
% self aligning torque (Mz).
% It was observed that when negative inclination angles where used as
% inputs, there was a discrepancy between the TNO solver and mfeval.
% This difference comes from the term QDZ3, that in the original equation
% is multiplied by abs(gamma_star). But in the paper the equation is
% different and the abs() term is not written. Equation (A60) from the
% paper resulted into a perfect match with TNO.
% Keep in mind that the equations from the paper don't include turn slip
% effects. The term zeta5 has been added although it doesn't appear in the
% paper.

% Paper definition:
Dt = (QDZ1 + QDZ2.*dfz).*(1 - PPZ1.*dpi).* (1 + QDZ3.*gamma + QDZ4.*gamma.^2).*Fz.*(Ro./Fzo_prime).*LTR.*zeta5; % (A60)

% Bt = (QBZ1 + QBZ2.*dfz + QBZ3.*dfz.^2).*(1 + QBZ5.*abs(gamma_star) + QBZ6.*gamma_star.^2).*LKY./LMUY_star; %(> 0)(4.E40)
%
% IMPORTANT NOTE: In the above original equation (4.E40) it is used the
% parameter QBZ6, which doesn't exist in the standard TIR files. Also note
% that on page 190 and 615 of the book a full set of parameters is given
% and QBZ6 doesn't appear.
% The equation has been replaced with equation (A58) from the paper.

% Paper definition:
Bt = (QBZ1 + QBZ2.*dfz + QBZ3.*dfz.^2).*(1 + QBZ4.*gamma + QBZ5.*abs(gamma) ).*LKY./LMUY_star; %(> 0) (A58)
Ct = QCZ1; % (> 0) (4.E41)
Et = (QEZ1 + QEZ2.*dfz + QEZ3.*dfz.^2).*(1+(QEZ4+QEZ5.*gamma_star).*(2./pi).*atan(Bt.*Ct.*alphat)); % (<=1) (4.E44)

% Limits check
if(useLimitsCheck)
    if(any(Et > 1))
        warning('Et over limit (>1), Eqn(4.E44)')
        Et(Et > 1) = 1;
    end
end

to = Dt.*cos(Ct.*atan(Bt.*alphat - Et.*(Bt.*alphat - atan(Bt.*alphat)))).*cos(alpha_prime); %t(aplhat)(4.E33)

% Create a structure of zeta with phit = 0
zetaStr_sub0.zeta0 = 1;
zetaStr_sub0.zeta1 = 1;
zetaStr_sub0.zeta2 = 1;
zetaStr_sub0.zeta3 = 1;
zetaStr_sub0.zeta4 = 1;
zetaStr_sub0.zeta5 = 1;
zetaStr_sub0.zeta6 = 1;
zetaStr_sub0.zeta7 = 1;
zetaStr_sub0.zeta8 = 1;

% Evaluate Fyo with gamma = 0 and phit = 0
[Fyo_sub0,~,~,~,~,~,~,~] = calculateFyo(TIRparam, useAlphaStar, useLimitsCheck, Fz, alpha_star, 0, Fzo_prime, dfz, dpi, Vcx, epsilony, epsilonk, LMUY_star, LMUY_prime, zetaStr_sub0);
Mzo_prime = -to.*Fyo_sub0; % gamma=phi=0 (4.E32)
Dr = Fz.*Ro.*((QDZ6 + QDZ7.*dfz).*LRES.*zeta2 + ((QDZ8 + QDZ9.*dfz).*(1 + PPZ2.*dpi) + (QDZ10 + QDZ11.*dfz).*abs(gamma_star)).*gamma_star.*LKZC.*zeta0).*LMUY_star.*sign(Vcx).*cos(alpha_star) + zeta8 - 1; % (4.E47)
Br = (QBZ9.*LKY./LMUY_star + QBZ10.*By.*Cy).*zeta6; % preferred: qBz9 = 0 (4.E45)
Cr = zeta7; % (4.E46)

% Evaluate Kya with gamma = 0
[~,~,Kya_sub0,~,~,~,~,~] = calculateFyo(TIRparam, useAlphaStar, useLimitsCheck, Fz, alpha_star, 0, Fzo_prime, dfz, dpi, Vcx, epsilony, epsilonk, LMUY_star, LMUY_prime, zetaStr);

Kzao = Dto.*Kya_sub0; %#ok<NASGU> % Kya for gamma=0 (= -dMzo./dAlphay at alphay=gamma=0)(=CMa) (4.E48)
Kzgo = Fz.*Ro.*(QDZ8 + QDZ9.*dfz).*(1 + PPZ2.*dpi).*LKZC.*LMUY_star - Dto.*Kygo; %#ok<NASGU> %(=dMzo./dgamma at alpha = gamma = 0) (= CMg)(4.E49)
Mzro = Dr.*cos(Cr.*atan(Br.*alphar)).*cos(alpha_prime); % =Mzr(alphar)(4.E36)
Mzo = Mzo_prime + Mzro; % (4.E31)
end

function [Fx] = calculateFx(TIRparam, useLimitsCheck, kappa, alpha_star, gamma_star, dfz, Fxo)%#codegen
% Declare extrinsic functions
coder.extrinsic('warning')

%[SCALING_COEFFICIENTS]
LXAL                    = TIRparam.LXAL             	; % Scale factor of alpha influence on Fx

%[LONGITUDINAL_COEFFICIENTS]
RBX1                    =  TIRparam.RBX1      			; %Slope factor for combined slip Fx reduction
RBX2                    =  TIRparam.RBX2      			; %Variation of slope Fx reduction with kappa
RBX3                    =  TIRparam.RBX3     			; %Influence of camber on stiffness for Fx combined
RCX1                    =  TIRparam.RCX1    			; %Shape factor for combined slip Fx reduction
REX1                    =  TIRparam.REX1     			; %Curvature factor of combined Fx
REX2                    =  TIRparam.REX2   				; %Curvature factor of combined Fx with load
RHX1                    =  TIRparam.RHX1  				; %Shift factor for combined slip Fx reduction

Cxa = RCX1; % (4.E55)
Exa = REX1 + REX2.*dfz; % (<= 1) (4.E56)

% Limits check
if(useLimitsCheck)
    if(any(Exa > 1))
        warning('Exa over limit (>1), Eqn(4.E56)')
        Exa(Exa > 1) = 1;
    end
end

SHxa = RHX1; % (4.E57)
Bxa = (RBX1 + RBX3.*gamma_star.^2).*cos(atan(RBX2.*kappa)).*LXAL; % (> 0) (4.E54)
alphas = alpha_star + SHxa; % (4.E53)
Gxao = cos(Cxa.*atan(Bxa.*SHxa-Exa.*(Bxa.*SHxa-atan(Bxa.*SHxa)))); % (4.E52)
Gxa = cos(Cxa.*atan(Bxa.*alphas-Exa.*(Bxa.*alphas-atan(Bxa.*alphas))))./Gxao;  % (> 0)(4.E51)
Fx = Gxa.*Fxo; % (4.E50)

end

function [Fy,Gyk] = calculateFy(TIRparam, useLimitsCheck, Fz, alpha_star, kappa, gamma_star, dfz, Fyo, muy, zetaStr, isLowSpeed, reduction_factor)%#codegen
% Declare extrinsic functions
coder.extrinsic('warning')

% Unpack zeta structure
zeta2= zetaStr.zeta2;

%[SCALING_COEFFICIENTS]
LYKA                    = TIRparam.LYKA             	; % Scale factor of alpha influence on Fx
LVYKA                   = TIRparam.LVYKA             	; % Scale factor of kappa induced Fy

%[LATERAL_COEFFICIENTS]
RBY1                    =  TIRparam.RBY1     			; %Slope factor for combined Fy reduction
RBY2                    =  TIRparam.RBY2    			; %Variation of slope Fy reduction with alpha
RBY3                    =  TIRparam.RBY3  				; %Shift term for alpha in slope Fy reduction
RBY4                    =  TIRparam.RBY4     			; %Influence of camber on stiffness of Fy combined
RCY1                    =  TIRparam.RCY1    			; %Shape factor for combined Fy reduction
REY1                    =  TIRparam.REY1    			; %Curvature factor of combined Fy
REY2                    =  TIRparam.REY2    			; %Curvature factor of combined Fy with load
RHY1                    =  TIRparam.RHY1  				; %Shift factor for combined Fy reduction
RHY2                    =  TIRparam.RHY2   				; %Shift factor for combined Fy reduction with load
RVY1                    =  TIRparam.RVY1  				; %Kappa induced side force Svyk./Muy.*Fz at Fznom
RVY2                    =  TIRparam.RVY2  				; %Variation of Svyk./Muy.*Fz with load
RVY3                    =  TIRparam.RVY3    			; %Variation of Svyk./Muy.*Fz with camber
RVY4                    =  TIRparam.RVY4   				; %Variation of Svyk./Muy.*Fz with alpha
RVY5                    =  TIRparam.RVY5     			; %Variation of Svyk./Muy.*Fz with kappa
RVY6                    =  TIRparam.RVY6   				; %Variation of Svyk./Muy.*Fz with atan(kappa)

DVyk = muy.*Fz.*(RVY1 + RVY2.*dfz + RVY3.*gamma_star).*cos(atan(RVY4.*alpha_star)).*zeta2; % (4.E67)
SVyk = DVyk.*sin(RVY5.*atan(RVY6.*kappa)).*LVYKA; % (4.E66)
SHyk = RHY1 + RHY2.*dfz; % (4.E65)
Eyk = REY1 + REY2.*dfz; % (<=1) (4.E64)

% Limits check
if(useLimitsCheck)
    if(any(Eyk > 1))
        warning('Eyk over limit (>1), Eqn(4.E64)')
        Eyk(Eyk > 1) = 1;
    end
end

Cyk = RCY1; % (4.E63)
Byk = (RBY1 + RBY4.*gamma_star.^2).*cos(atan(RBY2.*(alpha_star - RBY3))).*LYKA; % (> 0) (4.E62)
kappas = kappa + SHyk; % (4.E61)
Gyko = cos(Cyk.*atan(Byk.*SHyk - Eyk.*(Byk.*SHyk - atan(Byk.*SHyk)))); % (4.E60)
Gyk = cos(Cyk.*atan(Byk.*kappas - Eyk.*(Byk.*kappas-atan(Byk.*kappas))))./Gyko; % (> 0)(4.E59)
Fy = Gyk.*Fyo + SVyk; % (4.E58)

% Low speed model
if any(isLowSpeed)
    Fy(isLowSpeed) = Fy(isLowSpeed).* reduction_factor ;
end

end

function [Mx] = calculateMx(TIRparam, Fz, Fy, gamma, dpi)%#codegen

%[VERTICAL]
Fzo                     = TIRparam.FNOMIN         		; %Nominal wheel load

%[DIMENSION]
Ro                      = TIRparam.UNLOADED_RADIUS      ; %Free tyre radius

%[SCALING_COEFFICIENTS]
LVMX                    = TIRparam.LVMX             	; % Scale factor of Mx vertical shift
LMX                     = TIRparam.LMX             		; % Scale factor of overturning couple

%[OVERTURNING_COEFFICIENTS]
QSX1                    =  TIRparam.QSX1  				; %Vertical shift of overturning moment
QSX2                    =  TIRparam.QSX2    			; %Camber induced overturning couple
QSX3                    =  TIRparam.QSX3   				; %Fy induced overturning couple
QSX4                    =  TIRparam.QSX4				; %Mixed load lateral force and camber on Mx
QSX5                    =  TIRparam.QSX5	    		; %Load effect on Mx with lateral force and camber
QSX6                    =  TIRparam.QSX6	    		; %B-factor of load with Mx
QSX7                    =  TIRparam.QSX7	    		; %Camber with load on Mx
QSX8                    =  TIRparam.QSX8				; %Lateral force with load on Mx
QSX9                    =  TIRparam.QSX9				; %B-factor of lateral force with load on Mx
QSX10                   =  TIRparam.QSX10				; %Vertical force with camber on Mx
QSX11                   =  TIRparam.QSX11	    		; %B-factor of vertical force with camber on Mx
QSX12                   =  TIRparam.QSX12	  			; %Camber squared induced overturning moment
QSX13                   =  TIRparam.QSX13	  			; %Lateral force induced overturning moment
QSX14                   =  TIRparam.QSX14	  			; %Lateral force induced overturning moment with camber
PPMX1                   =  TIRparam.PPMX1				; %Influence of inflation pressure on overturning moment

% Mx = Ro.*Fz.*(QSX1.*LVMX - QSX2.*gamma.*(1 + PPMX1.*dpi) + QSX3.*((Fy)/Fzo)...
%     + QSX4.*cos(QSX5.*atan((QSX6.*(Fz./Fzo)).^2)).*sin(QSX7.*gamma + QSX8.*atan...
%     (QSX9.*((Fy)/Fzo))) + QSX10.*atan(QSX11.*(Fz./Fzo)).*gamma).*LMX; %(4.E69)
%
% IMPORTANT NOTE: The above original equation (4.E69) is not used
% because is not matching the results of the official TNO solver. Also, in
% the book definition and in the official paper, parameters QSX12 QSX13 and
% QSX14 are not used
% Instead of using equations (4.E69) from the book or (A47) from the
% official paper, it has been used the equation (49) described in the draft
% paper of Besselink (Not the official paper). This draft can be downloaded
% from:
% https://pure.tue.nl/ws/files/3139488/677330157969510.pdf
% purl.tue.nl/677330157969510.pdf

% Draft paper definition:
Mx = Ro .* Fz .* LMX .* (QSX1 .* LVMX - QSX2 .* gamma .* (1 + PPMX1 .* dpi) - QSX12 .* gamma .* abs(gamma) + QSX3 .* Fy ./ Fzo + ...
    QSX4 .* cos(QSX5 .* atan((QSX6 .* Fz ./ Fzo) .^ 2)) .* sin(QSX7 .* gamma + QSX8 .* atan(QSX9 .* Fy ./ Fzo)) + ...
    QSX10 .* atan(QSX11 .* Fz ./ Fzo) .* gamma) + Ro .* Fy .* LMX .* (QSX13 + QSX14 .* abs(gamma)); % (49)
end

function [My] = calculateMy(TIRparam, Fz_unlimited, Fx, Vcx, kappa, gamma, p, isLowSpeed, reduction_factor)%#codegen

%[OPERATING_CONDITIONS]
pio                     = TIRparam.NOMPRES        		; %Nominal tyre inflation pressure

%[MODEL]
Vo                      = TIRparam.LONGVL               ; %Nominal speed

%[VERTICAL]
Fzo                     = TIRparam.FNOMIN         		; %Nominal wheel load

%[DIMENSION]
Ro                      = TIRparam.UNLOADED_RADIUS      ; %Free tyre radius

%[SCALING_COEFFICIENTS]
LMY                     = TIRparam.LMY             		; % Scale factor of rolling resistance torque

%[ROLLING_COEFFICIENTS]
QSY1                    =  TIRparam.QSY1         		; %Rolling resistance torque coefficient
QSY2                    =  TIRparam.QSY2           		; %Rolling resistance torque depending on Fx
QSY3                    =  TIRparam.QSY3       			; %Rolling resistance torque depending on speed
QSY4                    =  TIRparam.QSY4       			; %Rolling resistance torque depending on speed .^4
QSY5                    =  TIRparam.QSY5            	; %Rolling resistance torque depending on camber squared
QSY6                    =  TIRparam.QSY6            	; %Rolling resistance torque depending on load and camber squared
QSY7                    =  TIRparam.QSY7         		; %Rolling resistance torque coefficient load dependency
QSY8                    =  TIRparam.QSY8         		; %Rolling resistance torque coefficient pressure dependency

% My = Fz.Ro*(QSY1 + QSY2.*(Fx./Fzo) + QSY3.*abs(Vcx./Vo) + QSY4.*(Vcx./Vo).^4 ...
%     +(QSY5 + QSY6.*(Fz./Fzo)).*gamma.^2).*((Fz./Fzo).^QSY7.*(p./pio).^QSY8).*LMY.; %(4.E70)
%
% IMPORTANT NOTE: The above equation from the book (4.E70) is not used
% because is not matching the results of the official TNO solver.
% This equation gives a positive output of rolling resistance, and in the
% ISO coordinate system, My should be negative. Furthermore, the equation
% from the book has an error, multiplying all the equation by Fz instead of
% Fzo (first term).
% Because of the previous issues it has been used the equation (A48) from
% the paper.

% Paper definition:
My = -Ro.*Fzo.*LMY.*(QSY1 + QSY2.*(Fx./Fzo) + QSY3.*abs(Vcx./Vo) + QSY4.*(Vcx./Vo).^4 ...
    +(QSY5 + QSY6.*(Fz_unlimited./Fzo)).*gamma.^2).*((Fz_unlimited./Fzo).^QSY7.*(p./pio).^QSY8); %(A48)

% Backward speed check
My(Vcx < 0) = -My(Vcx < 0);

% Negative SR check
My(kappa<-1) = -My(kappa<-1);

% Low speed model
if any(isLowSpeed)
    My(isLowSpeed) = My(isLowSpeed).*reduction_factor;
end

end

function [Mz,t] = calculateMz(TIRparam, useAlphaStar, useLimitsCheck, alphar, alphat, Kxk, Kya_prime, kappa, Fz, Fy, Fx, Fzo_prime, dfz, gamma, Dr, Cr, Br, Dt, Ct, Bt, Et, alpha_star, Fyo, muy, alpha_prime, Vcx, dpi, epsilony, epsilonk, LMUY_star, LMUY_prime, isLowSpeed, reduction_factor)%#codegen

%[VERTICAL]
Fzo                     = TIRparam.FNOMIN         		; %Nominal wheel load

%[DIMENSION]
Ro                      = TIRparam.UNLOADED_RADIUS      ; %Free tyre radius

%[SCALING_COEFFICIENTS]
LS                      = TIRparam.LS             		; % Scale factor of moment arm of Fx
LFZO                    = TIRparam.LFZO             	; % Scale factor of nominal (rated) load

%[ALIGNING_COEFFICIENTS]
SSZ1                    =  TIRparam.SSZ1   				; %Nominal value of s./R0: effect of Fx on Mz
SSZ2                    =  TIRparam.SSZ2     			; %Variation of distance s./R0 with Fy./Fznom
SSZ3                    =  TIRparam.SSZ3    			; %Variation of distance s./R0 with camber
SSZ4                    =  TIRparam.SSZ4   				; %Variation of distance s./R0 with load and camber

% alphar_eq = sqrt(alphar.^2+(Kxk./Kya_prime).^2.*kappa.^2).*sign(alphar); % (4.E78)
% alphat_eq = sqrt(alphat.^2+(Kxk./Kya_prime).^2.*kappa.^2).*sign(alphat); % (4.E77)
% s = Ro.*(SSZ1 + SSZ2.*(Fy./Fzo_prime) + (SSZ3 + SSZ4.*dfz).*gamma_star).*LS; % (4.E76)

% IMPORTANT NOTE: The equations 4.E78 and 4.E77 are not used due to small
% differences discovered at negative camber angles with the TNO solver.
% Instead equations A54 and A55 from the paper are used.
%
% IMPORTANT NOTE: The coefficient "s" (Equation 4.E76) determines the
% effect of Fx into Mz. The book uses "Fzo_prime" in the formulation,
% but the paper uses "Fzo". The equation (A56) from the paper has a better
% correlation with TNO.

alphar_eq = atan(sqrt(tan(alphar).^2+(Kxk./Kya_prime).^2.*kappa.^2)).*sign(alphar); % (A54)
alphat_eq = atan(sqrt(tan(alphat).^2+(Kxk./Kya_prime).^2.*kappa.^2)).*sign(alphat); % (A55)
s = Ro.*(SSZ1 + SSZ2.*(Fy./Fzo) + (SSZ3 + SSZ4.*dfz).*gamma).*LS; % (A56)
Mzr = Dr.*cos(Cr.*atan(Br.*alphar_eq)); % (4.E75)

% Create a structure of zeta with phit = 0
zetaStr_sub0.zeta0 = 1;
zetaStr_sub0.zeta1 = 1;
zetaStr_sub0.zeta2 = 1;
zetaStr_sub0.zeta3 = 1;
zetaStr_sub0.zeta4 = 1;
zetaStr_sub0.zeta5 = 1;
zetaStr_sub0.zeta6 = 1;
zetaStr_sub0.zeta7 = 1;
zetaStr_sub0.zeta8 = 1;

% Evaluate Gyk with gamma = 0 and phit  = 0
[~,Gyk_sub0] = calculateFy(TIRparam,useLimitsCheck,Fz,alpha_star,kappa,0,dfz,Fyo,muy,zetaStr_sub0, isLowSpeed, reduction_factor);

% Evaluate Fyo with gamma = 0 and phit  = 0
[Fyo_sub0,~,~,~,~,~,~,~] = calculateFyo(TIRparam, useAlphaStar, useLimitsCheck, Fz, alpha_star, 0, Fzo_prime, dfz, dpi, Vcx, epsilony, epsilonk, LMUY_star, LMUY_prime, zetaStr_sub0);

Fy_prime = Gyk_sub0.*Fyo_sub0; % (4.E74)
t = Dt.*cos(Ct.*atan(Bt.*alphat_eq - Et.*(Bt.*alphat_eq - atan(Bt.*alphat_eq)))).*cos(alpha_prime); % (4.E73)

t = t*LFZO;

% IMPORTANT NOTE: the above equation is not written in any source, but "t"
% is multiplied by LFZO in the TNO solver. This has been empirically
% discovered.

% Low speed model
if any(isLowSpeed)
    t(isLowSpeed) = t(isLowSpeed).*reduction_factor;
    Mzr(isLowSpeed) = Mzr(isLowSpeed).*reduction_factor;
end

Mz_prime = -t.*Fy_prime; % (4.E72)
Mz = Mz_prime + Mzr + s.*Fx; % (4.E71)
end

function [rho, Rl] = calculateRhoRl(TIRparam, omega, Fz_unlimited, Fx, Fy, dpi, Romega)%#codegen
% Declare extrinsic functions
coder.extrinsic('warning')

%[MODEL]
Vo                      = TIRparam.LONGVL               ; %Nominal speed

%[VERTICAL]
Fzo                     = TIRparam.FNOMIN         		; %Nominal wheel load
Cz0                     = TIRparam.VERTICAL_STIFFNESS	; %Tyre vertical stiffness
Q_V2                    = TIRparam.Q_V2    				; %Vertical stiffness increase with speed
Q_FZ2                   = TIRparam.Q_FZ2      			; %Quadratic term in load vs. deflection
Q_FCX                   = TIRparam.Q_FCX          		; %Longitudinal force influence on vertical stiffness
Q_FCY                   = TIRparam.Q_FCY           		; %Lateral force influence on vertical stiffness
PFZ1                    = TIRparam.PFZ1    				; %Pressure effect on vertical stiffness

%[DIMENSION]
Ro                      = TIRparam.UNLOADED_RADIUS      ; %Free tyre radius

% Model parameters as QFZ1 that normally aren't present in the TIR files
Q_FZ1 = sqrt((Cz0.*Ro./Fzo).^2 - 4.* Q_FZ2); % Rearranging [Eqn (4) Page 2 - Paper]

% Split Eqn (A3.3) Page 619 of the Book into different bits:
speed_effect    = Q_V2 .* (Ro ./ Vo) .*abs(omega);
fx_effect       = (Q_FCX .* Fx ./ Fzo).^2;
fy_effect       = (Q_FCY .* Fy ./ Fzo).^2;
pressure_effect   = (1+PFZ1.*dpi);

% Joining all the effects except tyre deflection terms:
external_effects = (1 + speed_effect - fx_effect - fy_effect) .* pressure_effect.*Fzo;

% Equation (A3.3) can be written as:
% Fz = (Q_FZ2*(rho/Ro)^2 + Q_FZ1*(rho/Ro)) * external_effects

% Rearranging all the terms we end up with a quadratic equation as:
% ax^2 + bx + c = 0
% Q_FZ2*(rho/Ro)^2 + Q_FZ1*(rho/Ro) -(Fz/(external_effects)) = 0

% Note: use of capital letters to avoid confusion with contact patch
% lengths "a" and "b"

A = Q_FZ2;
B = Q_FZ1;
C = -(Fz_unlimited ./ (external_effects));

if all((B^2 - 4*A*C) > 0)
    x = (-B + sqrt(B^2 - 4*A.*C)) ./ (2*A);
else
    warning('No positive solution for rho calculation')
    x = (-B + sqrt(B^2 - 4*A.*C)) ./ (2*A);
end

rho = x.*Ro;

% The loaded radius is the free-spinning radius minus the deflection
Rl = Romega - rho; % Eqn A3.2 Page 619 - Book assuming positive rho at all the time
end

function [a, b] = calculateContactPatch(TIRparam, Fz_unlimited, dpi)%#codegen

% Rename the TIR file variables in the Pacejka style
Ro  = TIRparam.UNLOADED_RADIUS;         % Free tyre radius
w   = TIRparam.WIDTH;                   % Nominal width of the tyre
Cz0 = TIRparam.VERTICAL_STIFFNESS;      % Vertical stiffness

PFZ1                    = TIRparam.PFZ1                 ; %Pressure effect on vertical stiffness

% Nominal stiffness (pressure corrected)
Cz = Cz0 .* (1 + PFZ1.*dpi); % [Eqn (5) Page 2 - Paper] - Vertical stiffness adapted for tyre inflation pressure


%[CONTACT_PATCH]
Q_RA1                    = TIRparam.Q_RA1    			; %Square root term in contact length equation
Q_RA2                    = TIRparam.Q_RA2     			; %Linear term in contact length equation
Q_RB1                    = TIRparam.Q_RB1     			; %Root term in contact width equation
Q_RB2                    = TIRparam.Q_RB2    			; %Linear term in contact width equation

% (12) Contact patch length
a = Ro.*(Q_RA2.*(Fz_unlimited./(Cz.*Ro)) + Q_RA1.*sqrt(Fz_unlimited./(Cz.*Ro))); %[Eqn (9) Page 3 - Paper]
b = w.*(Q_RB2.*(Fz_unlimited./(Cz.*Ro)) + Q_RB1.*(Fz_unlimited./(Cz.*Ro)).^(1/3)); %[Eqn (10) Page 3 - Paper] (Not used later in this code)
end

function [Re, Romega, omega] = calculateRe(TIRparam, ncolumns, INPUTS, dpi, Fz_unlimited, kappa_unlimited)%#codegen

% Excerpt from OpenTIRE MF6.1 implementation
% Date: 2016-12-01
% Prepared for Marco Furlan/JLR
% Questions: henning.olsson@calspan.com

% Forward velocity (m/s)
Vcx     = INPUTS(:,6);

% Rename the TIR file variables in the Pacejka style
Fzo     = TIRparam.FNOMIN;               % Nominal (rated) wheel load
Ro      = TIRparam.UNLOADED_RADIUS;       % Free tyre radius
Vo      = TIRparam.LONGVL;                % Nominal speed (LONGVL)
Cz0     = TIRparam.VERTICAL_STIFFNESS;   % Vertical stiffness

PFZ1                    = TIRparam.PFZ1                 ; %Pressure effect on vertical stiffness
BREFF                   = TIRparam.BREFF      			; %Low load stiffness effective rolling radius
DREFF                   = TIRparam.DREFF     			; %Peak value of effective rolling radius
FREFF                   = TIRparam.FREFF    			; %High load stiffness effective rolling radius
Q_RE0                   = TIRparam.Q_RE0       			; %Ratio of free tyre radius with nominal tyre radius
Q_V1                    = TIRparam.Q_V1  				; %Tyre radius increase with speed

% Nominal stiffness (pressure corrected)
Cz = Cz0 .* (1 + PFZ1.*dpi); % [Eqn (5) Page 2 - Paper] - Vertical stiffness adapted for tyre inflation pressure

% Check if omega is one of the inputs
% If it is, use it to calculate Re, otherwise it can be approximated with a
% short iteration

if ncolumns > 7
    % Omega is one of the inputs
    omega = INPUTS(:,8); % rotational speed         (rad/s)

    Romega = Ro .* (Q_RE0 + Q_V1.*((omega.*Ro)./Vo).^2); % [Eqn (1) Page 2 - Paper] - Centrifugal growth of the free tyre radius

    % Eff. Roll. Radius
    Re = Romega -(Fzo./Cz) .* ( DREFF .*atan( BREFF.*(Fz_unlimited./Fzo)) + FREFF.*(Fz_unlimited./Fzo)); % [Eqn (7) Page 2 - Paper]
else
    % Omega is not specified and is going to be approximated
    % Initial guess of Re based on something slightly less than R0
    Re = (Ro .* 0.965);

    for i=1:10
        % Use the most up to date Re to calculate an omega
        % omega = Vcx ./ Re; % Old definition of Henning without kappa, not valid for brake and drive
        omega = (kappa_unlimited.*Vcx+Vcx)./Re; % [Eqn (2.5) Page 65 - Book]

        % Then we calculate free-spinning radius
        Romega = Ro .* (Q_RE0 + Q_V1.*((omega.*Ro)./Vo).^2); % [Eqn (1) Page 2 - Paper] - Centrifugal growth of the free tyre radius

        % Effective Rolling Radius
        Re = Romega -(Fzo./Cz) .* ( DREFF .*atan( BREFF.*(Fz_unlimited./Fzo)) + FREFF.*(Fz_unlimited./Fzo)); % [Eqn (7) Page 2 - Paper]
    end
end

end