diff --git a/Physics/MovingParticles.mch b/Physics/MovingParticles.mch
index 5bd288488156fc80ea0ed15441264a4254cc9e6e..f1ccb7cc1041a3fbce737127178d38ed650e15f7 100644
--- a/Physics/MovingParticles.mch
+++ b/Physics/MovingParticles.mch
@@ -1,11 +1,12 @@
MACHINE MovingParticles
+// Encoding of an n-body problem in B using the REAL datatype
INCLUDES Vectors
DEFINITIONS
"LibraryReals.def"; "SORT.def";
MASS == FLOAT;
- POS == FLOAT;
- SPEED == FLOAT;
+
OBJ == 1..nobjs;
+
G == 6.67e-11;
SCALE == 1.25e11; // 100 % in visualisation
SCALE100 == (SCALE/100.0);
@@ -20,30 +21,36 @@ DEFINITIONS
unit_vec(delta_vec(pos(j),pos(i)))))
) ) ) // acceleration
ABSTRACT_CONSTANTS
- force_between
-CONSTANTS mass, nobjs
+ force_between /*@desc "gravitational force between two objects" */
+CONSTANTS
+ mass /*@desc "mass of each object OBJ --> MASS" */,
+ nobjs
PROPERTIES
- nobjs : 1..10 &
- mass : OBJ --> MASS &
+ nobjs : 1..10 &
+ mass : OBJ --> MASS &
- force_between = %(a,b,pa,pb).(a:OBJ & b:OBJ & a/=b | (G * mass(a) * mass(b)) / (distance(pa,pb)**2))
+ force_between = %(a,b,pa,pb).(a:OBJ & b:OBJ & a/=b | (G * mass(a) * mass(b)) / (distance(pa,pb)**2))
- &
- n = 2 & // two dimensional space
- nobjs = 3 & mass = [5.97e24, 1.989e30, 1.989e30]
+ & // specific configuration:
+ n = 2 & // two dimensional space
+ nobjs = 3 &
+ mass = [5.97e24, 1.989e30, 1.989e30]
VARIABLES pos, speed
INVARIANT
- 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] ]
+ 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
- r <-- ForceBetween(a,b) = PRE a:OBJ & b:OBJ & a/= b THEN
- r := force_between(a,b,pos(a),pos(b))
+ r <-- ForceBetween(a,b) =
+ PRE a:OBJ & b:OBJ & a/= b & bfalse THEN
+ r := force_between(a,b,pos(a),pos(b))
END;
- Move(delta) = PRE delta:REAL & delta = 10000.0 THEN
- speed := %i.(i:OBJ | add_vec(speed(i),times_vec(delta,gravitation(i)) )) ||
- pos := %i.(i:OBJ | add_vec(pos(i),times_vec(delta,speed(i))) )
+ Move(deltaT) = // move time by deltaT time units
+ PRE deltaT:REAL & deltaT = 10000.0 THEN
+ speed := %i.(i:OBJ | add_vec(speed(i),times_vec(deltaT,gravitation(i)) )) ||
+ pos := %i.(i:OBJ | add_vec(pos(i), times_vec(deltaT,speed(i))) ) // we could use new speed
END
END
diff --git a/Physics/Vectors.mch b/Physics/Vectors.mch
index 369ca3d09a8da78b633f9e3ce0dca615e2d321ec..843ad50dd783ee3bb93e4c202415dfdf3fa5b71f 100644
--- a/Physics/Vectors.mch
+++ b/Physics/Vectors.mch
@@ -4,41 +4,45 @@ DEFINITIONS
ABSTRACT_CONSTANTS
n,
VEC, DIM,
- add_vec,
- delta_vec,
- dot_product,
- magnitude,
- unit_vec,
- times_vec,
- distance,
+ add_vec /*@desc "sum of two vectors" */,
+ delta_vec /*@desc "difference of two vectors" */,
+ dot_product /*@desc "dot product of two vectors" */,
+ magnitude /*@desc "Magnitude of a vector" */,
+ unit_vec /*@desc "unit vector of a vector" */,
+ times_vec /*@desc "multiply scalar by a vector" */,
+ distance /*@desc "Eucledian distance between two vectors" */,
add_first,
- sigma_reals,
+ sigma_reals /*@desc "sum of a numbers of a vector" */,
add_firstv,
- sigma_vecs
+ sigma_vecs /*@desc "sum of a sequence of vectors" */
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_first = %v.(v:seq1(REAL) | (v(1)+v(2)) -> tail(tail(v))) & // auxiliary function
+ sigma_reals = %v.(v:VEC | (iterate(add_first,n-1)(v))(1)) & // sum the values in a vector
- add_firstv = %v.(v:seq1(VEC) | add_vec(v(1),v(2)) -> tail(tail(v))) &
- sigma_vecs = %v.(v:seq1(VEC) | (iterate(add_firstv,size(v)-1)(v))(1)) &
+ add_firstv = %v.(v:seq1(VEC) | add_vec(v(1),v(2)) -> tail(tail(v))) & // auxiliary function
+ sigma_vecs = %v.(v:seq1(VEC) | (iterate(add_firstv,size(v)-1)(v))(1)) & // sum a sequence of vectors
- 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))) &
+ add_vec = %(v1,v2).(v1:VEC & v2:VEC | %i.(i:DIM| v1(i)+v2(i)))
+ /*@desc "sum of two vectors" */ &
+ delta_vec = %(v1,v2).(v1:VEC & v2:VEC | %i.(i:DIM| v1(i)-v2(i)))
+ /*@desc "difference of two vectors" */ &
// 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)))) &
+ dot_product = %(v1,v2).(v1:VEC & v2:VEC | sigma_reals(%i.(i:DIM| v1(i)*v2(i))))
+ /*@desc "dot product of two vectors" */ &
- magnitude = %v1.(v1:VEC | RSQRT(dot_product(v1,v1))) &
+ magnitude = %v1.(v1:VEC | RSQRT(dot_product(v1,v1))) /*@desc "Magnitude of a vector" */ &
unit_vec = %v1.(v1:VEC | LET m BE m=magnitude(v1) IN
- %i.(i:DIM| v1(i)/m) END ) &
+ %i.(i:DIM| v1(i)/m) END ) /*@desc "unit vector of a vector" */ &
- times_vec = %(k,v1).(k:REAL & v1:VEC | %i.(i:DIM| k*v1(i))) &
+ times_vec = %(k,v1).(k:REAL & v1:VEC | %i.(i:DIM| k*v1(i))) /*@desc "multiply scalar by a vector" */ &
distance = %(v1,v2).(v1:VEC & v2:VEC | magnitude(delta_vec(v1,v2)))
+ /*@desc "Eucledian distance between two vectors" */
END
\ No newline at end of file