% Erich Bender
% AERO 463/464 -- Senior Project
close all
clear all
clc
% Mass Moments of Inertia and Stability Chart
% Reference frame: The moments of inertia will be calculated along a body
% frame with:
% - x pointing in the RAM direction
% - y pointing in the cross-track direction to the right of the RAM
% - z pointing in the Nadir direction
%Tmass = .40; % tip mass (kg)
Tmass = .20; % tip mass (kg)
Theight = .010; % height of cylindrical tip mass (m)
Tradius = .005; % radius of cylindrical tip mass (m)
Rmass = .010; % mass of boom (kg)
%Rlength = .34775; % length of boom (m)
Rlength = .3; % length of boom (m)
Bheight = .10; % height of rectangular body (m)
Bwidth = .3405; % width of rectangular body (m)
Bdepth = .10; % depth of rectangular body (m)
Rtip1 = [0; 0; -.3]; % vector describing position (x,y,z) of tip mass 1 from C.M. to component in m
Rtip2 = [0; 0; .3]; % vector describing position (x,y,z) of tip mass 2 from C.M. to component in m
Rrod1 = [0; 0; -.105125]; % vector describing position (x,y,z) of deployment device 1 from C.M. to component in m
Rrod2 = [0; 0; .105125]; % vector describing position (x,y,z) of deployment device 2 from C.M. to component in m
global Ix;
global Iy;
global Iz;
%caso A
%[Ix, Iy, Iz] = massmoment3U(Tmass, Tradius, Theight, Rmass, Rlength, Bheight, Bwidth, Bdepth, Rtip1, Rtip2, Rrod1, Rrod2)
%caso B
[Ix, Iy, Iz] = massmoment3U1boom(Tmass, Tradius, Theight, Rmass, Rlength, Bheight, Bwidth, Bdepth, Rtip1, 0, Rrod1, 0)
%caso C
%[Ix, Iy, Iz] = massmoment3Unoboom(0, Tradius, Theight, 0, Rlength, Bheight, Bwidth, Bdepth, 0, 0, 0, 0)
global K1
K1 = ((Ix) - (Iy))/(Iz);
% K1 = -.6;
global K2
K2 = ((Iy) - (Iz))/(Ix);
% K2 = .8824;
global K3
K3 = -((K1 + K2)/(1 +K1*K2));
% Stability Chart
x1=linspace(-1,1,1000);
for i=1:length(x1)
y1(i)=-x1(i);
end
figure(1)
plot(x1,y1, K1, K2, 'xr','MarkerSize',10)
axis([-1 1 -1 1])
grid off
hold on
line([-1,1],[0,0])
line([0,0],[-1,1])
hold off
title('Stability Chart')
xlabel('K1')
ylabel('K2')
% Orbital Parameters
R = 6378 + 500; %km
mu = 398600; %km^3/s^2
T = ((2*pi)/sqrt(mu))*(R^(3/2)); % period in seconds
omega = sqrt(mu/(R^3)); % spin of the earth m/s
% Initial Position (Euler Angles)
phi=0;
theta=0;
psi=0;
[e1i, e2i, e3i, e4i, a, E] = dcm2quat(phi, theta, psi);
% Time Interval
ti = 0;
tf = 16*T; %propagate over about 24 hours
% Inertia Properties of Gyro/Reaction Wheel
m = .120; % kg
r = .050; % m
global J
J = (m*r^2)/2;
% Gravity Gradient Simulation
x = [.01*omega, .01*omega, omega, e1i, e2i, e3i, e4i];
options=odeset('RelTol',1e-10);
[t,y]=ode45('simulation',[ti,tf],x,options);
for i = 1:length(t)
e1n = y(i,4);
e2n = y(i,5);
e3n = y(i,6);
e4n = y(i,7);
a(i) = sqrt(((y(i,4))^2)+((y(i,5))^2)+((y(i,6))^2)+((y(i,7)^2)));
E = [(1 - 2*(e2n^2+e3n^2)),2*(e1n*e2n + e3n*e4n),2*(e1n*e3n - e2n*e4n);2*(e1n*e2n - e3n*e4n),(1 - 2*(e1n^2+e3n^2)),2*(e2n*e3n + e1n*e4n);2*(e1n*e3n + e2n*e4n),2*(e2n*e3n - e1n*e4n),(1 - 2*(e1n^2+e2n^2))];
yaw(i,1)=acos(E(3,3))*180/pi;
pitch(i,1)=acos(E(2,2))*180/pi;
roll(i,1)=acos(E(1,1))*180/pi;
end
figure(2)
plot(t, yaw)
title('Yaw (RAM)')
xlabel('Time (s)')
ylabel('(deg)')
figure(3)
plot(t, pitch)
title('Pitch (Crosstrack)')
xlabel('Time (s)')
ylabel('(deg)')
figure(4)
plot(t, roll)
title('Roll (Nadir)')
xlabel('Time (s)')
ylabel('(deg)')
figure(5)
plot(t,a)
axis([0 tf .99 1.01])
title('Validity Check For Quaternions')
xlabel('Time (s)')
ylabel('Sum of the Squares')
评论0