% time-intation of the SLS-model given the deformation history
% Author: D. Baroudi, Dr.
% Version: 12.4.2017
% Calls:   function [DY] = Standard_Linear_Model_SLS(t, Y) 
% ---------------------------------------------------

%% global E_1 E_2 eta 
global E_0 E_infty tau
global t_TAB EPSILON_TAB dot_EPSILON_TAB 

tau = 0.8 % s , tau = eta / (E1 + E2) % (s) relaxation time

E_0     = 1500e6; % Pa , E_0 = G_0 = E1
E_infty =  500e6; % Pa, E_infty = E1*E2/(E1 + E2) 

sigma_0     = 0.0;      % alkujännitys (Pa);  
t_0         = 0;        % t = 0; initial time

%-- integration tiem
time          = [0 50*tau];     %  integrointiväli

% ------------------------------------------
% Given strain histories - keep active only example 1 or example 2 at the same time
% -------------------------------------------------------
% - Here you program the given deformation history epsilon(t)
% -------------------------------------------------------
% BEGIN Example 1: Heaviside ----
% ---------------------------
dt = tau/500;
times = 0:dt:max(time);
t_TAB = times;
k_0 = 5;
k = k_0  / 1;
H = 1 ./ ( 1 + exp(-2* k * times));

times = 0:dt:3*max(time);
NT     = length(times);

times(length(H)+1:length(times)) = [];

% epsilon history - given
epsilon_0       = 1e-2;
EPSILON_TAB     = epsilon_0 * H; 

% -------------------------------------------------------
% END Example 1: Heaviside ----
% ---------------------------

% -------------------------------------------------------
% BEGIN Example 2: periodic  ----
% ---------------------------
dt = tau/500;
times = 0:dt:max(time);
t_TAB = times;
NT     = length(times);

epsilon_0  = 5e-2; % (mm/mm)
omega      = 1; % (s)
epsilon = epsilon_0 * sin(omega * times);

% epsilon history - given
EPSILON_TAB     = epsilon; 

% -------------------------------------------------------
% END Example 2: periodic ----
% ---------------------------


dot_EPSILON_TAB = diff(EPSILON_TAB) ./ diff(times);
dot_EPSILON_TAB = [dot_EPSILON_TAB 0];


[Time, Y] = ode23s('Standard_Linear_Model_SLS',time,[sigma_0 t_0]);
Sigma = Y(:,1);
t     = Y(:,2);

epsilon_interp = interp1(times, epsilon, Time);

plot(Time/tau, Sigma/1e6);
xlabel('t/\tau')
ylabel('\sigma  (MPa)')
grid on

figure
plot(times/tau, 1000*EPSILON_TAB)
xlabel('t/\tau')
ylabel('\epsilon(t) \times 1000 mm/mm')
grid on

% --------------------------
figure
plot(times/tau, 1000*EPSILON_TAB)
hold on
plot(Time/tau, Sigma/1e6,'r-');
xlabel('t/\tau')
ylabel('\epsilon(t), \sigma(t)')
grid on
legend('1000\times \epsilon(t)','\sigma (MPa)')

% --------------------------



figure
plot(epsilon_interp, Sigma/1e6,'k.-')
grid on
xlabel('\epsilon')
ylabel('\sigma  (MPa)')
title('Response of SLS-Standard Linear Solid to given strain history \epsilon(t)')