%% tilt-free seismometer primary mass
% estimate of pitch moment of inertia 
% KLD Apr. 30, 2015

%  ===========       <-- cuboid mass
%  |         |
%  |         |
%  |         |
%  |         |
%  |         |       <--> pitch axis about which system rotates
%  |         |
%  |         |
%  |         |
%  |         |
%  ===========       <-- cuboid mass

% parameters
M = 5;      % [kg] cuboid mass
h = 0.025;  % [m]  cuboid height
d = 0.30;   % [m]  cuboid depth
w = 0.25;   % [m]  cuboid width

m = 0.5;    % [kg] rod mass
l = 0.82;   % [m]  rod length

%% cuboid mass (one each at top and bottom)

% MOI about its own center of mass (axis parallel to system pitch axis)
I_cub_CoM = 1/12*M*(h^2+d^2);

% MOI about system pitch axis (apply parallel axis theorem)
I_cub = I_cub_CoM + M*(l/2+h/2)^2;



%% rod (one at each corner)

% MOI about its own center of mass (axis parallel to system pitch axis)
% infinitely thin approximation
I_rod_CoM = 1/12*m*l^2;

% MOI about system pitch axis (apply parallel axis theorem)
I_rod = I_rod_CoM + m*(d/2)^2;


%% total system moment of inertia
% sum of parts

I_pitch = 2*I_cub + 4*I_rod;
sprintf('Moment of inertia (pitch) = %2.1f kg m^2' , I_pitch)


%% Eigenfrequencies
% analytic equations from Malik's T000134

% parameters
m_tot = 2*M + 4*m;
g = 9.8;                % [m/s^2]
L = l/2 + 0.15;         % [m]       wire length to suspension point
b = 0.001;              % [w]       distance from CoM to wire attachment

w_pitch = sqrt(m_tot*g*b*(L+b)/(I_pitch*L));
%w_roll = sqrt(m_tot*g*b*(L+b)/(I_roll*L));        (???) two wires!
%w_yaw = sqrt(m_tot*g*R1*R2/(I_yaw*L));

sprintf('Natural frequency (pitch) = %1.3f Hz' , w_pitch/(2*pi))