1 Goal In orbital mechanics, the position of an orbiting body at time t, e.g., a

Question

1 Goal In orbital mechanics, the position of an orbiting body at time t, e.g., a satellite orbiting the earth or a planet revolving around the sun, is determined by the Kepler's equation M(t) 3 E(t) esin E(t), (1) which relates the eccentric anomaly E, with the eccentricity of the orbit e and the average anomaly M. For validation purposes, astronomers want to build a MATLAB toolbox to find roots of the Kepler's equation. 2 Detailed description 2.1 Mathematical Model The position of the satellite at time t is given by the true anomaly v(t r(t) r cos v(t) y(t) v(t) T Sin where km 1+ecos, v' The true anomaly v relates to (1) through the equation e cos E (2) arccos e cos E 1 and the dependence on time is established through M by (3) 2T where T is the period of the orbit. To determine the v, take into account that VE 10, T] if E E 0, TJ, and Then, for all i E 10, n) find M(ti) from (3), you have to solve the Kepler's equation (1), and finally find v(ti) using (2). To solve the Kepler's equation, you shall combine Newton's method and bisection method: a possible way of doing that is described in Algorithm 1

Explanation / Answer

function [root, a, b] = bis(a,b,tolb,func,pars) %routine to find a root of a function using bisection method a=a+pars; %pars stores all constants in Kepler's equation b=b+pars; %func=@(x)(x-e*sin(x)-Mi); %Kepler’s equation to find root E(t) while abs(b-a)>tolb p=a+(b-a)/2; %compute midpoint p if (func(a)==0 || func(b)==0 || func(p)==0 || (b-a)/20 %same sign indicates root not in interval a=p; %consider second interval else b=p; %consider first interval end end root=p; end function [root, error, found, iter] = nwt(seed,toln,nmax,func,dfunc,pars) %routine to find root of a function using Newton's method %found is a logical indicating whether root is found i=1; p0=seed; %initial approximation while ipi; pars=pi; end toln=1e-16; %suppose we want the error in the coordinates of the orbit %to be smaller than 10^-6 nmax=100; %test to validate code func=@(x)(x-e*sin(x)-Mi); %Kepler’s equation to find root E(t) dfunc=@(x)1-e*cos(x); %derivative of f Ei=bisnwt(a,b,tolb,toln,nmax,func,dfunc,pars); orbit(i+1,2)=Ei; %use Ei to get true anomaly v vi=acos((e-cos(Ei))/(e*cos(Ei)-1)); %use v to get x,y position if Mi>pi; vi=abs(vi-2*pi); end orbit(i+1,3)=vi; u=3.986012e5; %constant mu to get a and r a=nthroot(u*((T*3600/(2*pi))^2), 3); ri=a*(1-e^2)/(1+e*cos(vi)); orbit(i+1, 4)=ri; xi=ri*cos(vi); yi=ri*sin(vi); orbit(i+1 ,5)=xi; orbit(i+1 ,6)=yi; end end

Related Questions

Navigate

• Browse (All)
• Browse 1
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.