% ESZ framework code to plot the stress-length relations that were fitted
% to experimental data (Figure 1).
% The family of work-loop ESSLRs is generated in Figure 2.
% 13 Sept 2018
% Author: Kenneth Tran

clear all;

% Define functions here
% Quadratic function for isometric work-loop stress
quadS = @(x,a,b,c) a*x.^2 + b*x + c;    

% Exponetial function for the passive stress
expS = @(x,a,b) a*exp(b*x);             

% Solving the quadratic function for the end-systolic length
quadL = @(a,b,c) (-b+sqrt(b*b-4*a.*c))./(2*a);

%% Parameters

clear Les_WL S_pass L_low S_WL Li

fs = 12;

% Parameters for isometric ESSLR
a_isom = 581.3;
b_isom = -852.63; 
c_isom = 314.41;

% Parameters for passive relation
a_pass = 3.289e-7;
b_pass = 17.202;

% Parameters for work-loop ESSLR at Lo 
a_WL1 = 451.9; 
b_WL1 = -489.96;
c_WL1 = 81.709;

% Parameters for work-loop ESSLR at 0.95 Lo 
a_WL2 = 730.79;
b_WL2 = -1066.87;
c_WL2 = 382.73;

% Parameters for work-loop ESSLR at 0.9 Lo 
a_WL3 = 215.28;
b_WL3 = -177.93;
c_WL3 = 4.3608;

% Parameters for normalised work-loop ESSLR
% Used to construct ESSLRs at any Li
a_WL = 15.0202;
b_WL = -19.8935;
c_WL = 5.8743;

% Normalised length array
Ln = linspace(0.74,1);

%% Isometric stress-length relation
S_isom = quadS(Ln,a_isom,b_isom,c_isom);

%% Passive-stress length relation
S_passive = expS(Ln,a_pass,b_pass);

%% Work-loop ESSLRs

% Work-looop ESSLR at Lo
Li = 1; % Initial muscle length at Lo
S_pass = expS(Li,a_pass,b_pass); % Passive force at Li
L_low = quadL(a_WL1,b_WL1,c_WL1-S_pass);
Les_WL1 = linspace(L_low, Li);
S_WL1 = quadS(Les_WL1,a_WL1,b_WL1,c_WL1);

% Work-looop ESSLR at 0.95 Lo
Li = 0.95; % Initial muscle length
S_pass = expS(Li,a_pass,b_pass); % Passive force at Li
L_low = quadL(a_WL2,b_WL2,c_WL2-S_pass);
Les_WL2 = linspace(L_low, Li);
S_WL2 = quadS(Les_WL2,a_WL2,b_WL2,c_WL2);

% Work-looop ESSLR at 0.9 Lo
Li = 0.9; % Initial muscle length
S_pass = expS(Li,a_pass,b_pass); % Passive force at Li
L_low = quadL(a_WL3,b_WL3,c_WL3-S_pass);
Les_WL3 = linspace(L_low, Li);
S_WL3 = quadS(Les_WL3,a_WL3,b_WL3,c_WL3);

figure(1)
subplot(2,1,1)
plot(Ln, S_isom, 'k-', 'linewidth', 2); hold on;
plot(Ln, S_passive, 'k-', 'linewidth', 2);
plot(Les_WL1, S_WL1, 'b.',  'linewidth', 1);
plot(Les_WL2, S_WL2, 'g.',  'linewidth', 1);
plot(Les_WL3, S_WL3, 'm.',  'linewidth', 1);
xlabel('L/L_o','fontangle','italic','fontsize',fs+1)
ylabel('Stress (kPa)','fontsize',fs+1)
axis([0.74 1.01 0 45])
set(gca,'xtick',[0.75 0.8 0.85 0.9 0.95 1],'fontsize',fs)
set(gca,'ytick',[0 20 40],'fontsize',fs)
box off

%% Using the double normalised work-loop ESSLR to interpolate family of curves

clear Les_WL S_pass L_low S_WL Li

n = 100; % Specify number of work-loop ESSLRs
Li = linspace(0.80, 1,n); % Specify range of initial muscle length

for i=1:length(Li)
    S_pass(i) = expS(Li(i),a_pass,b_pass)/quadS(Li(i),a_isom,b_isom,c_isom); % Normalising passive force at each Li
    L_low(i) = quadL(a_WL,b_WL,c_WL-S_pass(i))*Li(i); % Calculating the lower limit and scaling to L/Lo units
    Les_WL(i,:) = linspace(L_low(i), Li(i)); % Specifying range of end-systolic values 
    S_WL(i,:) = quadS(Les_WL(i,:)/Li(i),a_WL,b_WL,c_WL)*quadS(Li(i),a_isom,b_isom,c_isom); % Calculate afterload stress given Les range of values
end

subplot(2,1,2)
plot(Les_WL', S_WL','Color', 0.7.*[1 1 1],'linewidth', 1); hold on;
plot(Ln, S_isom, 'k-', 'linewidth', 2); hold on;
plot(Ln, S_passive, 'k-', 'linewidth', 2);
xlabel('L/L_o','fontangle','italic','fontsize',fs+1)
ylabel('Stress (kPa)','fontsize',fs+1)
axis([0.74 1.01 0 45])
set(gca,'xtick',[0.75 0.8 0.85 0.9 0.95 1],'fontsize',fs)
set(gca,'ytick',[0 20 40],'fontsize',fs)
box off


