
% some of the values that are set in this file are not used until later 
% parts of the problem

global r r2 rg lthg;
global iwk iwu1 iwu2 ium;
global w_data;
global gc g_by_gc; 
global m ai;

% figure counter
ifig = 0;

% some of geometry of linkage
%
% lengths of bars in linkage (meters)
r = 0.01*[10.0 40.0 30.0 40.0]';
r2 = r*r';
% angle (relative to local coordinates on link) to the center of mass;
% see m-file linkage_ag.m
lthg = [0.0 0.0 pi/6]';
% distance to centers of mass from starting point of link
rg = [r(1)*0.5 r(2)*0.5 r(3)/cos(lthg(3))]';

% masses, moments of inertia
m = [0.0 0.0 10.0]'; % in kg
% moment of inertia for uniform distribution of mass over links
ai = zeros(3,1);
ai(1) = m(1)*r2(1,1)/12.0;
ai(2) = m(2)*r2(2,2)/12.0;
diam_fan = 1.0;
ai(3) = m(3)*diam_fan*diam_fan/16.0;

% gravitation
g       = [0.0 9.8]'; % direction 2 is down, in m/s^2
gc      = 1.0; % metric, trust mass to be in kg
g_by_gc = g/gc;

th1_init = atan2(3,4);
th_init = fan_angles(r,r2,th1_init)

% data for initial angular velocities
w_data = [2.0*pi 0.0 0.0]';
% indices for known (iwk) and unknown (iwu*) angular velocities
% and index for unknown moment (ium)
iwk = 1; iwu1 = 2; iwu2 = 3; ium = 1;

tstart = 0.0; tstop = 2.0*pi/abs(w_data(iwk)); nt = 101;
tspan = linspace(tstart, tstop, nt);
dtime = (tstop - tstart)/(nt-1);


options = odeset('AbsTol',1e-10);
[t, th_num] = ode45('fan_solve',tspan,th_init,options);



% obtain accelerations over a full cycle

ag_fan = zeros(nt,2);   mag_ag_fan = zeros(nt,1);
thddot_fan = zeros(nt); mag_thddot_fan = zeros(nt,1);
for it = 1:nt;
   [thdot,thddot,ag] = fan_solve(t(it), th_num(it,:));
   % store results at time it
   ag_fan(it,:) = ag(:,3)';
   mag_ag_fan(it) = norm(ag(:,3));
   thddot_fan(it) = thddot(3);
   mag_thddot_fan(it) = abs(thddot(3));
end;

% plot angular acceleration of link 3

ifig = ifig+1; figure(ifig);
plot(t,thddot_fan);
xlabel('time (s)'); ylabel('angular acceleration (rad/s^2)');
title('angular acceleration of fan for a full cycle');

% plot componenets of linear acceleration of fan's center of mass over cycle

ifig = ifig+1; figure(ifig);
plot(t,ag_fan(:,1),t,ag_fan(:,2));
xlabel('time (s)'); ylabel('linear acceleration (m/s^2)');
title('components of acceleration of fan for a full cycle');

disp(sprintf('expect magnitudes to occur at same time, because'));
disp(sprintf('fan is pinned to ground at D and the linear acceleration'));
disp(sprintf('of the center of mass is a multiple of the angular'));
disp(sprintf('acceleration. \n\n'));

% find maximum magnitudes of angular and linear acceleration of fan
[max_ag_fan,i_max_ag_fan] = max(mag_ag_fan);
[max_thddot_fan,i_max_thddot_fan] = max(mag_thddot_fan);

disp(sprintf('approximate max magnitudes of linear and angular'));
disp(sprintf('  accelerations are, respectively:\n%6.2f m/s^2 \n%6.2f rad/s^2\n', ...
               max_ag_fan, max_thddot_fan));

if (i_max_ag_fan == i_max_thddot_fan);
 disp(sprintf('both maxima occur roughly %6.2f seconds into the cycle\n', ...
         dtime*(i_max_ag_fan-1)));
else;
 disp(sprintf('something is wrong, should occur at same time\n'));
end;

