clear
clc

% This script solves Exercice 2 of Week #7


%% Parameters of the material law
ke=400000;
thetaP=0.02;
thetaPc=0.05;
My=4000;
Mu=4500;
thetaY=My/ke;
thetaC=thetaP+thetaY;
thetaU=thetaPc+thetaC;
ks=(Mu-My)/thetaP;
kpc=-Mu/thetaPc;

%% Loading history
% peaksVector=[0,0.10]; % Monotonic
% peaksVector=-[0,0.06,0.05,0.08]; % Cyclic
peaksVector=[0,0.06,-0.06,0]; % Cyclic
nbSteps=100;
strain=generateInputStrainHistory(peaksVector, nbSteps);

nStepsTot=size(strain,2);

%% Initialize
theta_previous=0;
M_previous=0;
MmaxPos=My;
MmaxNeg=-My;
yieldFlag_Pos=0;
cappingFlag_Pos=0;
yieldFlag_Neg=0;
cappingFlag_Neg=0;
reversalFlag=0;
residualFlag=0;
theta_M0_currentPos=0;
theta_M0_currentNeg=0;
theta_M0_projected=0;
Di_previous=0;


%% Call the function
M_vector=zeros(size(strain));
for i=1:nStepsTot
    theta=strain(i);
    [M,k,MmaxPos,MmaxNeg,yieldFlag_Pos,cappingFlag_Pos,yieldFlag_Neg,cappingFlag_Neg,reversalFlag,residualFlag,theta_M0_currentPos,theta_M0_currentNeg,theta_M0_projected,Di]...
        =computeTrilinearMatLaw(theta,theta_previous,M_previous,theta_M0_currentPos,theta_M0_currentNeg,theta_M0_projected,ke,ks,kpc,thetaY,thetaC,thetaU,My,Mu,MmaxPos,MmaxNeg,yieldFlag_Pos,cappingFlag_Pos,yieldFlag_Neg,cappingFlag_Neg,reversalFlag,residualFlag,Di_previous);
    M_vector(i)=M;
    
    M_previous=M;
    theta_previous=theta;
    Di_previous=Di;
end


%% Plot the results
h1=figure;
plot(strain,M_vector/1000,'-k')
xlabel('\theta [rad]')
ylabel('M [kNm]')
grid on

plot_settings_ASCE(h1)



%% Function to return the moment and tangent
function [M,k,MmaxPos,MmaxNeg,yieldFlag_Pos,cappingFlag_Pos,yieldFlag_Neg,cappingFlag_Neg,reversalFlag,residualFlag,theta_M0_currentPos,theta_M0_currentNeg,theta_M0_projected,Di]...
    =computeTrilinearMatLaw(theta,theta_previous,M_previous,theta_M0_currentPos,theta_M0_currentNeg,theta_M0_projected,ke,ks,kpc,thetaY,thetaC,thetaU,My,Mu,MmaxPos,MmaxNeg,yieldFlag_Pos,cappingFlag_Pos,yieldFlag_Neg,cappingFlag_Neg,reversalFlag,residualFlag,Di_previous)
    % Find the direction of the current increment "Di": 1:if Ri is moving towards positive and -1: if Ri is moving negative
    if theta>=theta_previous
        Di=1;
	else
		Di=-1;
    end
    
    % Trial elastic moment
    M=M_previous+ke*(theta-theta_previous);
    
    % Quick return if failed
    if residualFlag==1
        M=0;
        k=0;
        return
    end
    
    %If inelastic and going to positive rotation
    if Di==1 && M>MmaxPos
        yieldFlag_Neg=0;
        cappingFlag_Neg=0;
        if yieldFlag_Pos==1 || cappingFlag_Pos==1
        else
            theta_M0_currentPos=theta_M0_projected;
        end
        
        if abs(theta_M0_currentPos-theta)>=thetaY && abs(theta_M0_currentPos-theta)<thetaC %hardening path
            yieldFlag_Pos=1;
            theta_M0_currentNeg=theta_M0_projected;
            M=Di*My+ks*(theta-theta_M0_currentPos-thetaY);
            MmaxNeg=-My;
        elseif abs(theta_M0_currentPos-theta)>=thetaC && abs(theta_M0_currentPos-theta)<thetaU %softening path
            yieldFlag_Pos=0;
            cappingFlag_Pos=1;
            M=Di*Mu+kpc*(theta-theta_M0_currentPos-thetaC);
        else %residual path
            residualFlag=1;
            M=0;
        end
    
        MmaxPos=M;
    end
    
    
    %If inelastic and going to negative rotation
    if Di==-1 && M<MmaxNeg
        yieldFlag_Pos=0;
        cappingFlag_Pos=0;
        if yieldFlag_Neg==1 || cappingFlag_Neg==1
        else
            theta_M0_currentNeg=theta_M0_projected;
        end
        
        if abs(theta_M0_currentNeg-theta)>=thetaY && abs(theta_M0_currentNeg-theta)<thetaC %hardening path
            yieldFlag_Neg=1;
            theta_M0_currentPos=theta_M0_projected;
            M=Di*My+ks*(theta-theta_M0_currentNeg+thetaY);
            MmaxPos=My;
        elseif abs(theta_M0_currentNeg-theta)>=thetaC && abs(theta_M0_currentNeg-theta)<thetaU
            yieldFlag_Neg=0;
            cappingFlag_Neg=1;
            M=Di*Mu+kpc*(theta-theta_M0_currentNeg+thetaC);
        else %failed
            residualFlag=1;
            M=0;
        end
    
        MmaxNeg=M;
    end
   
    
    % Check if unloading
    if Di_previous*Di<0 
       theta_M0_projected=theta_previous-M_previous/ke;
    end    
    
    % Compute the tangent stiffness
    k=(M-M_previous)/(theta-theta_previous);
    if theta==theta_previous
        k=0;
    elseif M==M_previous
        k=10^-6;
    end
end


%% Function to compute the input strain history
function strain=generateInputStrainHistory(peaksVector, nbSteps)
%nbSteps is used  in maxmim excursion, other excursions are scaled proportionaly

% General informations
peaksVector=[0,peaksVector];
nbLoadExcursions=size(peaksVector,2)-1;

% Find maximum load excursion
strainDiffLoadExcursions=abs([0,diff(peaksVector)]);
maxStrainDifference=max(strainDiffLoadExcursions);

% Generate the points for each load excursions
strain=[];
for i=1:nbLoadExcursions
    nbStepsTemp=round(strainDiffLoadExcursions(i+1)/maxStrainDifference*nbSteps);
    startStrain=peaksVector(i);
    endStrain=peaksVector(i+1);
    epsi11Temp=linspace(startStrain,endStrain,nbStepsTemp);
    if i<nbLoadExcursions
        epsi11Temp=epsi11Temp(1:end-1); %Remove duplicate strain
    end
    strain=[strain,epsi11Temp];
end
end