I'm trying to find out how to calculate the initial velocity needed to keep a spacecraft in a circular orbit around the Earth given its initial position and an initial angle theta in a Cartesian coordinate system. I know the formula for the magnitude of a circular orbit is $v_{\mathrm{circular}} = \sqrt{GM/r_{\mathrm{orbit}}}$ but this only works nicely as it is (i.e. in its magnitude form) if the orbital radius vector makes an angle of 0, 90 or 270 degrees in a Cartesian coordinate system. What I want to do is be able to use the formula for any angle $\theta$.

As an example, let's say I have a spacecraft at a radius of 300000 meters above the Earth's surface and its radius vector makes an angle of 45 degrees in an x-y Cartesian coordinate system. Then we have $r = r_{\mathrm{earth}} + 300000$, whereby $r_x = r \cos(45^\circ)$ and $r_y = r \sin(45^\circ)$. Here's a little picture I made to better understand the problem:

enter image description here

If theta = 45 degrees, then I get a nice circular orbit, but any other value of theta gives me very undesirable results (i.e. the spacecraft crashes into the planet or flies of into the cosmos)

Does anyone know what I'm doing wrong?

G = 6.672*10^-11;
m[1] = AstronomicalData["Earth", "Mass"];
tmax = 20000;
r[1] = AstronomicalData["Earth", "Radius"];
rx = (r[1] + 300000 ) Cos[45 Degree]
ry = (r[1] + 300000 ) Sin[45  Degree]
vx = Sqrt[(G  m[1])/(r[1] + 300000)] Cos[45  Degree]
vy = Sqrt[(G  m[1])/(r[1] + 300000)] Sin[45  Degree]

soln = NDSolve[{
x''[t] == -((G  m[1] x[t])/(x[t]^2 + y[t]^2)^(3/2)),
y''[t] == -((G  m[1] y[t])/(x[t]^2 + y[t]^2)^(3/2)), 

x[0] == rx, y[0] == ry, 
x'[0] == vx, y'[0] == -vy}, {x[t], 
y[t]}, {t, 0, tmax} , MaxSteps -> 1000000, 
Method -> "StiffnessSwitching"]

Show[ParametricPlot[Evaluate[{x[t], y[t]} /. soln], {t, 0, tmax}, 
AxesLabel -> {x, y}, PlotStyle -> Automatic, PlotRange -> Full, 
ImageSize -> Large], Graphics[{Green, Disk[{0, 0}, r[1]]}]]
  • $\begingroup$ What has your question to do with Mathematica? $\endgroup$
    – halirutan
    Nov 6, 2013 at 13:08
  • 2
    $\begingroup$ This question should probably be migrated to physics.stackexchange.com $\endgroup$
    – halirutan
    Nov 6, 2013 at 13:09
  • 3
    $\begingroup$ Hold on guys, the OP has asked more question on space physics and Mathematica. So this may quite well be Mathematica related. However, I'd like to ask the poster to provide his/her code so far to make this a real Mathematica question. Voting to close while awaiting additional information. $\endgroup$ Nov 6, 2013 at 13:32
  • $\begingroup$ @SjoerdC.deVries Yes, I have seen that too, but on the other hand he/she is not using Mathematica syntax. $\endgroup$
    – halirutan
    Nov 6, 2013 at 13:44
  • $\begingroup$ Hi there, sorry guys. I'm so used to posting Mathematica questions that I forgot to post in the correct section! I'll edit my above post with the code. $\endgroup$ Nov 6, 2013 at 13:49

2 Answers 2


More physics than Mathematica as pointed out in the comments..., but here's my shot at it.

G = 6.672*10^-11;
m[1] = AstronomicalData["Earth", "Mass"];
tmax = 20000;
r[1] = AstronomicalData["Earth", "Radius"];

sol = NDSolve[{θ'[t] == 1/(r[1] + 300000) Sqrt[(G m[1])/(r[1] + 300000)],  θ[0] == 50 °}, 
θ[t], {t, tmax}];

Show[ParametricPlot[{(r[1] + 300000) Cos[θ[t]], (r[1] + 300000) Sin[θ[t]]} /. sol[[1]], 
{t, 0, tmax}], Graphics[{Green, Disk[{0, 0}, r[1]]}]]

enter image description here

  • $\begingroup$ Thanks Suba Thomas, thankfully I managed to come up with that as well last night. I guess we use sin instead of cos for the x-component of velocity and cos instead of sin for the y-component of velocity because sin and cos are out of phase by 90 degrees, which would make the velocity vector perpendicular to the radial vector since the radial vector uses cos for its x-component and sin for its y-component? $\endgroup$ Nov 7, 2013 at 10:51
  • $\begingroup$ According to the figure the velocity at $\{(r(1)+300000) \cos (\theta ),(r(1)+300000) \sin (\theta )\}$ is $\left\{-v_c \sin (\theta),v_c \cos (\theta )\right\}$ $\endgroup$ Nov 7, 2013 at 15:23

There is a more general way to solve this. You say that the velocity is $\sqrt{\frac{GM}{r_\rm{orbit}}}$, but let's recognize that this is just a special case of the Vis-viva equation $v = \sqrt{GM(\frac{2}{r}-\frac{1}{a})}$ where $r$ is the distance between the large body and the small body, and $a$ is the semimajor axis. A body with the velocity given by this equation tangent to the orbit with semimajor axis $a$ will stay in orbit. An ellipse can be written as $(x, y) = (a \rm{cos}(\theta), b \rm{cos}(\theta))$, where b is the semiminor axis. Taking the derivative of each component with respect to $\theta$ yields $(-a\rm{sin}(\theta), b\rm{cos}(\theta))$. This is the tangent we were looking for.

Now we have the velocity and a direction for any point $(x,y)$, which is enough to solve for any orbit we might want. But since we're integrating numerically we need to know for how long we want to integrate, so it's interesting to know how much time it takes for the smaller body to travel around its orbit. This is called the period time and it is given by Kepler's third law as $T=2\pi\sqrt{\frac{a^3}{GM}}$.

For simplicity we'll be working in an earth-centric frame. Let's start by defining the gravitational constant, the mass of the earth and its radius. Let's also define a graphics object that will represent the earth in our plot:

G = 6.672*10^-11;
M = QuantityMagnitude@AstronomicalData["Earth", "Mass"];
R = QuantityMagnitude@AstronomicalData["Earth", "Radius"];
earth = Graphics[{
    Disk[{0.02, -0.02} R, R],
    RGBColor[0.6470588235294118`, 0.7450980392156863`, 
     0.9058823529411765`], Disk[{0, 0}, R]

The problem itself can be solved by NDSolve alone. In addition to the formulas for the velocity, the direction of the velocity and the period time described above we'll also use the law of gravitation to get the acceleration $\rm{acc} = -\frac{GM}{r^2}$.

orbit[{x0_, y0_}, a_, b_] := {
  2 Pi Sqrt[a^3/(G M)],
    {x''[t], y''[t]} == -G M {x[t], y[t]}/Norm[{x[t], y[t]}]^3,
    {x'[0], y'[0]} == 
     Sqrt[G M (2/Norm[{x[0], y[0]}] - 1/a)] Normalize[{-a Sin[
          ArcTan[x[0], y[0]]], b Cos[ArcTan[x[0], y[0]]]}],
    {x[0], y[0]} == {x0, y0}
   {x, y},
   {t, 0, 2 Pi Sqrt[a^3/(G M)]}

The function returns a list of the form {period time, solution}. Finally, to plot the solution we can define plotOrbit:

plotOrbit[{x0_, y0_}, a_, b_] := 
 With[{orbitsol = orbit[{x0, y0}, a, b]},
  ParametricPlot[{x[t], y[t]} /. Last@orbitsol, {t, 0, First@orbitsol}]

Which in turn can be used as:

 plotOrbit[{1.5 R, 0}, 5 R, R]

enter image description here

My example is of a satellite whose orbit around the earth is a an ellipse with a semimajor axis $5R$, semiminor axis $R$ and the periapsis - the closest distance to the earth - is $1.5 R$ from the center of the earth. $R$ is the radius of the earth. (Don't think that the x coordinate is always the periapsis, this is a special case.) However, by setting $a = r$ you can easily get a circular orbit thus proving that we've solved the problem in this question but in a more general way.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.