getPolynomialExpressions.m

% This function is part of the NMSM Pipeline, see file for full license.
%
% This function creates muscle-specific symbolic polynomial expressions for  
% muscle tendon length and moment arms.
%
% Inputs:
% jointAngles (numberFrames x degreesOfFreedom)
% polynomialDegree (value)
%
% (2D Number matrix, Number) -> (Symbol array, 2D Symbol matrix)
%
% returns muscle-specific polynomial expressions based on polynomial
% degree and the number of degrees of freedom.

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

function [polynomialExpressionMuscleTendonLength, ...
    polynomialExpressionMuscleTendonVelocities, ...
    polynomialExpressionMomentArms, theta] = ...
    getPolynomialExpressions(jointAngles, polynomialDegree)

% Initialize symbolic thetas
theta = sym('theta', [1 size(jointAngles, 2)]);
thetaDot = sym('thetaDot', [1 size(jointAngles, 2)]);
% Create polynomial basis function
basisFunction = 1;
for i = 1 : size(jointAngles, 2)
    basisFunction = theta(i) + basisFunction;
end
% Approx. muscle tendon length using polynomial basis equation to nth degree
polynomialExpressionMuscleTendonLength = basisFunction .^ polynomialDegree;
% Reformate polynomial expression
polynomialExpressionMuscleTendonLength = ...
    children(expand(polynomialExpressionMuscleTendonLength));
polynomialExpressionMuscleTendonLength = ...
    cat(2, polynomialExpressionMuscleTendonLength{:});
% Differentiate muscle tendon length w.r.t. time
for i = 1 : size(jointAngles, 2)
    for j = 1 : size(polynomialExpressionMuscleTendonLength, 2)
        syms(sprintf('theta%d(t)', i))
        theta_t(i) = symfun(eval(sprintf('theta%d(t)', i)), t);
        tempPolynomialExpressionMuscleTendonVelocity = subs(polynomialExpressionMuscleTendonLength(1, j), theta(i), theta_t(i));
        tempPolynomialExpressionMuscleTendonVelocity = diff(tempPolynomialExpressionMuscleTendonVelocity, t);
        tempPolynomialExpressionMuscleTendonVelocity = subs(tempPolynomialExpressionMuscleTendonVelocity, diff(theta_t(i), t), thetaDot(i));
        polynomialExpressionMuscleTendonVelocity(i, j) = subs(tempPolynomialExpressionMuscleTendonVelocity, theta_t(i), theta(i));
    end
end
polynomialExpressionMuscleTendonVelocities = polynomialExpressionMuscleTendonVelocity(1, :);
for i = 1 : size(jointAngles, 2)
    for j = 1 : size(polynomialExpressionMuscleTendonLength, 2)
        polynomialExpressionMuscleTendonVelocities(1, j) = polynomialExpressionMuscleTendonVelocity(i, j) + polynomialExpressionMuscleTendonVelocities(1, j);
    end
end
% Differentiate -muscle tendon length w.r.t. associated joint angle
for i = 1 : size(jointAngles, 2)
    polynomialExpressionMomentArms(i, :) = ...
        diff(-polynomialExpressionMuscleTendonLength, theta(i));
end
end