Ski boot walking attachment MATLAB code
From DDL Wiki
gmodel
This is the model itself. It calls a few functions as data inputs - namely, makelayers and flexdata_cd_dat.
function gmodel
clear all
close all
A = [1 1 2 2 5 5];
B = [3 4 3 4 3 4];
%heel is zero
xlen = 12;
xres = .25;
zlen = 3;
ddstep = [0 0];
compression_modifier = [0 0];
xpos = 0:xres:xlen;
y0 = spline([0 4 6 10 12],[1.2,1.6,1.8,1.4,1.2],xpos); %initial height
body_mass = 70; %kg
for j = 1:6
dd = zeros(size(y0)); %deflection
rubber_layers = makelayers(A(j),B(j));
%begin step!
[foot_theta foot_length ankle_x ankle_y hip_x hip_y_nod lint all_y] = flexdata_cd_dat;
hip_y = all_y(:,1)-mean(all_y(:,1));
step_time = .50;
theta = linspace(-pi/6,pi/6,round(length(lint)*step_time));
time = linspace(0,step_time,length(theta));
%make subsets of everything (we care about foot contact only here)
foot_theta = foot_theta(1:length(theta));
ankle_x = ankle_x(1:length(theta));
ankle_y = ankle_y(1:length(theta));
hip_x = hip_x(1:length(theta));
hip_y_nod = hip_y_nod(1:length(theta))-mean(hip_y_nod(1:length(theta)));
hip_y = hip_y(1:length(theta));
theta = foot_theta; %LOL
hip_y_adj = hip_y;
for i = 3:length(theta)
stepsecs = time(i)-time(i-1);
%find lengths
lens = -((xpos-mean(xpos))*sin(theta(i)))+abs(y0*cos(theta(i)))-abs(dd*cos(theta(i)));
%Figure out which layers are contacting the ground (or are very
%close to doing so)
[lens_sort lens_sort_ix] = sort(lens);
incontact = sum(lens>=lens_sort(end)-ddstep(i-1)*compression_modifier(i-1));
contact_indexes = lens_sort_ix(end-incontact+1:end);
numcontact(i) = incontact;
%Full dispalcement continues until more layers are brought into
%contact. Because of the lenient "which in contact" logic, this new contact
%theoretically happens towards the end of the timestep, making such
%an assumption more reasonable.
ddstep(i) = lens_sort(end) - lens_sort(end-incontact);
hip_dely_dumb = hip_y(i) - hip_y(i-1);
hip_dy_dumb(i) = hip_dely_dumb/stepsecs; %dumb calc of hip y velocity
hip_dy2_dumb(i) = (hip_dy_dumb(i) - ((hip_y(i-1) - hip_y(i-2))/stepsecs))/stepsecs;
%force applied
k = rubber_layers(1,contact_indexes);
c = rubber_layers(2,contact_indexes);
%force expected based on body CoM motion
f_expected = hip_dy2_dumb(i) * body_mass;
%verify full compression
w0 = (sum(k)/body_mass)^.5;
gamma = sum(c)/(2*(sum(k)*body_mass)^.5);
wd = w0*(1-gamma)^.5;
compression_modifier(i) = stepsecs/wd*1000;
if compression_modifier(i) > 1
compression_modifier(i) = 1;
elseif compression_modifier(i) <.4
compression_modifier(i) = .4;
end
%Update displacements of layers
dd(contact_indexes) = dd(contact_indexes) + ddstep(i)*compression_modifier(i);
%catch layers which will be in contact for most of the step
dd(dd>y0*.5) = y0(dd>y0*.5)*.5;
%Solve for the vertical velocity of the *device*
f_app_k = sum(k ./ body_mass .* dd(contact_indexes));
f_app_c = f_expected - f_app_k;
hip_dy_adj(i) = -f_app_c/sum(c/body_mass)/5000; %5000 is a made-up constant,
%designed to adjust for lack of accurate dynamic properties.
%Thus, results presented are *proportional* to the correct
%results.
%Keep things from flying off the handle
if hip_dy_adj(i) < -10
hip_dy_adj(i) = -10;
end
%Aggregate dy*time into y
hip_y_adj(i) = hip_y_adj(i-1) + hip_dy_adj(i) * stepsecs;
end %for i
figure(3)
plot(time,hip_y,'k')
plot(time(1:round(.9*end)+1),hip_y_nod(round(.1*end):end)*.7,'r')
hold on
plot(time,hip_y+hip_y_adj/10)
figure(2)
%plot(time,hip_y)
hold on
plot(time,hip_y_adj)
figure(5)
hold on
plot(time,smooth(hip_dy_adj,21))
end %for
makelayers
Makes the physical model of the device's springs and dampers.
function rubber_layers = makelayers(A,B) %multipliers k0 = 200; c0 = 10; %neoprene soft foam k(1) = .07*k0; c(1) = .6*c0; %r1 = [k1;c1]; %neoprene hard foam k(2) = .15*k0; c(2) = .7*c0; %r2 = [k2;c2]; %neoprene soft sheet k(3) = .6*k0; c(3) = .1*c0; %r3 = [k3;c3]; %neoprene hard sheet k(4) = k0; c(4) = .1*c0; %r4 = [k4;c4]; %silicone sheet k(5) = .12*k0; c(5) = .8*c0; %r5 = [k5;c5]; %which layers where? sel = [repmat([A B],1,24) A]; rubber_layers(1,01:10) = [k(sel(1)) k(sel(2)) k(sel(3)) k(sel(4)) k(sel(5)) k(sel(6)) k(sel(7)) k(sel(8)) k(sel(9)) k(sel(10))]; rubber_layers(1,11:20) = [k(sel(11)) k(sel(12)) k(sel(13)) k(sel(14)) k(sel(15)) k(sel(16)) k(sel(17)) k(sel(18)) k(sel(19)) k(sel(20))]; rubber_layers(1,21:30) = [k(sel(21)) k(sel(22)) k(sel(23)) k(sel(24)) k(sel(25)) k(sel(26)) k(sel(27)) k(sel(28)) k(sel(29)) k(sel(30))]; rubber_layers(1,31:40) = [k(sel(31)) k(sel(32)) k(sel(33)) k(sel(34)) k(sel(35)) k(sel(36)) k(sel(37)) k(sel(38)) k(sel(39)) k(sel(40))]; rubber_layers(1,41:49) = [k(sel(41)) k(sel(42)) k(sel(43)) k(sel(44)) k(sel(45)) k(sel(46)) k(sel(47)) k(sel(48)) k(sel(49))]; rubber_layers(2,01:10) = [c(sel(1)) c(sel(2)) c(sel(3)) c(sel(4)) c(sel(5)) c(sel(6)) c(sel(7)) c(sel(8)) c(sel(9)) c(sel(10))]; rubber_layers(2,11:20) = [c(sel(11)) c(sel(12)) c(sel(13)) c(sel(14)) c(sel(15)) c(sel(16)) c(sel(17)) c(sel(18)) c(sel(19)) c(sel(20))]; rubber_layers(2,21:30) = [c(sel(21)) c(sel(22)) c(sel(23)) c(sel(24)) c(sel(25)) c(sel(26)) c(sel(27)) c(sel(28)) c(sel(29)) c(sel(30))]; rubber_layers(2,31:40) = [c(sel(31)) c(sel(32)) c(sel(33)) c(sel(34)) c(sel(35)) c(sel(36)) c(sel(37)) c(sel(38)) c(sel(39)) c(sel(40))]; rubber_layers(2,41:49) = [c(sel(41)) c(sel(42)) c(sel(43)) c(sel(44)) c(sel(45)) c(sel(46)) c(sel(47)) c(sel(48)) c(sel(49))];
flexdata_cd_dat
Prepares the kinesiologic data; calculates foot motion with the device.
function [foot_theta foot_length ankle_x ankle_y hip_x hip_y lint all_y] = flexdata_cd_dat
lint = linspace(0,1000,1000);
hip_height = -10*cos(lint*2*3.1516/1000);
%hip_theta = [162 162 170 178 185 193 183 170 160 156 158 162];
hip_theta = [157 164 170 178 185 193 183 170 160 154 156 157];
hip_time = [0 100 200 300 400 500 600 700 800 850 900 1000];
hip_theta_cubic = interp1(hip_time,hip_theta,lint,'spline');
knee_theta = [10 5 15 5 7 45 70 20 10];
knee_time = [0 100 250 450 500 600 750 950 1000];
knee_theta_cubic = interp1(knee_time,knee_theta,lint,'spline');
%ankle_theta = [0 -5 0 22 0 -8 8 0];
%ankle_time = [0 25 50 250 440 600 900 1000];
ankle_theta = [15 15];
ankle_time = [0 1000];
ankle_theta_cubic = interp1(ankle_time,ankle_theta,lint,'spline');
foot_length = 28;
shin_length = 56;
thigh_length = 70;
%{
figure(1)
%plot(hip_time,hip_theta)
plot(lint,hip_theta_cubic)
figure(2)
%plot(knee_time,knee_theta)
plot(lint,knee_theta_cubic)
figure(3)
%plot(ankle_time,ankle_theta)
plot(lint,ankle_theta_cubic)
figure(5)
plot(lint,hip_height)
%}
%
%
% %% That's the data. Now for the motion.
%
%
pi = 3.1416;
twopi = 2*pi;
hip_theta_rad = (180-hip_theta_cubic)*pi/180;
knee_theta_rad = knee_theta_cubic*pi/180;
ankle_theta_rad = ankle_theta_cubic*pi/180;
%hip_x(1:1000) = 0;
hip_x = linspace(1,170,1000);
hip_y = hip_height;
knee_x = hip_x + thigh_length.*sin(hip_theta_rad);
knee_y = hip_y - thigh_length.*cos(hip_theta_rad);
ankle_x = knee_x + shin_length.*sin(hip_theta_rad-knee_theta_rad);
ankle_y = knee_y - shin_length.*cos(hip_theta_rad-knee_theta_rad);
toe_x = ankle_x + foot_length.*sin(hip_theta_rad-knee_theta_rad+(1.9-ankle_theta_rad));
toe_y = ankle_y - foot_length.*cos(hip_theta_rad-knee_theta_rad+(1.9-ankle_theta_rad));
device_y = [.3 2 .3];
device_x = [0 .5 1];
device_y_cubic = interp1([0 .3 .5 .8 1],[1,1.5,2,1.6,1.2],linspace(0,1,100),'spline');
%device_y_cubic = interp1(device_x, device_y,linspace(0,1,100),'spline');
foot_theta = (hip_theta_rad-knee_theta_rad+(1.9-ankle_theta_rad))'-pi/2;
[device_theta, device_phi, device_r] = cart2sph(linspace(0,foot_length,100),-device_y_cubic,zeros(1,100));
device_theta_rot = repmat(device_theta,1000,1) + repmat(foot_theta,1,100);
device_r_rot = repmat(device_r,1000,1);
[device_x_rot device_y_rot dummy2] = sph2cart(device_theta_rot,zeros(1000,100),device_r_rot);
device_y_contact = min(device_y_rot,[],2);
device_y_add = min([device_y_rot(:,1) device_y_rot(:,100)],[],2)-device_y_contact;
device_y_add(600:1000) = zeros(1,401);
%figure(9)
%plot(lint,device_y_add)
all_x = [hip_x' knee_x' ankle_x' toe_x'];
all_y_raw = [hip_y' knee_y' ankle_y' toe_y'];
ankle_low = ankle_y.*(ankle_y<toe_y);
toe_low = toe_y.*(toe_y<ankle_y);
all_y_adjustment = repmat(ankle_low' + toe_low',1,4);
all_y = all_y_raw - all_y_adjustment;
all_y(600:1000,:) = all_y_raw(600:1000,:) - repmat(linspace(all_y_adjustment(600,1),all_y_adjustment(1000,1),401)',1,4);
all_y_device = all_y + repmat(device_y_add,1,4);
figure(4)
clf
step = 40;
hold on
i = 1;
k = 1;
while k < 600
plot(all_x(k,:),all_y_device(k,:))
i = i + 1;
k = k + step;
end
%{
figure(6)
clf
step = 40;
hold on
i = 1;
k = 1;
while k < 1000
plot(all_x(k,:),all_y_raw(k,:))
i = i + 1;
k = k + step;
end
figure(8)
clf
step = 40;
hold on
i = 1;
k = 1;
while k < 1000
plot(all_x(k,:),all_y_device(k,:))
i = i + 1;
k = k + step;
end
figure(7)
clf
hold on
plot(all_x(:,1),all_y(:,1)-mean(all_y(:,1)),'-')
plot(all_x(:,1),all_y_raw(:,1),'-.')
plot(all_x(:,1),all_y_device(:,1)-mean(all_y_device(:,1)),'--')
%}
