ode45 problems -- not getting expected solution
9 views (last 30 days)
Show older comments
Rick Sellers
on 22 Jul 2020
Edited: David Goodmanson
on 24 Jul 2020
Hello,
I have a set of nonlinear satellite attitude equations to solve. I've implemented this as I've done in the past (at least I think I've done everything the same), and the solutions I get are linear. There should be an oscillation in each of the angles phi, tht, psi over time. My suspicion is that I didn't define the vector correctly being sent to ode45, my ode function doesn't update the contents of this vector correctly or there's some sort of scaling issue. I've asked other classmates already, the TA for my class, and stared at it for hours without seeing the error, so I'm turning to the experts for assistance.
Thank you,
Rick
% Derive Equations of Motion
clc; clear;
% PRY = [phi,tht,psi,dphi,dtht,dpsi]';
% Create sym variables
syms Sphi(t) Stht(t) Spsi(t) dSphi(t) dStht(t) dSpsi(t) ddSphi(t) ddStht(t)...
ddSpsi(t)
assume([Sphi(t) Stht(t) Spsi(t) dSphi(t) dStht(t) dSpsi(t) ddSphi(t) ...
ddStht(t) ddSpsi(t)],'real')
syms SOmega SR SG SmE SIxx SIyy SIzz real positive
% Form rotation matrices and cosine direction matrices
Rxphi = [1,0,0;0,cos(Sphi),sin(Sphi);0,-sin(Sphi),cos(Sphi)];
Rytht(t) = [cos(Stht),0,-sin(Stht);0,1,0;sin(Stht),0,cos(Stht)];
Rzpsi(t) = [cos(Spsi),sin(Spsi),0;-sin(Spsi),cos(Spsi),0;0,0,1];
BCA = Rxphi*Rytht*Rzpsi;
ACB = transpose(BCA);
% Compute O_omega_A in B-frame
OwA_B = BCA*[0;0;SOmega];
% Compute A_omega_B in B-frame
AwBcross = BCA*diff(ACB,t);
AwB_Bx = [0,0,1]*AwBcross*[0;1;0];
AwB_By = [1,0,0]*AwBcross*[0;0;1];
AwB_Bz = [0,1,0]*AwBcross*[1;0;0];
AwB_B = [ AwB_Bx;AwB_By;AwB_Bz ];
% Compute O_omega_B in B-frame
OwB_B = OwA_B + AwB_B;
% Decompose O_omega_B
OwB = eye(3)*OwB_B(t);
OwB = [OwB(1,:);OwB(2,:);OwB(3,:)];
% Moment of Inertia Matrix
I = [SIxx,0,0;0,SIyy,0;0,0,SIzz];
% Angular Momentum
h = I*OwB;
hdot = diff(h,t);
% Torques
K = 3*SG*SmE/SR^5;
Gx = K*(sin(2*Sphi)*(cos(Stht))^2*(SIzz-SIyy));
Gy = K*(sin(2*Stht)*cos(Sphi)*(SIzz-SIxx));
Gz = K*(sin(2*Stht)*sin(Sphi)*(SIxx-SIyy));
Tau = hdot + cross(OwB,h);
TauX = Tau(1,:);
TauY = Tau(2,:);
TauZ = Tau(3,:);
% Substitute
TauX = subs(TauX,[diff(Sphi(t), t, t),diff(Stht(t), t, t),diff(Spsi(t), t, t)],[ddSphi,ddStht,ddSpsi]);
TauY = subs(TauY,[diff(Sphi(t), t, t),diff(Stht(t), t, t),diff(Spsi(t), t, t)],[ddSphi,ddStht,ddSpsi]);
TauZ = subs(TauZ,[diff(Sphi(t), t, t),diff(Stht(t), t, t),diff(Spsi(t), t, t)],[ddSphi,ddStht,ddSpsi]);
TauX = subs(TauX,[diff(Sphi(t), t),diff(Stht(t), t),diff(Spsi(t), t)],[dSphi,dStht,dSpsi]);
TauY = subs(TauY,[diff(Sphi(t), t),diff(Stht(t), t),diff(Spsi(t), t)],[dSphi,dStht,dSpsi]);
TauZ = subs(TauZ,[diff(Sphi(t), t),diff(Stht(t), t),diff(Spsi(t), t)],[dSphi,dStht,dSpsi]);
TauX = subs(TauX,[Sphi(t),Stht(t),Spsi(t)],[Sphi,Stht,Spsi]);
TauY = subs(TauY,[Sphi(t),Stht(t),Spsi(t)],[Sphi,Stht,Spsi]);
TauZ = subs(TauZ,[Sphi(t),Stht(t),Spsi(t)],[Sphi,Stht,Spsi]);
% Tx = TauX == Gx;
% Ty = TauY == Gy;
% Tz = TauZ == Gz;
Tx = TauX == 0;
Ty = TauY == 0;
Tz = TauZ == 0;
isolate(Tx,ddSphi)
isolate(Ty,ddStht)
isolate(Tz,ddSpsi)
% Main
% Sets intitial values, creates PRY vector of unknowns and paramters; calls
% other functions and solver; sorts data and plots results
clc; clear; close all
% Set initial values
PRY = zeros(6,1);
phi = 0; tht = 0; psi = 0;
dphi = 0; dtht = 0; dpsi = 0;
RE = 6378.137e3 ; % m - radius of earth
RS = 500e3; % m - orbital altitude
R = (RE+RS); % m
mE = 5.972e24; % kg
G = 6.6743e-11; % m^3/kg s^2
Omega = sqrt(G*mE/R^3); % rad/s
Ixx = 6; Iyy = 8; Izz = 4; % kg m^2
mu = G*mE;
per = 2*pi*(sqrt(R^3/mu));
tend = 1.6*per;
tspan = [0 tend];
ic = [1;5;10;15;20;25];
for k= 1:6
phi =ic(k);
tht=ic(k);
psi=ic(k);
PRY = [phi,tht,psi,dphi,dtht,dpsi]';
[t,prydot] = ode45(@(t,pry) proj_ODE(t,PRY,Omega,R,mE,G,Ixx,Iyy,Izz),tspan,PRY);
% [t,prydot] = ode45(@(t,pry) proj_ODE_L(t,PRY,Omega,R,mE,G,Ixx,Iyy,Izz),tspan,PRY);
% Plotting results; pitch, roll, yaw for each initial condition
xticks = linspace(0,1.6,9);
figure(k)
% prydot = rad2deg(prydot);
subplot(3,1,1)
plot(t,prydot(:,1))
grid on
title('Pitch, \phi');
xlabel('Time, [s]'); ylabel('');
subplot(3,1,2)
plot(t,prydot(:,2));
grid on
title('Roll, \theta');
xlabel('Time, [s]');
subplot(3,1,3)
plot(t,prydot(:,3));
grid on
title('Yaw, \psi');
xlabel('Time, [s]'); ylabel('Attitude, [rad]');
end
% pry_ODE.m
% myODE function for solver
% PRY = [phi,tht,psi,dphi,dtht,dpsi];
function prydot = proj_ODE(t,pry,Omega,R,mE,G,Ixx,Iyy,Izz)
prydot = zeros(size(pry));
phi = pry(1);
tht = pry(2);
psi = pry(3);
dphi = pry(4);
dtht = pry(5);
dpsi = pry(6);
K = 3*G*mE/R^5;
w0 = Omega;
Gx = K*(sin(2*phi)*(cos(tht))^2*(Izz-Iyy));
Gy = K*(sin(2*tht)*cos(phi)*(Izz-Ixx));
Gz = K*(sin(2*tht)*sin(phi)*(Ixx-Iyy));
prydot(1) = dphi;
prydot(2) = dtht;
prydot(3) = dpsi;
prydot(4) = (Ixx*((sin(phi)*sin(psi) ...
+ cos(phi)*cos(psi)*sin(tht))*(cos(phi)*cos(psi)*prydot(6) ...
- cos(phi)*sin(psi)*dphi^2 - cos(phi)*sin(psi)*dpsi^2 ...
+ cos(psi)*sin(phi)*sin(tht)*dphi^2 + cos(psi)*sin(phi)*sin(tht)*dpsi^2 ...
+ cos(psi)*sin(phi)*sin(tht)*dtht^2 - 2*cos(psi)*sin(phi)*dphi*dpsi ...
- cos(psi)*cos(tht)*sin(phi)*prydot(5) + sin(phi)*sin(psi)*sin(tht)*prydot(6) ...
- 2*cos(phi)*cos(psi)*cos(tht)*dphi*dtht ...
+ 2*cos(phi)*sin(psi)*sin(tht)*dphi*dpsi ...
+ 2*cos(tht)*sin(phi)*sin(psi)*dpsi*dtht) - (cos(psi)*sin(phi) ...
- cos(phi)*sin(psi)*sin(tht))*(cos(phi)*sin(psi)*prydot(6) ...
+ cos(phi)*cos(psi)*dphi^2 + cos(phi)*cos(psi)*dpsi^2 ...
+ sin(phi)*sin(psi)*sin(tht)*dphi^2 + sin(phi)*sin(psi)*sin(tht)*dpsi^2 ...
+ sin(phi)*sin(psi)*sin(tht)*dtht^2 - 2*sin(phi)*sin(psi)*dphi*dpsi ...
- cos(psi)*sin(phi)*sin(tht)*prydot(6) - cos(tht)*sin(phi)*sin(psi)*prydot(5) ...
- 2*cos(phi)*cos(psi)*sin(tht)*dphi*dpsi - 2*cos(phi)*cos(tht)*sin(psi)*dphi*dtht ...
- 2*cos(psi)*cos(tht)*sin(phi)*dpsi*dtht) - (cos(phi)*sin(psi)*dphi ...
+ cos(psi)*sin(phi)*dpsi + cos(phi)*cos(psi)*cos(tht)*dtht ...
- cos(psi)*sin(phi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*sin(tht)*dpsi)*(sin(phi)*sin(psi)*dphi ...
- cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
+ cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
- (sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
+ cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
- sin(phi)*sin(psi)*sin(tht)*dphi)*(cos(phi)*sin(psi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
+ cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
+ cos(phi)*cos(tht)*(sin(phi)*sin(tht)*prydot(5) ...
+ cos(tht)*sin(phi)*dphi^2 + cos(tht)*sin(phi)*dtht^2 ...
+ 2*cos(phi)*sin(tht)*dphi*dtht) + w0*cos(tht)*dtht ...
+ cos(tht)*sin(phi)*dphi*(cos(phi)*cos(tht)*dphi - sin(phi)*sin(tht)*dtht) ...
+ cos(phi)*sin(tht)*dtht*(cos(phi)*cos(tht)*dphi - sin(phi)*sin(tht)*dtht)) ...
+ Iyy*(sin(tht)*(cos(tht)*sin(phi)*dphi + cos(phi)*sin(tht)*dtht) ...
+ w0*cos(tht)*sin(phi) + cos(psi)*cos(tht)*(cos(phi)*sin(psi)*dphi ...
+ cos(psi)*sin(phi)*dpsi + cos(phi)*cos(psi)*cos(tht)*dtht ...
- cos(psi)*sin(phi)*sin(tht)*dphi - cos(phi)*sin(psi)*sin(tht)*dpsi) ...
+ cos(tht)*sin(psi)*(sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
+ cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
- sin(phi)*sin(psi)*sin(tht)*dphi))*((cos(phi)*cos(psi) ...
+ sin(phi)*sin(psi)*sin(tht))*(cos(psi)*cos(tht)*dpsi ...
- sin(psi)*sin(tht)*dtht) + (cos(phi)*sin(psi) ...
- cos(psi)*sin(phi)*sin(tht))*(cos(tht)*sin(psi)*dpsi ...
+ cos(psi)*sin(tht)*dtht) + w0*cos(phi)*cos(tht) - cos(tht)^2*sin(phi)*dtht) ...
- Izz*(sin(tht)*(cos(tht)*sin(phi)*dphi + cos(phi)*sin(tht)*dtht) ...
+ w0*cos(tht)*sin(phi) + cos(psi)*cos(tht)*(cos(phi)*sin(psi)*dphi ...
+ cos(psi)*sin(phi)*dpsi + cos(phi)*cos(psi)*cos(tht)*dtht ...
- cos(psi)*sin(phi)*sin(tht)*dphi - cos(phi)*sin(psi)*sin(tht)*dpsi) ...
+ cos(tht)*sin(psi)*(sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
+ cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
- sin(phi)*sin(psi)*sin(tht)*dphi))*((cos(phi)*cos(psi) ...
+ sin(phi)*sin(psi)*sin(tht))*(cos(psi)*cos(tht)*dpsi ...
- sin(psi)*sin(tht)*dtht) + (cos(phi)*sin(psi) ...
- cos(psi)*sin(phi)*sin(tht))*(cos(tht)*sin(psi)*dpsi ...
+ cos(psi)*sin(tht)*dtht) + w0*cos(phi)*cos(tht) ...
- cos(tht)^2*sin(phi)*dtht))/(Ixx*(cos(phi)^2*cos(tht)^2 + (sin(phi)*sin(psi) ...
+ cos(phi)*cos(psi)*sin(tht))^2 + (cos(psi)*sin(phi) ...
- cos(phi)*sin(psi)*sin(tht))^2));
prydot(5) = (Iyy*(cos(tht)*sin(psi)*(cos(phi)*cos(psi)*prydot(4) ...
- sin(phi)*sin(psi)*prydot(6) - cos(psi)*sin(phi)*dphi^2 ...
- cos(psi)*sin(phi)*dpsi^2 + cos(phi)*sin(psi)*sin(tht)*dphi^2 ...
+ cos(phi)*sin(psi)*sin(tht)*dpsi^2 + cos(phi)*sin(psi)*sin(tht)*dtht^2 ...
- 2*cos(phi)*sin(psi)*dphi*dpsi - cos(phi)*cos(psi)*sin(tht)*prydot(6) ...
+ sin(phi)*sin(psi)*sin(tht)*prydot(4) ...
- 2*cos(phi)*cos(psi)*cos(tht)*dpsi*dtht ...
+ 2*cos(psi)*sin(phi)*sin(tht)*dphi*dpsi ...
+ 2*cos(tht)*sin(phi)*sin(psi)*dphi*dtht) ...
- sin(tht)*(cos(tht)*sin(phi)*prydot(4) + cos(phi)*cos(tht)*dphi^2 ...
+ cos(phi)*cos(tht)*dtht^2 - 2*sin(phi)*sin(tht)*dphi*dtht) ...
+ cos(psi)*cos(tht)*(sin(phi)*sin(psi)*dphi^2 + sin(phi)*sin(psi)*dpsi^2 ...
- cos(phi)*sin(psi)*prydot(4) - cos(psi)*sin(phi)*prydot(6) ...
+ cos(phi)*cos(psi)*sin(tht)*dphi^2 + cos(phi)*cos(psi)*sin(tht)*dpsi^2 ...
+ cos(phi)*cos(psi)*sin(tht)*dtht^2 - 2*cos(phi)*cos(psi)*dphi*dpsi ...
+ cos(psi)*sin(phi)*sin(tht)*prydot(4) ...
+ cos(phi)*sin(psi)*sin(tht)*prydot(6) ...
+ 2*cos(psi)*cos(tht)*sin(phi)*dphi*dtht ...
+ 2*cos(phi)*cos(tht)*sin(psi)*dpsi*dtht ...
- 2*sin(phi)*sin(psi)*sin(tht)*dphi*dpsi) ...
- cos(tht)*dtht*(cos(tht)*sin(phi)*dphi + cos(phi)*sin(tht)*dtht) ...
+ w0*sin(phi)*sin(tht)*dtht - cos(psi)*cos(tht)*dpsi*(sin(phi)*sin(psi)*dpsi ...
- cos(phi)*cos(psi)*dphi + cos(phi)*cos(psi)*sin(tht)*dpsi ...
+ cos(phi)*cos(tht)*sin(psi)*dtht - sin(phi)*sin(psi)*sin(tht)*dphi) ...
+ cos(tht)*sin(psi)*dpsi*(cos(phi)*sin(psi)*dphi + cos(psi)*sin(phi)*dpsi ...
+ cos(phi)*cos(psi)*cos(tht)*dtht - cos(psi)*sin(phi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*sin(tht)*dpsi) ...
+ cos(psi)*sin(tht)*dtht*(cos(phi)*sin(psi)*dphi ...
+ cos(psi)*sin(phi)*dpsi + cos(phi)*cos(psi)*cos(tht)*dtht ...
- cos(psi)*sin(phi)*sin(tht)*dphi - cos(phi)*sin(psi)*sin(tht)*dpsi) ...
+ sin(psi)*sin(tht)*dtht*(sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
+ cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
- sin(phi)*sin(psi)*sin(tht)*dphi) - w0*cos(phi)*cos(tht)*dphi) ...
- Ixx*((cos(phi)*cos(psi) + sin(phi)*sin(psi)*sin(tht))*(cos(psi)*cos(tht)*dpsi ...
- sin(psi)*sin(tht)*dtht) + (cos(phi)*sin(psi) ...
- cos(psi)*sin(phi)*sin(tht))*(cos(tht)*sin(psi)*dpsi ...
+ cos(psi)*sin(tht)*dtht) + w0*cos(phi)*cos(tht) ...
- cos(tht)^2*sin(phi)*dtht)*((sin(phi)*sin(psi) ...
+ cos(phi)*cos(psi)*sin(tht))*(sin(phi)*sin(psi)*dphi ...
- cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
+ cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
- (cos(psi)*sin(phi) ...
- cos(phi)*sin(psi)*sin(tht))*(cos(phi)*sin(psi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
+ cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
- w0*sin(tht) + cos(phi)*cos(tht)*(cos(phi)*cos(tht)*dphi ...
- sin(phi)*sin(tht)*dtht)) + Izz*((cos(phi)*cos(psi) ...
+ sin(phi)*sin(psi)*sin(tht))*(cos(psi)*cos(tht)*dpsi ...
- sin(psi)*sin(tht)*dtht) + (cos(phi)*sin(psi) ...
- cos(psi)*sin(phi)*sin(tht))*(cos(tht)*sin(psi)*dpsi ...
+ cos(psi)*sin(tht)*dtht) + w0*cos(phi)*cos(tht) ...
- cos(tht)^2*sin(phi)*dtht)*((sin(phi)*sin(psi) ...
+ cos(phi)*cos(psi)*sin(tht))*(sin(phi)*sin(psi)*dphi ...
- cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
+ cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
- (cos(psi)*sin(phi) ...
- cos(phi)*sin(psi)*sin(tht))*(cos(phi)*sin(psi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
+ cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
- w0*sin(tht) + cos(phi)*cos(tht)*(cos(phi)*cos(tht)*dphi ...
- sin(phi)*sin(tht)*dtht)))/(Iyy*(cos(phi)*cos(psi)^2*cos(tht)^2 ...
+ cos(phi)*cos(tht)^2*sin(psi)^2 + cos(phi)*sin(tht)^2));
prydot(6) = (Izz*((cos(phi)*cos(psi) ...
+ sin(phi)*sin(psi)*sin(tht))*(sin(psi)*sin(tht)*prydot(5) ...
+ cos(tht)*sin(psi)*dpsi^2 + cos(tht)*sin(psi)*dtht^2 ...
+ 2*cos(psi)*sin(tht)*dpsi*dtht) + (cos(tht)*sin(psi)*dpsi ...
+ cos(psi)*sin(tht)*dtht)*(sin(phi)*sin(psi)*dphi ...
- cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
+ cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
- (cos(psi)*cos(tht)*dpsi ...
- sin(psi)*sin(tht)*dtht)*(cos(phi)*sin(psi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
+ cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
- (cos(phi)*sin(psi) ...
- cos(psi)*sin(phi)*sin(tht))*(cos(psi)*sin(tht)*prydot(5) ...
+ cos(psi)*cos(tht)*dpsi^2 + cos(psi)*cos(tht)*dtht^2 ...
- 2*sin(psi)*sin(tht)*dpsi*dtht) + cos(tht)^2*sin(phi)*prydot(5) ...
- 2*cos(tht)*sin(phi)*sin(tht)*dtht^2 + cos(phi)*cos(tht)^2*dphi*dtht ...
+ w0*cos(tht)*sin(phi)*dphi + w0*cos(phi)*sin(tht)*dtht) ...
+ Ixx*(sin(tht)*(cos(tht)*sin(phi)*dphi + cos(phi)*sin(tht)*dtht) ...
+ w0*cos(tht)*sin(phi) + cos(psi)*cos(tht)*(cos(phi)*sin(psi)*dphi ...
+ cos(psi)*sin(phi)*dpsi + cos(phi)*cos(psi)*cos(tht)*dtht ...
- cos(psi)*sin(phi)*sin(tht)*dphi - cos(phi)*sin(psi)*sin(tht)*dpsi) ...
+ cos(tht)*sin(psi)*(sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
+ cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
- sin(phi)*sin(psi)*sin(tht)*dphi))*((sin(phi)*sin(psi) ...
+ cos(phi)*cos(psi)*sin(tht))*(sin(phi)*sin(psi)*dphi ...
- cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
+ cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
- (cos(psi)*sin(phi) ...
- cos(phi)*sin(psi)*sin(tht))*(cos(phi)*sin(psi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
+ cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
- w0*sin(tht) + cos(phi)*cos(tht)*(cos(phi)*cos(tht)*dphi ...
- sin(phi)*sin(tht)*dtht)) - Iyy*(sin(tht)*(cos(tht)*sin(phi)*dphi ...
+ cos(phi)*sin(tht)*dtht) + w0*cos(tht)*sin(phi) ...
+ cos(psi)*cos(tht)*(cos(phi)*sin(psi)*dphi + cos(psi)*sin(phi)*dpsi ...
+ cos(phi)*cos(psi)*cos(tht)*dtht - cos(psi)*sin(phi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*sin(tht)*dpsi) ...
+ cos(tht)*sin(psi)*(sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
+ cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
- sin(phi)*sin(psi)*sin(tht)*dphi))*((sin(phi)*sin(psi) ...
+ cos(phi)*cos(psi)*sin(tht))*(sin(phi)*sin(psi)*dphi ...
- cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
+ cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
- (cos(psi)*sin(phi) ...
- cos(phi)*sin(psi)*sin(tht))*(cos(phi)*sin(psi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
+ cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
- w0*sin(tht) + cos(phi)*cos(tht)*(cos(phi)*cos(tht)*dphi ...
- sin(phi)*sin(tht)*dtht)))/(Izz*(cos(psi)*cos(tht)*(cos(phi)*cos(psi) ...
+ sin(phi)*sin(psi)*sin(tht)) + cos(tht)*sin(psi)*(cos(phi)*sin(psi) ...
- cos(psi)*sin(phi)*sin(tht))));
end
15 Comments
David Goodmanson
on 24 Jul 2020
Edited: David Goodmanson
on 24 Jul 2020
Hi RIck,
It's good that you got it working. I feel that the code would be better (certainly way more compact) without using syms at all, but if it ain't broke ... One thing I didn't mention is that in your calculation of the prydots, prydot(4) is calculated in terms of the old prydots, but prydot(5) is calculated with the new prydot(4), and prydot(6) is calculated with new prydot(4) and new prydot(5). It's more consistent to calculate the new ones strictly in terms of the old ones. This might make only a small difference numerically, or it could have a measurable effect.
Accepted Answer
Steven Lord
on 22 Jul 2020
K = 3*G*mE/R^5;
It's been a while since I've studied any physics, but this looks odd to me.
RE = 6378.137e3 ; % m - radius of earth
RS = 500e3; % m - orbital altitude
R = (RE+RS); % m
mE = 5.972e24; % kg
G = 6.6743e-11; % m^3/kg s^2
By dimensional analysis, K has units of or times whatever the unit of 3 is. If that's supposed to be the gravitational force between a 3 kg satellite and the Earth, shouldn't that be over R^2 (which would give Newtons, ) instead of over R^5?
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!