diff --git a/Physics/MovingParticles.mch b/Physics/MovingParticles.mch
index a0aa6feaebaf1883c2601a367aa81c2af6bffaaf..60cfdf00dc9dcf3feed2b36bef8723a0a275cde6 100644
--- a/Physics/MovingParticles.mch
+++ b/Physics/MovingParticles.mch
@@ -1,32 +1,36 @@
 MACHINE MovingParticles
+INCLUDES Vectors
 DEFINITIONS
+  "LibraryReals.def";
   MASS == FLOAT;
   POS == FLOAT;
   SPEED == FLOAT;
-  OBJ == 1..n
-ABSTRACT_CONSTANTS sum_vec
-CONSTANTS mass, n
+  OBJ == 1..nobjs;
+  G == 6.67e-11;
+  SCALE == 1.25e11; // 100 % in visualisation
+  SCALE100 == (SCALE/100.0);
+  SPEED_MULTIPLIER == 1000000.0;
+  MAXMASS == 1.989e30
+
+CONSTANTS mass, nobjs
 PROPERTIES
-  n : 1..10 &
+  nobjs : 1..10 &
   mass : OBJ --> MASS &
 
- sum_vec = %(v1,v2).(v1:OBJ-->REAL & v2:OBJ-->REAL | %i.(i:OBJ| v1(i)+v2(i)))
-
-  & n = 3 & mass = [1.0, 100.0, 2.0]
-VARIABLES xpos, ypos, xspeed, yspeed
+  n = 2 & // two dimensional space
+  nobjs = 3 & mass = [5.97e24, 1.989e30, 1.989e30]
+VARIABLES pos, speed
 INVARIANT
- xpos : OBJ --> POS &
- ypos : OBJ --> POS &
- xspeed : OBJ --> SPEED &
- yspeed : OBJ --> SPEED
-INITIALISATION xpos :=   [ -1.0, 0.0 , 1.0 ] ||
-               ypos :=   [ 0.0,  0.0,  0.0 ] ||
-               xspeed := [ 0.0,  0.0,  0.0 ] ||
-               yspeed := [ 1.0,  0.0, -1.0 ]
+  pos : OBJ --> VEC &
+  speed : OBJ --> VEC
+INITIALISATION pos :=   [ [0.0,0.0]  ,  [0.0,4.5e10] ,   [0.0,-4.5e10 ] ] ||
+               speed := [ [0.05e04,0.0],  [3.0e04,0.0],  [-3.0e04,0.0]  ] 
 OPERATIONS
-  Move = BEGIN
-           xpos := sum_vec(xpos,xspeed) ||
-           ypos := sum_vec(ypos,yspeed)
+  r <-- ForceBetween(a,b) = PRE a:OBJ & b:OBJ & a/= b THEN
+            r := (G * mass(a) * mass(b)) / (distance(pos(a),pos(b))**2)
+         END;
+  Move(delta) = PRE delta:REAL & delta = 100000.0 THEN
+            pos := %i.(i:OBJ | add_vec(pos(i),times_vec(delta,speed(i))) )
          END
 END
 
diff --git a/Physics/Vectors.mch b/Physics/Vectors.mch
new file mode 100644
index 0000000000000000000000000000000000000000..53807d4d47a512067f6930351261b8bb7e755538
--- /dev/null
+++ b/Physics/Vectors.mch
@@ -0,0 +1,39 @@
+MACHINE Vectors
+DEFINITIONS
+  "LibraryReals.def"
+ABSTRACT_CONSTANTS
+  n,
+  VEC, DIM,
+  add_vec,
+  delta_vec,
+  dot_product,
+  magnitude,
+  unit_vec,
+  times_vec,
+  distance,
+  
+  add_first,
+  sigma_reals
+PROPERTIES
+  n:NATURAL1 &
+  DIM = (1..n) &
+  VEC = DIM --> REAL &
+  
+  add_first = %v.(v:seq1(REAL) | (v(1)+v(2)) -> tail(tail(v))) &
+  sigma_reals = %v.(v:VEC | (iterate(add_first,n-1)(v))(1)) &
+  
+  add_vec = %(v1,v2).(v1:VEC & v2:VEC | %i.(i:DIM| v1(i)+v2(i))) &
+  delta_vec = %(v1,v2).(v1:VEC & v2:VEC | %i.(i:DIM| v1(i)-v2(i))) &
+  // dot_product = %(v1,v2).(v1:VEC & v2:VEC | SIGMA(i).(i:DIM| v1(i)*v2(i))) &
+  dot_product = %(v1,v2).(v1:VEC & v2:VEC | sigma_reals(%i.(i:DIM| v1(i)*v2(i)))) &
+  
+  magnitude = %v1.(v1:VEC | RSQRT(dot_product(v1,v1))) &
+  
+  unit_vec = %v1.(v1:VEC | LET m BE m=magnitude(v1) IN
+                              %i.(i:DIM| v1(i)/m) END ) &
+                              
+  times_vec = %(k,v1).(k:REAL & v1:VEC | %i.(i:DIM| k*v1(i))) &
+                              
+  distance = %(v1,v2).(v1:VEC & v2:VEC | magnitude(delta_vec(v1,v2)))
+
+END
\ No newline at end of file
diff --git a/Physics/three_bodies.json b/Physics/three_bodies.json
index ee027a0083da301690087b30b37f1c0f5f4cf47f..221fce5415c857b3e02ee540ed1930049d5a2375 100644
--- a/Physics/three_bodies.json
+++ b/Physics/three_bodies.json
@@ -5,19 +5,45 @@
 	  "for": {"from":1, "to":3},
       "id": "body%0",
       "attr": "cx",
-      "value" : "10.0*xpos(%0)"
+      "value" : "pos(%0)(1) / SCALE100"
     },
 	{
 	  "for": {"from":1, "to":3},
       "id": "body%0",
       "attr": "cy",
-      "value" : "10.0*ypos(%0)"
+      "value" : "pos(%0)(2) / SCALE100"
     },
+    
 	{
 	  "for": {"from":1, "to":3},
       "id": "body%0",
       "attr": "r",
-      "value" : "2.0+mass(%0)/10.0"
+      "value" : "2.0+5.0*mass(%0)/MAXMASS"
+    },
+    
+	{
+	  "for": {"from":1, "to":3},
+      "id": "speed%0",
+      "attr": "x1",
+      "value" : "pos(%0)(1) / SCALE100"
+    },
+	{
+	  "for": {"from":1, "to":3},
+      "id": "speed%0",
+      "attr": "y1",
+      "value" : "pos(%0)(2) / SCALE100"
+    },
+	{
+	  "for": {"from":1, "to":3},
+      "id": "speed%0",
+      "attr": "x2",
+      "value" : "(pos(%0)(1)+ SPEED_MULTIPLIER*speed(%0)(1)) / SCALE100"
+    },
+	{
+	  "for": {"from":1, "to":3},
+      "id": "speed%0",
+      "attr": "y2",
+      "value" : "(pos(%0)(2)+ SPEED_MULTIPLIER*speed(%0)(2)) / SCALE100"
     }
   ],
   "events": [
diff --git a/Physics/three_bodies.svg b/Physics/three_bodies.svg
index 0127f1b114d53a982b2c11fbb9614cacea2d8cc7..3f84ed9b7e942a49f98c52b4eff4fe7979456eb3 100644
--- a/Physics/three_bodies.svg
+++ b/Physics/three_bodies.svg
@@ -5,13 +5,46 @@
    viewBox="-100 -100 200 200"
    >
   
-    <circle id="body1" cx="45" cy="15" r="25"
-    stroke="black" stroke-width="1" fill="green" />
+   
+   <defs>
+    <marker id="arrowhead" markerWidth="10" markerHeight="7" 
+    refX="0" refY="3.5" orient="auto">
+      <polygon points="0 0, 10 3.5, 0 7" />
+    </marker>
+    
+    <marker id="startarrow" markerWidth="10" markerHeight="7" 
+    refX="10" refY="3.5" orient="auto">
+      <polygon points="10 0, 10 7, 0 3.5" fill="blue" />
+    </marker>
+    
+    <marker id="endarrow" markerWidth="10" markerHeight="7" 
+    refX="0" refY="3.5" orient="auto">
+      <polygon points="0 0, 10 3.5, 0 7"  fill="blue" />
+    </marker>
+  </defs>
+  
+  
+    <circle id="body1" cx="45" cy="15" r="10"
+    stroke="black" stroke-width="0.5" fill="green" />
     
     <circle id="body2" cx="5" cy="15" r="5"
-    stroke="black" stroke-width="1" fill="red" />
+    stroke="black" stroke-width="0.5" fill="red" />
+    
+    <circle id="body3" cx="35" cy="45" r="5"
+    stroke="black" stroke-width="0.5" fill="blue" />
+    
     
-    <circle id="body3" cx="35" cy="45" r="20"
-    stroke="black" stroke-width="1" fill="blue" />
+  <line x1="45" y1="15" x2="55" y2="15" stroke="#000" 
+  stroke-width="0.3" marker-end="url(#arrowhead)"
+   id = "speed1"
+   />
+  <line x1="5" y1="15" x2="0" y2="0" stroke="#000" 
+  stroke-width="0.3" marker-end="url(#arrowhead)"
+   id = "speed2"
+   />
+  <line x1="35" y1="45" x2="35" y2="55" stroke="#000" 
+  stroke-width="0.3" marker-end="url(#arrowhead)"
+   id = "speed3"
+   />
     
 </svg>