- Author:
- Kenneth Tran <k.tran@auckland.ac.nz>
- Date:
- 2019-01-08 16:52:34+13:00
- Desc:
- Uploading Matlab files.
- Permanent Source URI:
- https://models.fieldml.org/workspace/567/rawfile/902a5740b3e6ac6b210dc436d5f25a82ac90070a/createFigure.m
% ESZ framework function to create the background ESZ figure
% This function is called by ESZ_PlotAllFigures.m
% 13 Sept 2018
% Author: Kenneth Tran
function figH = createFigure(isom,a,b,c,d)
% Input arguments to control aspects of plotting
% isom = 1 to plot the isometric ESSLR, isom = 0 to hide it
% 'a' sets the xticks and 'b' sets the yticks
% 'c' sets the x-axis upper lim, 'd' sets y-axis upper lim
if nargin == 5
xtick = a;
ytick = b;
xlim = c;
ylim = d;
else
xtick = [0.75 0.8 0.85 0.9 0.95 1];
ytick = [0 20 40];
xlim = [0.73,1.01];
ylim = [0,45];
end
% 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
% 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.83, 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
% Interpolating to define the lower boundary from the work-loop ESSLRs
Ln_lower = linspace(min(Les_WL(:,1)), max(Les_WL(:,1)));
S_lower = interp1(Les_WL(:,1),S_WL(:,1),Ln_lower,'spline');
figH = figure; hold on;
% Shading in the ESZ
cc = 0.9.*[1 1 1];
xarea1 = Ln_lower;
yarea1 = quadS(Ln_lower,a_isom,b_isom,c_isom);
yarea2 = S_lower;
xarea2 = [xarea1,fliplr(xarea1)];
inBetween = [yarea1,fliplr(yarea2)];
fill(xarea2,inBetween,cc,'LineStyle','none');
xarea1 = linspace(max(Les_WL(:,1)),1);
yarea1 = quadS(xarea1,a_isom,b_isom,c_isom);
yarea2 = quadS(Les_WL(size(Les_WL,1),:)/Li(size(Les_WL,1)),a_WL,b_WL,c_WL)*quadS(Li(size(Les_WL,1)),a_isom,b_isom,c_isom);
xarea2 = [xarea1,fliplr(xarea1)];
inBetween = [yarea1,fliplr(yarea2)];
fill(xarea2,inBetween,cc,'LineStyle','none');
if isom
plot(Ln, S_isom, 'k-', 'linewidth', 1);
end
plot(Ln, S_passive, 'k-', 'linewidth', 1);
xlabel('L/L_o','fontangle','italic','fontsize',fs+1)
ylabel('Stress (kPa)','fontsize',fs+1)
axis([xlim ylim])
set(gca,'xtick',xtick,'fontsize',fs)
set(gca,'ytick',ytick,'fontsize',fs)
box off
