ode45 problems -- not getting expected solution

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

[Alan Stevens]

You calculate Gx, Gy and Gz inside your function but they don't seem to be used anywhere. Should they be?

[Roger J]

Are you sure you want to plot these:

plot(t,prydot(:,1))
plot(t,prydot(:,2))
plot(t,prydot(:,3))

or would you rather see:

plot(t,prydot(:,4))
plot(t,prydot(:,5))
plot(t,prydot(:,6))

[Rick Sellers]

Roger,

I'm thinking the output of ode45 will be the integration of my components in prydot. e.g. prydot(1) = dphi/dt, but after the solver this should be phi(t), which is what I'm looking for. prydot(4),(5),(6), should be the first derivatives. However, it's quite possible I'm mistaken in the output of the solver.

Thanks,

Rick

[James Tursa]

I don't think anyone is going to be able to decipher your derivative equations to help you. You might back up a step and give us the differential equations you are solving and the symbolic code you used to generate your derivatives.

[Rick Sellers]

James,

I understand. I guess I was looking for another set of eyes to spot anything obviously wrong. I'll upload the symbolic portion I did. The basic process is:

Find angular velocity of satellite body reference frame, O_omega_B

Compute angular momentum -- h = I*O_omega_B

The gravity gradient terms were derived outside of matlab and I just typed those in, but the supply the external torques on the satellite

These torques are equal to dh/dt + O_omega_B X h

This last step is how I obtained the mess of equation in the above proj_ODE function, then isolated the second derivatives, and copied and pasted over into my ODE function.

Hopefully that makes some sense.

THank you

Rick

[Alan Stevens]

In your function proj_ODE prydot(4) contains a few mentions of prydot(6) and prydot(5). Since you zero these on entry to the function they will be zero. Similarly comments apply for prydot(5).

[Rick Sellers]

Alan,

I agree this is the problem. My understanding is the ode45 solver will update these values as it computes new approximations, and those values will get populated after the first iteration step. However, I'm not 100% sure about that, and I still may have coded it incorrectly. If you have any ideas about how to incorporate the fact that my other second deriatives appear in each equation, that would be great.

If I have this, how do I tell it to use the last approximation for in computing (as an example)

= cos(ϕ) + *sin(θ)

= *cos(ψ) - cos(θ)sin(θ)*

I thought I could say, if my variables are y = [phi,tht,psi, phi', tht', psi'], then my ode function would be something like

dydt = myode(t,y)

dydt(1) = y(4)

dydt(2) = y(5)

...

dydt(4) = cos(y(1)) + dydy(6)*sin(y(2))

dydt(5) = dydt(4)*cos(y(3)) - cos(y(3))*sin(y(3))*dydt(6)

Is this the right idea? Will this work?

Thank you,

Rick

[Steven Lord]

If I may make a suggestion related to James Tursa's comment, those expressions for prydot(4), prydot(5), and prydot(6) are a bit wall of text-like. [It's a common problem; whenever I write up PowerPoint slides I have to force myself not to write a novel on each one. I also sometimes have that problem in Answers.] Consider putting comments after the ellipsis (...) on lines indicating the purpose of that term and/or defining intermediate variables (with descriptive names) for specific parts of those calculations and then combining those variables.

Splitting your equations for prydot(4), prydot(5), and prydot(6) would require some work (perhaps a lot of work) but the act of splitting them and determining the correct names for each term may help expose where a particular term isn't quite right. It will also help with dimensional analysis; if one term is supposed to be a torque but instead it turns out to be a force, it should be easier to detect that in the smaller piece of code.

[Rick Sellers]

Steven,

Good tips. I understnad it's nearly impossible to read on here. I may try defining some other "helper" expressions like you mentions. e.g. let A = cos(phi)*cos(tht)*(cos(phi)*cos(tht), then I would have a bunch of A*dphi - B*C*dtht, and so on. This should help my troubleshooting.

It seems like the way I've written everything it evaluates to phi'' = 0; etc.

It's pretty late for me to keep tinkering with it, so I'll probably have to submit what I have, highlight that it's incorrect, point out why, and hope for the best.

Thank you,

Rick

[Alan Stevens]

With regard to your comment: "My understanding is the ode45 solver will update these values as it computes new approximations, and those values will get populated after the first iteration step.". This won't happen here because you explicitly zero prydot at the start of each entry into the function proj_ODE.

[Rick Sellers]

Alan,

Aha. I just thought of the way forward, I think. I should have been sending it a separate set of values for the initial conditions.

THank you

[David Goodmanson]

Hi Rick,

I don't know if this is what you meant, but you need to change

[t,prydot] = ode45(@(t,pry) proj_ODE(t,PRY,Omega,R,mE,G,Ixx,Iyy,Izz),tspan,PRY);
to
[t,prydot] = ode45(@(t,pry) proj_ODE(t,pry,Omega,R,mE,G,Ixx,Iyy,Izz),tspan,PRY);
                                        ^

because initial conditions don't belong in the second argument. This runs, but unfortunately the system goes singular after a few hundred seconds.

[David Goodmanson]

Hi Rick,

It's not so clear how all these equations were arrived at. For example, the 60 lines of "derive equations of motion" code at the top can be commented out and make no difference. Apparently those are used in some way to provide the lines of code for prydot(4) etc. in the proj_ODE function. The prydot lines are typical symbolic variable bloat and undecipherable as they stand, which is all right as long as their source is clear.

Could you post all the equations that are the starting point? Something like (if torque is N)

I1*w1dot -w2*w3*(I2-I3) = N1     etc.

or similar, and the expressions for the torques, the motion of the satellite and whatever else is needed.

[Rick Sellers]

David,

I ended up getting it working. It doesn't quite match the plots from the paper we were given to re-create, but that paper was basically garbage, and so I'm happy with what I have now. As someone else pointed out, my return vector from the ODE function was zeroized at the beginning of the function. This, along with the problem you noted were the source of my issues.

Another comment had asked about the derivation of these equations, adn that's what the EOM code is about. You're correct, that was my back-end symbolic work to get these equations. They're very simply derived from 1) getting the angular velociy components of the satellite body frame wrt an inertial reference frame; 2) computing the angular momentum as h = [I] * omega, where I is the moment of inertia tensor; 3) the torques on the satellite are then given by {T} = d/dt(h) + omega X h; The LHS of the equations comes from gravity gradient moments (the 3 G m/R^3 terms).

I obviously have a very rudimentary understanding of this, but hope this shed some light on what the problem was.

'

Thank you for the help,

Rick

[David Goodmanson]

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.

Answer 1

[Steven Lord]

    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?

[Rick Sellers]

Steven,

Those terms are multiplied by the moments of inertia [kg * m^2], but that still leaves a unit of length^2 missing to make it Newton-m (these are torques). So, I'm thinking it should be over R^3. But, I just made this change, ran it, and got the same results.

Thank you

Rick