Error in Orbital Motion Code

[Nick Sorabella]

Hello!

I've put together a little orbital motion code but it seems to have trouble running in glowscript.  This is strange as I can run it in Jupyter notebook (albeit slowly) so I'm not really sure what the error is and, or, how to fix it.  Any thoughts?  The code is attached below:

Web VPython 3.2
m1 = 1.988e30
m2 = 1.988e30
mu = (m1*m2)/(m1+m2)
G = 6.67e-11
a = 1.5e11
ecc = 0
aph = a*(1+ecc)
x0 = (m2*aph)/(m1+m2)
X0 = -(m1*aph)/(m1+m2)
y0 = 0
Y0 = 0
vel1 = ((G*mu/a)*(m2/m1)*((1-ecc)/(1+ecc)))**0.5
vel2 = -((G*mu/a)*(m2/m1)*((1-ecc)/(1+ecc)))**0.5 * (m1/m2)
canvas(background=color.white,height=800, width=1520)
bh = sphere(pos=vector(x0,0,0), radius=0.1e10, color=color.black,
             make_trail=True, trail_type='curve', interval=1000)
bh.p = vector(0, vel1, 0) * m1

star = sphere(pos=vector(X0,0,0), radius=1e10, color=vector(0.98039,0.3921568,0.12549),
              make_trail=True, trail_type='curve', interval=1000)
star.p = vector(0,vel2,0) * m2
dt = 10
while True:
    rate(200000)
    r = star.pos-bh.pos
    F = 6.67e-11 * m1 * m2 * r.hat/mag(r)**2
    bh.p = bh.p + F*dt
    star.p = star.p - F*dt
    bh.pos = bh.pos + (bh.p/m1) * dt
    star.pos = star.pos + (star.p/m2) * dt

[Bruce Sherwood]

Ouch. Change

vel2 = -((G*mu/a)*(m2/m1)*((1-ecc)/(1+ecc)))**0.5 * (m1/m2)

to this:

vel2 = (-1)*((G*mu/a)*(m2/m1)*((1-ecc)/(1+ecc)))**0.5 * (m1/m2)

Python supports "operator overloading" but JavaScript does not. Web VPython uses a library that converts x*y to x["*"](y). When x["*"](y) is executed, the arguments x and y are inspected, and if (for example) x is scalar and y is a vector it multiplies the vector y by the scalar x. Unfortunately, that library apparently makes a mistake in how the initial minus sign is treated. If you choose "Share or export this program" you'll find that lines 13 and 14 have been transpiled to JavaScript like this, and line 14 is just wrong, alas:

    vel1 = Math.pow((G["*"](mu)["/"](a)["*"](m2["/"](m1))["*"](1["-"](ecc)["/"](1["+"](ecc)))), .5);
    vel2 = Math.pow(G["*"](mu)["/"](a)["*"](m2["/"](m1))["*"](1["-"](ecc)["/"](1["+"](ecc)))["-u"](), .5)["*"](m1["/"](m2));

When running from installed Python, these transformation are not needed, because Python supports operator overloading.

Thanks for the clear example. This came up once before, and at the time I couldn't figure out how to fix the error in the way the minus is handled. It's possible that, for example, converting an initial -( to (-1)* would fix the problem. I should try that.

Bruce

[Nick Sorabella]

Thank you very much for the help Bruce!  It is amazing how the smallest of mistakes can cause the most annoying of issues!  It can be like finding a needle in a haystack.  I also appreciate the insight into how the issue was working behinds the scenes.

-Nick

[Bruce Sherwood]

I should share how I identified the problem. At many places in the program I used print() to display significant quantities at that point in the program. When I saw some variable having the value "NaN" ("not a number", a JavaScript entity) it was easy to identify in how the program failed.

Bruce

[Bruce Sherwood]

The bug has been fixed, and the original form of your program now runs correctly. Thanks again for the clear report.

Bruce