I'm having trouble trying to have Mathematica tell me when the period of the orbit of my system is. Namely, I'm trying to use Mathematica to find the orbital period of the Earth around the Sun. Here's my code:
eq = {r''[t] == b/(r[t].r[t]) r[t]/Norm[r[t]]};
inits = {r[0] == {Ro, 0}, r'[0] == {0, Vo}};
event = WhenEvent[r[t] == {Ro, -0.001}, Period = t];
SEperiod = NDSolve[{eq,inits,event},r,{t,0,9.4*10^7}]
It seemed to run fine, but Period
is never given a value from t
. I tried it again by putting it in NDSolve
explicitly:
SEperiod =
NDSolve[
{
r''[t] == b/(r[t].r[t]) r[t]/Norm[r[t]],
r[0] == {Ro, 0},
r'[0] == {0, Vo},
WhenEvent[r[t] == {Ro, -0.001}, Period = t]
},
r,
{t, 0, 9.4*10^7}
]
Which came with no error messages; but once again, Period
wasn't saved. I chose the condition r[t]=={Ro,-0.001}
as the triggering for the event because I expect it to be at the position approximately right before the Earth completes a full orbit. However changing around the values still sets no value for Period
. How do I get it to be assigned a value like I want?
Edit: Just to make things clear, Ro 1 AU in meters, b is the product of the Gravitational Constant, the mass of the sun, and Vo is the initial velocity of the Earth.
b
,Ro
andVo
. $\endgroup$r[t] == {Ro, -0.001}
. Try using something likeNorm[r[t]] < 10^-3
. Some experimentation may be necessary to find the right number for comparison. $\endgroup$