I have attempted to use the source code to write a MATLAB script for determining the CLF at a prescribed angle (attached below).
I'm uncertain whether I have correclty implemented the simmspline function and the CLF transitions and would very much appreciate any feedback.
Code: Select all
import org.opensim.modeling.*;
%CLF parameters for the ARMS Lab Hand and Wrist model (McFarland et al. 2023)
%FLEXION-EXTENSION DOF
%Units: deg and Nm/deg
upperLimit = 60;
lowerLimit = -60;
upperStiffness = 200;
lowerStiffness = 135.3232;
transitionRange = 92.916;
% Calculate limit force at angle = -10 degrees
limitForce = calcLimitForce(upperLimit, lowerLimit, upperStiffness, lowerStiffness, transitionRange);
disp(['CLF at angle = -10: ', num2str(limitForce)]);
%GRAPH
angleRange = -200:0.5:200;
limitForces = zeros(size(angleRange));
%Simmsplinbe in transition regions for lower and upper stiffness
%upper
upperSpline = SimmSpline();
upperSpline.addPoint(upperLimit, 0);
upperSpline.addPoint(upperLimit + transitionRange, upperStiffness);
%lower
lowerSpline = SimmSpline();
lowerSpline.addPoint(lowerLimit - transitionRange, lowerStiffness);
lowerSpline.addPoint(lowerLimit, 0);
% Calculate limit forces for each angle in the range
for i = 1:length(angleRange)
angle = angleRange(i);
% Determine the current stiffness using SimmSpline functions
upperK = calculateUpperStiffness(angle, upperSpline, upperLimit, upperStiffness, transitionRange);
lowerK = calculateLowerStiffness(angle, lowerSpline, lowerLimit, lowerStiffness, transitionRange);
% Calculate the total limit force
upperForce = -upperK * (angle - upperLimit);
lowerForce = lowerK * (lowerLimit - angle);
limitForces(i) = upperForce + lowerForce;
end
% Plot the results
figure;
plot(angleRange, limitForces, 'k', 'LineWidth', 2);
grid on;
xlabel('Coordinate Position (deg)');
ylabel('Coordinate Limit Force (Nm)');
title('Coordinate Limit Force for Extension-Flexion DOF');
% Mark the important lines
% ROM for flexon-extension
xline(-70, '--', 'LineWidth', 1);
xline(70, '--', 'LineWidth', 1);
% show upper and lower limit
xline(lowerLimit, 'r', 'LineWidth', 1);
xline(upperLimit, 'r', 'LineWidth', 1);
% show angle at which upper and lower stiffness is reached
xline(lowerLimit - transitionRange, 'b', 'LineWidth', 1);
xline(upperLimit + transitionRange, 'b', 'LineWidth', 1);
%Limits in flexion-extension = 70-70
limitForceAtNeg70 = interp1(angleRange, limitForces, -70);
limitForceAt70 = interp1(angleRange, limitForces, 70);
text(-70, limitForceAtNeg70, sprintf('%.2f', limitForceAtNeg70), ...
'HorizontalAlignment', 'right', 'VerticalAlignment', 'bottom', 'FontSize', 12);
text(70, limitForceAt70, sprintf('%.2f', limitForceAt70), ...
'HorizontalAlignment', 'left', 'VerticalAlignment', 'bottom', 'FontSize', 12);
% Calculate CLF at a specific angle
function limitForce = calcLimitForce(upperLimit, lowerLimit, upperStiffness, lowerStiffness, transitionRange)
import org.opensim.modeling.*;
angle = -10;
upperSpline = SimmSpline();
upperSpline.addPoint(upperLimit, 0);
upperSpline.addPoint(upperLimit + transitionRange, upperStiffness);
lowerSpline = SimmSpline();
lowerSpline.addPoint(lowerLimit - transitionRange, lowerStiffness);
lowerSpline.addPoint(lowerLimit, 0);
%Stiffness
upperK = calculateUpperStiffness(angle, upperSpline, upperLimit, upperStiffness, transitionRange);
lowerK = calculateLowerStiffness(angle, lowerSpline, lowerLimit, lowerStiffness, transitionRange);
%Force
upperForce = -upperK * (angle - upperLimit);
lowerForce = lowerK * (lowerLimit - angle);
limitForce = upperForce + lowerForce;
end
% Calculate upper stiffness using SimmSpline with constant behavior beyond the transition zone
function upperK = calculateUpperStiffness(angle, upperSpline, upperLimit, upperStiffness, transitionRange)
import org.opensim.modeling.*;
if angle < upperLimit
upperK = 0; % Zero stiffness below the upper limit
elseif angle > (upperLimit + transitionRange)
upperK = upperStiffness; % Constant stiffness beyond the transition range
else
upperK = upperSpline.calcValue(Vector(1, angle));
end
end
% Calculate lower stiffness using SimmSpline with constant behavior beyond the transition zone
function lowerK = calculateLowerStiffness(angle, lowerSpline, lowerLimit, lowerStiffness, transitionRange)
import org.opensim.modeling.*;
if angle > lowerLimit
lowerK = 0; % Zero stiffness above the lower limit
elseif angle < (lowerLimit - transitionRange)
lowerK = lowerStiffness; % Constant stiffness beyond the transition range
else
lowerK = lowerSpline.calcValue(Vector(1, angle));
end
end