#N-body gravity simulation prediction isn't accurate

4 messages · Page 1 of 1 (latest)

north lichen
#

N-body gravity simulation prediction isn't accurate

#

gravity code:


var dt = 1 / global.SUBSTEPS;

// --- Compute acceleration at current position ---
var accel  = scr_get_accel(x, y, id);
var ax     = accel[0];
var ay     = accel[1];

// --- Velocity Verlet Integration ---
// half velocity update
vx += 0.5 * ax * dt;
vy += 0.5 * ay * dt;

// position update
x += vx * dt;
y += vy * dt;

// recompute acceleration at new position
var accel1 = scr_get_accel(x, y, id);
var ax1    = accel1[0];
var ay1    = accel1[1];

// finish velocity update
vx += 0.5 * ax1 * dt;
vy += 0.5 * ay1 * dt;```
#

function scr_get_accel:

/// @function scr_get_accel(x0, y0, self_id)
/// @desc Compute gravitational acceleration at (x0, y0) from all gravity_mass objects.
/// @param {real} x0
/// @param {real} y0
/// @param {id}   self_id (the calling instance, so it doesn’t attract itself)
/// @return [ax, ay]

function scr_get_accel(x0, y0, self_id) {
    var G    = global.G;
    var soft = global.SOFTENING;
    var ax   = 0;
    var ay   = 0;

    with (all) {
        if (id != self_id &&
            object_index != global.OBJ_Controller &&
            object_index != global.OBJ_Protagonist) 
        {
            if (variable_instance_exists(id, "gravity_mass")) {
                var dx = x - x0;
                var dy = y - y0;
                var r2 = dx*dx + dy*dy + soft;
                var r  = sqrt(r2);

                if (r > 0.0001) {
                    var a = (G * gravity_mass) / r2;
                    ax += a * (dx / r);
                    ay += a * (dy / r);
                }
            }
        }
    }

    return [ax, ay];
}