File: ef230-2023-08/www/examples/curvature.m Download
%| Demonstrate symbolic math in MATLAB
%| Calculate the radius of curvature for trajectory equation
%| Calculate normal and tangential acceleration components
function curvature(~,~)
% demonstrate symbolic math in MATLAB
% calculate the radius of curvature for trajectory equation
% calculate normal and tangential acceleration components
% from radius of curvature, determine v based on slope and constant vx
% calculate normal (v^2/r) and tengential (atot - at) acceleration components
% plot results
% Will Schleter, UTK, Fall, 2012

clc;
format compact;

txt = {'This app allows parameters for a trajectory to be calculated'
       'and displayed. It demonstrates the use of the symbolic toolbox.'
       'It is also an example of how to combine a plot and a number table'}

waitfor(msgbox(txt));

figure(1)

% The trajectory equation from EF 151
ysym = sym('y0+tan(th)*(x-x0)-g/(2*v0^2)*(1+(tan(th))^2)*(x-x0)^2')

% first derivative (slope)
dydx = diff(ysym)

% 2nd derivative
d2ydx2 = diff(ysym,2)

% general radius of curvature formula (from ef151 lecture 4.2)
rsym = (1+dydx^2)^(3/2)/abs(d2ydx2)


g=32.2; % gravity in ft/s/s
x0=0; % starting position in ft
y0=0;
v0=15; % launch velocity in ft/s
thd=0; % launch angle in degrees
xf=5; % final x in ft

prompts={'Gravity'
    'Initial X'
    'Initial Y'
    'Initial Velocity'
    'Initial Angle (degrees)',
    'Final X'}


while 1
defs = {num2str(g),num2str(x0),num2str(y0),num2str(v0),num2str(thd),num2str(xf)};

a = inputdlg(prompts,'Trajectory Curvature Values',1,defs);
if isempty(a)
    close(1);
    return
end
pause(0.1); % not sure why this is needed, but it helps prevent hanging on multiple runs
% now apply some numbers
g=str2double(a{1});
x0=str2double(a{2});
y0=str2double(a{3});
v0=str2double(a{4});
thd=str2double(a{5});
xf=str2double(a{6});

x = linspace(x0,xf,20); % a list of x values
th = thd*pi/180;

% calculated values
y = double(subs(ysym)) % calculate the corresponding y
r = double(subs(rsym)) % calculate the corresponding radius of curvature
slope = double(subs(dydx)) % the slope at each point

vx = v0.*cos(th).*ones(size(x)) % vx is constant - use the ones command to make it a list
vy = slope.*vx % v is tangent to the path, so vy/vx = slope
v =  sqrt(sum(([vx;vy]).^2)) % determine total v based on components
a_tan = g.*vy./v;
a_normal = g.*vx./v;
%a_normal = v.^2./r % definition of normal acceleration magnitude
r1 = v.^2./a_normal;

% Plotting
figure(1)
clf
% the graph on top
subplot(2,1,1);
plot(x,y,'-')
axis equal
title('Trajectory Plot')
% the table on bottom
subplot(2,1,2)
cnames = {'X','Y','R','R1','A_n','A_t','V','Vx','Vy','dy/dx'}; % column names
dat = [x;y;r;r1;a_normal;a_tan;v;vx;vy;slope]'; % put data in one big 2d matrix
t = uitable('Parent',gcf,'Data',dat,'ColumnName',cnames,...
            'units','normalized','Position',[0 0 1 .5]);

end

end % curvature function