
% 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
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);

p_num = zeros(nt,2,4);
for it = 1:nt;
	p_num(it,:,:) = linkage_pos(th_num(it,:),r);
end;

ifig = ifig+1; figure(ifig);
plot(...
	p_num(:,1,1),p_num(:,2,1),'*', ...
	p_num(:,1,2),p_num(:,2,2), ...
	p_num(:,1,3),p_num(:,2,3), ...
	p_num(:,1,4),p_num(:,2,4),'*')
axis('equal');
xlabel('horizontal position'); ylabel('vertical position');
title('positions of pins, from numerical integration of angular velocities');

% get approximate maxiumum angle of tilt of fan blades from vertical
max_th3_num = max(th_num(:,3));
max_tilt_num = max_th3_num + pi/2;
disp(sprintf('max tilt of blades from vertical degrees is about %3.0f degrees\n', ...
	     max_tilt_num*180/pi));
% 

