# Gravity problem solver in Factor

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.

• 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
;
: >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.

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.