Implement N-body simulation in Modelica

diff --git a/nbody.mo b/nbody.mo
new file mode 100644 (file)
index 0000000..ef1ef5e
--- /dev/null
+++ b/nbody.mo
@@ -0,0 +1,60 @@
+package nbody
+
+       type Mass = Real;
+       type Acceleration = Real[3];
+       type Speed = Real[3];
+       type RelPosition = Real[3];
+       type Position = Real[3];
+
+       model Planet
+               parameter Mass  mass = 1e9;
+               Acceleration    accel;
+               Speed                   speed;
+               Position                pos;
+
+       equation
+               accel = der(speed);
+               speed = der(pos);
+       end Planet;
+
+       function norm
+               input Real[:]   vec;
+               output Real             res;
+       algorithm
+               res := sqrt(vec * vec);
+       end norm;
+
+       function unitvec
+               input Real[:]   vec;
+               output Real[:]  res;
+       algorithm
+               res := vec / norm(vec);
+       end unitvec;
+
+       model System
+               parameter Integer       N = 4;
+               Planet                          p[N];
+               RelPosition                     pos[N, N];
+
+       initial equation
+               for i in 1:N loop
+                       p[i].pos[1] = (i / N) ^ 2;
+                       p[i].pos[2] = ((i - 1) / N) ^ 2;
+                       p[i].pos[3] = 0;
+                       p[i].speed[1] = 0;
+                       p[i].speed[2] = 0;
+                       p[i].speed[3] = 0;
+               end for;
+
+       equation
+               for i in 1:N loop
+                       for j in 1:N loop
+                               pos[i, j] = p[i].pos - p[j].pos;
+                       end for;
+                       p[i].accel = Modelica.Constants.G * sum(
+                               p[j].mass / (pos[i, j] * pos[i, j]) * (-unitvec(pos[i, j]))
+                               for j in vector([1:(i-1); (i+1):N]));
+               end for;
+       end System;
+
+end nbody;

diff --git a/script.mos b/script.mos
new file mode 100644 (file)
index 0000000..6be1b36
--- /dev/null
+++ b/script.mos
@@ -0,0 +1,3 @@
+loadFile("nbody.mo");
+instantiateModel(System);
+simulate(nbody.System, numberOfIntervals = 5000);