Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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]]}]]
share|improve this question
    
What has your question to do with Mathematica? –  halirutan Nov 6 '13 at 13:08
2  
This question should probably be migrated to physics.stackexchange.com –  halirutan Nov 6 '13 at 13:09
1  
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

up vote 1 down vote accepted

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

share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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