Gravity problem solver in Factor

Question

After answering a question about Factor, I decided that it looked interesting. To get the hang of the syntax and such, I decided to write some basic physics code to "solve" the three-body problem. Then, I realized that math is harder than I thought it would be, so instead I just made a vocabulary to solve some basic problems about classical gravitational attraction.

I'm looking for advice on:

• Idiomatic-ness: The functions feel pretty long for Factor, but they seem about as short as they'll get.
• Naming: Some of my function names are subpar, but I can't find better names for them.
• Separation: I feel like the vector arithmetic should be in its own vocabulary, but I'm not sure. Also, I have no idea how to do that.
• Removing rounding errors: The program is fairly accurate, but some rounding error leads to ever-so-slightly wrong numbers (i.e. \49ドル.99999999999999\$ instead of \50ドル\$)

USING: locals math.functions math.vectors math.trig math.constants ;
IN: q.physics.gravitational-attraction
TUPLE: body
 { mass initial: 0 }
 { position initial: { 0 0 } }
 { velocity initial: { 0 0 } }
;
:: array>body ( array -- obj )
 body new
 array first >>mass
 array second 2 head >>position
 array third 2 head >>velocity
;
: >body ( mass position velocity -- obj )
 3array array>body
;
! Constants
SYMBOLS: G ;
! Throws an error if it's 6.674e−11 so let's use 1
1 G set
:: vector-between ( from to -- v )
 to first from first -
 to second from second -
 2array
;
:: dot ( u v -- rad )
 u first v first * u second v second * +
 u first sq u second sq + sqrt v first sq v second sq + sqrt *
 /
;
: angle-between ( u v -- rad ) dot ;
:: v+ ( u v -- sum )
 u first v first +
 u second v second +
 2array
;
:: s*v ( k v -- scaled )
 v first k *
 v second k *
 2array
;
:: -1*v ( v -- neg )
 v first -1 *
 v second -1 *
 2array 
;
:: magnitude ( v -- num )
 v first sq v second sq + sqrt
;
:: unit ( v -- unit )
 v magnitude :> len
 v first len /
 v second len /
 2array
;
:: force-on ( obj from -- force )
 obj position>> from position>> vector-between :> vec
 G get obj mass>> * from mass>> *
 vec magnitude
 /
 vec unit s*v
;
:: forces ( bodies -- forces )
 bodies [
 :> body
 body bodies remove [
 body swap force-on
 ] map
 ] map
;

I don't know the math well enough to test the edge cases, but I did check with this code:

SYMBOLS: b1 b2 b3 ;
10 { 0 0 } { 0 0 } >body b1 set
10 { 1 0 } { 0 0 } >body b2 set
10 { 0 1 } { 0 0 } >body b3 set
b1 get b2 get b3 get 3array forces [| f | f first f second v+ ] map .

which gave me the result:

{
 { 100.0 100.0 }
 { -150.0 49.99999999999999 }
 { 49.99999999999999 -150.0 }
}

And according to my math, that's correct. Note that the results are in terms of \$G\,ドル not multiplying by it, since it gives me an error when I try to put in 6.674e−11.

Answer 1

I'm not quite skilled enough at Factor yet to review your code in full, but I would just like to note as a response to

I don't know the math well enough to test the edge cases, but I did check with this code

Factor has really slick builtin unit-testing. Here are the docs; I'll pull an example from my factor/work directory:

USING: tools.test mk-initialism ;
IN: mk-initialism.tests
{ "LASER" } [ "Light Ablah by S of Ed Radiof" >initialism ] unit-test
{ "USA" } [ "United S OF Americaof" >initialism ] unit-test
{ "USA" } [ "united states and america" >initialism ] unit-test
{ "JTW" } [ "Jordan Of the World" >initialism ] unit-test 

To create the test file, open the listener and do:

USE: tools.scaffold
"q.physics.gravitational-attraction" scaffold-tests

To run the tests, just run "q.physics.gravitational-attraction" test.