Circular Orbit Velocity in Vector Components

Question

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:

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]]}]]

Answer 1

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]]}]]

Answer 2

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[{
    Opacity[0.6],
    Disk[{0.02, -0.02} R, R],
    Opacity[1],
    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)],
  NDSolve[{
    {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:

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

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.