clc
clear all
close all

% structural data
m = 1;
k1=1;
k2=0;
d1=1;
d2=1000; %useless, choose it large enough

omega=sqrt(k1/m); T1=2*pi/omega;
R = k1*d1;

listTd=10.^[-1:0.05:2];
listR=[0.2 0.4 0.6 0.8 1 1.02 1.05 1.1 1.2 1.3 1.4 1.6 1.8 2];

% rectangular pulse
for i2=1:length(listR)
    for i1=1:length(listTd)
        Td = listTd(i1)*T1; DT=min([Td,T1])/100;
        P = R/listR(i2);

        nL=floor(Td/DT);   Nstep=10*nL;
        time = (0:Nstep-1)*DT;
        
        p= P * [ones(nL,1);zeros(Nstep-nL,1)];
        % p(1:nL)=P*(1-(0:nL-1)*DT/Td);
        
        [dep2,fint2,regime2] = DiffCentraleElastoPlastique(m,k1,k2,p,DT,Nstep,d1);
        ductility(i1,i2) = max(abs(dep2/d1))
    end
end

%%
subplot(121)
loglog(listTd,ductility); xlabel ('t_d/T'); ylabel ('\mu')
ylim ([0.1 100])

subplot(122)
semilogy(listR,ductility); xlabel ('R/P'); ylabel ('\mu')
ylim ([0.1 100])