function [alpha] = calculateYld2000Parameters(sigma, r, a);

% Input:
% sigma = [s0, s45, s90, sb] - Flow stress ratios with respect to hardening
%   law, usually assumed that hardening follows behaviour along RD
% r = [r0, r45, r90, rb] - Lankford coefficients and strain ratio under
%   equibiaxial load
% a - exponent used in Yld2000-2D

hardStress = 1;

% Initial guess for parameter set
initGuess = ones(1,8);

% Least-square matching using Levenberg-Marquardt algorithm
alpha = lsqnonlin(@(x)costFunction(x, sigma, r, a, hardStress), initGuess, [], [], ...
    optimoptions(@lsqnonlin, 'Algorithm', 'Levenberg-Marquardt', 'MaxIter', 100, 'TypicalX', initGuess, 'Disp', 'Iter'));


end


%% Supplementary functions
% Cost function for calculating alphas
function [res] = costFunction(alpha, sigma, r, a, hardStress)

% Preallocate output
res = zeros(1,8);

%% 0 deg uniaxial
ang = 0; %Angle between the load and the material csys
sig = sigma(1); % Fetch experimental flow ratio

% Stress state in the experiment (in material csys)
stress_in = [cosd(ang)^2*sig; sind(ang)^2*sig; cosd(ang)*sind(ang)*sig];

% Predicted values (with model)
    % Equivalent stress, comparable to hardening
eqStress = calcEqStressYld2000(stress_in, alpha, a);
    % Lankford coefficient 
rModelled = calcLankYld2000(ang, stress_in, alpha, a, 1);

% Compare the predicted values with experimental measurements
res(1) = eqStress-hardStress;
res(2) = rModelled-r(1);

%% 45 deg uniaxial
ang = 45; %Angle between the load and material csys
sig = sigma(2); % Fetch experimental flow ratio

% Stress state in the experiment (in material csys)
stress_in = [cosd(ang)^2*sig; sind(ang)^2*sig; cosd(ang)*sind(ang)*sig];

% Predicted values (with model)
    % Equivalent stress, comparable to hardening
eqStress = calcEqStressYld2000(stress_in, alpha, a);
    % Lankford coefficient 
rModelled = calcLankYld2000(ang, stress_in, alpha, a, 1);

% Compare the predicted values with experimental measurements
res(3) = eqStress-hardStress;
res(4) = rModelled-r(2);

%% 90 deg uniaxial
ang = 90; %Angle between the load and material csys
sig = sigma(3);

% Stress state in the experiment (in material csys)
stress_in = [cosd(ang)^2*sig; sind(ang)^2*sig; cosd(ang)*sind(ang)*sig];

% Predicted values (with model)
    % Equivalent stress, comparable to hardening
eqStress = calcEqStressYld2000(stress_in, alpha, a);
    % Lankford coefficient 
rModelled = calcLankYld2000(ang, stress_in, alpha, a, 1);

% Compare the predicted values with experimental measurements
res(5) = eqStress-hardStress;
res(6) = rModelled-r(3);

%% Biaxial test
ang = 0;
sig = sigma(4);
stress_in = [sig; sig; 0];

% Predicted values (with model)
    % Equivalent stress, comparable to hardening
eqStress = calcEqStressYld2000(stress_in, alpha, a);
    % Lankford coefficient 
rModelled = calcLankYld2000(ang, stress_in, alpha, a, 2);

% Compare the predicted values with experimental measurements
res(7) = eqStress-hardStress;
res(8) = rModelled-r(4);

end

function s_eq = calcEqStressYld2000(stress, alpha, a);
% This function calculates equivalent stress using Yld2000-2d model given
% stress state and model parameters. The annotations refer to equation
% numbers in [1]
% 
% [1]F. Barlat, J. C. Brem, J. W. Yoon, K. Chung, R. E. Dick, D. J.
%   Lege, F. Pourboghrat, S.-H. Choi, and E. Chu. Plane stress yield
%   function for aluminum alloy sheets—part 1: theory. International
%   Journal of Plasticity, 19(9):1297–1319, 2003.


% Eq (16)
L_1 = [2*alpha(1)/3, -alpha(1)/3 0;
    -alpha(2)/3, 2*alpha(2)/3, 0;
    0 0 alpha(7)];

% Eq (16)
L_2 = [(8*alpha(5)-2*alpha(3)-2*alpha(6)+2*alpha(4))/9,...
    (4*alpha(6)-4*alpha(4)-4*alpha(5)+alpha(3))/9, 0;
    (4*alpha(3)-4*alpha(5)-4*alpha(4)+alpha(6))/9,...
    (8*alpha(4)-2*alpha(6)-2*alpha(3)+2*alpha(5))/9, 0;
    0 0 alpha(8)];

% Eq(14)
X1 = L_1*stress;
% Eq(14)
X2 = L_2*stress;

% Eq(17)
eta1 = [0.5*(X1(1)+X1(2) + sqrt((X1(1)-X1(2)).^2+4.0*X1(3)^2));
    0.5*(X1(1)+X1(2) - sqrt((X1(1)-X1(2)).^2+4.0*X1(3)^2))];
% Eq(17)
eta2 = [0.5*(X2(1)+X2(2) + sqrt((X2(1)-X2(2)).^2+4.0*X2(3)^2));
    0.5*(X2(1)+X2(2) - sqrt((X2(1)-X2(2)).^2+4.0*X2(3)^2))];

% Eq(10)+Eq(11)
s_eq = (0.5*(abs(eta1(1)-eta1(2)).^a+...
    abs(2*eta2(2)+eta2(1)).^a+...
    abs(2*eta2(1)+eta2(2)).^a))^(1/a);

end

function rCalc = calcLankYld2000(ang, stress, alpha, a, mode)
% This function calculates plastic flow ratios (Lankford coefficients and
% rb) based on desired stress state and model parameters using Yld2000-2D


% Calculate plastic flow rates from associated flow rule
[e_norm] = calcFlowRateYld2000(stress, alpha, a);
% Define rotation matrix (to rotate strains from the material to global
% csys)
Rmat = [cosd(ang) -sind(ang); sind(ang) cosd(ang)];
% Write flow as tensorial strain
e_norm = [e_norm(1) 0.5*e_norm(3); 0.5*e_norm(3) e_norm(2)];
% Rotate the flow rates to the global csys
e_norm = Rmat'*e_norm*Rmat;
% Revert back to engineering components
e_norm = [e_norm(1,1); e_norm(2,2); 2*e_norm(1,2)];

% Calculate plastic flow ratios
switch mode
    case 1 % uniaxial
        rCalc = e_norm(2)/(-e_norm(2)-e_norm(1));
    case 2 % biaxial
        rCalc = e_norm(2)/e_norm(1);
end
end

function [e_norm] = calcFlowRateYld2000(stress, alpha, a);
% This function calculates plastic flow rates using Yld2000-2d model given
% stress state and model parameters. The annotations refer to equation
% numbers in [1]
% 
% [1]F. Barlat, J. C. Brem, J. W. Yoon, K. Chung, R. E. Dick, D. J.
%   Lege, F. Pourboghrat, S.-H. Choi, and E. Chu. Plane stress yield
%   function for aluminum alloy sheets—part 1: theory. International
%   Journal of Plasticity, 19(9):1297–1319, 2003.

% Eq (16)
L_1 = [2*alpha(1)/3, -alpha(1)/3 0;
    -alpha(2)/3, 2*alpha(2)/3, 0;
    0 0 alpha(7)];
% Eq (16)
L_2 = [(8*alpha(5)-2*alpha(3)-2*alpha(6)+2*alpha(4))/9,...
    (4*alpha(6)-4*alpha(4)-4*alpha(5)+alpha(3))/9, 0;
    (4*alpha(3)-4*alpha(5)-4*alpha(4)+alpha(6))/9,...
    (8*alpha(4)-2*alpha(6)-2*alpha(3)+2*alpha(5))/9, 0;
    0 0 alpha(8)];

% Eq(14)
X1 = L_1*stress;
% Eq(14)
X2 = L_2*stress;

% Eq(17)
eta1 = [0.5*((X1(1)+X1(2)) + sqrt((X1(1)-X1(2)).^2+4.0*X1(3)^2));
    0.5*((X1(1)+X1(2)) - sqrt((X1(1)-X1(2)).^2+4.0*X1(3)^2))];
% Eq(17)
eta2 = [0.5*((X2(1)+X2(2)) + sqrt((X2(1)-X2(2)).^2+4.0*X2(3)^2));
    0.5*((X2(1)+X2(2)) - sqrt((X2(1)-X2(2)).^2+4.0*X2(3)^2))];

% Eq(10)+Eq(11)
s_eq = (0.5*(abs(eta1(1)-eta1(2)).^a+...
    abs(2*eta2(2)+eta2(1)).^a+...
    abs(2*eta2(1)+eta2(2)).^a))^(1/a);

% First term in Eq(18)
dPhi_eta1 = zeros(2,1);

% Eq(A1.4)
dPhi_eta1(1) = a*abs(eta1(1)-eta1(2))^(a-1)*sign(eta1(1)-eta1(2));
dPhi_eta1(2) = -a*abs(eta1(1)-eta1(2))^(a-1)*sign(eta1(1)-eta1(2));

% Second term in Eq(18)
dPhi_eta2 = zeros(2,1);

% Eq(A1.8)
dPhi_eta2(1) = a*abs(2*eta2(2)+eta2(1))^(a-1)*sign(2*eta2(2)+eta2(1))+...
               2*a*abs(2*eta2(1)+eta2(2))^(a-1)*sign(2*eta2(1)+eta2(2));
dPhi_eta2(2) = 2*a*abs(2*eta2(2)+eta2(1))^(a-1)*sign(2*eta2(2)+eta2(1))+...
               a*abs(2*eta2(1)+eta2(2))^(a-1)*sign(2*eta2(1)+eta2(2));

% Deltas
r = zeros(2,1);
% Eq(A1.2)
r(1) = sqrt(((X1(1)-X1(2))/2.0)^2.0+X1(3)^2.0);
% Eq(A1.2) for Phi''
r(2) = sqrt(((X2(1)-X2(2))/2.0)^2.0+X2(3)^2.0);

% Check for singular case
dPhi1_dXx = zeros(1,3);
% Corrections for singular cases
if r(1) == 0 || eta1(1)==eta1(2)
    dPhi1_dXx(1) = 0;
    dPhi1_dXx(2) = 0;
    dPhi1_dXx(3) = 0;
else
    % Eq(A1.5)
    dEta1_X1 = zeros(2,3);
    dEta1_X1(1,1) = 1.0/2.0 + (X1(1) - X1(2))/(4.0*r(1));
    dEta1_X1(1,2) = 1.0/2.0 - (X1(1) - X1(2))/(4.0*r(1));
    dEta1_X1(1,3) = X1(3)/r(1);
    dEta1_X1(2,1) = 1.0/2.0 - (X1(1) - X1(2))/(4.0*r(1));
    dEta1_X1(2,2) = 1.0/2.0 + (X1(1) - X1(2))/(4.0*r(1));
    dEta1_X1(2,3) = -X1(3)/r(1);
    
    % Eq(A1.3)
    dPhi1_dXx = dPhi_eta1'*dEta1_X1;
end

% Check for singular case
dPhi2_dXx = zeros(1,3);
% Corrections for singular cases
if r(2) == 0 || eta2(1) == eta2(2)
    dPhi2_dXx(1) = dPhi_eta2(1);
    dPhi2_dXx(2) = dPhi_eta2(1);
    dPhi2_dXx(3) = 0;
else
    % Eq(A1.9)
    dEta2_X2 = zeros(2,3);
    dEta2_X2(1,1) = 1.0/2.0 + (X2(1) - X2(2))/(4.0*r(2));
    dEta2_X2(1,2) = 1.0/2.0 - (X2(1) - X2(2))/(4.0*r(2));
    dEta2_X2(1,3) = X2(3)/r(2);
    dEta2_X2(2,1) = 1.0/2.0 - (X2(1) - X2(2))/(4.0*r(2));
    dEta2_X2(2,2) = 1.0/2.0 + (X2(1) - X2(2))/(4.0*r(2));
    dEta2_X2(2,3) = -X2(3)/r(2);
    
    % Eq(1.7)
    dPhi2_dXx = dPhi_eta2'*dEta2_X2;
end

% Eq(18)
e_norm = transpose((dPhi1_dXx*L_1 + dPhi2_dXx*L_2));

end