%%
% ========================================================================
% © Ciara McGrath 20/11/2018
% 
% Equations to calculate change in RAAN and AoL after 3 phase maneuver  
% using low-thrust propulsion.
% From appendix A of PhD Thesis: Analytical Methods for Satellite
% Constellation Reconfiguration and Reconnaissance using Low-Thrust
% Manoeuvres
% 
% Note: The semi-major axis used in these equations is assumed to be the 
%       mean semi-major axis.
% 
% Note: For lowering maneuvers, a negative propulsive acceleration, A, 
%       should be used. For raising maneuvers, a positive acceleration, B, 
%       should be used.
%
% * DISCLAIMER *
%   The S/W remains property of Ciara McGrath and shall not be modified,
%   nor shall it be forwarded to third parties without prior written consent
%
% =======================================================================
%%
close all
clear all
clc

% ---------------------------
% Input Values
% ---------------------------

% Constants
mu =  3.986004418*10^14;    % standard gravitational parameter, m^3/s^2
Re = 6371000;               % mean Earth radius, m
J2 = 1082.7*(10^-6);        % coefficient of the Earth's gravitational zonal harmonic of the 2nd degree, -

% Variables
RAAN0 = 0*pi/180;       % RAAN at maneuver start, rads
u0 = 0*pi/180;          % AoL at maneuver start, rads
A = -0.00011667;        % propulsive acceleration, m/s^2 (negative for lowering case) 
B = -A;                 % propulsive acceleration, m/s^2 (positive for raising case) 
a0 = 6.96627*10^6;      % initial mean semi-major axis, m
a3 = a0;                % final mean semi-major axis, m
aref = a3;              % mean semi-major axis of non-maneuvering reference satellite for calculating RAANdif, m
incl = 45*pi/180;       % orbit inclination, rads
deltaVtotal = 30;       % deltaV used for maneuver, m/s
tmin = abs(deltaVtotal/A);   % time required to achieve deltaV, secs 
ttotal = 7*24*60*60;    % total maneuver time, sec

% Note ttotal must be larger than tmin to have enough time to achieve required deltaV
if ttotal < tmin
    disp('Maneuver time too short to achieve required deltaV')
    return
end

%%
% ---------------------------
% Lowering Maneuver Equations
% ---------------------------
% RAAN
% ---------------------------
deltaRAAN1_low = (3/134217728).*A.^(-1).*a0.^(-16).*J2.*mu.^(-7).*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^8.*Re.^2.*cos(incl).*(32768.*a0.^8.* ...
  mu.^4.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4).*(a0.^4+( ...
  -256).*a0.^4.*mu.^4.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).* ...
  (deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^( ...
  -4))+1024.*a0.^4.*J2.*mu.^2.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu) ...
  .^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^( ...
  1/2))).^(-2).*(a0.^6+(-4096).*a0.^6.*mu.^6.*(mu+a0.*(deltaVtotal+( ...
  a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+( ...
  a3.^(-1).*mu).^(1/2))).^(-6)).*Re.^2.*(1+3.*cos(2.*incl))+9.*J2.^2.*( ...
  a0.^8+(-65536).*a0.^8.*mu.^8.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu) ...
  .^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^( ...
  1/2))).^(-8)).*Re.^4.*(1+3.*cos(2.*incl)).^2);

deltaRAAN2_low = (-3/256).*a0.^(-2).*J2.*mu.^(-2).*(mu+a0.*(deltaVtotal+(a3.^(-1).* ...
  mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu) ...
  .^(1/2))).^2.*(a0.^(-3).*mu.^(-2).*(mu+a0.*(deltaVtotal+(a3.^(-1) ...
  .*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).* ...
  mu).^(1/2))).^3).^(1/2).*Re.^2.*cos(incl).*(1+(3/32).*a0.^(-2).*J2.* ...
  mu.^(-2).*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^2.* ...
  Re.^2.*(1+(-3/2).*sin(incl).^2)).*(ttotal+(-1/640).*A.^(-1).*a3.^( ...
  -5/2).*mu.^(1/2).*(a0.*mu.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^( ...
  1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2) ...
  )).^(-1)).^(-5/2).*((-320).*a0.^2.*a3.^(5/2).*mu.^2.*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2)+3.*a3.^(5/2).*J2.*Re.^2.*(( ...
  -2)+3.*sin(incl).^2)+32.*(a0.*mu.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu) ...
  .^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^( ...
  1/2))).^(-1)).^(5/2).*(20.*a3.^2+3.*J2.*Re.^2.*(2+(-3).*sin(incl).^2) ...
  ))+(1/640).*A.^(-1).*a0.^(-5/2).*mu.^(1/2).*(a0.*mu.*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(-5/2).*((-640).*a0.^2.*( ...
  a0.*mu.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(5/2)+ ...
  96.*J2.*(a0.*mu.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1) ...
  ).^(5/2).*Re.^2.*((-2)+3.*sin(incl).^2)+a0.^(5/2).*(320.*a0.^2.* ...
  mu.^2.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2)+3.*J2.* ...
  Re.^2.*(2+(-3).*sin(incl).^2))));

deltaRAAN3_low = (-3/134217728).*A.^(-1).*a0.^(-8).*a3.^(-8).*J2.*mu.^(-7).*(mu+ ...
  a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1) ...
  .*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^8.*Re.^2.*cos(incl).*(32768.* ...
  a0.^4.*a3.^4.*mu.^4.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).* ...
  (deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^( ...
  -4).*((-1).*a3.^4+256.*a0.^4.*mu.^4.*(mu+a0.*(deltaVtotal+(a3.^( ...
  -1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1) ...
  .*mu).^(1/2))).^(-4))+1024.*a0.^2.*a3.^2.*J2.*mu.^2.*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2).*((-1).*a3.^6+4096.*a0.^6.* ...
  mu.^6.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-6)).*Re.^2.*( ...
  1+3.*cos(2.*incl))+9.*J2.^2.*((-1).*a3.^8+65536.*a0.^8.*mu.^8.*(mu+ ...
  a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1) ...
  .*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-8)).*Re.^4.*(1+3.*cos(2.* ...
  incl)).^2);

RAANdif_low = (3/134217728).*A.^(-1).*a0.^(-16).*J2.*mu.^(-7).*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^8.*Re.^2.*cos(incl).*(32768.*a0.^8.* ...
  mu.^4.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4).*(a0.^4+( ...
  -256).*a0.^4.*mu.^4.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).* ...
  (deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^( ...
  -4))+1024.*a0.^4.*J2.*mu.^2.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu) ...
  .^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^( ...
  1/2))).^(-2).*(a0.^6+(-4096).*a0.^6.*mu.^6.*(mu+a0.*(deltaVtotal+( ...
  a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+( ...
  a3.^(-1).*mu).^(1/2))).^(-6)).*Re.^2.*(1+3.*cos(2.*incl))+9.*J2.^2.*( ...
  a0.^8+(-65536).*a0.^8.*mu.^8.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu) ...
  .^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^( ...
  1/2))).^(-8)).*Re.^4.*(1+3.*cos(2.*incl)).^2)+(-3/134217728).*A.^(-1) ...
  .*a0.^(-8).*a3.^(-8).*J2.*mu.^(-7).*(mu+a0.*(deltaVtotal+(a3.^(-1) ...
  .*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).* ...
  mu).^(1/2))).^8.*Re.^2.*cos(incl).*(32768.*a0.^4.*a3.^4.*mu.^4.*(mu+ ...
  a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1) ...
  .*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4).*((-1).*a3.^4+256.* ...
  a0.^4.*mu.^4.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4) ...
  )+1024.*a0.^2.*a3.^2.*J2.*mu.^2.*(mu+a0.*(deltaVtotal+(a3.^(-1).* ...
  mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu) ...
  .^(1/2))).^(-2).*((-1).*a3.^6+4096.*a0.^6.*mu.^6.*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-6)).*Re.^2.*(1+3.*cos(2.*incl))+ ...
  9.*J2.^2.*((-1).*a3.^8+65536.*a0.^8.*mu.^8.*(mu+a0.*(deltaVtotal+( ...
  a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+( ...
  a3.^(-1).*mu).^(1/2))).^(-8)).*Re.^4.*(1+3.*cos(2.*incl)).^2)+(3/2).* ...
  aref.^(-2).*J2.*(aref.^(-3).*mu).^(1/2).*Re.^2.*ttotal.*cos(incl).*( ...
  1+(3/2).*aref.^(-2).*J2.*Re.^2.*(1+(-3/2).*sin(incl).^2))+(-3/256).* ...
  a0.^(-2).*J2.*mu.^(-2).*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2) ...
  ).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))) ...
  .^2.*(a0.^(-3).*mu.^(-2).*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^( ...
  1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2) ...
  )).^3).^(1/2).*Re.^2.*cos(incl).*(1+(3/32).*a0.^(-2).*J2.*mu.^(-2).*( ...
  mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^2.*Re.^2.*(1+(-3/2).*sin( ...
  incl).^2)).*(ttotal+(-1/640).*A.^(-1).*a3.^(-5/2).*mu.^(1/2).*(a0.* ...
  mu.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*( ...
  a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(-5/2).*(( ...
  -320).*a0.^2.*a3.^(5/2).*mu.^2.*(mu+a0.*(deltaVtotal+(a3.^(-1).* ...
  mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu) ...
  .^(1/2))).^(-2)+3.*a3.^(5/2).*J2.*Re.^2.*((-2)+3.*sin(incl).^2)+32.*( ...
  a0.*mu.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(5/2).*( ...
  20.*a3.^2+3.*J2.*Re.^2.*(2+(-3).*sin(incl).^2)))+(1/640).*A.^(-1).* ...
  a0.^(-5/2).*mu.^(1/2).*(a0.*mu.*(mu+a0.*(deltaVtotal+(a3.^(-1).* ...
  mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu) ...
  .^(1/2))).^(-1)).^(-5/2).*((-640).*a0.^2.*(a0.*mu.*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(5/2)+96.*J2.*(a0.*mu.*( ...
  mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(5/2).*Re.^2.*(( ...
  -2)+3.*sin(incl).^2)+a0.^(5/2).*(320.*a0.^2.*mu.^2.*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2)+3.*J2.*Re.^2.*(2+(-3).*sin( ...
  incl).^2)))); 


RAANtotal_low = RAAN0 + deltaRAAN1_low + deltaRAAN2_low + deltaRAAN3_low;

% ---------------------------
% AOL
% ---------------------------

deltau1_low = (-1/536870912).*A.^(-1).*a0.^(-16).*mu.^(-7).*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^8.*(8388608.*a0.^12.*mu.^6.*(mu+ ...
  a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1) ...
  .*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-6).*(a0.^2+(-16).*a0.^2.* ...
  mu.^2.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2))+27.* ...
  J2.^3.*(a0.^8+(-65536).*a0.^8.*mu.^8.*(mu+a0.*(deltaVtotal+(a3.^( ...
  -1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1) ...
  .*mu).^(1/2))).^(-8)).*Re.^6.*(1+3.*cos(2.*incl)).^2.*(3+5.*cos(2.*incl) ...
  )+98304.*a0.^8.*J2.*mu.^4.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^( ...
  1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2) ...
  )).^(-4).*(a0.^4+(-256).*a0.^4.*mu.^4.*(mu+a0.*(deltaVtotal+(a3.^( ...
  -1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1) ...
  .*mu).^(1/2))).^(-4)).*Re.^2.*(5+11.*cos(2.*incl))+768.*a0.^4.* ...
  J2.^2.*mu.^2.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2) ...
  .*(a0.^6+(-4096).*a0.^6.*mu.^6.*(mu+a0.*(deltaVtotal+(a3.^(-1).* ...
  mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu) ...
  .^(1/2))).^(-6)).*Re.^4.*(53+68.*cos(2.*incl)+39.*cos(4.*incl)));

deltau2_low = (ttotal+(1/2).*A.^(-1).*((-2).*a3.^(-1/2).*mu.^(1/2)+mu.^(1/2).*( ...
  a0.*mu.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(-1/2)+( ...
  -3/5).*a3.^(-5/2).*J2.*mu.^(1/2).*Re.^2+(3/160).*J2.*mu.^(1/2).*( ...
  a0.*mu.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(-5/2).* ...
  Re.^2+(9/10).*a3.^(-5/2).*J2.*mu.^(1/2).*Re.^2.*sin(incl).^2+(-9/320) ...
  .*J2.*mu.^(1/2).*(a0.*mu.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^( ...
  1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2) ...
  )).^(-1)).^(-5/2).*Re.^2.*sin(incl).^2)+(-1/2).*A.^(-1).*(2.*a0.^( ...
  -1/2).*mu.^(1/2)+(-1).*mu.^(1/2).*(a0.*mu.*(mu+a0.*(deltaVtotal+( ...
  a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+( ...
  a3.^(-1).*mu).^(1/2))).^(-1)).^(-1/2)+(3/5).*a0.^(-5/2).*J2.*mu.^( ...
  1/2).*Re.^2+(-3/160).*J2.*mu.^(1/2).*(a0.*mu.*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(-5/2).*Re.^2+(-9/10).* ...
  a0.^(-5/2).*J2.*mu.^(1/2).*Re.^2.*sin(incl).^2+(9/320).*J2.*mu.^(1/2) ...
  .*(a0.*mu.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1) ...
  ).^(-5/2).*Re.^2.*sin(incl).^2)).*((1/8).*(a0.^(-3).*mu.^(-2).*(mu+ ...
  a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1) ...
  .*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^3).^(1/2).*(1+(3/32).*a0.^( ...
  -2).*J2.*mu.^(-2).*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^2.* ...
  Re.^2.*(1+(-3/2).*sin(incl).^2))+(3/256).*a0.^(-2).*J2.*mu.^(-2).*( ...
  mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^2.*(a0.^(-3).*mu.^(-2).*( ...
  mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^3).^(1/2).*Re.^2.*(2+( ...
  -5/2).*sin(incl).^2).*(1+(3/32).*a0.^(-2).*J2.*mu.^(-2).*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^2.*Re.^2.*(1+(-3/2).*sin(incl).^2))) ...
  ;

deltau3_low = (1/536870912).*A.^(-1).*a0.^(-8).*a3.^(-8).*mu.^(-7).*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^8.*(8388608.*a0.^6.*a3.^6.* ...
  mu.^6.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-6).*((-1).* ...
  a3.^2+16.*a0.^2.*mu.^2.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2) ...
  ).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))) ...
  .^(-2))+27.*J2.^3.*((-1).*a3.^8+65536.*a0.^8.*mu.^8.*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-8)).*Re.^6.*(1+3.*cos(2.*incl)) ...
  .^2.*(3+5.*cos(2.*incl))+98304.*a0.^4.*a3.^4.*J2.*mu.^4.*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4).*((-1).*a3.^4+256.*a0.^4.* ...
  mu.^4.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4)).*Re.^2.*( ...
  5+11.*cos(2.*incl))+768.*a0.^2.*a3.^2.*J2.^2.*mu.^2.*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2).*((-1).*a3.^6+4096.*a0.^6.* ...
  mu.^6.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-6)).*Re.^4.*( ...
  53+68.*cos(2.*incl)+39.*cos(4.*incl)));

udif_low = (-1/536870912).*A.^(-1).*a0.^(-16).*mu.^(-7).*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^8.*(8388608.*a0.^12.*mu.^6.*(mu+ ...
  a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1) ...
  .*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-6).*(a0.^2+(-16).*a0.^2.* ...
  mu.^2.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2))+27.* ...
  J2.^3.*(a0.^8+(-65536).*a0.^8.*mu.^8.*(mu+a0.*(deltaVtotal+(a3.^( ...
  -1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1) ...
  .*mu).^(1/2))).^(-8)).*Re.^6.*(1+3.*cos(2.*incl)).^2.*(3+5.*cos(2.*incl) ...
  )+98304.*a0.^8.*J2.*mu.^4.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^( ...
  1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2) ...
  )).^(-4).*(a0.^4+(-256).*a0.^4.*mu.^4.*(mu+a0.*(deltaVtotal+(a3.^( ...
  -1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1) ...
  .*mu).^(1/2))).^(-4)).*Re.^2.*(5+11.*cos(2.*incl))+768.*a0.^4.* ...
  J2.^2.*mu.^2.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2) ...
  .*(a0.^6+(-4096).*a0.^6.*mu.^6.*(mu+a0.*(deltaVtotal+(a3.^(-1).* ...
  mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu) ...
  .^(1/2))).^(-6)).*Re.^4.*(53+68.*cos(2.*incl)+39.*cos(4.*incl)))+( ...
  1/536870912).*A.^(-1).*a0.^(-8).*a3.^(-8).*mu.^(-7).*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^8.*(8388608.*a0.^6.*a3.^6.* ...
  mu.^6.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-6).*((-1).* ...
  a3.^2+16.*a0.^2.*mu.^2.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2) ...
  ).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))) ...
  .^(-2))+27.*J2.^3.*((-1).*a3.^8+65536.*a0.^8.*mu.^8.*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-8)).*Re.^6.*(1+3.*cos(2.*incl)) ...
  .^2.*(3+5.*cos(2.*incl))+98304.*a0.^4.*a3.^4.*J2.*mu.^4.*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4).*((-1).*a3.^4+256.*a0.^4.* ...
  mu.^4.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4)).*Re.^2.*( ...
  5+11.*cos(2.*incl))+768.*a0.^2.*a3.^2.*J2.^2.*mu.^2.*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2).*((-1).*a3.^6+4096.*a0.^6.* ...
  mu.^6.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-6)).*Re.^4.*( ...
  53+68.*cos(2.*incl)+39.*cos(4.*incl)))+ttotal.*((-1).*(aref.^(-3).*mu) ...
  .^(1/2).*(1+(3/2).*aref.^(-2).*J2.*Re.^2.*(1+(-3/2).*sin(incl).^2))+( ...
  -3/2).*aref.^(-2).*J2.*(aref.^(-3).*mu).^(1/2).*Re.^2.*(2+(-5/2).* ...
  sin(incl).^2).*(1+(3/2).*aref.^(-2).*J2.*Re.^2.*(1+(-3/2).*sin(incl).^2) ...
  ))+(ttotal+(1/2).*A.^(-1).*((-2).*a3.^(-1/2).*mu.^(1/2)+mu.^(1/2) ...
  .*(a0.*mu.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1) ...
  ).^(-1/2)+(-3/5).*a3.^(-5/2).*J2.*mu.^(1/2).*Re.^2+(3/160).*J2.* ...
  mu.^(1/2).*(a0.*mu.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1) ...
  ).^(-5/2).*Re.^2+(9/10).*a3.^(-5/2).*J2.*mu.^(1/2).*Re.^2.*sin(incl) ...
  .^2+(-9/320).*J2.*mu.^(1/2).*(a0.*mu.*(mu+a0.*(deltaVtotal+(a3.^( ...
  -1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1) ...
  .*mu).^(1/2))).^(-1)).^(-5/2).*Re.^2.*sin(incl).^2)+(-1/2).*A.^(-1).* ...
  (2.*a0.^(-1/2).*mu.^(1/2)+(-1).*mu.^(1/2).*(a0.*mu.*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(-1/2)+(3/5).*a0.^(-5/2).* ...
  J2.*mu.^(1/2).*Re.^2+(-3/160).*J2.*mu.^(1/2).*(a0.*mu.*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(-5/2).*Re.^2+(-9/10).* ...
  a0.^(-5/2).*J2.*mu.^(1/2).*Re.^2.*sin(incl).^2+(9/320).*J2.*mu.^(1/2) ...
  .*(a0.*mu.*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1) ...
  ).^(-5/2).*Re.^2.*sin(incl).^2)).*((1/8).*(a0.^(-3).*mu.^(-2).*(mu+ ...
  a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1) ...
  .*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^3).^(1/2).*(1+(3/32).*a0.^( ...
  -2).*J2.*mu.^(-2).*(mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^2.* ...
  Re.^2.*(1+(-3/2).*sin(incl).^2))+(3/256).*a0.^(-2).*J2.*mu.^(-2).*( ...
  mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^2.*(a0.^(-3).*mu.^(-2).*( ...
  mu+a0.*(deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^3).^(1/2).*Re.^2.*(2+( ...
  -5/2).*sin(incl).^2).*(1+(3/32).*a0.^(-2).*J2.*mu.^(-2).*(mu+a0.*( ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*(deltaVtotal+2.*(a0.^(-1).*mu) ...
  .^(1/2)+(a3.^(-1).*mu).^(1/2))).^2.*Re.^2.*(1+(-3/2).*sin(incl).^2))) ...
  ; 

utotal_low = u0 + deltau1_low + deltau2_low + deltau3_low; 


% --------------
% Print Results
% --------------
disp('deltaRAAN_low')
disp(deltaRAAN1_low)
disp('deltaRAAN2_low')
disp(deltaRAAN2_low)
disp('deltaRAAN3_low')
disp(deltaRAAN3_low)
disp('RAANdif_low')
disp(RAANdif_low)
disp('RAANtotal_low')
disp(RAANtotal_low)
disp('deltau1_low')
disp(deltau1_low)
disp('deltau2_low')
disp(deltau2_low)
disp('deltau3_low')
disp(deltau3_low)
disp('udif_low')
disp(udif_low)
disp('utotal_low')
disp(utotal_low)



%%
% ---------------------------
% Raising Maneuver Equations
% ---------------------------
% RAAN
% ---------------------------
deltaRAAN1_raise = (3/134217728).*B.^(-1).*a0.^(-16).*J2.*mu.^(-7).*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^8.*Re.^2.*cos(incl).*( ...
  32768.*a0.^8.*mu.^4.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^( ...
  1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu) ...
  .^(1/2))).^(-4).*(a0.^4+(-256).*a0.^4.*mu.^4.*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4))+1024.*a0.^4.*J2.* ...
  mu.^2.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).* ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2) ...
  .*(a0.^6+(-4096).*a0.^6.*mu.^6.*(mu+a0.*((-1).*deltaVtotal+(a3.^( ...
  -1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+( ...
  a3.^(-1).*mu).^(1/2))).^(-6)).*Re.^2.*(1+3.*cos(2.*incl))+9.*J2.^2.*( ...
  a0.^8+(-65536).*a0.^8.*mu.^8.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1) ...
  .*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^( ...
  -1).*mu).^(1/2))).^(-8)).*Re.^4.*(1+3.*cos(2.*incl)).^2);

deltaRAAN2_raise = (-3/256).*a0.^(-2).*J2.*mu.^(-2).*(mu+a0.*((-1).*deltaVtotal+( ...
  a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2) ...
  +(a3.^(-1).*mu).^(1/2))).^2.*(a0.^(-3).*mu.^(-2).*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^3).^(1/2).*Re.^2.*cos(incl) ...
  .*(1+(3/32).*a0.^(-2).*J2.*mu.^(-2).*(mu+a0.*((-1).*deltaVtotal+( ...
  a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2) ...
  +(a3.^(-1).*mu).^(1/2))).^2.*Re.^2.*(1+(-3/2).*sin(incl).^2)).*( ...
  ttotal+(-1/640).*B.^(-1).*a3.^(-5/2).*mu.^(1/2).*(a0.*mu.*(mu+a0.* ...
  ((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*( ...
  a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(-5/2).*(( ...
  -320).*a0.^2.*a3.^(5/2).*mu.^2.*(mu+a0.*((-1).*deltaVtotal+(a3.^( ...
  -1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+( ...
  a3.^(-1).*mu).^(1/2))).^(-2)+3.*a3.^(5/2).*J2.*Re.^2.*((-2)+3.* ...
  sin(incl).^2)+32.*(a0.*mu.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu) ...
  .^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).* ...
  mu).^(1/2))).^(-1)).^(5/2).*(20.*a3.^2+3.*J2.*Re.^2.*(2+(-3).*sin( ...
  incl).^2)))+(1/640).*B.^(-1).*a0.^(-5/2).*mu.^(1/2).*(a0.*mu.*(mu+ ...
  a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(-5/2).*( ...
  (-640).*a0.^2.*(a0.*mu.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu) ...
  .^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).* ...
  mu).^(1/2))).^(-1)).^(5/2)+96.*J2.*(a0.*mu.*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(5/2).*Re.^2.*(( ...
  -2)+3.*sin(incl).^2)+a0.^(5/2).*(320.*a0.^2.*mu.^2.*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2)+3.*J2.*Re.^2.*(2+( ...
  -3).*sin(incl).^2))));

deltaRAAN3_raise = (-3/134217728).*B.^(-1).*a0.^(-8).*a3.^(-8).*J2.*mu.^(-7).*(mu+ ...
  a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^8.*Re.^2.*cos(incl) ...
  .*(32768.*a0.^4.*a3.^4.*mu.^4.*(mu+a0.*((-1).*deltaVtotal+(a3.^( ...
  -1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+( ...
  a3.^(-1).*mu).^(1/2))).^(-4).*((-1).*a3.^4+256.*a0.^4.*mu.^4.*(mu+ ...
  a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4))+1024.* ...
  a0.^2.*a3.^2.*J2.*mu.^2.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu) ...
  .^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).* ...
  mu).^(1/2))).^(-2).*((-1).*a3.^6+4096.*a0.^6.*mu.^6.*(mu+a0.*((-1) ...
  .*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-6)).*Re.^2.*(1+3.*cos( ...
  2.*incl))+9.*J2.^2.*((-1).*a3.^8+65536.*a0.^8.*mu.^8.*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-8)).*Re.^4.*(1+3.*cos( ...
  2.*incl)).^2);

RAANdif_raise = (3/134217728).*B.^(-1).*a0.^(-16).*J2.*mu.^(-7).*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^8.*Re.^2.*cos(incl).*( ...
  32768.*a0.^8.*mu.^4.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^( ...
  1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu) ...
  .^(1/2))).^(-4).*(a0.^4+(-256).*a0.^4.*mu.^4.*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4))+1024.*a0.^4.*J2.* ...
  mu.^2.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).* ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2) ...
  .*(a0.^6+(-4096).*a0.^6.*mu.^6.*(mu+a0.*((-1).*deltaVtotal+(a3.^( ...
  -1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+( ...
  a3.^(-1).*mu).^(1/2))).^(-6)).*Re.^2.*(1+3.*cos(2.*incl))+9.*J2.^2.*( ...
  a0.^8+(-65536).*a0.^8.*mu.^8.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1) ...
  .*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^( ...
  -1).*mu).^(1/2))).^(-8)).*Re.^4.*(1+3.*cos(2.*incl)).^2)+( ...
  -3/134217728).*B.^(-1).*a0.^(-8).*a3.^(-8).*J2.*mu.^(-7).*(mu+a0.* ...
  ((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*( ...
  a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^8.*Re.^2.*cos(incl).*( ...
  32768.*a0.^4.*a3.^4.*mu.^4.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).* ...
  mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1) ...
  .*mu).^(1/2))).^(-4).*((-1).*a3.^4+256.*a0.^4.*mu.^4.*(mu+a0.*(( ...
  -1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*( ...
  a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4))+1024.*a0.^2.* ...
  a3.^2.*J2.*mu.^2.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2) ...
  ).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^( ...
  1/2))).^(-2).*((-1).*a3.^6+4096.*a0.^6.*mu.^6.*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-6)).*Re.^2.*(1+3.*cos( ...
  2.*incl))+9.*J2.^2.*((-1).*a3.^8+65536.*a0.^8.*mu.^8.*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-8)).*Re.^4.*(1+3.*cos( ...
  2.*incl)).^2)+(3/2).*aref.^(-2).*J2.*(aref.^(-3).*mu).^(1/2).*Re.^2.* ...
  ttotal.*cos(incl).*(1+(3/2).*aref.^(-2).*J2.*Re.^2.*(1+(-3/2).*sin(incl) ...
  .^2))+(-3/256).*a0.^(-2).*J2.*mu.^(-2).*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^2.*(a0.^(-3).*mu.^(-2).*( ...
  mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).* ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^3) ...
  .^(1/2).*Re.^2.*cos(incl).*(1+(3/32).*a0.^(-2).*J2.*mu.^(-2).*(mu+ ...
  a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^2.*Re.^2.*(1+( ...
  -3/2).*sin(incl).^2)).*(ttotal+(-1/640).*B.^(-1).*a3.^(-5/2).*mu.^( ...
  1/2).*(a0.*mu.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).* ...
  ((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2)) ...
  ).^(-1)).^(-5/2).*((-320).*a0.^2.*a3.^(5/2).*mu.^2.*(mu+a0.*((-1) ...
  .*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2)+3.*a3.^(5/2).*J2.* ...
  Re.^2.*((-2)+3.*sin(incl).^2)+32.*(a0.*mu.*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(5/2).*(20.*a3.^2+ ...
  3.*J2.*Re.^2.*(2+(-3).*sin(incl).^2)))+(1/640).*B.^(-1).*a0.^(-5/2).* ...
  mu.^(1/2).*(a0.*mu.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^( ...
  1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu) ...
  .^(1/2))).^(-1)).^(-5/2).*((-640).*a0.^2.*(a0.*mu.*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(5/2)+96.*J2.*( ...
  a0.*mu.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).* ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1) ...
  ).^(5/2).*Re.^2.*((-2)+3.*sin(incl).^2)+a0.^(5/2).*(320.*a0.^2.* ...
  mu.^2.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).* ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2) ...
  +3.*J2.*Re.^2.*(2+(-3).*sin(incl).^2))));


RAANtotal_raise = RAAN0 + deltaRAAN1_raise + deltaRAAN2_raise + deltaRAAN3_raise;

% ---------------------------
% AOL
% ---------------------------

deltau1_raise = (-1/536870912).*B.^(-1).*a0.^(-16).*mu.^(-7).*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^8.*(8388608.*a0.^12.* ...
  mu.^6.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).* ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-6) ...
  .*(a0.^2+(-16).*a0.^2.*mu.^2.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1) ...
  .*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^( ...
  -1).*mu).^(1/2))).^(-2))+27.*J2.^3.*(a0.^8+(-65536).*a0.^8.* ...
  mu.^8.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).* ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-8) ...
  ).*Re.^6.*(1+3.*cos(2.*incl)).^2.*(3+5.*cos(2.*incl))+98304.*a0.^8.*J2.* ...
  mu.^4.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).* ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4) ...
  .*(a0.^4+(-256).*a0.^4.*mu.^4.*(mu+a0.*((-1).*deltaVtotal+(a3.^( ...
  -1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+( ...
  a3.^(-1).*mu).^(1/2))).^(-4)).*Re.^2.*(5+11.*cos(2.*incl))+768.* ...
  a0.^4.*J2.^2.*mu.^2.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^( ...
  1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu) ...
  .^(1/2))).^(-2).*(a0.^6+(-4096).*a0.^6.*mu.^6.*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-6)).*Re.^4.*(53+68.* ...
  cos(2.*incl)+39.*cos(4.*incl)));

deltau2_raise = (ttotal+(1/2).*B.^(-1).*((-2).*a3.^(-1/2).*mu.^(1/2)+mu.^(1/2).*( ...
  a0.*mu.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).* ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1) ...
  ).^(-1/2)+(-3/5).*a3.^(-5/2).*J2.*mu.^(1/2).*Re.^2+(3/160).*J2.* ...
  mu.^(1/2).*(a0.*mu.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^( ...
  1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu) ...
  .^(1/2))).^(-1)).^(-5/2).*Re.^2+(9/10).*a3.^(-5/2).*J2.*mu.^(1/2) ...
  .*Re.^2.*sin(incl).^2+(-9/320).*J2.*mu.^(1/2).*(a0.*mu.*(mu+a0.*((-1) ...
  .*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(-5/2).*Re.^2.* ...
  sin(incl).^2)+(-1/2).*B.^(-1).*(2.*a0.^(-1/2).*mu.^(1/2)+(-1).*mu.^( ...
  1/2).*(a0.*mu.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).* ...
  ((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2)) ...
  ).^(-1)).^(-1/2)+(3/5).*a0.^(-5/2).*J2.*mu.^(1/2).*Re.^2+(-3/160) ...
  .*J2.*mu.^(1/2).*(a0.*mu.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).* ...
  mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1) ...
  .*mu).^(1/2))).^(-1)).^(-5/2).*Re.^2+(-9/10).*a0.^(-5/2).*J2.* ...
  mu.^(1/2).*Re.^2.*sin(incl).^2+(9/320).*J2.*mu.^(1/2).*(a0.*mu.*(mu+ ...
  a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(-5/2).* ...
  Re.^2.*sin(incl).^2)).*((1/8).*(a0.^(-3).*mu.^(-2).*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^3).^(1/2).*(1+(3/32).* ...
  a0.^(-2).*J2.*mu.^(-2).*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu) ...
  .^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).* ...
  mu).^(1/2))).^2.*Re.^2.*(1+(-3/2).*sin(incl).^2))+(3/256).*a0.^(-2).* ...
  J2.*mu.^(-2).*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  (-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))) ...
  .^2.*(a0.^(-3).*mu.^(-2).*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).* ...
  mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1) ...
  .*mu).^(1/2))).^3).^(1/2).*Re.^2.*(2+(-5/2).*sin(incl).^2).*(1+(3/32) ...
  .*a0.^(-2).*J2.*mu.^(-2).*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).* ...
  mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1) ...
  .*mu).^(1/2))).^2.*Re.^2.*(1+(-3/2).*sin(incl).^2)));

deltau3_raise = (1/536870912).*B.^(-1).*a0.^(-8).*a3.^(-8).*mu.^(-7).*(mu+a0.*(( ...
  -1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*( ...
  a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^8.*(8388608.*a0.^6.* ...
  a3.^6.*mu.^6.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  (-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))) ...
  .^(-6).*((-1).*a3.^2+16.*a0.^2.*mu.^2.*(mu+a0.*((-1).*deltaVtotal+ ...
  (a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^( ...
  1/2)+(a3.^(-1).*mu).^(1/2))).^(-2))+27.*J2.^3.*((-1).*a3.^8+ ...
  65536.*a0.^8.*mu.^8.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^( ...
  1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu) ...
  .^(1/2))).^(-8)).*Re.^6.*(1+3.*cos(2.*incl)).^2.*(3+5.*cos(2.*incl))+ ...
  98304.*a0.^4.*a3.^4.*J2.*mu.^4.*(mu+a0.*((-1).*deltaVtotal+(a3.^( ...
  -1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+( ...
  a3.^(-1).*mu).^(1/2))).^(-4).*((-1).*a3.^4+256.*a0.^4.*mu.^4.*(mu+ ...
  a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4)).*Re.^2.*( ...
  5+11.*cos(2.*incl))+768.*a0.^2.*a3.^2.*J2.^2.*mu.^2.*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2).*((-1).*a3.^6+4096.* ...
  a0.^6.*mu.^6.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  (-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))) ...
  .^(-6)).*Re.^4.*(53+68.*cos(2.*incl)+39.*cos(4.*incl)));

udif_raise = (-1/536870912).*B.^(-1).*a0.^(-16).*mu.^(-7).*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^8.*(8388608.*a0.^12.* ...
  mu.^6.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).* ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-6) ...
  .*(a0.^2+(-16).*a0.^2.*mu.^2.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1) ...
  .*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^( ...
  -1).*mu).^(1/2))).^(-2))+27.*J2.^3.*(a0.^8+(-65536).*a0.^8.* ...
  mu.^8.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).* ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-8) ...
  ).*Re.^6.*(1+3.*cos(2.*incl)).^2.*(3+5.*cos(2.*incl))+98304.*a0.^8.*J2.* ...
  mu.^4.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).* ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4) ...
  .*(a0.^4+(-256).*a0.^4.*mu.^4.*(mu+a0.*((-1).*deltaVtotal+(a3.^( ...
  -1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+( ...
  a3.^(-1).*mu).^(1/2))).^(-4)).*Re.^2.*(5+11.*cos(2.*incl))+768.* ...
  a0.^4.*J2.^2.*mu.^2.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^( ...
  1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu) ...
  .^(1/2))).^(-2).*(a0.^6+(-4096).*a0.^6.*mu.^6.*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-6)).*Re.^4.*(53+68.* ...
  cos(2.*incl)+39.*cos(4.*incl)))+(1/536870912).*B.^(-1).*a0.^(-8).*a3.^( ...
  -8).*mu.^(-7).*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).* ...
  ((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2)) ...
  ).^8.*(8388608.*a0.^6.*a3.^6.*mu.^6.*(mu+a0.*((-1).*deltaVtotal+( ...
  a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2) ...
  +(a3.^(-1).*mu).^(1/2))).^(-6).*((-1).*a3.^2+16.*a0.^2.*mu.^2.*( ...
  mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).* ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-2) ...
  )+27.*J2.^3.*((-1).*a3.^8+65536.*a0.^8.*mu.^8.*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-8)).*Re.^6.*(1+3.*cos( ...
  2.*incl)).^2.*(3+5.*cos(2.*incl))+98304.*a0.^4.*a3.^4.*J2.*mu.^4.*(mu+ ...
  a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-4).*((-1).* ...
  a3.^4+256.*a0.^4.*mu.^4.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu) ...
  .^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).* ...
  mu).^(1/2))).^(-4)).*Re.^2.*(5+11.*cos(2.*incl))+768.*a0.^2.*a3.^2.* ...
  J2.^2.*mu.^2.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  (-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))) ...
  .^(-2).*((-1).*a3.^6+4096.*a0.^6.*mu.^6.*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-6)).*Re.^4.*(53+68.* ...
  cos(2.*incl)+39.*cos(4.*incl)))+ttotal.*((-1).*(aref.^(-3).*mu).^(1/2).* ...
  (1+(3/2).*aref.^(-2).*J2.*Re.^2.*(1+(-3/2).*sin(incl).^2))+(-3/2).* ...
  aref.^(-2).*J2.*(aref.^(-3).*mu).^(1/2).*Re.^2.*(2+(-5/2).*sin(incl) ...
  .^2).*(1+(3/2).*aref.^(-2).*J2.*Re.^2.*(1+(-3/2).*sin(incl).^2)))+( ...
  ttotal+(1/2).*B.^(-1).*((-2).*a3.^(-1/2).*mu.^(1/2)+mu.^(1/2).*( ...
  a0.*mu.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).* ...
  deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1) ...
  ).^(-1/2)+(-3/5).*a3.^(-5/2).*J2.*mu.^(1/2).*Re.^2+(3/160).*J2.* ...
  mu.^(1/2).*(a0.*mu.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^( ...
  1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu) ...
  .^(1/2))).^(-1)).^(-5/2).*Re.^2+(9/10).*a3.^(-5/2).*J2.*mu.^(1/2) ...
  .*Re.^2.*sin(incl).^2+(-9/320).*J2.*mu.^(1/2).*(a0.*mu.*(mu+a0.*((-1) ...
  .*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(-5/2).*Re.^2.* ...
  sin(incl).^2)+(-1/2).*B.^(-1).*(2.*a0.^(-1/2).*mu.^(1/2)+(-1).*mu.^( ...
  1/2).*(a0.*mu.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).* ...
  ((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2)) ...
  ).^(-1)).^(-1/2)+(3/5).*a0.^(-5/2).*J2.*mu.^(1/2).*Re.^2+(-3/160) ...
  .*J2.*mu.^(1/2).*(a0.*mu.*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).* ...
  mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1) ...
  .*mu).^(1/2))).^(-1)).^(-5/2).*Re.^2+(-9/10).*a0.^(-5/2).*J2.* ...
  mu.^(1/2).*Re.^2.*sin(incl).^2+(9/320).*J2.*mu.^(1/2).*(a0.*mu.*(mu+ ...
  a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+ ...
  2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^(-1)).^(-5/2).* ...
  Re.^2.*sin(incl).^2)).*((1/8).*(a0.^(-3).*mu.^(-2).*(mu+a0.*((-1).* ...
  deltaVtotal+(a3.^(-1).*mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^( ...
  -1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))).^3).^(1/2).*(1+(3/32).* ...
  a0.^(-2).*J2.*mu.^(-2).*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu) ...
  .^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).* ...
  mu).^(1/2))).^2.*Re.^2.*(1+(-3/2).*sin(incl).^2))+(3/256).*a0.^(-2).* ...
  J2.*mu.^(-2).*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).*mu).^(1/2)).*( ...
  (-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1).*mu).^(1/2))) ...
  .^2.*(a0.^(-3).*mu.^(-2).*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).* ...
  mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1) ...
  .*mu).^(1/2))).^3).^(1/2).*Re.^2.*(2+(-5/2).*sin(incl).^2).*(1+(3/32) ...
  .*a0.^(-2).*J2.*mu.^(-2).*(mu+a0.*((-1).*deltaVtotal+(a3.^(-1).* ...
  mu).^(1/2)).*((-1).*deltaVtotal+2.*(a0.^(-1).*mu).^(1/2)+(a3.^(-1) ...
  .*mu).^(1/2))).^2.*Re.^2.*(1+(-3/2).*sin(incl).^2)));

utotal_raise = u0 + deltau1_raise + deltau2_raise + deltau3_raise; 


% --------------
% Print Results
% -------------
disp('deltaRAAN1_raise')
disp(deltaRAAN1_raise)
disp('deltaRAAN2_raise')
disp(deltaRAAN2_raise)
disp('deltaRAAN3_raise')
disp(deltaRAAN3_raise)
disp('RAANdif_raise')
disp(RAANdif_raise)
disp('RAANtotal_raise')
disp(RAANtotal_raise)
disp('deltau1_raise')
disp(deltau1_raise)
disp('deltau2_raise')
disp(deltau2_raise)
disp('deltau3_raise')
disp(deltau3_raise)
disp('udif_raise')
disp(udif_raise)
disp('utotal_raise')
disp(utotal_raise)