Home » Posts tagged 'Kepler’s equation'
Tag Archives: Kepler’s equation
Kepler Equation Solver Without Transcendental Function Evaluations
clc;clear all;
tic; j = 1; for M = 0:0.005:pi i = 1; for e = 0:0.005:1 c3 = 5/2 + 560*e; a = 15/c3*(1-e); b = -M/c3; y = (b^2/4+ a^3/27)^0.5; x = (-b/2 + y)^(1/3) - (b/2 + y)^(1/3); c15 = 3003/14336 + 16384*e; c13 = 3465/13312 - 61440*e; c11 = 945/2816 + 92160*e; c9 = 175/384 -70400*e; c7 = 75/112 + 28800*e; c5 = 9/8-6048*e; x2 = x^2;x3 = x2*x;x4 = x3*x; x5 = x4*x;x6 = x5*x;x7 = x6*x; x8 = x7*x;x9 = x8*x; x10 =x9*x; x11 = x10*x; x12 = x11*x; x13 = x12*x; x14 = x13*x; x15 = x14*x; f = c15*x15 + c13*x13 + c11*x11 + c9*x9 + c7*x7 + c5*x5 + c3*x3 +... 15*(1 - e)*x -M; f1 = 15* c15 *x14 + 13* c13 *x12 + 11* c11 *x10 + 9* c9 *x8 + 7* c7 *x6 + ... 5 *c5 *x4 + 3 *c3 *x2 + 15*(1 - e); f2 = 210* c15 *x13 + 156* c13 *x11 + 110* c11 *x9 + 72* c9 *x7 + ... 42* c7 *x5 + 20* c5 *x3 + 6* c3 *x; f3 = 2730*c15 *x12 + 1716* c13 *x10 + 990* c11 *x8 + 504* c9 *x6 +... 210* c7 *x4 + 60* c5 *x2 + 6* c3; f4 = 32760* c15 *x11 + 17160* c13 *x9 + 7920* c11 *x7 + 3024* c9 *x5 +... 840* c7 *x3 + 120* c5 *x; f5 = 360360* c15 *x10 + 154440* c13 *x8 + 55440* c11 *x6 + 15120* c9 *x4 ... + 2520* c7 *x2 + 120* c5; f6 = 3603600* c15 *x9 + 1235520* c13 *x7 + 332640* c11 *x5 + 60480* c9 *x3... + 5040* c7 *x; f7 = 32432400* c15 *x8 + 8648640* c13 *x6 + 1663200* c11 *x4 + 181440* c9 *x2 ... + 5040* c7; f8 = 25945920* c15 *x7 + 51891840* c13 *x5 + 6652800* c11 *x3 + ... 362880* c9 *x; f9 = 1.8162144e9* c15 *x6 + 259459200* c13 *x4 + 19958400* c11 *x2 + ... 362880* c9; f10 = 1.08972864e10* c15 *x5 + 1.0378368e9*c13 *x3 + 39916800* c11 *x; f11 = 5.4486432e10* c15 *x4 + 3.1135104e9* c13 *x2 + 39916800* c11; f12 = 2.17945728e11* c15 *x3 + 6.2270208e9* c13 *x; f13 = 6.53837184e11* c15 *x2 + 6.2270208e9* c13; f14 = 1.307674368e13* c15 *x; f15 = 1.307674368e13* c15; g1 = 1/2; g2 = 1/6; g3 = 1/24; g4 = 1/120; g5 = 1/720; g6 = 1/5040; g7 = 1/40320; g8 = 1/362880; g9 = 1/3628800; g10 = 1/39916800; g11 = 1/479001600; g12 =1/6.2270208e9 ; g13 = 1/8.71782912e10; g14 = 1/1.307674368e12; u1 = -f/f1; h2 = f1 + g1* u1* f2; u2 = -f/h2; h3 = f1 + g1* u2* f2 + g2* u2^2*f3; u3 = -f/h3; h4 = f1 + g1* u3* f2 + g2* u3^2*f3 + g3* u3^3*f4; u4 = -f/h4; h5 = f1 + g1* u4* f2 + g2* u4^2*f3 + g3* u4^3*f4 + g4* u4^4*f5; u5 = -f/h5; h6 = f1 + g1* u5* f2 + g2* u5^2*f3 + g3* u5^3*f4 + g4* u5^4*f5 +... g5* u5^5*f6; u6 = -f/h6; h7 = f1 + g1* u6* f2 + g2* u6^2*f3 + g3* u6^3*f4 + g4* u6^4*f5 +... g5* u6^5*f6 + g6* u6^6*f7; u7 = -f/h7; h8 = f1 + g1* u7* f2 + g2* u7^2*f3 + g3* u7^3*f4 + g4* u7^4*f5 +... g5* u7^5*f6 + g6* u7^6*f7 + g7* u7^7*f8; u8 = -f/h8; h9 = f1 + g1* u8* f2 + g2* u8^2*f3 + g3* u8^3*f4 + g4* u8^4*f5 +... g5* u8^5*f6 + g6* u8^6*f7+ g7* u8^7*f8 + g8* u8^8*f9; u9 = - f/h9; h10 = f1 + g1* u9* f2 + g2* u9^2*f3 + g3* u9^3*f4 + g4* u9^4*f5 +... g5* u9^5*f6 + g6* u9^6*f7+ g7* u9^7*f8 + g8* u9^8*f9 + ... g9* u9^9*f10; u10 = -f/h10; h11 = f1 + g1* u10* f2 + g2* u10^2* f3 + g3* u10^3* f4 + g4* u10^4*f5... + g5* u10^5* f6 + g6* u10^6*f7+ g7* u10^7* f8 + g8* u10^8* f9... + g9* u10^9* f10 + g10* u10^10* f11; u11 = -f/h11; h12 = f1 + g1* u11* f2 + g2* u11^2* f3 + g3* u11^3* f4 + ... g4* u11^4* f5 + g5* u11^5* f6 + g6* u11^6* f7... + g7* u11^7* f8 + g8* u11^8* f9 + g9* u11^9* f10 + g10* u11^10* f11... + g11* u11^11* f12; u12 = -f/h12; h13 = f1 + g1* u12* f2 + g2* u12^2* f3 + g3* u12^3* f4 +... g4* u12^4* f5 + g5* u12^5* f6 + g6* u12^6* f7... + g7* u12^7* f8 + g8* u12^8* f9 + g9* u12^9* f10 +... g10* u12^10* f11 + g11* u12^11* f12 + g12* u12^12* f13; u13 = -f/h13; h14 = f1 + g1* u13* f2 + g2* u13^2* f3 + g3* u13^3* f4 + ... g4* u13^4* f5 + g5* u13^5* f6 + g6* u13^6* f7... + g7* u13^7* f8 + g8* u13^8* f9 + g9* u13^9* f10 + g10* u13^10*f11... + g11* u13^11* f12+ g12* u13^12* f13 + g13* u13^13* f14; u14 = -f/h14; h15 = f1 + g1* u14* f2 + g2* u14^2* f3 + g3* u14^3* f4 + ... g4* u14^4* f5 + g5* u14^5* f6 + g6* u14^6* f7... + g7* u14^7* f8 + g8* u14^8* f9 + g9* u14^9* f10 +... g10* u14^10* f11 + g11*u14^11* f12+ g12* u14^12* f13 + ... g13* u14^13* f14 + g14* u14^14* f15; u15 = -f/h15; x = x + u15; w = x - 0.01171875*x^17/(1 + e); E = M + e*(-16384*w^15 + 61440*w^13 - 92160*w^11+ 70400*w^9 -... 28800*w^7 + 6048*w^5 - 560*w^3 + 15*w); E_nit(i,j) = E; i = i + 1; end j = j + 1; end t_nit = toc; fprintf('Calculation time = %4.6f [s] \n',t_nit);
Calculation time = 2.591335 [s]
Newton Iterative Solver
tic; j = 1; eps = 1E-15; % Tolerance for M = 0:0.005:pi i = 1; for e = 0:0.005:1 En = M; Ens = En - (En-e*sin(En)- M)/(1 - e*cos(En)); while ( abs(Ens-En) > eps ) En = Ens; Ens = En - (En - e*sin(En) - M)/(1 - e*cos(En)); end; E_it(i,j) = Ens; i = i + 1; end j = j + 1; end t_it = toc; fprintf('Calculation time = %4.6f [s] \n',t_it);
Calculation time = 3.520123 [s]
Error = log10(abs(E_it - E_nit)); e = 0:0.005:1; M = (0:0.005:pi); contourf (M,e,Error); xlabel('M - Mean anomaly [rad]'); ylabel('e - Eccentricity'); title('Error in E calculation on Log scale [rad] '); colorbar; REFERENCES Adonis Pimienta-Pe˜nalver and John L. Crassidis†,ACCURATE KEPLER EQUATION SOLVER WITHOUT TRANSCENDENTAL FUNCTION EVALUATIONS, University at Buffalo, State University of New York, Amherst, NY, 14260-4400
Kepler’s equation, Iterative and Non-Iterative Solver Comparison
Kepler’s equation M = E-e*sin(E) M – Mean anomaly [0..pi] rad e – Eccentricity [0..1] E – Eccentric anomaly [rad]
clear all;
clc;
Non-Iterative Solver
tic; j = 1; for M = 0:0.01:pi i = 1; for e = 0:0.01:1 a = (1 - e)*3/(4*e+0.5); b = -M/(4*e+0.5); y = (b^2/4 + a^3/27)^0.5; x = (-b/2 + y)^(1/3) - (b/2 + y)^(1/3); w = x - 0.078*x^5/(1 + e); E = M + e*(3*w - 4*w^3); % Two times applying 5th-order Newton correction for k = 1:2 f = (E - e*sin(E)- M); fd = 1 - e*cos(E); f2d = e*sin(E); f3d = -e*cos(E); f4d = e*sin(E); E = E -f/fd*(1 + f*f2d/(2*fd^2) + f^2*(3*f2d^2 - fd*f3d)/(6*fd^4)+... (10*fd*f2d*f3d - 15*f2d^3-fd^2*f4d)*f^3/(24*fd^6)); end E_nit(i,j) = E; i = i + 1; end j = j + 1; end t_nit = toc; fprintf('Calculation time = %4.6f [s] \n',t_nit); Calculation time = 0.149860 [s]
Newton Iterative Solver
tic; j = 1; eps = 1E-15; % Tolerance for M = 0:0.01:pi i = 1; for e = 0:0.01:1 En = M; Ens = En - (En-e*sin(En)- M)/(1 - e*cos(En)); while ( abs(Ens-En) > eps ) En = Ens; Ens = En - (En - e*sin(En) - M)/(1 - e*cos(En)); end; E_it(i,j) = Ens; i = i + 1; end j = j + 1; end t_it = toc; fprintf('Calculation time = %4.6f [s] \n',t_it);
Calculation time = 0.893673 [s]
Error = log10(abs(E_it - E_nit)); e = 0:0.01:1; M = (0:0.01:pi); contourf (M,e,Error); xlabel('M - Mean anomaly [rad]'); ylabel('e - Eccentricity'); title('Non-Iterative Solver,Error in E calculation on Log scale [rad] '); colorbar;[6] S. Mikkola, “A Cubic Approximation for Kepler’s Equation,” Celestial Mechanics and Dynamical Astronomy,Vol. 40, 1987, pp. 329–334.
Kepler’s equation Solver
In this example we will show how to solve Kepler’s equation M = E-e*sin(E)
clear all; clc; M = 0.3458; % Mean anomaly [rad] eps = 1E-9; % Tolerance e = 0.1; % Eccentricity En = M; Ens = En - (En-e*sin(En)- M)/(1 - e*cos(En)); while ( abs(Ens-En) > eps ) En = Ens; Ens = En - (En - e*sin(En) - M)/(1 - e*cos(En)); end; E = Ens; % Eccentric anomaly E [rad]. fprintf('Eccentric anomaly E = %4.4f [rad]\n',E);
Eccentric anomaly E = 0.3832 [rad]
Small Satellites 