2
\$\begingroup\$

So I am pretty new to JavaScript and I thought I might make an animation demonstrating a planet orbiting a sun in an elliptical orbit using Kepler's laws.

I am struggling to get the speed right. So I want the speed of the planet to be slower further away from the sun, and faster closer to the sun.

Within this method I use Keplers laws to calculate the x and y coordinate. I have also calculated the (real) velocity at each point in its trajectory (this may need to be scaled appropriately as I hope to use real astronomical data).

I tried to use a scaled version of the velocity. The velocity values produced were in the range of 10,000 to 80,000 - so I thought of using that to update theta using the following:

  this.speed = function(){
    this.theta += this.velocity/10000 ;
  }

However this didn't work and the planet disappears.

I've got full code is below:

HTML:

<canvas id="mybox"></canvas>

CSS:

canvas {
  border: 2px solid #000;
  background-color: #102233;
}

JS:

var canvas = document.getElementById('mybox');
canvas.width = window.innerWidth ;
canvas.height = window.innerHeight ;

var ctx = canvas.getContext('2d');

function sun(){
// Draws Sun at center
  ctx.beginPath() ;
  ctx.fillStyle = "orange";
  ctx.arc(canvas.width/2, canvas.height/2, 40, 0, Math.PI*2);
  ctx.fill();
};

function planet(M, G, T, e) {
  this.x = 0 ;
  this.y = 0 ;
  this.radius = 20 ;
  this.theta = 0 ;
  this.dtheta = 0.5 ;
  this.vel = 0 ;
  
  this.updateTheta = function() {
    if (this.theta > Math.PI){
      this.theta += this.dtheta ;
      this.dtheta += 0.001 ;

    }
    else if (this.theta < Math.PI){
      this.theta +=this.dtheta ;
      this.dtheta -= 0.001 ;
    }
    if(this.theta>=2*Math.PI || this.theta < 0) {
      this.theta = 0;
    }
  }
  
  
  this.position = function(theta) {
// Calculates position
    var R = Math.pow(((Math.pow(T,2))/((4*Math.pow(Math.PI, 2))/(G*M))),(1/3)) ;
    var a = R/(1-e) ;
    var radius = a*(1-Math.pow(e,2))/(1+e*Math.cos(theta)) ;
    this.x = radius*Math.cos(theta)*10e-10 + canvas.width/2 ;
    this.y = radius*Math.sin(theta)*10e-10 + canvas.height/2;
    this.vel = Math.sqrt(G*M*((2/radius)-(1/a))) ;
  }
  
  this.draw = function() {
    ctx.beginPath();
    ctx.fillStyle = 'green';
    ctx.arc(this.x, this.y, this.radius, 0, Math.PI*2);
    ctx.fill();
    
    this.position(this.theta) ;
    this.updateTheta() ;
  }
  
}

// Values to display elliptical orbit
var G = 6.e-11;
var M = 5e+30;
var T = 1e+7;
var e = 0.6;

function update(){
  ctx.clearRect(0, 0, canvas.width, canvas.height); 
  sun();
  earth.draw();
}

let earth = new planet(M, G, T, e)
setInterval(update, 50);

Any help or comments will be appreciated. Thanks!

\$\endgroup\$
5
  • \$\begingroup\$ Other than looking at applied math, can you be a little more specific about what you want from a review of the code snippet? \$\endgroup\$ – Drew Reese Apr 30 at 3:15
  • \$\begingroup\$ I was hoping for help with ensuring the animation displayed a scaled version of the velocity. Apologies if this is not the correct way of asking a question. This has just been bugging me for a while and was unsure of who to ask \$\endgroup\$ – Yasser Apr 30 at 3:32
  • 2
    \$\begingroup\$ Your total program isn't that big (as in, we've seen way worse). Could you include the full code into the question? We tend to only review what's in the question itself and for licensing (and other) reasons we encourage the author of the post to do that themself. \$\endgroup\$ – Mast Apr 30 at 5:28
  • \$\begingroup\$ Welcome to Code Review! Can you confirm that the code is complete and that it functions correctly? If so, I recommend that you edit to add a summary of the testing (ideally as reproducible unit-test code). If it's not working, it isn't ready for review (see help center) and the question may be deleted. \$\endgroup\$ – Toby Speight Apr 30 at 10:49
  • \$\begingroup\$ @Mast Thanks for letting me know. I've updated the question to include everything \$\endgroup\$ – Yasser Apr 30 at 13:52
3
\$\begingroup\$

Review

Orbit

You are incorrectly calculating the position of the planet.

The eccentric anomaly which you attempt to calculate with a linear, integrated 2nd derivative, in the function this.updateTheta will never work as the solution in non linear. That said the example below uses a linear approximation.

There is no easy way to change your code such that time steps can be scaled while maintaining the same orbit. It is best to use the standard orbital solution.

Animation

Never use 'setInterval' for animation

For the best result when animating via JS in the browser use the animation callback requestAnimationFrame (rAF)

It will call the main render loop function at 60FPS (If your render function takes no more than ~80% of 1/60th second (12 - 13ms)) Please read the documentation for rAF as there are many caveats.

Example of basic animation loop

    requestAnimationFrame(renderLoop);     // request first frame

    function renderLoop() {
        ctx.clearRect(0, 0, W, H);         // clear canvas ready for next frame render

        ... // call render function

        requestAnimationFrame(renderLoop); // Request Next frame
    }

Readable

Your calculations are very hard to read and you are using some outdated calls.

  • Use the ** operator rather than Math.pow or Math.sqrt. Eg 2 ** 2 is two squared and 2 ** 0.5 is the square root of two.

  • The main reason that your calcs are hard to read is because you don't space out the operators. Always use spaces between operators.

    The line...

    var R = Math.pow(((Math.pow(T,2))/((4*Math.pow(Math.PI, 2))/(G*M))),(1/3)) ;
    

    is far easier to read as...

     const R = (T ** 2 / (4 * Math.PI ** 2 / (G * M))) ** (1 / 3);
    

    Note that variables that do not change should be constants.

Orbits and controlling their animation periods.

Note I use radians

Kepler's laws modify a circular orbit into an ellipse such that the line from the planet (Earth) to the parent (Sun) sweeps equal areas in equal time.

We define the ellipse by the semi-major axis, the eccentricity, and either the closest (periapsis) or furthermost (apoapsis) points in the orbit. Both points are on the semi-major axis

We want to crate a function that gives the planets position at a given time. We can then scale that time to control the rate that the animation plays out.

  • For a circular orbit the object must travel 2 * MATH.PI.

  • Using ms we can define the period of the orbit. eg 1 second is 1000ms

  • We can then scale the orbit position as a function of time.

    Example

    • The orbit position scale is scale = (period) => 2 * MATH.PI / period
    • With an orbit of 1 second the position is pos = (ms) => scale(1000) * ms where ms is the current time. The pos returned is the mean anomaly M

Kepler's law

  • The orbital position is called the mean anomaly M in radian (range -PI to PI)

  • When traveling on a non circular orbit the mean anomaly M is advanced or retarded depending on the orbital distance. That value is the eccentric anomaly E in radians

  • The ellipse is defined by eccentricity e range 0 <= e < 1

  • Kepler states that M = E - e * Math.sin(E)

    • We know M and e, but not E so we must solve the above equation in terms of M. This is non trivial so we must approximate.

      To do this we use Newtons method to solve M = E - e * Math.sin(E) for M with a linear approximation of the non linear part e * sin(E) see function eccentricAnomaly in Demo

  • With E we can then calculate the 2D position of the body using e and the distance to apoapsis r

      const x = r * (Math.cos(E) - e);
      const y = r * Math.sin(E) * (1 - e ** 2) ** 0.5; 
    
  • This 2D position can then easily be aligned to the semi-major axis, and if needed rotated into 3D to match the inclination, position of ascending node and argument of periapsis to fully define any 2 body orbit (Newtonian).

DEMO

All that above was clear as mud I am sure, thus a Demo.

The code below animates 3 objects (two 2 body solutions) at a specified orbital period of 1 and 2 seconds. The time is given by rAF as first argument given to the callback renderLoop.

Time 0 is set to apoapsis which is defined by the start positions of the planet relative to the parent body. The orbit is rotated in 2D to align to the apoapsis.

The eccentricity of planets e is 0.5

requestAnimationFrame(renderLoop);
const ctx = canvas.getContext("2d");
ctx.lineCap = "round";
var frameCount = 0, startTime, prevTimeScale = timeScale.value, lastTime;
const W = canvas.width, H = canvas.height, CX = W * 0.25, CY = H * 0.5;

const SUN_RADIUS = 20, SR = SUN_RADIUS;               // in pixels
const PLANET_SCALE = 0.4, PS = PLANET_SCALE;          // in sun radi
const FIRST_PLANET_DIST = 4, FPD = FIRST_PLANET_DIST; // in sun radi
const ECCENTRICITY = 0.5, E = ECCENTRICITY;           // 0 circlular
const KEPLER_PERIOD = 1000, KP = KEPLER_PERIOD;       // ms per orbit

const MBLUR_LINE_GRAD = 8; // Steps to grade blur line
const FRAME_RATE = 60; // Max 60 Frames per second. Only use valid values
                       // Valid vals 60, 30, 20, 15, 12, 10, 6, 5, 4, 3, 2, 1
const FRAME_SKIP = 60 / FRAME_RATE;
Math.TAU = Math.PI * 2;
Math.R90 = Math.PI * 0.5;


const calcVectors = (A, B) => {
    const dx = B.x - A.x;
    const dy = B.y - A.y;
    const rSqr = dx * dx + dy * dy, r = rSqr ** 0.5;    
    const nx = dx / r;
    const ny = dy / r;
    return {nx, ny, rSqr, r, axis: Math.atan2(ny, nx)};
}
const eccentricAnomaly = (M, e) => {  // newtons method to solve M = E-e sin(E) for E
    var d, E = M; // guess E
    do {
        d = (E - e * Math.sin(E) - M) / (1 - e * Math.cos(E));
        E -= d;
    } while (d > 1e-6);
    return E;
}


const Planet = {
    x:  0, y:  0, vx: 0, vy: 0, ox: 0, oy: 0,
    init() { return this }, 
    update() { },
    draw(ctx, lastSubFrame) {
        const A = this;
        const dx = A.ox - A.x;
        const dy = A.oy - A.y;
        const dist = (dx * dx + dy * dy) ** 0.5;    
        ctx.setTransform(1, 0, 0, 1, CX + A.x, CY + A.y);
        ctx.fillStyle = ctx.strokeStyle = A.col;
        if (dist > 0.1 && MBLUR_LINE_GRAD > 0) {
            ctx.globalAlpha = (1 / (dist - 0.1)) ** 1.4;
            ctx.beginPath();
            ctx.strokeStyle = A.col;
            const r = A.radius * 2;
            ctx.moveTo(dx, dy);
            ctx.lineTo(0, 0);
            const MBLG = MBLUR_LINE_GRAD;
            let i = MBLG;
            while (i--) {
                ctx.lineWidth = Math.cos(((i - 1) / MBLG) * Math.R90) * r;
                ctx.stroke();
            }
        } else if (lastSubFrame) {
            ctx.beginPath();  
            ctx.arc(dx / 2, dy / 2, A.radius, 0, Math.TAU);
            ctx.fill();
        }
        ctx.globalAlpha = 1;
        A.ox = A.x;
        A.oy = A.y;
    }, 
};
const KeplerPlanet = {
    period: 1,
    init(e, B) {  // B at relative position at time 0. A is at apoapsis        
        const A = this;
        A.period =  Math.TAU / A.period;
        A.e = e;
        Object.assign(A, calcVectors(A, B)); 
        A.px = Math.cos(A.axis); // Transform to Rotate semi major axis
        A.py = Math.sin(A.axis);
        return A;
    },
    update(time) {
        const A = this;
        const E = eccentricAnomaly(A.period * time, A.e);
        const x = A.r * (Math.cos(E) - A.e);
        const y = A.r * Math.sin(E) * (1 - A.e ** 2) ** 0.5;   
        A.x = x * A.px - y * A.py;
        A.y = x * A.py + y * A.px;
    },    
};
const FixedPlanet = {
    radius: 0,
    col: "#EC8",
};

function createOrbitObj(type, x, y, radius, col, period, e, orbits) {
    return  ({...Planet, ...type, x, y, ox: x, oy: y, period, radius, col }).init(e, orbits);
}
const sun = createOrbitObj(FixedPlanet , 0, 0, SR, "#FF8");
const p1 = createOrbitObj(KeplerPlanet,  FPD * SR, 0, SR * PS, "#6C1", KP, E, sun);
const p9 = createOrbitObj(KeplerPlanet,  FPD * SR * 2, 0, SR * PS, "#C61", KP * 2, E, sun);
const system = Object.assign([sun ,p1, p9], {      
    update(time) { 
        var i = this.length - 1;
        while (i) { this[i--].update(time) }
    },      
    draw(ctx, lastSubFrame) {
        for (const A of this) { A.draw(ctx, lastSubFrame) }
    },
});

function renderLoop(time) {
    var tT, i, pTS = prevTimeScale; // tT is totalTime
    if (frameCount++ % FRAME_SKIP === 0) {
        startTime = startTime ?? time;
        tT = time - startTime;

        ctx.setTransform(1, 0, 0, 1, 0, 0);
        ctx.clearRect(0, 0, W, H);

        const inputTimeScale = timeScale.value;
        if (inputTimeScale !== pTS) {
            startTime = time - tT * pTS / inputTimeScale;
            pTS = prevTimeScale = inputTimeScale;
            tT = time - startTime;
        }
        const dT = tT * pTS - lastTime; // delta time
        i = 0;
        while (i < MBLUR_LINE_GRAD && lastTime) {
             const current = lastTime + dT * (i++ / MBLUR_LINE_GRAD);
             system.update(current);
             system.draw(ctx, false);
        }
        
        system.update(tT * pTS);
        system.draw(ctx, true);
        lastTime = tT * pTS;
    }
    requestAnimationFrame(renderLoop);
}
canvas { background: black; }
input { position: fixed; top: 10px; left: 10px; color: white; width: 340px }
<canvas id="canvas" width="350" height="300"></canvas>
<input type="range" id="timeScale" min="0.01" max="4" step ="0.01" value="1">

Why retro? The view is up from under the ecliptic (to matche canvas coords), orbits are thus retrograde to the solar system's direct (prograde)

Extras

I have added some additional features to demo.

  • Time Scale. The range bar changes the time scale that lets you change the rate of time. Time is a value based on a start time. When scaling time that start time must be change so that the planet does not skip time.

  • Frame Rate. There is a value called FRAME_RATE which sets the frame rate. Max 60. This rate is a requested rate, the actual rate will vary depending on device and CPU/GPU load, but over time it will average the requested rate.

  • Motion Blur. To improve the look when the time scale is large there is a simple motion blur applied to the moving planets.

    The blur is just a linear interpolation between sub time steps and thus only works when movement per sub step is less than about 0.5 radians.

    MBLUR_LINE_GRAD controls how the blur is drawn and the number of sub time steps are used. A value of 0 will turn off motion blur. Values 1 or greater will increase the persistence. Current value is 8 and is tuned for the current frame rate of 60FPS.

\$\endgroup\$
3
  • \$\begingroup\$ Thank you for this thorough review and explanation! It really helped! \$\endgroup\$ – Yasser May 1 at 16:03
  • \$\begingroup\$ just happened to have this post-it note on your refrigerator, eh? \$\endgroup\$ – radarbob May 1 at 23:33
  • \$\begingroup\$ Read Kepler's Witch. "Witch" refers to Kepler's mother who at age 70 was convicted of witchcraft. His fame helped her avoid physical torture or death. Kepler's life was hellacious: born with physical abnormalities, insanity in the family, a wanderlust absent father, pestilence killing children and wife, frequent near destitution, a pious Lutherian perpetually harangued by the church. and that odd-ball personality of a god-like genius. The final insult was a skirmish of The 30 Years War destroyed his grave marker. Today no one knows where (within the cemetery) Johann Kepler is buried. \$\endgroup\$ – radarbob May 2 at 0:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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