t=[100 101 102 103 104 105]
theta=[0.75 0.72 0.7 0.68 0.67 0.66]
radius=[1560 1636 1694 1767 1838 1902]
%Delta=1s
Delta=1
figure(1)
plot(t,theta,’ko-‘)
xlabel(‘t(s)’);ylabel(‘\theta’)
figure(2)
plot(t,radius,’bo-‘)
xlabel(‘t(s)’);ylabel(‘r(m)’)
%Index 4 correspond to t=103
rdot4=(radius(5)-radius(3))/(2*Delta);
thetadot4=(theta(5)-theta(3))/(2*Delta);
rdoubledot4=(radius(5)-2*radius(4)+radius(3))/(Delta.^2);
thetadoubledot4=(theta(5)-2*theta(4)+theta(3))/(Delta.^2);
ur4=rdot4;
utheta4=radius(4)*thetadot4;
ar4=rdoubledot4-radius(4)*thetadot4.^2;
atheta4=radius(4)*thetadoubledot4+2*rdot4*thetadot4;
sprintf(‘The radial and azimuthal velocity at t=103s is %f %f ‘,ur4,utheta4)
sprintf(‘The radial and azimuthal acceleration at t=103s is %f %f ‘,ar4,atheta4)
%Get data to plot magnitude of velocity and acceleration
ur=zeros(size(t));
utheta=zeros(size(t));
ar=zeros(size(t));
atheta=zeros(size(t));
rdot=zeros(size(t));
thetadot=zeros(size(t));
rdoubledot=zeros(size(t));
thetadoubledot=zeros(size(t));
rdot(1)=(radius(2)-radius(1))/Delta;
thetadot(1)=(theta(2)-theta(1))/Delta;
rdoubledot(1)=(radius(3)-2*radius(2)+radius(1))/(Delta.^2);
thetadoubledot(1)=(theta(3)-2*theta(2)+theta(1))/(Delta.^2);
N=length(t) %number of data points
for i=2:N-1
rdot(i)=(radius(i+1)-radius(i-1))/(2*Delta);
thetadot(i)=(theta(i+1)-theta(i-1))/(2*Delta);
rdoubledot(i)=(radius(i+1)-2*radius(i)+radius(i-1))/(Delta.^2);
thetadoubledot(i)=(theta(i+1)-2*theta(i)+theta(i-1))/(Delta.^2);
end
rdot(N)=(radius(N)-radius(N-1))/Delta;
thetadot(N)=(theta(N)-theta(N-1))/Delta;
rdoubledot(N)=(radius(N-2)-2*radius(N-1)+radius(N))/(Delta.^2);
thetadoubledot(N)=(theta(N-2)-2*theta(N-1)+theta(N))/(Delta.^2);
ur=rdot;
utheta=radius.*thetadot;
ar=rdoubledot-radius.*thetadot.^2;
atheta=radius.*thetadoubledot+2.*rdot.*thetadot;
figure(3)
subplot(221)
plot(t,ur,’ko-‘)
xlabel(‘t(s)’);ylabel(‘U_r(m/s)’)
subplot(222)
plot(t,utheta,’ko-‘)
xlabel(‘t(s)’);ylabel(‘U_{\theta}(m/s)’)
subplot(223)
plot(t,ar,’ko-‘)
xlabel(‘t(s)’);ylabel(‘A_r(m/s^2)’)
subplot(224)
plot(t,atheta,’ko-‘)
xlabel(‘t(s)’);ylabel(‘A_{\theta}(m/s^2)’)