ex2_P2

% exam 2
% P2

clear all; clc; close all;

syms q dq ddq L m g F beta

% d(q)/dt = dq
% d^2(q)/dt^2 = ddq

c = cos(q);
s = sin(q);

omega = [0 0 dq ];
alpha = [0 0 ddq];

rC = [0 0 0];
xB = L*c;
yB = L*s;
rB = [xB yB 0];
rA = -rB;

aC = [0 0 0];

IC  = m*(2*L)^2/12;

G = [0 -m*g 0];
FA = [F*cos(beta) F*sin(beta) 0];

fprintf('IC*alpha = OA x FA = 0 => \n')

eqM  = -IC*alpha + cross(rA, FA);
eqMz = eqM(3);

fprintf('%s = 0 \n',char(vpa(eqMz)))

slist = {L, m, g, F, beta};
nlist = {1, 1, 10, 100, pi/6};

eqMn = subs(eqMz,slist,nlist);
fprintf('%s = 0 \n',char(vpa(eqMn)))

fprintf('\n')
fprintf('m*aC = G + FA + F01 = 0 => \n')
F01 = -(G+FA);
F01n = subs(F01,slist,nlist);
fprintf('F01 = [%6.3f,%6.3f,%g] (N)\n',F01n)

fprintf('\n')
alphas = solve(eqMn,'ddq');
alpha0 = subs(alphas,q,pi/4);
fprintf('alpha0 = %6.3f (rad/s^2) \n',alpha0)

% end of program

% IC*alpha = OA x FA = 0 =>
% F*L*cos(beta)*sin(q) - 0.333*L^2*ddq*m - 1.0*F*L*cos(q)*sin(beta) = 0
% 86.6*sin(q) - 50.0*cos(q) - 0.333*ddq = 0
%
% m*aC = G + FA + F01 = 0 =>
% F01 = [-86.603,-40.000,0] (N)
%
% alpha0 = 77.646 (rad/s^2)

IC*alpha = OA x FA = 0 => 
F*L*cos(beta)*sin(q) - 0.33333333333333333333333333333333*L^2*ddq*m - 1.0*F*L*cos(q)*sin(beta) = 0 
86.602540378443864676372317075294*sin(q) - 50.0*cos(q) - 0.33333333333333333333333333333333*ddq = 0 

m*aC = G + FA + F01 = 0 => 
F01 = [-86.603,-40.000,0] (N)

alpha0 = 77.646 (rad/s^2)