N-body with Yampa FRP, haskell

Question

I am trying to make an n-body solver (a bunch of objects gravitationally attracted to each-other). The problem is that it looks like the gravity1 function does not feed back the return objects, resulting in linear motion of the objects:

The code looks like this:

    updateGame :: Game -> SF AppInput Game
    updateGame game = 
      proc input -> do
        ...
        objs        <- updateObjects' $ _foreground (Game._objects game) -< ()
    
        ...
    updateObjects' :: [Object] -> SF () [Object]
    updateObjects' objs =  parB . fmap (updateObject1 objs ) $ objs
    
    updateObject1 :: [Object] -> Object -> SF () Object
    updateObject1 objs0 obj0 =
      proc () -> do
        obj  <- gravity1 (objs0, obj0) -< ()
        returnA -< obj
    
    g = 6.673**(-11.0) :: Double
    
    gravity1 :: ([Object], Object) -> SF () (Object)    
    gravity1 (objs0, obj0) =
      proc () -> do
        let
          m0     =  _mass obj0               :: Double
          xform0 = (head . _transforms) obj0 :: M44 Double
          p0     = ( view (_w._xyz)) xform0  :: V3 Double
          
          ms     = foldr1 (+) $ fmap (_mass) objs0              :: Double
          xforms = fmap (head . _transforms) objs0              :: [M44 Double]
          ps     = foldr1 (^+^) $ fmap ( view (_w._xyz)) xforms :: V3 Double
    
          dir  = ps ^-^ p0                 :: V3 Double
          dist = norm dir                  :: Double
          f    = g * m0 * ms / dist**2.0   :: Double
          acc  = (f / ms) *^ (dir ^/ dist) :: V3 Double
          s    = 1000000000000000.0
          
        --vel <- ((_velocity obj0) ^+^) ^<< integral -< (s *^ acc)
        vel <- ((_velocity obj0) ^+^) ^<< integral -< (s *^ (DT.trace ("acc :" ++ show (s *^ acc)) $ acc))
    
        let mtx =
              mkTransformationMat
              rot
              tr
              where
                rot = (view _m33 xform0)
                tr  = vel + p0
          
        returnA -< obj0 { _transforms = [mtx]
                        , _velocity   = vel }

Running the code, there's the output I read in the console:

acc :V3 12105.49700148636 12105.49700148636 0.0
acc :V3 NaN NaN NaN
acc :V3 (-12105.49700148636) 12105.49700148636 0.0
acc :V3 12105.49700148636 12105.49700148636 0.0
acc :V3 NaN NaN NaN
acc :V3 (-12105.49700148636) 12105.49700148636 0.0
acc :V3 12105.49700148636 12105.49700148636 0.0
acc :V3 NaN NaN NaN
acc :V3 (-12105.49700148636) 12105.49700148636 0.0
acc :V3 12105.49700148636 12105.49700148636 0.0
acc :V3 NaN NaN NaN
acc :V3 (-12105.49700148636) 12105.49700148636 0.0
acc :V3 12105.49700148636 12105.49700148636 0.0
acc :V3 NaN NaN NaN
acc :V3 (-12105.49700148636) 12105.49700148636 0.0
acc :V3 12105.49700148636 12105.49700148636 0.0
acc :V3 NaN NaN NaN
acc :V3 (-12105.49700148636) 12105.49700148636 0.0
acc :V3 12105.49700148636 12105.49700148636 0.0
acc :V3 NaN NaN NaN

Basically the values are all the same (the NaN's are due to object computed against itself, which I should fix, but that's not the issue here), it looks like the gravity1 function does not feed back the return object, despite return value being:

    returnA -< obj0 { _transforms = [mtx]
                    , _velocity   = vel }

The result is the linear motion of the objects, since acc seems to be a constant.

My expectation here is that aftergravity1 :: ([Object], Object) -> SF () (Object) has updated and returned the object, the updateObjects' objs = parB . fmap (updateObject1 objs ) $ objs and updateObject1 objs0 obj0 = ... returnA -< obj should result in all objects being updated and next iteration cycle should feed the gravity1 :: ([Object], Object) -> SF () (Object) with an updated set of objects, so that the value of acc is different every frame...

Am I misunderstanding the logic of how things are supposed to work here?

Answer 1

Let's consider the type of gravity1:

gravity1 :: ([Object], Object) -> SF () (Object)

This is a function that, when provided with a list of Object and a particular Object will produce a signal function. This signal function is a so-called generator, meaning that it can be provided with empty input and continually produce a stream of Object output.

Thus, in essence, gravity1 as two different kinds of inputs:

• It has static input ([Object], Object), which it gets once and is forever unchanging, and
• It has streaming input () (which is basically equivalent to no streaming input).

With this in mind, it does make sense that, once supplied with its static arguments, gravity1 will produce a constant stream of Object — after all, it is never receiving updated data about where any Objects are!

In order to make the output stream dynamically react to changes in the positions of the objects, those positions need to be streaming, not static. In particular, the [Object] input should be streaming:

gravity1 :: Object -> SF [Object] Object

The other input argument, which specifies the starting position of the Object we care about, probably doesn't need to be streaming, but it does need to be dealt with carefully (it is only a starting position after all).

But if gravity1 takes the positions of the Objects as a streaming argument, how will you run it from updateObject? You will likely need to use some form of delay, such as:

updateObjects :: [Object] -> SF () [Object]
updateObjects objs0 = proc () -> do
  rec objs  <- iPre objs0 -< objs'
      objs' <- parB (fmap gravity1 objs0) -< objs
  returnA -< objs'

Incidentally, this strategy of using delay (i.e., the rec keyword along with iPre or something similar) is exactly what you'll need to use in gravity1 also to keep track of the current position of the particular object that you're calculating the gravity for.