EqZonalMeanProp.m

%% EqZonalMeanProp
% Desc: Propagates given differential mean equinoctial elements for posigrade elliptic orbit
% that includes the effects of any zonal greater than J2 up to
% second-order. The zonals higher than J2 are considered of second-order.
% All the zonals up to given degree are included for computing mean rates.
% 
% Equinoctial set:
% (a,l+AOP+RAAN,tan(i/2)*cos(RAAN),tan(i/2)*sin(RAAN),e*cos(AOP+RAAN),e*sin(AOP+RAAN))
% 
% Reference: 
% B. Mahajan, S. R. Vadali and K. T. Alfriend, �Analytic Solution For Satellite 
% Relative Motion: The Complete Zonal Gravitational Problem,� 26th AAS/AIAA 
% Space Flight Mechanics Meeting, Napa, CA, February 2016.
%
% Author: Bharat Mahajan (https://github.com/princemahajan)


function [Xt, phi] = EqZonalMeanProp(dt, X0, degree, mu, Re, Jcoeff) 


a = X0(1);
p1 = X0(3);
p2 = X0(4);
q1 = X0(5);
q2 = X0(6);


e = sqrt(q1^2 + q2^2);
inc = 2*atan(sqrt(p1^2 + p2^2));
eta = sqrt(1 - e^2);
L = sqrt(mu*a);
G = L*eta;
H = G*cos(inc);


% for each zonal
narr = 3:1:degree;
J2 = Jcoeff(2);

% if two-body propagation
if degree < 2
   J2 = 0; 
end


% Rates due to J2 up to order 3 (from Maple code)
R__e = Re;
ci = cos(inc);
lghdot2 = -(J2 ^ 3 * (708750 * ci ^ 12 * eta ^ 9 + 12656250 * eta ^ 8 * ci ^ 12 + 73883750 * eta ^ 7 * ci ^ 12 - 607500 * ci ^ 11 * eta ^ 8 - 1567125 * ci ^ 10 * eta ^ 9 + 188101250 * eta ^ 6 * ci ^ 12 - 10935000 * ci ^ 11 * eta ^ 7 - 28112625 * eta ^ 8 * ci ^ 10 + 154953750 * eta ^ 5 * ci ^ 12 - 50266500 * ci ^ 11 * eta ^ 6 - 161286000 * ci ^ 10 * eta ^ 7 + 1017000 * ci ^ 9 * eta ^ 8 + 1188075 * ci ^ 8 * eta ^ 9 - 101418750 * eta ^ 4 * ci ^ 12 - 64041000 * ci ^ 11 * eta ^ 5 - 398143500 * ci ^ 10 * eta ^ 6 + 18609000 * ci ^ 9 * eta ^ 7 + 22653975 * eta ^ 8 * ci ^ 8 - 328983750 * eta ^ 3 * ci ^ 12 + 10528500 * ci ^ 11 * eta ^ 4 - 314078250 * ci ^ 10 * eta ^ 5 + 86214900 * ci ^ 9 * eta ^ 6 + 130452350 * ci ^ 8 * eta ^ 7 - 583800 * ci ^ 7 * eta ^ 8 - 381150 * ci ^ 6 * eta ^ 9 - 339521250 * eta ^ 2 * ci ^ 12 + 97299000 * ci ^ 11 * eta ^ 3 + 173327250 * ci ^ 10 * eta ^ 4 + 108264600 * ci ^ 9 * eta ^ 5 + 315242950 * ci ^ 8 * eta ^ 6 - 11223600 * ci ^ 7 * eta ^ 7 - 8231430 * eta ^ 8 * ci ^ 6 - 185562500 * eta * ci ^ 12 + 110368500 * ci ^ 11 * eta ^ 2 + 533043000 * ci ^ 10 * eta ^ 3 - 11457100 * ci ^ 9 * eta ^ 4 + 236289800 * ci ^ 8 * eta ^ 5 - 52909880 * ci ^ 7 * eta ^ 6 - 48611520 * ci ^ 6 * eta ^ 7 + 154800 * ci ^ 5 * eta ^ 8 + 54600 * ci ^ 4 * eta ^ 9 - 43137500 * ci ^ 12 + 62925000 * ci ^ 11 * eta + 518722500 * ci ^ 10 * eta ^ 2 - 133779400 * ci ^ 9 * eta ^ 3 - 127554800 * ci ^ 8 * eta ^ 4 - 65797520 * ci ^ 7 * eta ^ 5 - 116077800 * ci ^ 6 * eta ^ 6 + 3142800 * ci ^ 5 * eta ^ 7 + 1459836 * eta ^ 8 * ci ^ 4 + 15225000 * ci ^ 11 + 281760375 * ci ^ 10 * eta - 143396100 * ci ^ 9 * eta ^ 2 - 355070950 * ci ^ 8 * eta ^ 3 + 5620920 * ci ^ 7 * eta ^ 4 - 81407580 * ci ^ 6 * eta ^ 5 + 15075672 * ci ^ 5 * eta ^ 6 + 9006066 * ci ^ 4 * eta ^ 7 - 18780 * ci ^ 3 * eta ^ 8 - 3213 * ci ^ 2 * eta ^ 9 + 67374375 * ci ^ 10 - 80915000 * ci ^ 9 * eta - 324449150 * ci ^ 8 * eta ^ 2 + 71882480 * ci ^ 7 * eta ^ 3 + 48557580 * ci ^ 6 * eta ^ 4 + 18618288 * ci ^ 5 * eta ^ 5 + 21316590 * ci ^ 4 * eta ^ 6 - 415224 * ci ^ 3 * eta ^ 7 - 121833 * eta ^ 8 * ci ^ 2 - 20199500 * ci ^ 9 - 174470475 * ci ^ 8 * eta + 72956120 * ci ^ 7 * eta ^ 2 + 115207920 * ci ^ 6 * eta ^ 3 - 1465128 * ci ^ 5 * eta ^ 4 + 13488222 * ci ^ 4 * eta ^ 5 - 2035524 * ci ^ 3 * eta ^ 6 - 797392 * ci ^ 2 * eta ^ 7 + 648 * ci * eta ^ 8 + 63 * eta ^ 9 - 42732975 * ci ^ 8 + 40635600 * ci ^ 7 * eta + 96283320 * ci ^ 6 * eta ^ 2 - 18806352 * ci ^ 5 * eta ^ 3 - 10006902 * ci ^ 4 * eta ^ 4 - 2509896 * ci ^ 3 * eta ^ 5 - 1861388 * ci ^ 2 * eta ^ 6 + 20232 * ci * eta ^ 7 + 3795 * eta ^ 8 + 10413600 * ci ^ 7 + 50677770 * ci ^ 6 * eta - 18188088 * ci ^ 5 * eta ^ 2 - 18652098 * ci ^ 4 * eta ^ 3 + 182420 * ci ^ 3 * eta ^ 4 - 968314 * ci ^ 2 * eta ^ 5 + 104884 * ci * eta ^ 6 + 26410 * eta ^ 7 + 12749610 * ci ^ 6 - 9995760 * ci ^ 5 * eta - 13208862 * ci ^ 4 * eta ^ 2 + 2396120 * ci ^ 3 * eta ^ 3 + 1071010 * ci ^ 2 * eta ^ 4 + 131512 * ci * eta ^ 5 + 59338 * eta ^ 6 - 2618280 * ci ^ 5 - 6591126 * ci ^ 4 * eta + 2233140 * ci ^ 3 * eta ^ 2 + 1366232 * ci ^ 2 * eta ^ 3 - 6348 * ci * eta ^ 4 + 17796 * eta ^ 5 - 1741590 * ci ^ 4 + 1214056 * ci ^ 3 * eta + 610948 * ci ^ 2 * eta ^ 2 - 118504 * ci * eta ^ 3 - 47580 * eta ^ 4 + 323800 * ci ^ 3 + 237951 * ci ^ 2 * eta - 109252 * ci * eta ^ 2 - 30354 * eta ^ 3 + 77391 * ci ^ 2 - 59352 * ci * eta + 10878 * eta ^ 2 - 16044 * ci + 11621 * eta + 1793) / L ^ 15 / eta ^ 12 / (1 + eta) ^ 4 / (5 * ci ^ 2 - 1) ^ 3 * R__e ^ 6 * mu ^ 8) / 0.2048e4 + 0.3e1 / 0.128e3 / (eta ^ 8) / (L ^ 11) * (J2 ^ 2) * (25 * ci ^ 4 * eta ^ 3 + 189 * ci ^ 4 * eta ^ 2 + 465 * ci ^ 4 * eta - 20 * ci ^ 3 * eta ^ 2 - 90 * ci ^ 2 * eta ^ 3 + 385 * ci ^ 4 - 144 * ci ^ 3 * eta - 222 * ci ^ 2 * eta ^ 2 - 140 * ci ^ 3 - 162 * ci ^ 2 * eta + 36 * ci * eta ^ 2 + 25 * eta ^ 3 + 90 * ci ^ 2 + 48 * ci * eta + 41 * eta ^ 2 - 20 * ci + 9 * eta - 35) * (R__e ^ 4) * (mu ^ 6) + 0.3e1 / 0.4e1 / (eta ^ 4) / (L ^ 7) * J2 * (3 * ci ^ 2 * eta + 5 * ci ^ 2 - 2 * ci - eta - 1) * (R__e ^ 2) * (mu ^ 4) + (mu ^ 2 / L ^ 3);
ghdot2 = -(J2 ^ 3 * (1316250 * eta ^ 8 * ci ^ 12 + 25110000 * eta ^ 7 * ci ^ 12 - 607500 * ci ^ 11 * eta ^ 8 + 118781250 * eta ^ 6 * ci ^ 12 - 10935000 * ci ^ 11 * eta ^ 7 - 2584125 * eta ^ 8 * ci ^ 10 + 144150000 * eta ^ 5 * ci ^ 12 - 50266500 * ci ^ 11 * eta ^ 6 - 50557500 * ci ^ 10 * eta ^ 7 + 1017000 * ci ^ 9 * eta ^ 8 - 45898750 * eta ^ 4 * ci ^ 12 - 64041000 * ci ^ 11 * eta ^ 5 - 242052000 * ci ^ 10 * eta ^ 6 + 18609000 * ci ^ 9 * eta ^ 7 + 1771875 * eta ^ 8 * ci ^ 8 - 268810000 * eta ^ 3 * ci ^ 12 + 10528500 * ci ^ 11 * eta ^ 4 - 286270500 * ci ^ 10 * eta ^ 5 + 86214900 * ci ^ 9 * eta ^ 6 + 37482300 * ci ^ 8 * eta ^ 7 - 583800 * ci ^ 7 * eta ^ 8 - 307721250 * eta ^ 2 * ci ^ 12 + 97299000 * ci ^ 11 * eta ^ 3 + 71994750 * ci ^ 10 * eta ^ 4 + 108264600 * ci ^ 9 * eta ^ 5 + 184444250 * ci ^ 8 * eta ^ 6 - 11223600 * ci ^ 7 * eta ^ 7 - 535950 * eta ^ 8 * ci ^ 6 - 177950000 * eta * ci ^ 12 + 110368500 * ci ^ 11 * eta ^ 2 + 432565500 * ci ^ 10 * eta ^ 3 - 11457100 * ci ^ 9 * eta ^ 4 + 213054500 * ci ^ 8 * eta ^ 5 - 52909880 * ci ^ 7 * eta ^ 6 - 12895240 * ci ^ 6 * eta ^ 7 + 154800 * ci ^ 5 * eta ^ 8 - 43137500 * ci ^ 12 + 62925000 * ci ^ 11 * eta + 467547000 * ci ^ 10 * eta ^ 2 - 133779400 * ci ^ 9 * eta ^ 3 - 51512700 * ci ^ 8 * eta ^ 4 - 65797520 * ci ^ 7 * eta ^ 5 - 66016480 * ci ^ 6 * eta ^ 6 + 3142800 * ci ^ 5 * eta ^ 7 + 73380 * eta ^ 8 * ci ^ 4 + 15225000 * ci ^ 11 + 268894500 * ci ^ 10 * eta - 143396100 * ci ^ 9 * eta ^ 2 - 286190300 * ci ^ 8 * eta ^ 3 + 5620920 * ci ^ 7 * eta ^ 4 - 73800120 * ci ^ 6 * eta ^ 5 + 15075672 * ci ^ 5 * eta ^ 6 + 2191032 * ci ^ 4 * eta ^ 7 - 18780 * ci ^ 3 * eta ^ 8 + 67374375 * ci ^ 10 - 80915000 * ci ^ 9 * eta - 290894450 * ci ^ 8 * eta ^ 2 + 71882480 * ci ^ 7 * eta ^ 3 + 20600020 * ci ^ 6 * eta ^ 4 + 18618288 * ci ^ 5 * eta ^ 5 + 11860134 * ci ^ 4 * eta ^ 6 - 415224 * ci ^ 3 * eta ^ 7 - 3861 * eta ^ 8 * ci ^ 2 - 20199500 * ci ^ 9 - 165656100 * ci ^ 8 * eta + 72956120 * ci ^ 7 * eta ^ 2 + 92953480 * ci ^ 6 * eta ^ 3 - 1465128 * ci ^ 5 * eta ^ 4 + 12564456 * ci ^ 4 * eta ^ 5 - 2035524 * ci ^ 3 * eta ^ 6 - 173244 * ci ^ 2 * eta ^ 7 + 648 * ci * eta ^ 8 - 42732975 * ci ^ 8 + 40635600 * ci ^ 7 * eta + 86273520 * ci ^ 6 * eta ^ 2 - 18806352 * ci ^ 5 * eta ^ 3 - 4694526 * ci ^ 4 * eta ^ 4 - 2509896 * ci ^ 3 * eta ^ 5 - 1016704 * ci ^ 2 * eta ^ 6 + 20232 * ci * eta ^ 7 + 63 * eta ^ 8 + 10413600 * ci ^ 7 + 47914680 * ci ^ 6 * eta - 18188088 * ci ^ 5 * eta ^ 2 - 15244824 * ci ^ 4 * eta ^ 3 + 182420 * ci ^ 3 * eta ^ 4 - 971812 * ci ^ 2 * eta ^ 5 + 104884 * ci * eta ^ 6 + 4892 * eta ^ 7 + 12749610 * ci ^ 6 - 9995760 * ci ^ 5 * eta - 11934486 * ci ^ 4 * eta ^ 2 + 2396120 * ci ^ 3 * eta ^ 3 + 570606 * ci ^ 2 * eta ^ 4 + 131512 * ci * eta ^ 5 + 31934 * eta ^ 6 - 2618280 * ci ^ 5 - 6204456 * ci ^ 4 * eta + 2233140 * ci ^ 3 * eta ^ 2 + 1158268 * ci ^ 2 * eta ^ 3 - 6348 * ci * eta ^ 4 + 23556 * eta ^ 5 - 1741590 * ci ^ 4 + 1214056 * ci ^ 3 * eta + 577432 * ci ^ 2 * eta ^ 2 - 118504 * ci * eta ^ 3 - 29128 * eta ^ 4 + 323800 * ci ^ 3 + 221220 * ci ^ 2 * eta - 109252 * ci * eta ^ 2 - 28924 * eta ^ 3 + 77391 * ci ^ 2 - 59352 * ci * eta + 7338 * eta ^ 2 - 16044 * ci + 11132 * eta + 1793) / L ^ 15 / eta ^ 12 / (1 + eta) ^ 4 / (5 * ci ^ 2 - 1) ^ 3 * R__e ^ 6 * mu ^ 8) / 0.2048e4 + 0.3e1 / 0.128e3 / (eta ^ 8) / (L ^ 11) * (J2 ^ 2) * (45 * ci ^ 4 * eta ^ 2 + 360 * ci ^ 4 * eta - 20 * ci ^ 3 * eta ^ 2 + 385 * ci ^ 4 - 144 * ci ^ 3 * eta - 126 * ci ^ 2 * eta ^ 2 - 140 * ci ^ 3 - 192 * ci ^ 2 * eta + 36 * ci * eta ^ 2 + 90 * ci ^ 2 + 48 * ci * eta + 25 * eta ^ 2 - 20 * ci + 24 * eta - 35) * (R__e ^ 4) * (mu ^ 6) + 0.3e1 / 0.4e1 / (eta ^ 4) / (L ^ 7) * J2 * (5 * ci ^ 2 - 2 * ci - 1) * (R__e ^ 2) * (mu ^ 4);
hdot2 = (ci * J2 ^ 3 * (151875 * ci ^ 10 * eta ^ 7 + 2581875 * ci ^ 10 * eta ^ 6 + 9984750 * ci ^ 10 * eta ^ 5 - 254250 * ci ^ 8 * eta ^ 7 + 6025500 * ci ^ 10 * eta ^ 4 - 4398000 * ci ^ 8 * eta ^ 6 - 8657625 * ci ^ 10 * eta ^ 3 - 17155725 * ci ^ 8 * eta ^ 5 + 145950 * ci ^ 6 * eta ^ 7 - 15667125 * ci ^ 10 * eta ^ 2 - 9910425 * ci ^ 8 * eta ^ 4 + 2659950 * ci ^ 6 * eta ^ 6 - 11925000 * ci ^ 10 * eta + 12774700 * ci ^ 8 * eta ^ 3 + 10567520 * ci ^ 6 * eta ^ 5 - 38700 * ci ^ 4 * eta ^ 7 - 3806250 * ci ^ 10 + 20670150 * ci ^ 8 * eta ^ 2 + 5881860 * ci ^ 6 * eta ^ 4 - 747000 * ci ^ 4 * eta ^ 6 + 15178875 * ci ^ 8 * eta - 7287090 * ci ^ 6 * eta ^ 3 - 3021918 * ci ^ 4 * eta ^ 5 + 4695 * ci ^ 2 * eta ^ 7 + 5049875 * ci ^ 8 - 10683530 * ci ^ 6 * eta ^ 2 - 1632654 * ci ^ 4 * eta ^ 4 + 99111 * ci ^ 2 * eta ^ 6 - 7555500 * ci ^ 6 * eta + 1998936 * ci ^ 4 * eta ^ 3 + 409770 * ci ^ 2 * eta ^ 5 - 162 * eta ^ 7 - 2603400 * ci ^ 6 + 2702652 * ci ^ 4 * eta ^ 2 + 217704 * ci ^ 2 * eta ^ 4 - 4896 * eta ^ 6 + 1844370 * ci ^ 4 * eta - 263309 * ci ^ 2 * eta ^ 3 - 21325 * eta ^ 5 + 654570 * ci ^ 4 - 335721 * ci ^ 2 * eta ^ 2 - 11553 * eta ^ 4 - 222564 * ci ^ 2 * eta + 13140 * eta ^ 3 - 80950 * ci ^ 2 + 16486 * eta ^ 2 + 10827 * eta + 4011) / L ^ 15 / eta ^ 12 / (1 + eta) ^ 3 / (5 * ci ^ 2 - 1) ^ 3 * R__e ^ 6 * mu ^ 8) / 0.512e3 - 0.3e1 / 0.32e2 * ci * (J2 ^ 2) / (L ^ 11) / (eta ^ 8) * (5 * ci ^ 2 * eta ^ 2 + 36 * ci ^ 2 * eta + 35 * ci ^ 2 - 9 * eta ^ 2 - 12 * eta + 5) * (R__e ^ 4) * (mu ^ 6) - 0.3e1 / 0.2e1 * J2 * (mu ^ 4) * (R__e ^ 2) / (L ^ 7) / (eta ^ 4) * ci;

dlghdotda_J2 = 0.5315625e7 / 0.2048e4 * sqrt(mu) * R__e ^ 6 * a ^ (-0.17e2 / 0.2e1) * ((-0.4930e4 / 0.81e2 + 0.5903e4 / 0.27e2 * eta ^ 5 - 0.3005e4 / 0.21e2 * eta ^ 4 - 0.29243e5 / 0.63e2 * eta ^ 3 - 0.90539e5 / 0.189e3 * eta ^ 2 - 0.148450e6 / 0.567e3 * eta + eta ^ 9 + 0.125e3 / 0.7e1 * eta ^ 8 + 0.59107e5 / 0.567e3 * eta ^ 7 + 0.150481e6 / 0.567e3 * eta ^ 6) * ci ^ 12 - 0.6e1 / 0.7e1 * (eta ^ 7 + 0.17e2 * eta ^ 6 + 0.26626e5 / 0.405e3 * eta ^ 5 + 0.5356e4 / 0.135e3 * eta ^ 4 - 0.23087e5 / 0.405e3 * eta ^ 3 - 0.41779e5 / 0.405e3 * eta ^ 2 - 0.2120e4 / 0.27e2 * eta - 0.2030e4 / 0.81e2) * (0.1e1 + eta) * ci ^ 11 + (0.35933e5 / 0.378e3 - 0.530858e6 / 0.945e3 * eta ^ 6 - 0.418771e6 / 0.945e3 * eta ^ 5 + 0.231103e6 / 0.945e3 * eta ^ 4 + 0.33844e5 / 0.45e2 * eta ^ 3 + 0.138326e6 / 0.189e3 * eta ^ 2 + 0.751361e6 / 0.1890e4 * eta - 0.199e3 / 0.90e2 * eta ^ 9 - 0.24989e5 / 0.630e3 * eta ^ 8 - 0.215048e6 / 0.945e3 * eta ^ 7) * ci ^ 10 + 0.452e3 / 0.315e3 * (eta ^ 7 + 0.5864e4 / 0.339e3 * eta ^ 6 + 0.228743e6 / 0.3390e4 * eta ^ 5 + 0.132139e6 / 0.3390e4 * eta ^ 4 - 0.255494e6 / 0.5085e4 * eta ^ 3 - 0.137801e6 / 0.1695e4 * eta ^ 2 - 0.40477e5 / 0.678e3 * eta - 0.40399e5 / 0.2034e4) * (0.1e1 + eta) * ci ^ 9 + (-0.569773e6 / 0.9450e4 + 0.302053e6 / 0.9450e4 * eta ^ 8 + 0.372721e6 / 0.2025e4 * eta ^ 7 + 0.6304859e7 / 0.14175e5 * eta ^ 6 + 0.4725796e7 / 0.14175e5 * eta ^ 5 - 0.2551096e7 / 0.14175e5 * eta ^ 4 - 0.7101419e7 / 0.14175e5 * eta ^ 3 - 0.6488983e7 / 0.14175e5 * eta ^ 2 - 0.2326273e7 / 0.9450e4 * eta + 0.2263e4 / 0.1350e4 * eta ^ 9) * ci ^ 8 - 0.556e3 / 0.675e3 * (eta ^ 7 + 0.17733e5 / 0.973e3 * eta ^ 6 + 0.1056752e7 / 0.14595e5 * eta ^ 5 + 0.196062e6 / 0.4865e4 * eta ^ 4 - 0.242903e6 / 0.4865e4 * eta ^ 3 - 0.1068353e7 / 0.14595e5 * eta ^ 2 - 0.50370e5 / 0.973e3 * eta - 0.17356e5 / 0.973e3) * (0.1e1 + eta) * ci ^ 7 + (0.424987e6 / 0.23625e5 - 0.121e3 / 0.225e3 * eta ^ 9 - 0.274381e6 / 0.23625e5 * eta ^ 8 - 0.540128e6 / 0.7875e4 * eta ^ 7 - 0.773852e6 / 0.4725e4 * eta ^ 6 - 0.2713586e7 / 0.23625e5 * eta ^ 5 + 0.1618586e7 / 0.23625e5 * eta ^ 4 + 0.142232e6 / 0.875e3 * eta ^ 3 + 0.458492e6 / 0.3375e4 * eta ^ 2 + 0.1689259e7 / 0.23625e5 * eta) * ci ^ 6 + 0.344e3 / 0.1575e4 * (eta ^ 7 + 0.830e3 / 0.43e2 * eta ^ 6 + 0.503653e6 / 0.6450e4 * eta ^ 5 + 0.90703e5 / 0.2150e4 * eta ^ 4 - 0.55526e5 / 0.1075e4 * eta ^ 3 - 0.225221e6 / 0.3225e4 * eta ^ 2 - 0.20493e5 / 0.430e3 * eta - 0.7273e4 / 0.430e3) * (0.1e1 + eta) * ci ^ 5 + (-0.19351e5 / 0.7875e4 + 0.52e2 / 0.675e3 * eta ^ 9 + 0.3862e4 / 0.1875e4 * eta ^ 8 + 0.55593e5 / 0.4375e4 * eta ^ 7 + 0.236851e6 / 0.7875e4 * eta ^ 6 + 0.2248037e7 / 0.118125e6 * eta ^ 5 - 0.61771e5 / 0.4375e4 * eta ^ 4 - 0.3108683e7 / 0.118125e6 * eta ^ 3 - 0.2201477e7 / 0.118125e6 * eta ^ 2 - 0.1098521e7 / 0.118125e6 * eta) * ci ^ 4 - 0.626e3 / 0.23625e5 * (eta ^ 7 + 0.33037e5 / 0.1565e4 * eta ^ 6 + 0.27318e5 / 0.313e3 * eta ^ 5 + 0.72568e5 / 0.1565e4 * eta ^ 4 - 0.263309e6 / 0.4695e4 * eta ^ 3 - 0.111907e6 / 0.1565e4 * eta ^ 2 - 0.74188e5 / 0.1565e4 * eta - 0.16190e5 / 0.939e3) * (0.1e1 + eta) * ci ^ 3 + (0.8599e4 / 0.78750e5 - 0.17e2 / 0.3750e4 * eta ^ 9 - 0.13537e5 / 0.78750e5 * eta ^ 8 - 0.398696e6 / 0.354375e6 * eta ^ 7 - 0.930694e6 / 0.354375e6 * eta ^ 6 - 0.484157e6 / 0.354375e6 * eta ^ 5 + 0.107101e6 / 0.70875e5 * eta ^ 4 + 0.97588e5 / 0.50625e5 * eta ^ 3 + 0.305474e6 / 0.354375e6 * eta ^ 2 + 0.1259e4 / 0.3750e4 * eta) * ci ^ 2 + 0.4e1 / 0.4375e4 * (eta ^ 7 + 0.272e3 / 0.9e1 * eta ^ 6 + 0.21325e5 / 0.162e3 * eta ^ 5 + 0.3851e4 / 0.54e2 * eta ^ 4 - 0.730e3 / 0.9e1 * eta ^ 3 - 0.8243e4 / 0.81e2 * eta ^ 2 - 0.401e3 / 0.6e1 * eta - 0.1337e4 / 0.54e2) * (0.1e1 + eta) * ci + 0.29669e5 / 0.354375e6 * eta ^ 6 + 0.253e3 / 0.47250e5 * eta ^ 8 + 0.11621e5 / 0.708750e6 * eta - 0.1586e4 / 0.23625e5 * eta ^ 4 + eta ^ 9 / 0.11250e5 + 0.2641e4 / 0.70875e5 * eta ^ 7 + 0.2966e4 / 0.118125e6 * eta ^ 5 - 0.5059e4 / 0.118125e6 * eta ^ 3 + 0.259e3 / 0.16875e5 * eta ^ 2 + 0.1793e4 / 0.708750e6) * J2 ^ 3 / eta ^ 12 / (0.5e1 * ci ^ 2 - 0.1e1) ^ 3 / (0.1e1 + eta) ^ 4 - 0.825e3 / 0.256e3 * sqrt(mu) * ((eta ^ 3 + 0.189e3 / 0.25e2 * eta ^ 2 + 0.93e2 / 0.5e1 * eta + 0.77e2 / 0.5e1) * ci ^ 4 + (-0.4e1 / 0.5e1 * eta ^ 2 - 0.144e3 / 0.25e2 * eta - 0.28e2 / 0.5e1) * ci ^ 3 + (-0.18e2 / 0.5e1 * eta ^ 3 - 0.222e3 / 0.25e2 * eta ^ 2 - 0.162e3 / 0.25e2 * eta + 0.18e2 / 0.5e1) * ci ^ 2 + (0.36e2 / 0.25e2 * eta ^ 2 + 0.48e2 / 0.25e2 * eta - 0.4e1 / 0.5e1) * ci + eta ^ 3 + 0.41e2 / 0.25e2 * eta ^ 2 + 0.9e1 / 0.25e2 * eta - 0.7e1 / 0.5e1) * R__e ^ 4 * a ^ (-0.13e2 / 0.2e1) * J2 ^ 2 / eta ^ 8 - 0.21e2 / 0.8e1 * (0.3e1 * ci ^ 2 * eta + 0.5e1 * ci ^ 2 - 0.2e1 * ci - eta - 0.1e1) * a ^ (-0.9e1 / 0.2e1) / eta ^ 4 * J2 * sqrt(mu) * R__e ^ 2 - 0.3e1 / 0.2e1 * a ^ (-0.5e1 / 0.2e1) * sqrt(mu);

dlghdotdp12byp12 = 0.5315625e7 / 0.512e3 * ((ci + 1) ^ 2) * a ^ (-0.15e2 / 0.2e1) * ((0.5903e4 / 0.27e2 * eta ^ 5 - 0.29243e5 / 0.63e2 * eta ^ 3 - 0.90539e5 / 0.189e3 * eta ^ 2 + eta ^ 9 + 0.125e3 / 0.7e1 * eta ^ 8 + 0.59107e5 / 0.567e3 * eta ^ 7 - 0.3005e4 / 0.21e2 * eta ^ 4 - 0.148450e6 / 0.567e3 * eta + 0.150481e6 / 0.567e3 * eta ^ 6 - 0.4930e4 / 0.81e2) * (ci ^ 13) - 0.5e1 / 0.7e1 * (0.1e1 + eta) * (eta ^ 7 + 0.17e2 * eta ^ 6 + 0.26626e5 / 0.405e3 * eta ^ 5 + 0.5356e4 / 0.135e3 * eta ^ 4 - 0.23087e5 / 0.405e3 * eta ^ 3 - 0.41779e5 / 0.405e3 * eta ^ 2 - 0.2120e4 / 0.27e2 * eta - 0.2030e4 / 0.81e2) * (ci ^ 12) + (0.16579e5 / 0.189e3 - 0.5222e4 / 0.27e2 * eta ^ 7 - 0.1085468e7 / 0.2835e4 * eta ^ 5 + 0.624476e6 / 0.2835e4 * eta ^ 4 + 0.649274e6 / 0.945e3 * eta ^ 3 - 0.253e3 / 0.135e3 * eta ^ 9 - 0.31739e5 / 0.945e3 * eta ^ 8 + 0.1048261e7 / 0.2835e4 * eta - 0.454226e6 / 0.945e3 * eta ^ 6 + 0.1926494e7 / 0.2835e4 * eta ^ 2) * (ci ^ 11) + 0.65e2 / 0.63e2 * (eta ^ 7 + 0.16777e5 / 0.975e3 * eta ^ 6 + 0.195823e6 / 0.2925e4 * eta ^ 5 + 0.12737e5 / 0.325e3 * eta ^ 4 - 0.152989e6 / 0.2925e4 * eta ^ 3 - 0.10291e5 / 0.117e3 * eta ^ 2 - 0.63797e5 / 0.975e3 * eta - 0.62729e5 / 0.2925e4) * (0.1e1 + eta) * (ci ^ 10) + (-0.17762279e8 / 0.42525e5 * eta ^ 3 - 0.3041539e7 / 0.14175e5 * eta + 0.2624e4 / 0.2025e4 * eta ^ 9 - 0.27187e5 / 0.525e3 + 0.338444e6 / 0.14175e5 * eta ^ 8 + 0.2038247e7 / 0.6075e4 * eta ^ 6 - 0.859663e6 / 0.6075e4 * eta ^ 4 - 0.16863433e8 / 0.42525e5 * eta ^ 2 + 0.5834767e7 / 0.42525e5 * eta ^ 7 + 0.11007361e8 / 0.42525e5 * eta ^ 5) * (ci ^ 9) - 0.8048e4 / 0.14175e5 * (eta ^ 7 + 0.70509e5 / 0.4024e4 * eta ^ 6 + 0.1657913e7 / 0.24144e5 * eta ^ 5 + 0.316275e6 / 0.8048e4 * eta ^ 4 - 0.201877e6 / 0.4024e4 * eta ^ 3 - 0.239449e6 / 0.3018e4 * eta ^ 2 - 0.465033e6 / 0.8048e4 * eta - 0.155909e6 / 0.8048e4) * (0.1e1 + eta) * (ci ^ 8) + (0.28405676e8 / 0.212625e6 * eta ^ 3 + 0.25955932e8 / 0.212625e6 * eta ^ 2 - 0.18903184e8 / 0.212625e6 * eta ^ 5 + 0.4652546e7 / 0.70875e5 * eta - 0.4526e4 / 0.10125e5 * eta ^ 9 - 0.604106e6 / 0.70875e5 * eta ^ 8 - 0.25219436e8 / 0.212625e6 * eta ^ 6 + 0.10204384e8 / 0.212625e6 * eta ^ 4 - 0.1490884e7 / 0.30375e5 * eta ^ 7 + 0.1139546e7 / 0.70875e5) * (ci ^ 7) + 0.11042e5 / 0.70875e5 * (0.1e1 + eta) * (eta ^ 7 + 0.99231e5 / 0.5521e4 * eta ^ 6 + 0.1177261e7 / 0.16563e5 * eta ^ 5 + 0.220065e6 / 0.5521e4 * eta ^ 4 - 0.273433e6 / 0.5521e4 * eta ^ 3 - 0.1225429e7 / 0.16563e5 * eta ^ 2 - 0.291111e6 / 0.5521e4 * eta - 0.99673e5 / 0.5521e4) * (ci ^ 6) + (-0.1062647e7 / 0.118125e6 * eta ^ 4 + 0.537001e6 / 0.23625e5 * eta ^ 6 - 0.328232e6 / 0.118125e6 + 0.373349e6 / 0.39375e5 * eta ^ 7 + 0.5892721e7 / 0.354375e6 * eta ^ 5 - 0.8412109e7 / 0.354375e6 * eta ^ 3 - 0.3969256e7 / 0.354375e6 * eta + 0.829e3 / 0.10125e5 * eta ^ 9 + 0.193279e6 / 0.118125e6 * eta ^ 8 - 0.1485371e7 / 0.70875e5 * eta ^ 2) * (ci ^ 5) - 0.547e3 / 0.23625e5 * (0.1e1 + eta) * (eta ^ 7 + 0.49963e5 / 0.2735e4 * eta ^ 6 + 0.597536e6 / 0.8205e4 * eta ^ 5 + 0.108838e6 / 0.2735e4 * eta ^ 4 - 0.403003e6 / 0.8205e4 * eta ^ 3 - 0.565163e6 / 0.8205e4 * eta ^ 2 - 0.130742e6 / 0.2735e4 * eta - 0.27448e5 / 0.1641e4) * (ci ^ 4) + (-0.94519e5 / 0.590625e6 * eta ^ 8 - 0.5019106e7 / 0.5315625e7 * eta ^ 7 - 0.367e3 / 0.50625e5 * eta ^ 9 + 0.4651852e7 / 0.5315625e7 * eta ^ 4 - 0.480386e6 / 0.212625e6 * eta ^ 6 + 0.11820938e8 / 0.5315625e7 * eta ^ 3 + 0.1800457e7 / 0.1771875e7 * eta - 0.1235236e7 / 0.759375e6 * eta ^ 5 + 0.10154122e8 / 0.5315625e7 * eta ^ 2 + 0.30103e5 / 0.118125e6) * (ci ^ 3) + 0.223e3 / 0.118125e6 * (0.1e1 + eta) * (eta ^ 7 + 0.19437e5 / 0.1115e4 * eta ^ 6 + 0.139237e6 / 0.2007e4 * eta ^ 5 + 0.121429e6 / 0.3345e4 * eta ^ 4 - 0.153809e6 / 0.3345e4 * eta ^ 3 - 0.595013e6 / 0.10035e5 * eta ^ 2 - 0.44113e5 / 0.1115e4 * eta - 0.9505e4 / 0.669e3) * (ci ^ 2) + (-0.2483e4 / 0.253125e6 + 0.1202e4 / 0.196875e6 * eta ^ 8 + 0.200621e6 / 0.5315625e7 * eta ^ 7 + 0.2e1 / 0.9375e4 * eta ^ 9 - 0.35731e5 / 0.1063125e7 * eta ^ 4 + 0.485659e6 / 0.5315625e7 * eta ^ 6 - 0.68711e5 / 0.1771875e7 * eta + 0.350687e6 / 0.5315625e7 * eta ^ 5 - 0.455461e6 / 0.5315625e7 * eta ^ 3 - 0.387059e6 / 0.5315625e7 * eta ^ 2) * ci - 0.2e1 / 0.65625e5 * (0.1e1 + eta) * (eta ^ 7 + 0.272e3 / 0.9e1 * eta ^ 6 + 0.21325e5 / 0.162e3 * eta ^ 5 + 0.3851e4 / 0.54e2 * eta ^ 4 - 0.730e3 / 0.9e1 * eta ^ 3 - 0.8243e4 / 0.81e2 * eta ^ 2 - 0.401e3 / 0.6e1 * eta - 0.1337e4 / 0.54e2)) * J2 ^ 3 * R__e ^ 6 * sqrt(mu) / eta ^ 12 / ((5 * ci ^ 2 - 1) ^ 4) / (0.1e1 + eta) ^ 4 - 0.75e2 / 0.32e2 * ((ci + 1) ^ 2) * a ^ (-0.11e2 / 0.2e1) * J2 ^ 2 * R__e ^ 4 * sqrt(mu) * (((ci ^ 3) - 0.9e1 / 0.5e1 * ci) * eta ^ 3 + (0.189e3 / 0.25e2 * (ci ^ 3) - 0.3e1 / 0.5e1 * (ci ^ 2) - 0.111e3 / 0.25e2 * ci + 0.9e1 / 0.25e2) * eta ^ 2 + (0.93e2 / 0.5e1 * (ci ^ 3) - 0.108e3 / 0.25e2 * (ci ^ 2) - 0.81e2 / 0.25e2 * ci + 0.12e2 / 0.25e2) * eta + 0.77e2 / 0.5e1 * (ci ^ 3) - 0.21e2 / 0.5e1 * (ci ^ 2) + 0.9e1 / 0.5e1 * ci - 0.1e1 / 0.5e1) / eta ^ 8 - 0.3e1 / 0.2e1 * ((ci + 1) ^ 2) * (0.3e1 * ci * eta + (5 * ci) - 0.1e1) * a ^ (-0.7e1 / 0.2e1) * J2 * sqrt(mu) * R__e ^ 2 / eta ^ 4;
dlghdotdp1_J2 = dlghdotdp12byp12*p1;
dlghdotdp2_J2 = dlghdotdp12byp12*p2;

dlghdotdq12byq12 = -0.2480625e7 / 0.1024e4 * a ^ (-0.15e2 / 0.2e1) * ((eta ^ 10 + 0.1021e4 / 0.49e2 * eta ^ 9 + 0.63607e5 / 0.441e3 * eta ^ 8 + 0.600115e6 / 0.1323e4 * eta ^ 7 + 0.251831e6 / 0.441e3 * eta ^ 6 - 0.35293e5 / 0.1323e4 * eta ^ 5 - 0.452279e6 / 0.441e3 * eta ^ 4 - 0.2057107e7 / 0.1323e4 * eta ^ 3 - 0.1647640e7 / 0.1323e4 * eta ^ 2 - 0.80930e5 / 0.147e3 * eta - 0.19720e5 / 0.189e3) * ci ^ 12 - 0.48e2 / 0.49e2 * (0.1e1 + eta) * (eta ^ 8 + 0.79e2 / 0.4e1 * eta ^ 7 + 0.30757e5 / 0.324e3 * eta ^ 6 + 0.36313e5 / 0.324e3 * eta ^ 5 - 0.275e3 / 0.6e1 * eta ^ 4 - 0.75091e5 / 0.324e3 * eta ^ 3 - 0.86299e5 / 0.324e3 * eta ^ 2 - 0.16735e5 / 0.108e3 * eta - 0.1015e4 / 0.27e2) * ci ^ 11 + (-0.199e3 / 0.90e2 * eta ^ 10 - 0.204091e6 / 0.4410e4 * eta ^ 9 + 0.227339e6 / 0.270e3 * eta - 0.695122e6 / 0.2205e4 * eta ^ 8 - 0.425588e6 / 0.441e3 * eta ^ 7 - 0.7791629e7 / 0.6615e4 * eta ^ 6 - 0.158161e6 / 0.6615e4 * eta ^ 5 + 0.11088236e8 / 0.6615e4 * eta ^ 4 + 0.2297048e7 / 0.945e3 * eta ^ 3 + 0.717229e6 / 0.378e3 * eta ^ 2 + 0.71866e5 / 0.441e3) * ci ^ 10 + 0.3616e4 / 0.2205e4 * (eta ^ 8 + 0.18157e5 / 0.904e3 * eta ^ 7 + 0.263927e6 / 0.2712e4 * eta ^ 6 + 0.305473e6 / 0.2712e4 * eta ^ 5 - 0.24671e5 / 0.678e3 * eta ^ 4 - 0.511579e6 / 0.2712e4 * eta ^ 3 - 0.558941e6 / 0.2712e4 * eta ^ 2 - 0.323621e6 / 0.2712e4 * eta - 0.40399e5 / 0.1356e4) * (0.1e1 + eta) * ci ^ 9 + (0.2263e4 / 0.1350e4 * eta ^ 10 + 0.2463947e7 / 0.66150e5 * eta ^ 9 - 0.154758533e9 / 0.99225e5 * eta ^ 3 - 0.46892389e8 / 0.39690e5 * eta ^ 2 - 0.11568457e8 / 0.22050e5 * eta + 0.8431247e7 / 0.33075e5 * eta ^ 8 + 0.3043753e7 / 0.3969e4 * eta ^ 7 + 0.17962582e8 / 0.19845e5 * eta ^ 6 + 0.493484e6 / 0.19845e5 * eta ^ 5 - 0.22545443e8 / 0.19845e5 * eta ^ 4 - 0.1139546e7 / 0.11025e5) * ci ^ 8 - 0.4448e4 / 0.4725e4 * (0.1e1 + eta) * (eta ^ 8 + 0.82231e5 / 0.3892e4 * eta ^ 7 + 0.1216349e7 / 0.11676e5 * eta ^ 6 + 0.1386731e7 / 0.11676e5 * eta ^ 5 - 0.67317e5 / 0.1946e4 * eta ^ 4 - 0.2044697e7 / 0.11676e5 * eta ^ 3 - 0.2126123e7 / 0.11676e5 * eta ^ 2 - 0.407205e6 / 0.3892e4 * eta - 0.26034e5 / 0.973e3) * ci ^ 7 + (-0.121e3 / 0.225e3 * eta ^ 10 - 0.2233163e7 / 0.165375e6 * eta ^ 9 + 0.85586e5 / 0.33075e5 * eta ^ 5 + 0.12574424e8 / 0.33075e5 * eta ^ 4 + 0.79494592e8 / 0.165375e6 * eta ^ 3 + 0.2297333e7 / 0.6615e4 * eta ^ 2 + 0.8460547e7 / 0.55125e5 * eta - 0.64004e5 / 0.675e3 * eta ^ 8 - 0.9358904e7 / 0.33075e5 * eta ^ 7 - 0.53065006e8 / 0.165375e6 * eta ^ 6 + 0.1699948e7 / 0.55125e5) * ci ^ 6 + 0.2752e4 / 0.11025e5 * (0.1e1 + eta) * (eta ^ 8 + 0.7685e4 / 0.344e3 * eta ^ 7 + 0.578353e6 / 0.5160e4 * eta ^ 6 + 0.651877e6 / 0.5160e4 * eta ^ 5 - 0.3035e4 / 0.86e2 * eta ^ 4 - 0.177083e6 / 0.1032e4 * eta ^ 3 - 0.176159e6 / 0.1032e4 * eta ^ 2 - 0.167259e6 / 0.1720e4 * eta - 0.21819e5 / 0.860e3) * ci ^ 5 + (0.52e2 / 0.675e3 * eta ^ 10 + 0.93988e5 / 0.39375e5 * eta ^ 9 + 0.46044997e8 / 0.826875e6 * eta ^ 6 - 0.855509e6 / 0.165375e6 * eta ^ 5 - 0.1535869e7 / 0.23625e5 * eta ^ 4 - 0.2351953e7 / 0.33075e5 * eta ^ 3 - 0.7698517e7 / 0.165375e6 * eta ^ 2 - 0.16727971e8 / 0.826875e6 * eta + 0.1609147e7 / 0.91875e5 * eta ^ 8 + 0.2868847e7 / 0.55125e5 * eta ^ 7 - 0.77404e5 / 0.18375e5) * ci ^ 4 - 0.5008e4 / 0.165375e6 * (eta ^ 8 + 0.152579e6 / 0.6260e4 * eta ^ 7 + 0.782061e6 / 0.6260e4 * eta ^ 6 + 0.877189e6 / 0.6260e4 * eta ^ 5 - 0.118173e6 / 0.3130e4 * eta ^ 4 - 0.1122359e7 / 0.6260e4 * eta ^ 3 - 0.1078851e7 / 0.6260e4 * eta ^ 2 - 0.610409e6 / 0.6260e4 * eta - 0.8095e4 / 0.313e3) * (0.1e1 + eta) * ci ^ 3 + (-0.17e2 / 0.3750e4 * eta ^ 10 - 0.109367e6 / 0.551250e6 * eta ^ 9 + 0.17198e5 / 0.91875e5 - 0.85154e5 / 0.55125e5 * eta ^ 8 - 0.2260084e7 / 0.496125e6 * eta ^ 7 - 0.10909891e8 / 0.2480625e7 * eta ^ 6 + 0.3036961e7 / 0.2480625e7 * eta ^ 5 + 0.13164548e8 / 0.2480625e7 * eta ^ 4 + 0.297848e6 / 0.70875e5 * eta ^ 3 + 0.1935749e7 / 0.992250e6 * eta ^ 2 + 0.428413e6 / 0.551250e6 * eta) * ci ^ 2 + 0.32e2 / 0.30625e5 * (eta ^ 8 + 0.277e3 / 0.8e1 * eta ^ 7 + 0.121313e6 / 0.648e3 * eta ^ 6 + 0.138179e6 / 0.648e3 * eta ^ 5 - 0.2719e4 / 0.54e2 * eta ^ 4 - 0.166289e6 / 0.648e3 * eta ^ 3 - 0.158219e6 / 0.648e3 * eta ^ 2 - 0.9959e4 / 0.72e2 * eta - 0.1337e4 / 0.36e2) * (0.1e1 + eta) * ci + 0.8429e4 / 0.165375e6 * eta ^ 8 + 0.5797e4 / 0.183750e6 * eta + 0.3586e4 / 0.826875e6 + eta ^ 10 / 0.11250e5 + 0.10183e5 / 0.1653750e7 * eta ^ 9 + 0.24181e5 / 0.165375e6 * eta ^ 7 + 0.91964e5 / 0.826875e6 * eta ^ 6 - 0.74398e5 / 0.826875e6 * eta ^ 5 - 0.43069e5 / 0.275625e6 * eta ^ 4 - 0.20149e5 / 0.826875e6 * eta ^ 3 + 0.699e3 / 0.12250e5 * eta ^ 2) * J2 ^ 3 * R__e ^ 6 * sqrt(mu) / eta ^ 14 / (0.5e1 * ci ^ 2 - 0.1e1) ^ 3 / (0.1e1 + eta) ^ 5 + (0.3e1 * (0.125e3 * eta ^ 3 + 0.1134e4 * eta ^ 2 + 0.3255e4 * eta + 0.3080e4) * J2 ^ 2 * R__e ^ 4 * ci ^ 4 - 0.360e3 * (eta ^ 2 + 0.42e2 / 0.5e1 * eta + 0.28e2 / 0.3e1) * J2 ^ 2 * R__e ^ 4 * ci ^ 3 + 0.3e1 * (-0.450e3 * eta ^ 3 - 0.1332e4 * eta ^ 2 - 0.1134e4 * eta + 0.720e3) * J2 ^ 2 * R__e ^ 4 * ci ^ 2 + 0.3e1 * (0.216e3 * eta ^ 2 + 0.336e3 * eta - 0.160e3) * J2 ^ 2 * R__e ^ 4 * ci + 0.3e1 * (0.125e3 * eta ^ 3 + 0.246e3 * eta ^ 2 + 0.63e2 * eta - 0.280e3) * J2 ^ 2 * R__e ^ 4) * sqrt(mu) * a ^ (-0.11e2 / 0.2e1) / eta ^ 10 / 0.128e3 + 0.3e1 / 0.4e1 * (0.9e1 * ci ^ 2 * eta + 0.20e2 * ci ^ 2 - 0.8e1 * ci - 0.3e1 * eta - 0.4e1) * a ^ (-0.7e1 / 0.2e1) / eta ^ 6 * J2 * sqrt(mu) * R__e ^ 2;
dlghdotdq1_J2 = dlghdotdq12byq12*q1;
dlghdotdq2_J2 = dlghdotdq12byq12*q2;

dghdotda_J2 = 0.9871875e7 / 0.2048e4 * sqrt(mu) * R__e ^ 6 * a ^ (-0.17e2 / 0.2e1) * ((-0.34510e5 / 0.1053e4 + 0.38440e5 / 0.351e3 * eta ^ 5 - 0.36719e5 / 0.1053e4 * eta ^ 4 - 0.215048e6 / 0.1053e4 * eta ^ 3 - 0.27353e5 / 0.117e3 * eta ^ 2 - 0.142360e6 / 0.1053e4 * eta + eta ^ 8 + 0.248e3 / 0.13e2 * eta ^ 7 + 0.31675e5 / 0.351e3 * eta ^ 6) * ci ^ 12 - 0.6e1 / 0.13e2 * (eta ^ 7 + 0.17e2 * eta ^ 6 + 0.26626e5 / 0.405e3 * eta ^ 5 + 0.5356e4 / 0.135e3 * eta ^ 4 - 0.23087e5 / 0.405e3 * eta ^ 3 - 0.41779e5 / 0.405e3 * eta ^ 2 - 0.2120e4 / 0.27e2 * eta - 0.2030e4 / 0.81e2) * (0.1e1 + eta) * ci ^ 11 + (0.35933e5 / 0.702e3 - 0.322736e6 / 0.1755e4 * eta ^ 6 - 0.381694e6 / 0.1755e4 * eta ^ 5 + 0.95993e5 / 0.1755e4 * eta ^ 4 + 0.576754e6 / 0.1755e4 * eta ^ 3 + 0.623396e6 / 0.1755e4 * eta ^ 2 + 0.358526e6 / 0.1755e4 * eta - 0.2297e4 / 0.1170e4 * eta ^ 8 - 0.1498e4 / 0.39e2 * eta ^ 7) * ci ^ 10 + 0.452e3 / 0.585e3 * (eta ^ 7 + 0.5864e4 / 0.339e3 * eta ^ 6 + 0.228743e6 / 0.3390e4 * eta ^ 5 + 0.132139e6 / 0.3390e4 * eta ^ 4 - 0.255494e6 / 0.5085e4 * eta ^ 3 - 0.137801e6 / 0.1695e4 * eta ^ 2 - 0.40477e5 / 0.678e3 * eta - 0.40399e5 / 0.2034e4) * (0.1e1 + eta) * ci ^ 9 + (-0.569773e6 / 0.17550e5 + 0.35e2 / 0.26e2 * eta ^ 8 + 0.83294e5 / 0.2925e4 * eta ^ 7 + 0.737777e6 / 0.5265e4 * eta ^ 6 + 0.852218e6 / 0.5265e4 * eta ^ 5 - 0.343418e6 / 0.8775e4 * eta ^ 4 - 0.5723806e7 / 0.26325e5 * eta ^ 3 - 0.5817889e7 / 0.26325e5 * eta ^ 2 - 0.1104374e7 / 0.8775e4 * eta) * ci ^ 8 - 0.3892e4 / 0.8775e4 * (eta ^ 7 + 0.17733e5 / 0.973e3 * eta ^ 6 + 0.1056752e7 / 0.14595e5 * eta ^ 5 + 0.196062e6 / 0.4865e4 * eta ^ 4 - 0.242903e6 / 0.4865e4 * eta ^ 3 - 0.1068353e7 / 0.14595e5 * eta ^ 2 - 0.50370e5 / 0.973e3 * eta - 0.17356e5 / 0.973e3) * (0.1e1 + eta) * ci ^ 7 + (0.424987e6 / 0.43875e5 - 0.397e3 / 0.975e3 * eta ^ 8 - 0.1289524e7 / 0.131625e6 * eta ^ 7 - 0.6601648e7 / 0.131625e6 * eta ^ 6 - 0.2460004e7 / 0.43875e5 * eta ^ 5 + 0.2060002e7 / 0.131625e6 * eta ^ 4 + 0.9295348e7 / 0.131625e6 * eta ^ 3 + 0.2875784e7 / 0.43875e5 * eta ^ 2 + 0.1597156e7 / 0.43875e5 * eta) * ci ^ 6 + 0.344e3 / 0.2925e4 * (eta ^ 7 + 0.830e3 / 0.43e2 * eta ^ 6 + 0.503653e6 / 0.6450e4 * eta ^ 5 + 0.90703e5 / 0.2150e4 * eta ^ 4 - 0.55526e5 / 0.1075e4 * eta ^ 3 - 0.225221e6 / 0.3225e4 * eta ^ 2 - 0.20493e5 / 0.430e3 * eta - 0.7273e4 / 0.430e3) * (0.1e1 + eta) * ci ^ 5 + (-0.19351e5 / 0.14625e5 + 0.2446e4 / 0.43875e5 * eta ^ 8 + 0.121724e6 / 0.73125e5 * eta ^ 7 + 0.152053e6 / 0.16875e5 * eta ^ 6 + 0.2094076e7 / 0.219375e6 * eta ^ 5 - 0.260807e6 / 0.73125e5 * eta ^ 4 - 0.2540804e7 / 0.219375e6 * eta ^ 3 - 0.221009e6 / 0.24375e5 * eta ^ 2 - 0.344692e6 / 0.73125e5 * eta) * ci ^ 4 - 0.626e3 / 0.43875e5 * (eta ^ 7 + 0.33037e5 / 0.1565e4 * eta ^ 6 + 0.27318e5 / 0.313e3 * eta ^ 5 + 0.72568e5 / 0.1565e4 * eta ^ 4 - 0.263309e6 / 0.4695e4 * eta ^ 3 - 0.111907e6 / 0.1565e4 * eta ^ 2 - 0.74188e5 / 0.1565e4 * eta - 0.16190e5 / 0.939e3) * (0.1e1 + eta) * ci ^ 3 + (0.8599e4 / 0.146250e6 - 0.11e2 / 0.3750e4 * eta ^ 8 - 0.28874e5 / 0.219375e6 * eta ^ 7 - 0.39104e5 / 0.50625e5 * eta ^ 6 - 0.485906e6 / 0.658125e6 * eta ^ 5 + 0.95101e5 / 0.219375e6 * eta ^ 4 + 0.579134e6 / 0.658125e6 * eta ^ 3 + 0.288716e6 / 0.658125e6 * eta ^ 2 + 0.2458e4 / 0.14625e5 * eta) * ci ^ 2 + 0.4e1 / 0.8125e4 * (eta ^ 7 + 0.272e3 / 0.9e1 * eta ^ 6 + 0.21325e5 / 0.162e3 * eta ^ 5 + 0.3851e4 / 0.54e2 * eta ^ 4 - 0.730e3 / 0.9e1 * eta ^ 3 - 0.8243e4 / 0.81e2 * eta ^ 2 - 0.401e3 / 0.6e1 * eta - 0.1337e4 / 0.54e2) * (0.1e1 + eta) * ci + 0.15967e5 / 0.658125e6 * eta ^ 6 + 0.7e1 / 0.146250e6 * eta ^ 8 + 0.5566e4 / 0.658125e6 * eta - 0.14564e5 / 0.658125e6 * eta ^ 4 + 0.2446e4 / 0.658125e6 * eta ^ 7 + 0.302e3 / 0.16875e5 * eta ^ 5 - 0.14462e5 / 0.658125e6 * eta ^ 3 + 0.1223e4 / 0.219375e6 * eta ^ 2 + 0.1793e4 / 0.1316250e7) * J2 ^ 3 / eta ^ 12 / (0.5e1 * ci ^ 2 - 0.1e1) ^ 3 / (0.1e1 + eta) ^ 4 - 0.1485e4 / 0.256e3 * sqrt(mu) * ((eta ^ 2 + 0.8e1 * eta + 0.77e2 / 0.9e1) * ci ^ 4 + (-0.4e1 / 0.9e1 * eta ^ 2 - 0.16e2 / 0.5e1 * eta - 0.28e2 / 0.9e1) * ci ^ 3 + (-0.14e2 / 0.5e1 * eta ^ 2 - 0.64e2 / 0.15e2 * eta + 0.2e1) * ci ^ 2 + (0.4e1 / 0.5e1 * eta ^ 2 + 0.16e2 / 0.15e2 * eta - 0.4e1 / 0.9e1) * ci + 0.5e1 / 0.9e1 * eta ^ 2 + 0.8e1 / 0.15e2 * eta - 0.7e1 / 0.9e1) * R__e ^ 4 * a ^ (-0.13e2 / 0.2e1) * J2 ^ 2 / eta ^ 8 - 0.21e2 / 0.8e1 * (0.5e1 * ci ^ 2 - 0.2e1 * ci - 0.1e1) * a ^ (-0.9e1 / 0.2e1) / eta ^ 4 * J2 * sqrt(mu) * R__e ^ 2;

dghdotdp1byp12 = 0.9871875e7 / 0.512e3 * ((ci + 1) ^ 2) * a ^ (-0.15e2 / 0.2e1) * ((eta ^ 8 + 0.248e3 / 0.13e2 * eta ^ 7 + 0.31675e5 / 0.351e3 * eta ^ 6 + 0.38440e5 / 0.351e3 * eta ^ 5 - 0.36719e5 / 0.1053e4 * eta ^ 4 - 0.215048e6 / 0.1053e4 * eta ^ 3 - 0.27353e5 / 0.117e3 * eta ^ 2 - 0.142360e6 / 0.1053e4 * eta - 0.34510e5 / 0.1053e4) * (ci ^ 13) - 0.5e1 / 0.13e2 * (0.1e1 + eta) * (eta ^ 7 + 0.17e2 * eta ^ 6 + 0.26626e5 / 0.405e3 * eta ^ 5 + 0.5356e4 / 0.135e3 * eta ^ 4 - 0.23087e5 / 0.405e3 * eta ^ 3 - 0.41779e5 / 0.405e3 * eta ^ 2 - 0.2120e4 / 0.27e2 * eta - 0.2030e4 / 0.81e2) * (ci ^ 12) + (-0.2999e4 / 0.1755e4 * eta ^ 8 - 0.19444e5 / 0.585e3 * eta ^ 7 - 0.835522e6 / 0.5265e4 * eta ^ 6 - 0.994028e6 / 0.5265e4 * eta ^ 5 + 0.265424e6 / 0.5265e4 * eta ^ 4 + 0.58652e5 / 0.195e3 * eta ^ 3 + 0.1739146e7 / 0.5265e4 * eta ^ 2 + 0.111308e6 / 0.585e3 * eta + 0.16579e5 / 0.351e3) * (ci ^ 11) + 0.5e1 / 0.9e1 * (eta ^ 7 + 0.16777e5 / 0.975e3 * eta ^ 6 + 0.195823e6 / 0.2925e4 * eta ^ 5 + 0.12737e5 / 0.325e3 * eta ^ 4 - 0.152989e6 / 0.2925e4 * eta ^ 3 - 0.10291e5 / 0.117e3 * eta ^ 2 - 0.63797e5 / 0.975e3 * eta - 0.62729e5 / 0.2925e4) * (0.1e1 + eta) * (ci ^ 10) + (-0.1166833e7 / 0.6075e4 * eta ^ 2 - 0.965668e6 / 0.8775e4 * eta + 0.1936e4 / 0.1755e4 * eta ^ 8 + 0.195644e6 / 0.8775e4 * eta ^ 7 + 0.341197e6 / 0.3159e4 * eta ^ 6 + 0.399460e6 / 0.3159e4 * eta ^ 5 - 0.91487e5 / 0.2925e4 * eta ^ 4 - 0.14375116e8 / 0.78975e5 * eta ^ 3 - 0.27187e5 / 0.975e3) * (ci ^ 9) - 0.8048e4 / 0.26325e5 * (0.1e1 + eta) * (eta ^ 7 + 0.70509e5 / 0.4024e4 * eta ^ 6 + 0.1657913e7 / 0.24144e5 * eta ^ 5 + 0.316275e6 / 0.8048e4 * eta ^ 4 - 0.201877e6 / 0.4024e4 * eta ^ 3 - 0.239449e6 / 0.3018e4 * eta ^ 2 - 0.465033e6 / 0.8048e4 * eta - 0.155909e6 / 0.8048e4) * (ci ^ 8) + (0.1373672e7 / 0.131625e6 * eta ^ 4 + 0.22895224e8 / 0.394875e6 * eta ^ 3 + 0.23271556e8 / 0.394875e6 * eta ^ 2 + 0.4417496e7 / 0.131625e6 * eta - 0.14e2 / 0.39e2 * eta ^ 8 - 0.333176e6 / 0.43875e5 * eta ^ 7 - 0.2951108e7 / 0.78975e5 * eta ^ 6 - 0.3408872e7 / 0.78975e5 * eta ^ 5 + 0.1139546e7 / 0.131625e6) * (ci ^ 7) + 0.11042e5 / 0.131625e6 * (0.1e1 + eta) * (eta ^ 7 + 0.99231e5 / 0.5521e4 * eta ^ 6 + 0.1177261e7 / 0.16563e5 * eta ^ 5 + 0.220065e6 / 0.5521e4 * eta ^ 4 - 0.273433e6 / 0.5521e4 * eta ^ 3 - 0.1225429e7 / 0.16563e5 * eta ^ 2 - 0.291111e6 / 0.5521e4 * eta - 0.99673e5 / 0.5521e4) * (ci ^ 6) + (0.5285936e7 / 0.658125e6 * eta ^ 5 - 0.1277581e7 / 0.658125e6 * eta ^ 4 - 0.6754544e7 / 0.658125e6 * eta ^ 3 - 0.2212757e7 / 0.219375e6 * eta ^ 2 - 0.417488e6 / 0.73125e5 * eta + 0.8273e4 / 0.131625e6 * eta ^ 8 + 0.71104e5 / 0.50625e5 * eta ^ 7 + 0.1541653e7 / 0.219375e6 * eta ^ 6 - 0.328232e6 / 0.219375e6) * (ci ^ 5) - 0.547e3 / 0.43875e5 * (eta ^ 7 + 0.49963e5 / 0.2735e4 * eta ^ 6 + 0.597536e6 / 0.8205e4 * eta ^ 5 + 0.108838e6 / 0.2735e4 * eta ^ 4 - 0.403003e6 / 0.8205e4 * eta ^ 3 - 0.565163e6 / 0.8205e4 * eta ^ 2 - 0.130742e6 / 0.2735e4 * eta - 0.27448e5 / 0.1641e4) * (0.1e1 + eta) * (ci ^ 4) + (0.30103e5 / 0.219375e6 - 0.441604e6 / 0.3290625e7 * eta ^ 7 - 0.521278e6 / 0.759375e6 * eta ^ 6 - 0.7705396e7 / 0.9871875e7 * eta ^ 5 + 0.613832e6 / 0.3290625e7 * eta ^ 4 + 0.9453484e7 / 0.9871875e7 * eta ^ 3 + 0.9047326e7 / 0.9871875e7 * eta ^ 2 + 0.188828e6 / 0.365625e6 * eta - 0.721e3 / 0.131625e6 * eta ^ 8) * (ci ^ 3) + 0.223e3 / 0.219375e6 * (0.1e1 + eta) * (eta ^ 7 + 0.19437e5 / 0.1115e4 * eta ^ 6 + 0.139237e6 / 0.2007e4 * eta ^ 5 + 0.121429e6 / 0.3345e4 * eta ^ 4 - 0.153809e6 / 0.3345e4 * eta ^ 3 - 0.595013e6 / 0.10035e5 * eta ^ 2 - 0.44113e5 / 0.1115e4 * eta - 0.9505e4 / 0.669e3) * (ci ^ 2) + (-0.1337e4 / 0.253125e6 + 0.6e1 / 0.40625e5 * eta ^ 8 + 0.5548e4 / 0.1096875e7 * eta ^ 7 + 0.268847e6 / 0.9871875e7 * eta ^ 6 + 0.309236e6 / 0.9871875e7 * eta ^ 5 - 0.7427e4 / 0.1096875e7 * eta ^ 4 - 0.362204e6 / 0.9871875e7 * eta ^ 3 - 0.343751e6 / 0.9871875e7 * eta ^ 2 - 0.2588e4 / 0.131625e6 * eta) * ci - 0.2e1 / 0.121875e6 * (0.1e1 + eta) * (eta ^ 7 + 0.272e3 / 0.9e1 * eta ^ 6 + 0.21325e5 / 0.162e3 * eta ^ 5 + 0.3851e4 / 0.54e2 * eta ^ 4 - 0.730e3 / 0.9e1 * eta ^ 3 - 0.8243e4 / 0.81e2 * eta ^ 2 - 0.401e3 / 0.6e1 * eta - 0.1337e4 / 0.54e2)) * J2 ^ 3 * R__e ^ 6 * sqrt(mu) / eta ^ 12 / ((5 * ci ^ 2 - 1) ^ 4) / (0.1e1 + eta) ^ 4 - 0.135e3 / 0.32e2 * ((ci + 1) ^ 2) * a ^ (-0.11e2 / 0.2e1) * (((ci ^ 3) - (ci ^ 2) / 0.3e1 - 0.7e1 / 0.5e1 * ci + 0.1e1 / 0.5e1) * eta ^ 2 + ((8 * ci ^ 3) - 0.12e2 / 0.5e1 * (ci ^ 2) - 0.32e2 / 0.15e2 * ci + 0.4e1 / 0.15e2) * eta + 0.77e2 / 0.9e1 * (ci ^ 3) - 0.7e1 / 0.3e1 * (ci ^ 2) + ci - 0.1e1 / 0.9e1) * J2 ^ 2 * R__e ^ 4 * sqrt(mu) / eta ^ 8 - 0.15e2 / 0.2e1 * ((ci + 1) ^ 2) * a ^ (-0.7e1 / 0.2e1) * (ci - 0.1e1 / 0.5e1) * J2 * R__e ^ 2 * sqrt(mu) / eta ^ 4;
dghdotdp1_J2 = dghdotdp1byp12*p1;
dghdotdp2_J2 = dghdotdp1byp12*p2;

dghdotdq1byq12 = -0.658125e6 / 0.128e3 * a ^ (-0.15e2 / 0.2e1) * J2 ^ 3 * R__e ^ 6 * ((ci ^ 11 + 0.7e1 / 0.13e2 * ci ^ 10 - 0.1667e4 / 0.1170e4 * ci ^ 9 - 0.763e3 / 0.1170e4 * ci ^ 8 + 0.406e3 / 0.585e3 * ci ^ 7 + 0.2198e4 / 0.8775e4 * ci ^ 6 - 0.55e2 / 0.351e3 * ci ^ 5 - 0.343e3 / 0.8775e4 * ci ^ 4 + 0.731e3 / 0.43875e5 * ci ^ 3 + 0.7e1 / 0.2925e4 * ci ^ 2 - 0.79e2 / 0.146250e6 * ci - 0.7e1 / 0.146250e6) * (ci - 0.1e1) * eta ^ 9 + (0.41e2 / 0.9750e4 - 0.243e3 / 0.1625e4 * ci ^ 2 + 0.57e2 / 0.3250e4 * ci + 0.571e3 / 0.26e2 * ci ^ 12 - 0.249e3 / 0.26e2 * ci ^ 11 - 0.25853e5 / 0.585e3 * ci ^ 10 + 0.19061e5 / 0.1170e4 * ci ^ 9 + 0.21261e5 / 0.650e3 * ci ^ 8 - 0.86123e5 / 0.8775e4 * ci ^ 7 - 0.164168e6 / 0.14625e5 * ci ^ 6 + 0.8029e4 / 0.2925e4 * ci ^ 5 + 0.833867e6 / 0.438750e6 * ci ^ 4 - 0.158839e6 / 0.438750e6 * ci ^ 3) * eta ^ 8 + (-0.93464e5 / 0.43875e5 * ci ^ 3 - 0.551663e6 / 0.526500e6 * ci ^ 2 + 0.115e3 / 0.1053e4 * ci + 0.175115e6 / 0.1404e4 * ci ^ 12 - 0.18578e5 / 0.351e3 * ci ^ 11 - 0.356441e6 / 0.1404e4 * ci ^ 10 + 0.159199e6 / 0.1755e4 * ci ^ 9 + 0.1015927e7 / 0.5265e4 * ci ^ 8 - 0.1463042e7 / 0.26325e5 * ci ^ 7 - 0.241547e6 / 0.3510e4 * ci ^ 6 + 0.53356e5 / 0.3375e4 * ci ^ 5 + 0.86371e5 / 0.7020e4 * ci ^ 4 + 0.191e3 / 0.5850e4) * eta ^ 7 + (0.3489497e7 / 0.175500e6 * ci ^ 4 - 0.6637e4 / 0.1755e4 * ci ^ 3 - 0.4197539e7 / 0.2632500e7 * ci ^ 2 + 0.129746e6 / 0.658125e6 * ci + 0.306445e6 / 0.1404e4 * ci ^ 12 - 0.33535e5 / 0.351e3 * ci ^ 11 - 0.613505e6 / 0.1404e4 * ci ^ 10 + 0.1460e4 / 0.9e1 * ci ^ 9 + 0.53081e5 / 0.162e3 * ci ^ 8 - 0.520616e6 / 0.5265e4 * ci ^ 7 - 0.154859e6 / 0.1350e4 * ci ^ 6 + 0.246046e6 / 0.8775e4 * ci ^ 5 + 0.626e3 / 0.14625e5) * eta ^ 6 + (-0.49889e5 / 0.1950e4 * ci ^ 6 + 0.469777e6 / 0.43875e5 * ci ^ 5 + 0.131737e6 / 0.43875e5 * ci ^ 4 - 0.640843e6 / 0.438750e6 * ci ^ 3 + 0.11147e5 / 0.2632500e7 * ci ^ 2 + 0.105551e6 / 0.1316250e7 * ci + 0.30551e5 / 0.702e3 * ci ^ 12 - 0.1651e4 / 0.54e2 * ci ^ 11 - 0.759971e6 / 0.7020e4 * ci ^ 10 + 0.206789e6 / 0.3510e4 * ci ^ 9 + 0.8732291e7 / 0.105300e6 * ci ^ 8 - 0.982829e6 / 0.26325e5 * ci ^ 7 - 0.5129e4 / 0.292500e6) * eta ^ 5 + (0.2448599e7 / 0.26325e5 * ci ^ 7 + 0.2288659e7 / 0.17550e5 * ci ^ 6 - 0.213503e6 / 0.8775e4 * ci ^ 5 - 0.1964491e7 / 0.87750e5 * ci ^ 4 + 0.271741e6 / 0.87750e5 * ci ^ 3 + 0.4905583e7 / 0.2632500e7 * ci ^ 2 - 0.198917e6 / 0.1316250e7 * ci - 0.42908e5 / 0.117e3 * ci ^ 12 + 0.89941e5 / 0.702e3 * ci ^ 11 + 0.4132873e7 / 0.7020e4 * ci ^ 10 - 0.67807e5 / 0.390e3 * ci ^ 9 - 0.8265151e7 / 0.21060e5 * ci ^ 8 - 0.50753e5 / 0.877500e6) * eta ^ 4 + (-0.1329647e7 / 0.2106e4 * ci ^ 8 + 0.834164e6 / 0.5265e4 * ci ^ 7 + 0.3407351e7 / 0.17550e5 * ci ^ 6 - 0.353242e6 / 0.8775e4 * ci ^ 5 - 0.187831e6 / 0.6500e4 * ci ^ 4 + 0.220121e6 / 0.43875e5 * ci ^ 3 + 0.925423e6 / 0.526500e6 * ci ^ 2 - 0.162254e6 / 0.658125e6 * ci - 0.99665e5 / 0.156e3 * ci ^ 12 + 0.80695e5 / 0.351e3 * ci ^ 11 + 0.1391833e7 / 0.1404e4 * ci ^ 10 - 0.35684e5 / 0.117e3 * ci ^ 9 - 0.3283e4 / 0.219375e6) * eta ^ 3 + (-0.441281e6 / 0.1755e4 * ci ^ 9 - 0.2696893e7 / 0.5265e4 * ci ^ 8 + 0.3347738e7 / 0.26325e5 * ci ^ 7 + 0.2635759e7 / 0.17550e5 * ci ^ 6 - 0.1382572e7 / 0.43875e5 * ci ^ 5 - 0.78671e5 / 0.3900e4 * ci ^ 4 + 0.168926e6 / 0.43875e5 * ci ^ 3 + 0.454631e6 / 0.526500e6 * ci ^ 2 - 0.4957e4 / 0.26325e5 * ci - 0.766195e6 / 0.1404e4 * ci ^ 12 + 0.68252e5 / 0.351e3 * ci ^ 11 + 0.1161185e7 / 0.1404e4 * ci ^ 10 + 0.2003e4 / 0.87750e5) * eta ^ 2 + (0.2519e4 / 0.175500e6 + 0.896851e6 / 0.2340e4 * ci ^ 10 - 0.404419e6 / 0.3510e4 * ci ^ 9 - 0.2784383e7 / 0.11700e5 * ci ^ 8 + 0.170447e6 / 0.2925e4 * ci ^ 7 + 0.676903e6 / 0.9750e4 * ci ^ 6 - 0.23433e5 / 0.1625e4 * ci ^ 5 - 0.1334923e7 / 0.146250e6 * ci ^ 4 + 0.772309e6 / 0.438750e6 * ci ^ 3 + 0.33997e5 / 0.97500e5 * ci ^ 2 - 0.4211e4 / 0.48750e5 * ci - 0.88255e5 / 0.351e3 * ci ^ 12 + 0.20795e5 / 0.234e3 * ci ^ 11) * eta + 0.8599e4 / 0.97500e5 * ci ^ 2 + 0.3238e4 / 0.8775e4 * ci ^ 3 - 0.1337e4 / 0.73125e5 * ci + 0.1793e4 / 0.877500e6 - 0.17255e5 / 0.351e3 * ci ^ 12 + 0.2030e4 / 0.117e3 * ci ^ 11 + 0.35933e5 / 0.468e3 * ci ^ 10 - 0.40399e5 / 0.1755e4 * ci ^ 9 - 0.569773e6 / 0.11700e5 * ci ^ 8 + 0.34712e5 / 0.2925e4 * ci ^ 7 + 0.424987e6 / 0.29250e5 * ci ^ 6 - 0.14546e5 / 0.4875e4 * ci ^ 5 - 0.19351e5 / 0.9750e4 * ci ^ 4) * sqrt(mu) / eta ^ 14 / (0.5e1 * ci ^ 2 - 0.1e1) ^ 3 / (0.1e1 + eta) ^ 5 + 0.405e3 / 0.64e2 * ((eta ^ 2 + 0.28e2 / 0.3e1 * eta + 0.308e3 / 0.27e2) * ci ^ 4 + (-0.4e1 / 0.9e1 * eta ^ 2 - 0.56e2 / 0.15e2 * eta - 0.112e3 / 0.27e2) * ci ^ 3 + (-0.14e2 / 0.5e1 * eta ^ 2 - 0.224e3 / 0.45e2 * eta + 0.8e1 / 0.3e1) * ci ^ 2 + (0.4e1 / 0.5e1 * eta ^ 2 + 0.56e2 / 0.45e2 * eta - 0.16e2 / 0.27e2) * ci + 0.5e1 / 0.9e1 * eta ^ 2 + 0.28e2 / 0.45e2 * eta - 0.28e2 / 0.27e2) * a ^ (-0.11e2 / 0.2e1) * J2 ^ 2 * R__e ^ 4 * sqrt(mu) / eta ^ 10 + 0.3e1 * (0.5e1 * ci ^ 2 - 0.2e1 * ci - 0.1e1) * a ^ (-0.7e1 / 0.2e1) / eta ^ 6 * J2 * sqrt(mu) * R__e ^ 2;
dghdotdq1_J2 = dghdotdq1byq12*q1;
dghdotdq2_J2 = dghdotdq1byq12*q2;

dhdotda_J2 = -0.15e2 / 0.1024e4 * ci * (151875 * ci ^ 10 * eta ^ 7 + 2581875 * ci ^ 10 * eta ^ 6 + 9984750 * ci ^ 10 * eta ^ 5 - 254250 * ci ^ 8 * eta ^ 7 + 6025500 * ci ^ 10 * eta ^ 4 - 4398000 * ci ^ 8 * eta ^ 6 - 8657625 * ci ^ 10 * eta ^ 3 - 17155725 * ci ^ 8 * eta ^ 5 + 145950 * ci ^ 6 * eta ^ 7 - 15667125 * ci ^ 10 * eta ^ 2 - 9910425 * ci ^ 8 * eta ^ 4 + 2659950 * ci ^ 6 * eta ^ 6 - 11925000 * ci ^ 10 * eta + 12774700 * ci ^ 8 * eta ^ 3 + 10567520 * ci ^ 6 * eta ^ 5 - 38700 * ci ^ 4 * eta ^ 7 - 3806250 * ci ^ 10 + 20670150 * ci ^ 8 * eta ^ 2 + 5881860 * ci ^ 6 * eta ^ 4 - 747000 * ci ^ 4 * eta ^ 6 + 15178875 * ci ^ 8 * eta - 7287090 * ci ^ 6 * eta ^ 3 - 3021918 * ci ^ 4 * eta ^ 5 + 4695 * ci ^ 2 * eta ^ 7 + 5049875 * ci ^ 8 - 10683530 * ci ^ 6 * eta ^ 2 - 1632654 * ci ^ 4 * eta ^ 4 + 99111 * ci ^ 2 * eta ^ 6 - 7555500 * ci ^ 6 * eta + 1998936 * ci ^ 4 * eta ^ 3 + 409770 * ci ^ 2 * eta ^ 5 - 162 * eta ^ 7 - 2603400 * ci ^ 6 + 2702652 * ci ^ 4 * eta ^ 2 + 217704 * ci ^ 2 * eta ^ 4 - 4896 * eta ^ 6 + 1844370 * ci ^ 4 * eta - 263309 * ci ^ 2 * eta ^ 3 - 21325 * eta ^ 5 + 654570 * ci ^ 4 - 335721 * ci ^ 2 * eta ^ 2 - 11553 * eta ^ 4 - 222564 * ci ^ 2 * eta + 13140 * eta ^ 3 - 80950 * ci ^ 2 + 16486 * eta ^ 2 + 10827 * eta + 4011) * a ^ (-0.17e2 / 0.2e1) * J2 ^ 3 * sqrt(mu) * R__e ^ 6 / (eta ^ 12) / ((5 * ci ^ 2 - 1) ^ 3) / ((1 + eta) ^ 3) + 0.165e3 / 0.64e2 * sqrt(mu) * R__e ^ 4 * (((eta ^ 2) + 0.36e2 / 0.5e1 * eta + 0.7e1) * (ci ^ 2) - 0.9e1 / 0.5e1 * (eta ^ 2) - 0.12e2 / 0.5e1 * eta + 0.1e1) * a ^ (-0.13e2 / 0.2e1) * ci * J2 ^ 2 / (eta ^ 8) + 0.21e2 / 0.4e1 * a ^ (-0.9e1 / 0.2e1) * ci / (eta ^ 4) * J2 * sqrt(mu) * R__e ^ 2;

dhdotdp1byp12 = 0.1875e4 / 0.2e1 * a ^ (-0.15e2 / 0.2e1) * ((ci + 1) ^ 2) * (R__e ^ 2) * (a ^ 4 * ((ci ^ 2) - 0.1e1 / 0.5e1) ^ 4 * eta ^ 11 + 0.3e1 * a ^ 4 * ((ci ^ 2) - 0.1e1 / 0.5e1) ^ 4 * eta ^ 10 + 0.3e1 * (0.5e1 / 0.16e2 * J2 * (ci ^ 2) * (R__e ^ 2) - 0.3e1 / 0.16e2 * (R__e ^ 2) * J2 + a ^ 2) * a ^ 2 * ((ci ^ 2) - 0.1e1 / 0.5e1) ^ 4 * eta ^ 9 + (0.153e3 / 0.16e2 * J2 * (ci ^ 2) * (R__e ^ 2) - 0.39e2 / 0.16e2 * (R__e ^ 2) * J2 + a ^ 2) * a ^ 2 * ((ci ^ 2) - 0.1e1 / 0.5e1) ^ 4 * eta ^ 8 + 0.237e3 / 0.8e1 * (-0.675e3 / 0.2528e4 * J2 * (ci ^ 12) * (R__e ^ 2) + (0.975e3 / 0.2528e4 * (R__e ^ 2) * J2 + a ^ 2) * (ci ^ 10) + (-0.503e3 / 0.2370e4 * (R__e ^ 2) * J2 - 0.1093e4 / 0.1185e4 * a ^ 2) * (ci ^ 8) + (0.5521e4 / 0.94800e5 * (R__e ^ 2) * J2 + 0.2002e4 / 0.5925e4 * a ^ 2) * (ci ^ 6) + (-0.547e3 / 0.63200e5 * (R__e ^ 2) * J2 - 0.606e3 / 0.9875e4 * a ^ 2) * (ci ^ 4) + (0.223e3 / 0.316000e6 * (R__e ^ 2) * J2 + 0.817e3 / 0.148125e6 * a ^ 2) * (ci ^ 2) - 0.9e1 / 0.790000e6 * (R__e ^ 2) * J2 - 0.29e2 / 0.148125e6 * a ^ 2) * (R__e ^ 2) * J2 * eta ^ 7 + 0.327e3 / 0.8e1 * (R__e ^ 2) * J2 * (-0.11475e5 / 0.3488e4 * J2 * (ci ^ 12) * (R__e ^ 2) + (0.16777e5 / 0.3488e4 * (R__e ^ 2) * J2 + a ^ 2) * (ci ^ 10) + (-0.23503e5 / 0.8720e4 * (R__e ^ 2) * J2 - 0.461e3 / 0.545e3 * a ^ 2) * (ci ^ 8) + (0.33077e5 / 0.43600e5 * (R__e ^ 2) * J2 + 0.754e3 / 0.2725e4 * a ^ 2) * (ci ^ 6) + (-0.49963e5 / 0.436000e6 * (R__e ^ 2) * J2 - 0.586e3 / 0.13625e5 * a ^ 2) * (ci ^ 4) + (0.19437e5 / 0.2180000e7 * (R__e ^ 2) * J2 + 0.209e3 / 0.68125e5 * a ^ 2) * (ci ^ 2) - 0.17e2 / 0.68125e5 * (R__e ^ 2) * J2 - a ^ 2 / 0.13625e5) * eta ^ 6 + 0.423e3 / 0.16e2 * (R__e ^ 2) * (-0.66565e5 / 0.3384e4 * J2 * (ci ^ 12) * (R__e ^ 2) + (0.195823e6 / 0.6768e4 * (R__e ^ 2) * J2 + a ^ 2) * (ci ^ 10) + (-0.1657913e7 / 0.101520e6 * (R__e ^ 2) * J2 - 0.559e3 / 0.705e3 * a ^ 2) * (ci ^ 8) + (0.1177261e7 / 0.253800e6 * (R__e ^ 2) * J2 + 0.826e3 / 0.3525e4 * a ^ 2) * (ci ^ 6) + (-0.37346e5 / 0.52875e5 * (R__e ^ 2) * J2 - 0.178e3 / 0.5875e4 * a ^ 2) * (ci ^ 4) + (0.139237e6 / 0.2538000e7 * (R__e ^ 2) * J2 + 0.121e3 / 0.88125e5 * a ^ 2) * (ci ^ 2) - 0.853e3 / 0.507600e6 * (R__e ^ 2) * J2 + a ^ 2 / 0.88125e5) * J2 * eta ^ 5 + 0.105e3 / 0.16e2 * (R__e ^ 2) * J2 * (-0.1339e4 / 0.28e2 * J2 * (ci ^ 12) * (R__e ^ 2) + (0.38211e5 / 0.560e3 * (R__e ^ 2) * J2 + a ^ 2) * (ci ^ 10) + (-0.4217e4 / 0.112e3 * (R__e ^ 2) * J2 - 0.79e2 / 0.105e3 * a ^ 2) * (ci ^ 8) + (0.14671e5 / 0.1400e4 * (R__e ^ 2) * J2 + 0.106e3 / 0.525e3 * a ^ 2) * (ci ^ 6) + (-0.54419e5 / 0.35000e5 * (R__e ^ 2) * J2 - 0.18e2 / 0.875e3 * a ^ 2) * (ci ^ 4) + (0.17347e5 / 0.150000e6 * (R__e ^ 2) * J2 + a ^ 2 / 0.13125e5) * (ci ^ 2) - 0.3851e4 / 0.1050000e7 * (R__e ^ 2) * J2 + a ^ 2 / 0.13125e5) * eta ^ 4 + 0.115435e6 / 0.256e3 * (R__e ^ 4) * ((ci ^ 12) - 0.152989e6 / 0.115435e6 * (ci ^ 10) + 0.403754e6 / 0.577175e6 * (ci ^ 8) - 0.546866e6 / 0.2885875e7 * (ci ^ 6) + 0.403003e6 / 0.14429375e8 * (ci ^ 4) - 0.153809e6 / 0.72146875e8 * (ci ^ 2) + 0.876e3 / 0.14429375e8) * (J2 ^ 2) * eta ^ 3 + 0.208895e6 / 0.256e3 * (R__e ^ 4) * (J2 ^ 2) * ((ci ^ 12) - 0.1255e4 / 0.1019e4 * (ci ^ 10) + 0.1915592e7 / 0.3133425e7 * (ci ^ 8) - 0.2450858e7 / 0.15667125e8 * (ci ^ 6) + 0.565163e6 / 0.26111875e8 * (ci ^ 4) - 0.595013e6 / 0.391678125e9 * (ci ^ 2) + 0.16486e5 / 0.391678125e9) * eta ^ 2 + 0.19875e5 / 0.32e2 * (R__e ^ 4) * (J2 ^ 2) * ((ci ^ 12) - 0.63797e5 / 0.53000e5 * (ci ^ 10) + 0.155011e6 / 0.265000e6 * (ci ^ 8) - 0.97037e5 / 0.662500e6 * (ci ^ 6) + 0.65371e5 / 0.3312500e7 * (ci ^ 4) - 0.44113e5 / 0.33125000e8 * (ci ^ 2) + 0.1203e4 / 0.33125000e8) * eta + 0.25375e5 / 0.128e3 * ((ci ^ 12) - 0.62729e5 / 0.50750e5 * (ci ^ 10) + 0.155909e6 / 0.253750e6 * (ci ^ 8) - 0.491e3 / 0.3125e4 * (ci ^ 6) + 0.13724e5 / 0.634375e6 * (ci ^ 4) - 0.1901e4 / 0.1268750e7 * (ci ^ 2) + 0.191e3 / 0.4531250e7) * (R__e ^ 4) * (J2 ^ 2)) * J2 * sqrt(mu) / eta ^ 12 / (0.1e1 + eta) ^ 3 / ((5 * ci ^ 2 - 1) ^ 4);
dhdotdp1_J2 = dhdotdp1byp12*p1;
dhdotdp2_J2 = dhdotdp1byp12*p2;

dhdotdq1byq12 = (607500 * ci ^ 10 * eta ^ 8 + 11998125 * ci ^ 10 * eta ^ 7 + 57669375 * ci ^ 10 * eta ^ 6 - 1017000 * ci ^ 8 * eta ^ 8 + 68086875 * ci ^ 10 * eta ^ 5 - 20426625 * ci ^ 8 * eta ^ 7 - 27843750 * ci ^ 10 * eta ^ 4 - 98972625 * ci ^ 8 * eta ^ 6 + 583800 * ci ^ 6 * eta ^ 8 - 140795625 * ci ^ 10 * eta ^ 3 - 114552375 * ci ^ 8 * eta ^ 5 + 12334650 * ci ^ 6 * eta ^ 7 - 161810625 * ci ^ 10 * eta ^ 2 + 37006500 * ci ^ 8 * eta ^ 4 + 60817450 * ci ^ 6 * eta ^ 6 - 154800 * eta ^ 8 * ci ^ 4 - 94134375 * ci ^ 10 * eta + 191842125 * ci ^ 8 * eta ^ 3 + 69336550 * ci ^ 6 * eta ^ 5 - 3458250 * ci ^ 4 * eta ^ 7 - 22837500 * ci ^ 10 + 209602875 * ci ^ 8 * eta ^ 2 - 20195100 * ci ^ 6 * eta ^ 4 - 17350590 * ci ^ 4 * eta ^ 6 + 18780 * ci ^ 2 * eta ^ 8 + 121357875 * ci ^ 8 * eta - 102234850 * ci ^ 6 * eta ^ 3 - 19556310 * ci ^ 4 * eta ^ 5 + 457737 * ci ^ 2 * eta ^ 7 + 30299250 * ci ^ 8 - 106306150 * ci ^ 6 * eta ^ 2 + 5463000 * ci ^ 4 * eta ^ 4 + 2346183 * ci ^ 2 * eta ^ 6 - 648 * eta ^ 8 - 61080750 * ci ^ 6 * eta + 26562450 * ci ^ 4 * eta ^ 3 + 2631567 * ci ^ 2 * eta ^ 5 - 22437 * eta ^ 7 - 15620400 * ci ^ 6 + 26423850 * ci ^ 4 * eta ^ 2 - 709038 * ci ^ 2 * eta ^ 4 - 121313 * eta ^ 6 + 15053310 * ci ^ 4 * eta - 3367077 * ci ^ 2 * eta ^ 3 - 138179 * eta ^ 5 + 3927420 * ci ^ 4 - 3236553 * ci ^ 2 * eta ^ 2 + 32628 * eta ^ 4 - 1831227 * ci ^ 2 * eta + 166289 * eta ^ 3 - 485700 * ci ^ 2 + 158219 * eta ^ 2 + 89631 * eta + 24066) * ci * a ^ (-0.15e2 / 0.2e1) * J2 ^ 3 * sqrt(mu) * R__e ^ 6 / (eta ^ 14) / ((5 * ci ^ 2 - 1) ^ 3) / ((1 + eta) ^ 4) / 0.256e3 - 0.45e2 / 0.16e2 * a ^ (-0.11e2 / 0.2e1) * ci * (((eta ^ 2) + 0.42e2 / 0.5e1 * eta + 0.28e2 / 0.3e1) * (ci ^ 2) - 0.9e1 / 0.5e1 * (eta ^ 2) - 0.14e2 / 0.5e1 * eta + 0.4e1 / 0.3e1) * J2 ^ 2 * R__e ^ 4 * sqrt(mu) / (eta ^ 10) - 0.6e1 * a ^ (-0.7e1 / 0.2e1) * ci / (eta ^ 6) * J2 * sqrt(mu) * R__e ^ 2;
dhdotdq1_J2 = dhdotdq1byq12*q1;
dhdotdq2_J2 = dhdotdq1byq12*q2;


% Higher order Zonals
lghdotn = 0;
hdotn = 0;
ghdotn = 0;

dlghdotdan = 0;
dlghdotdp1n = 0;
dlghdotdp2n = 0;
dlghdotdq1n = 0;
dlghdotdq2n = 0;

dghdotdan = 0;
dghdotdp1n = 0;
dghdotdp2n = 0;
dghdotdq1n = 0;
dghdotdq2n = 0;

dhdotdan = 0;
dhdotdp1n = 0;
dhdotdp2n = 0;
dhdotdq1n = 0;
dhdotdq2n = 0;

% No second order contributions from odd harmonics
for ctr = 2:2:length(narr)
    
    % current zonal
    n = narr(ctr);
    Jn = Jcoeff(n);
    
    % Get second order rates
    [ lghdot, ghdot, hdot ] = ZonalRates(a,e,inc, n, Jn, J2, mu, Re );
    lghdotn = lghdotn + lghdot;
    ghdotn = ghdotn + ghdot;
    hdotn = hdotn + hdot;
    
    % Get partials of second order rates
    [dlghdotda,dlghdotdp1,dlghdotdp2,dlghdotdq1,dlghdotdq2] = getMLdotRates(lghdotn,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re );
    [dghdotda,dghdotdp1,dghdotdp2,dghdotdq1,dghdotdq2] = getGHdotRates(ghdotn,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re );
    [dhdotda,dhdotdp1,dhdotdp2,dhdotdq1,dhdotdq2] = getHdotRates(hdotn,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re );
  
    dlghdotdan = dlghdotdan + dlghdotda;
    dlghdotdp1n = dlghdotdp1n + dlghdotdp1;
    dlghdotdp2n = dlghdotdp2n + dlghdotdp2;
    dlghdotdq1n = dlghdotdq1n + dlghdotdq1;
    dlghdotdq2n = dlghdotdq2n + dlghdotdq2;

    dghdotdan = dghdotdan + dghdotda;
    dghdotdp1n = dghdotdp1n + dghdotdp1;
    dghdotdp2n = dghdotdp2n + dghdotdp2;
    dghdotdq1n = dghdotdq1n + dghdotdq1;
    dghdotdq2n = dghdotdq2n + dghdotdq2;

    dhdotdan = dhdotdan + dhdotda;
    dhdotdp1n = dhdotdp1n + dhdotdp1;
    dhdotdp2n = dhdotdp2n + dhdotdp2;
    dhdotdq1n = dhdotdq1n + dhdotdq1;
    dhdotdq2n = dhdotdq2n + dhdotdq2;

end

% Total rates due to all harmonics
lghdot = lghdot2 + J2^2/2*lghdotn;
ghdot = ghdot2 + J2^2/2*ghdotn;
hdot = hdot2 + J2^2/2*hdotn;

dlghdotda = dlghdotda_J2 + J2^2/2*dlghdotdan;
dlghdotdp1 = dlghdotdp1_J2 + J2^2/2*dlghdotdp1n;
dlghdotdp2 = dlghdotdp2_J2 + J2^2/2*dlghdotdp2n;
dlghdotdq1 = dlghdotdq1_J2 + J2^2/2*dlghdotdq1n;
dlghdotdq2 = dlghdotdq2_J2 + J2^2/2*dlghdotdq2n;

dghdotda = dghdotda_J2 + J2^2/2*dghdotdan;
dghdotdp1 = dghdotdp1_J2 + J2^2/2*dghdotdp1n;
dghdotdp2 = dghdotdp2_J2 + J2^2/2*dghdotdp2n;
dghdotdq1 = dghdotdq1_J2 + J2^2/2*dghdotdq1n;
dghdotdq2 = dghdotdq2_J2 + J2^2/2*dghdotdq2n;

dhdotda = dhdotda_J2 + J2^2/2*dhdotdan;
dhdotdp1 = dhdotdp1_J2 + J2^2/2*dhdotdp1n;
dhdotdp2 = dhdotdp2_J2 + J2^2/2*dhdotdp2n;
dhdotdq1 = dhdotdq1_J2 + J2^2/2*dhdotdq1n;
dhdotdq2 = dhdotdq2_J2 + J2^2/2*dhdotdq2n;

% propagate Equinoctial elements
at = X0(1);
Lambdat = X0(2) + (lghdot)*dt;

p1t = X0(3)*cos(hdot*dt) - X0(4)*sin(hdot*dt);
p2t = X0(4)*cos(hdot*dt) + X0(3)*sin(hdot*dt);

q1t = X0(5)*cos(ghdot*dt) - X0(6)*sin(ghdot*dt);
q2t = X0(6)*cos(ghdot*dt) + X0(5)*sin(ghdot*dt);

Xt = [at Lambdat p1t p2t q1t q2t]';

% Mean STM

% Mean SMA is constant
dadx = [1,0,0,0,0,0];

% Mean Longitude
dLda0 = dlghdotda*dt;

dLdL0 = 1;

dLdp10 = dlghdotdp1*dt;

dLdp20 = dlghdotdp2*dt;

dLdq10 = dlghdotdq1*dt;

dLdq20 = dlghdotdq2*dt;

% p1
dp1da0 = -p2t*dhdotda*dt; 

dp1dL0 = 0;

dp1dp10 = cos(hdot*dt) - p2t*dhdotdp1*dt;
dp1dp20 = -sin(hdot*dt) - p2t*dhdotdp2*dt;

dp1dq10 = -p2t*dhdotdq1*dt;
dp1dq20 = -p2t*dhdotdq2*dt;

% p2
dp2da0 = p1t*dhdotda*dt; 

dp2dL0 = 0;

dp2dp10 = sin(hdot*dt) + p1t*dhdotdp1*dt;

dp2dp20 = cos(hdot*dt) + p1t*dhdotdp2*dt;

dp2dq10 = p1t*dhdotdq1*dt;

dp2dq20 = p1t*dhdotdq2*dt;

% q1
dq1da0 = -q2t*dghdotda*dt;

dq1dL0 = 0;

dq1dp10 = -q2t*dghdotdp1*dt;
dq1dp20 = -q2t*dghdotdp2*dt;

dq1dq10 = cos(ghdot*dt) - q2t*dghdotdq1*dt;
dq1dq20 = -sin(ghdot*dt) - q2t*dghdotdq2*dt;

% q2
dq2da0 = q1t*dghdotda*dt;

dq2dL0 = 0;

dq2dp10 = q1t*dghdotdp1*dt;
dq2dp20 = q1t*dghdotdp2*dt;

dq2dq10 = sin(ghdot*dt) + q1t*dghdotdq1*dt;
dq2dq20 = cos(ghdot*dt) + q1t*dghdotdq2*dt;


phi(1,1:6) = dadx;

phi(2,1:6) = [dLda0, dLdL0, dLdp10, dLdp20, dLdq10, dLdq20];

phi(3,1:6) = [dp1da0, dp1dL0, dp1dp10, dp1dp20, dp1dq10, dp1dq20];

phi(4,1:6) = [dp2da0, dp2dL0, dp2dp10, dp2dp20, dp2dq10, dp2dq20];

phi(5,1:6) = [dq1da0, dq1dL0, dq1dp10, dq1dp20, dq1dq10, dq1dq20];

phi(6,1:6) = [dq2da0, dq2dL0, dq2dp10, dq2dp20, dq2dq10, dq2dq20];


end

function [dlghdotda,dlghdotdp1,dlghdotdp2,dlghdotdq1,dlghdotdq2] = getMLdotRates(lghdot,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re )

dlghdotda = -(2*n + 3)/2/a*lghdot;

dlghdotdp12 = dabF(@getmldotp12,2,0,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re );

dlghdotdp1 = dlghdotdp12*p1;
dlghdotdp2 = dlghdotdp12*p2;

dLdotdq101 = dabF(@getmldotq121,0,0,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re );
dLdotdq102 = dabF(@getmldotq122,0,2,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re );

dlghdotdq1 = (dLdotdq101 + dLdotdq102)*q1;
dlghdotdq2 = (dLdotdq101 + dLdotdq102)*q2;

end

function [dghdotda,dghdotdp1,dghdotdp2,dghdotdq1,dghdotdq2] = getGHdotRates(ghdot,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re )

eta = sqrt(1 - e^2);

dghdotda = (-(n+1)/a-1/2/a)*ghdot;

dghdotdp12 = dabF(@getghdotp12,2,0,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re );
dghdotdp1 = dghdotdp12*p1;
dghdotdp2 = dghdotdp12*p2;

dq1dq10t2 = dabF(@getdq1dq1j0e0,0,0,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re );
dq1dq10t3 = dabF(@getdq1dq1j0e2,0,2,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re );
dq1dq10t4 = dabF(@getdq1dq1j0e4,0,4,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re );

dghdotdq1 = (2*n-1)/eta^2*ghdot*q1 + (dq1dq10t2 + dq1dq10t3 + dq1dq10t4)*q1;
dghdotdq2 = (2*n-1)/eta^2*ghdot*q2 + (dq1dq10t2 + dq1dq10t3 + dq1dq10t4)*q2;

end

function [dhdotda,dhdotdp1,dhdotdp2,dhdotdq1,dhdotdq2] = getHdotRates(hdot,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re )

dhdotda = hdot*(-(2*n+3)/(2*a));

dhdotdp12 = dabF(@gethdotdp12,2,0,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re );

dhdotdp1 = dhdotdp12*p1;
dhdotdp2 = dhdotdp12*p2;

dhdotdq01 = dabF(@gethdotq121,0,0,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re );
dhdotdq02 = dabF(@gethdotq122,0,2,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re );

dhdotdq1 = (dhdotdq01 + dhdotdq02)*q1;
dhdotdq2 = (dhdotdq01 + dhdotdq02)*q2;

end




function p12term = getmldotp12(n,j,k,a,e,i,eta,ci,p1,p2,q1,q2,mu)
si = sin(i);
p12term = (n-2*j)*((-3 - (2*n-1)/eta - k*eta/(1+eta) - (n-2*j)*ci/eta/(1+ci))*ci*(1+ci)^2 + si^2);

end

function q1term = getmldotq121(n,j,k,a,e,i,eta,ci,p1,p2,q1,q2,mu)

q1term = (((-3 - (2*n-1)/eta - k*eta/(1+eta) - (n-2*j)*ci/eta/(1+ci)))*((2*n-1)/eta^2 + 0) ...
                 + (-(2*n-1)/eta^3 + k/eta/(1+eta)^2 - (n-2*j)*ci/(1+ci)/eta^4));

end

function q1term = getmldotq122(n,j,k,a,e,i,eta,ci,p1,p2,q1,q2,mu)

q1term = (((-3 - (2*n-1)/eta - k*eta/(1+eta) - (n-2*j)*ci/eta/(1+ci)))*(k));

end

function hdotterm = gethdot(n,j,k,a,e,i,eta,ci,p1,p2,q1,q2,mu)

ci = cos(i);

hdotterm = -1/eta/ci*(n-2*j)*ci^2;

end

function hdotda = getdhdotda(n,j,k,a,e,i,eta,ci,p1,p2,q1,q2,mu)

ci = cos(i);
H = sqrt(mu*a)*eta*ci;

hdotda = (-1/eta/ci*(n-2*j)*ci^2)*(-(n+1)/a -1/(2*a*H));

end

function hdotp12term = gethdotdp12(n,j,k,a,e,i,eta,ci,p1,p2,q1,q2,mu)

ci = cos(i);

hdotp12term = 1/eta*(n-2*j)*(1+ci)^2*(1 - (n-2*j)*ci^2);

end



function hq1term = gethdotq121(n,j,k,a,e,i,eta,ci,p1,p2,q1,q2,mu)
ci = cos(i);

hq1term = 1/eta*(-(n-2*j)*ci)*(2*n-1)/eta^2;

end

function hq2term = gethdotq122(n,j,k,a,e,i,eta,ci,p1,p2,q1,q2,mu)
ci = cos(i);

hq2term = 1/eta*(-(n-2*j)*ci)*(k);

end

function ghdotp1term = getghdotp12(n,j,k,a,e,i,eta,ci,p1,p2,q1,q2,mu)
ci = cos(i);

ghdotp1term = 1/eta*(n-2*j)*(ci*(1+ci)^2 + 1);

end



function dq1dq1term1 = getdq1dq1j0e0(n,j,k,a,e,i,eta,ci,p1,p2,q1,q2,mu)

ci = cos(i);

dq1dq1term1 = 1/eta*((2*n-1)/eta^2 + (n-2*j)*ci/(1+ci)/eta^2);

end

function dq1dq1term2 = getdq1dq1j0e2(n,j,k,a,e,i,eta,ci,p1,p2,q1,q2,mu)

ci = cos(i);

dq1dq1term2 = 1/eta*((2*n-1)*(-k) + (n-2*j)*ci/(1+ci)*(-k) - k + k^2);

end

function dq1dq1term3 = getdq1dq1j0e4(n,j,k,a,e,i,eta,ci,p1,p2,q1,q2,mu)

ci = cos(i);

dq1dq1term3 = 1/eta*((2-k)*k);

end



function dab = dabF(func,jf,k0,a,e,inc,L,G,H,p1,p2,q1,q2, n, Jn, J2, mu, Re )

dab = 0;ci = cos(inc);

% No contributions from odd zonal harmonics
if mod(n,2) ~= 0
    return;
end

% intermediate terms
eta = sqrt(1 - e^2);

% outer j loop
sumblgh = 0;
for j = 0:(n/2 - jf/2)
    
    bn2 = nchoosek((n-2*j),(n/2-j));
    
    sumAlgh = 0;
    for k1 = (k0/2):(floor(n/2)-1)
        
        if mod(n,2) == 0
            % even harmonic
            k = 2*k1;
        else
            % odd harmonic
            k = 2*k1 + 1;
        end
        
        Ckk2 = nchoosek((k),(k/2));
        t1 = alphak0(k,n,e,k0)/2^k*Ckk2;
        
        % lghdot
        sumAlgh = sumAlgh + t1/L*func(n,j,k,a,e,inc,eta,ci,p1,p2,q1,q2,mu);
        
    end
    
    bt2 = betajfn(j,n,inc,jf);
    sumblgh = sumblgh + bt2*bn2*sumAlgh;
    
end

del = delta(n, a, eta, mu, Re, J2, Jn);
dab = del*sumblgh;

end


function [ lghdot, ghdot, hdot ] = ZonalRates(a,e,i, n, Jn, J2, mu, Re )

lghdot = 0;
ghdot = 0;
hdot = 0;

% No contributions from odd zonal harmonics
if mod(n,2) ~= 0
    return;
end

% intermediate terms
eta = sqrt(1 - e^2);
L = sqrt(mu*a);
G = L*eta;
H = G*cos(i);
incns = cos(i)/(2*(cos(i/2))^2);

% outer j loop
%     sumbl = 0;
%     sumbg = 0;
sumbh = 0;
sumbgh = 0;
sumblgh = 0;
for j = 0:floor(n/2)
    
    bn2 = nchoosek((n-2*j),(n/2-j));
    
    % inner k loop
    %         sumAl = 0;
    %         sumAg = 0;
    sumAh = 0;
    sumAgh1 = 0;
    sumAgh2 = 0;
    sumAlgh = 0;
    for k1 = 0:(floor(n/2)-1)
        
        if mod(n,2) == 0
            % even harmonic
            k = 2*k1;
        else
            % odd harmonic
            k = 2*k1 + 1;
        end
        
        Ckk2 = nchoosek((k),(k/2));
        t1 = alpha(k,n,e)/2^k*Ckk2;
        %             % ldot
        %             sumAl = sumAl + t1*(-3 + k*eta^2/e^2)/L;
        %             % gdot
        %             sumAg = sumAg + t1*(-(2*n - 1) - k*eta^2/e^2 + (n - 2*j)/(tan(i))^2)/G;
        
        % hdot
        if j < floor(n/2)
            sumAh = sumAh + t1*(-(n - 2*j)*cos(i)^2)/H;
        end
        
        % ghdot
        sumAgh1 = sumAgh1 + t1/G*(-(2*n - 1) - (n - 2*j)*incns);
        if k >= 2
            sumAgh2 = sumAgh2 + (1/2^k*Ckk2*e^(k-2)*(nchoosek((n-1),(k))))/G*(-k*eta^2);
        end
        
        % lghdot
        sumAlgh = sumAlgh + t1/L*(-3 -(2*n - 1)/eta - k*eta/(1 + eta) ...
            - (n - 2*j)/eta*incns);
        
    end
    
    bt = betajn(j,n,i);
    if j < floor(n/2)
        bt2 = betajn2(j,n,i);
        sumbh = sumbh + bt2*bn2*sumAh;
    end

    sumbgh = sumbgh + bt*bn2*(sumAgh1 + sumAgh2);
    sumblgh = sumblgh + bt*bn2*sumAlgh;
    
end

del = delta(n, a, eta, mu, Re, J2, Jn);
% ldot = del*sumbl;
% gdot = del*sumbg;
hdot = del*sumbh;
ghdot = del*sumbgh;
lghdot = del*sumblgh;


end


function y = delta(n, a, eta, mu, Re, J2, Jn)

y = 2*Jn/J2^2*mu*Re^n/(2^n*eta^(2*n-1)*a^(n+1));

end


function y = betajn(j,n,i)

n2j = n - 2*j;
y = (-1)^j*factorial(2*n-2*j)*((sin(i))^(n2j))/(factorial(j)*factorial(n-j)*factorial(n2j)*2^(n2j));

end

function y = betajn2(j,n,i)

n2j = n - 2*j;
y = (-1)^j*factorial(2*n-2*j)*((sin(i))^(n2j-2))/(factorial(j)*factorial(n-j)*factorial(n2j)*2^(n2j));

end


function y = betajfn(j,n,i,jf)

n2j = n - 2*j;
y = (-1)^j*factorial(2*n-2*j-jf)*((sin(i))^(n2j))/(factorial(j)*factorial(n-j)*factorial(n2j)*2^(n2j));

end

function y = alpha(k,n,e)

y = e^k*(nchoosek((n-1), (k)));

end

function y = alpha2(k,n,e)

y = e^(k-2)*(nchoosek((n-1), (k)));

end

function y = alphak0(k,n,e, k0)

y = e^(k-k0)*(nchoosek((n-1), (k)));

end