function [thdot,thddot,ag] = fan_solve(t, th)
%
% given angles of components, find:
% thdot  -- angular velocities
% thddot -- angular accelerations
% ag     -- linear acceleration of mass centers
%
global r r2 rg lthg;
global iwk iwu1 iwu2 ium;
global w_data;
global gc g_by_gc; 
global m ai;
%
nlink = 3; % number of moving links
if (length(th) ~= nlink);
      disp('error, th of wrong length');
      return;
end;
%
% get unit vectors:
% vectors nvec are along the links and vectors tvec are perpendicular to
% the links and in the direction of motion of the end of the link when 
% the angular velocity is positive
%
[nvec, tvec] = linkage_vec(th);
%
% solve for angular velocitiess
A = [r(iwu1)*tvec(:,iwu1) r(iwu2)*tvec(:,iwu2)];
y = [w_data(iwk)*r(iwk)*tvec(:,iwk)];
x = -A\y;
%
% set values in thdot
thdot = zeros(3,1);
thdot(iwk)  = w_data(iwk);
thdot(iwu1) = x(1);
thdot(iwu2) = x(2);
%
% 
if (nargout > 1);
% find angular accelerations
%
% get squares of angular velocities
thdot2 = thdot.*thdot;
%
A = [r(iwu1)*tvec(:,iwu1) r(iwu2)*tvec(:,iwu2)];
y = [-thdot2(iwk)* r(iwk)* nvec(:,iwk)  + ...
     -thdot2(iwu1)*r(iwu1)*nvec(:,iwu1) + ...
     -thdot2(iwu2)*r(iwu2)*nvec(:,iwu2)];
x = -A\y;
%
% set values in thddot
thddot = zeros(3,1);
thddot(iwk)  = 0.0;
thddot(iwu1) = x(1);
thddot(iwu2) = x(2);
%
end; % end of calculation of angular accelerations
%
%
if (nargout > 2);
% find linear acceleration of mass centers
%
[ag,ng] = linkage_ag(thdot2,thddot,r,rg,lthg,nvec,tvec);
%
end; % end calculation of linear acceleration of mass centers
%
%
