% 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



