}
+// returns <1.37136364,1.37136734> for 1 mass at <0,0>
+/*
Vector<float> World::ForceAt(Vector<float> p1, float m1) const {
Vector<float> force(0, 0);
return force;
}
+*/
+
+// returns <1.37136674,1.37136662> for 1 mass at <0,0>
+// (should be less prone to rounding loss, i.e. more accurate)
+Vector<float> World::ForceAt(Vector<float> p1, float m1) const {
+ Vector<float> force(0, 0);
+
+ for (Vector<int> p2(0, 0); p2.y < size.y; ++p2.y) {
+ Vector<float> rowForce(0, 0);
+ for (p2.x = 0; p2.x < size.x; ++p2.x) {
+ const int i = Index(p2);
+ if (p1 == Vector<float>(p2)) continue;
+
+ const Vector<float> diff(Vector<float>(p2) - p1);
+ const Vector<float> dir(Norm(diff));
+
+ const float m2 = masses[i];
+ const float r2 = Dot(diff, diff);
+
+ const float mag = G * ((m1 * m2) / r2) * 2; // double because m2 is our reference frame
+
+ rowForce += dir * mag;
+ }
+ force += rowForce;
+ }
+
+ return force;
+}
}