LinearSinkForce doesn't correctly sim gravity.

Awhile ago I experimented with the LinearSinkForce as a possible way of implementing gravity for a Solar System simulation. The documentation states “Attractor force. Think black hole.” which would imply that it was intended to accurately simulate gravity. However, I found that it did not actually create a true gravity well but instead generated a bowl shaped attracting force, meaning that the strength of the attraction increased as the object moved further away from the center instead of falling off.

The most obvious consequence of this is that each object (when starting at rest) will fall to the center point in the same amount of time, arriving at the origin simultaneously. This would be inverse to the effect of a true gravity well.

Though this did not fit the description of “think black hole”, I wasn’t sure if this was a bug with Panda or just an error in the documentation. The force includes fall_off_type, radius, and amplitude as parameters, but no combination seemed to accurately recreate gravity.

Below is a demo that includes the LinearSinkForce and a simple implementation of a true gravity well to show the difference. Sorry for the bloated sphere/planet creation code, but it was easier for me to just scale down from my actual code than re-implement using models and this should work by just pasting it into a file. There are some settings that can be changed to show different effects. If this is a bug I’ll submit a proper bug report, but I’d thought I’d get some opinions first.

# ====================================
# LinearSinkForce vs Gravity Well demo
# ====================================

# Panda3d v1.8

# Demonstrates that the LinearSinkForce creates a bowl shaped attraction
# rather than a realistic gravity well, the major symptom of a bowl shaped
# force is that all objects, regardless of their initial distance (when 
# starting at rest) will arrive at the center at the same time:
#
#                Bowl                         Well
#        .                   .      .                       .
#        .                   .           .             .
#         .                 .               .       .
#           .             .                   .   .
#              .       .                       . .
#                  .                            .
#
#
# "G_MODE" controls which which implentation of gravity the demo uses.
#
G_MODE = "sink"    # use Panda3d's LinearSinkForce; results in "bowl" shape.
#G_MODE = "well"   # use _apply_gravity_well_ to demonstrate a "well" shape.

ORBIT = False      # "True" gives planets an init vector that results in roughly 
                   # spherical orbits. This applies only in "well" mode.

# In "sink" mode the fall_off_type, radius, and amplitude of the LinearSinkForce
# can be set below (No values I tried would give the sink force a well shape): 

LSF_fo_type = 2             # 0 = by r, 1 = by r**2, 2 = by r**3
LSF_radius = 10000          # int - radius of field from its origin
LSF_amplitude = 1000000000  # int - strength of field


# Planet spawn positions, colours, and sizes.
sys_init_list = [{'init_pos':(-100,0,0),'color':(.7,0,0,1),'scale':6},
                 {'init_pos':(0,200,0),'color':(0,.7,0,1),'scale':3},
                 {'init_pos':(300,-150,0),'color':(0,0,.7,1),'scale':8}]

# Imports.
import direct.directbase.DirectStart
from panda3d.core import GeomVertexFormat, GeomVertexData, GeomVertexWriter
from panda3d.core import GeomTriangles, Geom, GeomNode, NodePath, PandaNode
from panda3d.core import AmbientLight, PointLight, Material, ClockObject, Vec3
from direct.actor.Actor import ActorNode, ForceNode, LinearSinkForce
from math import sqrt, cos, sin, degrees

# Basic sphere object to act as planet.
class Ico_Sphere:

    def __init__(self, sphere_init):
        pts, tris = self.__generate_Points(sphere_init)
        actor_node = ActorNode("sphere_actor")
        self.ACTOR_NP = NodePath(actor_node)
        self.ACTOR_NP.setPos(sphere_init['init_pos'])
        self.ACTOR_NP.setZ(0)  # Sim does not handle z axis.
        self.NP = self.__build_Sphere(pts, tris)
        self.NP.reparentTo(self.ACTOR_NP)
        self.NP.setScale(sphere_init['scale'])
        self.c_vec = Vec3(0, 0, 0)
        
        if ORBIT:
            dist = self.NP.getDistance(render)
            if dist:
                init_vec = self.__calc_Init_Vec(dist, sphere_init)
                self.c_vec = init_vec
        
        if G_MODE == "sink":
            base.physicsMgr.attachPhysicalNode(actor_node)

    def __generate_Points(self, sphere_init):
        # Set positions of 12 initial points.
        s = 0.5257311121191336
        t = s * 1.6180339887498948
        ipts = [[0, s, t], [0, -s, t], [t, 0, s], [s, -t, 0],
                [t, 0, -s], [0, -s, -t], [s, t, 0], [0, s, -t],
                [-s, t, 0], [-t, 0, -s], [-t, 0, s], [-s, -t, 0]]
        
        # Intitialize attrs of 12 points.
        pts = []
        col = sphere_init['color']
        for p in ipts:
            pt = {}
            pt['pos'] = [p[0], p[1], p[2]]
            pt['color'] = col
            pts.append(pt)
        
        # Map pts into list of 20 triangles that make of panels of sphere.
        tris = [[0, 1, 2], [2, 1, 3], [2, 3, 4], [4, 3, 5],
                [0, 2, 6], [6, 2, 4], [6, 4, 7], [7, 4, 5],
                [0, 6, 8], [8, 6, 7], [8, 7, 9], [9, 7, 5],
                [0, 8, 10], [10, 8, 9], [10, 9, 11], [11, 9, 5],
                [0, 10, 1], [1, 10, 11], [1, 11, 3], [3, 11, 5]]

        return pts, tris

    def __build_Sphere(self, pts, tris):
        sphere_np = NodePath(PandaNode("sphere"))

        # Initialize vertex writers.
        format =GeomVertexFormat.getV3n3c4()
        vdata = GeomVertexData("Data", format, Geom.UHStatic) 
        vertices = GeomVertexWriter(vdata, "vertex")
        colours = GeomVertexWriter(vdata, "color")
        normals = GeomVertexWriter(vdata, "normal")

        # Pass pts to Vertex Writers.
        for pt in pts:
            vertices.addData3f(*pt['pos'])
            colours.addData4f(*pt['color'])
            normals.addData3f(*pt['pos'])

        # Build Primitive.
        triangles = GeomTriangles(Geom.UHStatic)
        for tri in tris:
            triangles.addVertices(*tri)
        triangles.closePrimitive()
        
        # Create Geom and node.
        geom = Geom(vdata)
        geom.addPrimitive(triangles)
        node = GeomNode("sphere")
        node.addGeom(geom)
        node_path = NodePath(node)
        node_path.reparentTo(sphere_np)
        return sphere_np

    def __calc_Init_Vec(self, dist, sphere_init):
        vel = sqrt(2/dist-1/dist) * 4
        base_vec = Vec3(1, 0, 0)
        radial_vec = Vec3(*sphere_init['init_pos'])
        radial_vec.normalize()
        angle = radial_vec.angleRad(base_vec)
        if sphere_init['init_pos'][1] < 0:
            angle = -angle
        x = vel * sin(angle)
        y = vel * -cos(angle)
        return -Vec3(x, y, 0)

# Solar System.
class Solar_System:
    
    def __init__(self, system_init):
        self._prev_dt = 0
        planets = self.__build_System(system_init)
        if G_MODE == "well":
            taskMgr.add(self._apply_gravity_well_, "alt_gravity", extraArgs=[planets], appendTask=True)
        
    def _apply_gravity_well_(self, planets, task, clock=ClockObject()):
        new_dt = clock.getRealTime()
        dt = new_dt - self._prev_dt
        for planet in planets:
            dist = planet.NP.getDistance(render)
            pos = planet.NP.getPos(render)
            planet.c_vec -= pos / (dist**3) * dt * 1000
            pos += planet.c_vec
            planet.NP.setPos(render, pos)
        self._prev_dt = new_dt
        return task.cont
        
    def __build_System(self, sys_init_list):
        base.enableParticles()
        base.setBackgroundColor(0, 0, 0)
        
        # Ambient Light.
        ambient = AmbientLight("ambient")
        ambient.setColor((0.2, 0.2, 0.2, 1))
        ambient_np = render.attachNewNode(ambient)
        render.setLight(ambient_np)
        
        # Add star in center of gravity field.
        star = Ico_Sphere({'init_pos':(0, 0, 0), 'color':(1, 1, 1, 1), 'scale':10})
        star.ACTOR_NP.reparentTo(render)
        
        # Point light for sun.
        sun_mat = Material()
        sun_mat.setEmission((100, 100, 100, 1))
        star.NP.setMaterial(sun_mat)
        p_light = PointLight("sun")
        p_light.setColor((1, 1, 1, 1))
        p_light_np = NodePath(p_light)
        p_light_np.reparentTo(star.NP)
        
        # Init LinearSinkForce if chosen.
        if G_MODE == "sink":
            force_node = ForceNode("gravity")
            gravity = LinearSinkForce((0, 0, 0), LSF_fo_type, LSF_radius, LSF_amplitude)
            force_node.addForce(gravity)
            base.physicsMgr.addLinearForce(gravity)
            star.NP.attachNewNode(force_node)
        
        # Add "planets".
        planets = []
        for p_init in sys_init_list:
            planet = Ico_Sphere(p_init)
            planet.ACTOR_NP.reparentTo(render)
            planet.NP.setLight(p_light_np)
            planets.append(planet)
            
        # Set camera to look down from above.
        base.camera.setPos(0, 0, 2000)
        base.camera.lookAt(star.NP)
        return planets

# Start.
import sys
base.accept("escape", sys.exit)
base.disableMouse()
Solar_System(sys_init_list)
run()

Ok, I can see the problem now, within the source code. The falloff (e.g. one-over-r) uses the user-defined value for radius, and not the distance between the force center and the the physics object.
See:

LVector3 LinearSinkForce::get_child_vector(const PhysicsObject *po)
INLINE PN_stdfloat LinearDistanceForce::get_scalar_term() const

In my eyes this renders LinearSinkForce and LinearSourceForce pretty useless, since setting falloff and radius is only a complicated way of setting a fixed value scalar factor.

I could try to fix it, but since this code has been written way before I joined the community and the aothor of these classes is unknown to me I would like to have feeback first. I don’t want to break existing games working with these classes (LinearDistanceForce, LinearSinkforce, LinearSourceForce).

I’m not sure how the dev process works, but you could possibly add 2 new derivatives of LinearDistanceForce, something like LinearGravityForce and LinearRepulsionForce, which would implement spherical gravity (and anti-gravity) correctly. This wouldn’t break old code and would give users the option to still use the sink or source forces if that’s the effect they’re going for.

Is this going to be fixed in a new release anytime soon? I was planning to do a simple game using the panda physics engine, but I need actual 1/r^2 gravity wells. Or is there a simple workaround to make these forces behave correctly?

The demo posted above includes a sample of how to implement a gravity well not using built in physics. The relevant code is in the method “apply_gravity_well” of the Solar_System object; it’s called by a task that runs once every frame (it also assumes that render(0,0,0) is the center of the solar system.)

I do appreciate the demo code, I already modified it to use 2 different sources and the interaction seems just fine. Won’t this have a lot more work to tie into the collision detection though? That was the main goal of continuing to use the built in physics over your gravity well function.
edit: Got it, thanks for the quick response!

Collisions are separate from physics in Panda3d so your choice of implementation for your gravity well won’t have any effect on how you implement collisions. The amount of work involved in implementing collisions depends on what you want to have happen for collision events: objects collide and bounce off each other or objects collide and explode, etc.

Were you able to solve the problem? It seems like big issue…

Hi Craig, see the code in the “apply_gravity_well” method of the Solar_System object in the original demo above for an example of the work around I used.

Were you able to find the right code for calculating gravity? I do need the code so can you provide it? Waiting for reply thanks in advance:)

Hi Issac, in the demo posted above you can find the following method which shows my gravity solution:

# method called by panda 'task'.
def _apply_gravity_well_(self, planets, task, clock=ClockObject()):  
        new_dt = clock.getRealTime()
        dt = new_dt - self._prev_dt
        for planet in planets:
            dist = planet.NP.getDistance(render)
            pos = planet.NP.getPos(render)
            planet.c_vec -= pos / (dist**3) * dt * 1000
            pos += planet.c_vec
            planet.NP.setPos(render, pos)
        self._prev_dt = new_dt
        return task.cont

If you want to see it working just copy and paste the demo into a file and run it (make sure you set the variable ‘G_MODE’ to “well” and not “sink”). This code is just an approximation of gravity and will not produce exact results (small errors creep in over time) so you can’t use it for a completely real simulator though it’s fine for games.