Biturbowned
Triple Kill!
Joined: 29 Nov 2006
Posts: 266
Location: here 'n there
|
 Posted: Fri Feb 16, 2007 10:21 pm Post subject: why the fvck to Mechanical Engineers need to know this crap? |
 |
|
%develop Matlab program which determines classical orbital elements from
%(i) position and velocity (i.e. R and V) and (ii) 3 position measurements
%R1, R2, R3 (Gibbs method). Use algorithm to determine orbital elements
%for cases listed:
%Prob 2.2
%Prob 3.2a
%Prob 3.2b
%Prob 3.4
%Prob 2.15a
%--------------------------------------------------------------------------
function[]=main()
clc; clear all;
disp('**Astrodynamics HW5, problem 4-8**')
disp('Coded by Chad Hodgkins')
A=input('If input is in form r and v, press "i", if input is in form R1, R2, R3, press "ii" ','s')
if strcmp(A,'i')==1
astro_i()
% disp('future program here')
elseif strcmp(A,'ii')==1
astro_ii()
else
disp('Error: input not valid')
% main()
end
end
%--------------------------------------------------------------------------
%part (i)
function[]=astro_i()
clear all;
disp('You selected part i')
r=input('Input value for r in form of "[x,y,z]"')
v=input('Input value for v in form of "[x,y,z]"')
mu=input('Input value of mu') %request user preference for units
i=[1,0,0]; %definition of inertial coordinate frame
k=[0,0,1];
%for i:
h=cross(r,v); %definition of angular momentum
i=acosd(h(1,3)/norm(h))
%for e:
e=(cross(v,h)-(mu*r/norm(r)))/mu;
e_norm=norm(e)
%for w:
n=(cross(k,h));
w=acosd(dot(n,e)/(norm(n)*e_norm))
%for f0:
f0=acosd(dot(r,e)/(norm(r)*e_norm))
%for omega
omega=acosd(dot(n,i)/norm(n))
%for P:
P=norm(h)^2/mu
%for a:
a=P/(1-(e_norm)^2)
end
%--------------------------------------------------------------------------
function [] = astro_ii()
clear all;
disp('You selected part ii')
mu=input('Input value of mu ') %request user preference for units
% R1=input('Input vector R1 in form of "[x,y,z]"')
% R2=input('Input vector R2 in form of "[x,y,z]"')
% R3=input('Input vector R3 in form of "[x,y,z]"')
R1=[1.41422511,0,1.414202];
R2=[1.81065659,1.06066883,.3106515];
R3=[1.35353995,1.41422511,-0.6464495];
pass=true;
k=[0,0,1];
%check for coplanarity:
R4=dot(cross(R2,R3),R1)
N=cross(R1,R2);
n=cross(k,N);
if (R4>-.0001 && R4<.0001) %Vectors lie on the same plane
C1=dot(cross(R3,R2),N)/(norm(N))^2;
C2=dot(cross(R1,R3),N)/(norm(N))^2;
R_1=R1*C1;
R_2=R2*C2;
P=(norm(R3)-norm(R_1)-norm(R_2))/(1-C1-C2)
N_cross_e=-(R1*(P-norm(R2))-R2*(P-norm(R1)));
e=cross(N_cross_e,N)/(norm(N))^2;
e_norm=norm(e)
a=P/(1-e_norm^2);
i=acosd(dot(N,k)/norm(N))
w=dot(n,e)/(norm(n)*e_norm)
omega=atand(norm(N)/N(1,2))
else %Vectors do not lie on the same plane
pass=false;
end
%error handling:
if pass==false;
disp('ERROR: R1, R2, R3 are not coplanar, you will return to menu')
main()
end
end
omghey
_________________
PS3::XB360 |
|
