# Circular Orbit Velocity in Vector Components

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

-
What has your question to do with Mathematica? –  halirutan Nov 6 '13 at 13:08
This question should probably be migrated to physics.stackexchange.com –  halirutan Nov 6 '13 at 13:09
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. –  Sjoerd C. de Vries Nov 6 '13 at 13:32
@SjoerdC.deVries Yes, I have seen that too, but on the other hand he/she is not using Mathematica syntax. –  halirutan Nov 6 '13 at 13:44
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. –  user7388 Nov 6 '13 at 13:49
show 3 more comments

## 1 Answer

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


-
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? –  user7388 Nov 7 '13 at 10:51
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\}$ –  Suba Thomas Nov 7 '13 at 15:23
add comment