I am only thinking of this in terms of physics because the code reproduces Newtonian gravity (e.g. D = 3.0), where the acceleration is relatively small. In terms of geometry, I'm trying everything that I can think of to hammer out the widening of the bands. The bands, if not widening, are obviously the result of discrete steps, which is totally expected. However, the gradient reduces as distance increases, at an accelerating rate (the bands are curves in general). For fun, I set epsilon to 1 and it makes no difference.

const long long signed int collision_count_plus = get_intersecting_line_count(receiver_pos_plus, 1.0, D, true);

const long long signed int collision_count = get_intersecting_line_count(receiver_pos, 1.0, D, true);

vector_3 gradient;
gradient.x = static_cast<MyBig>(collision_count_plus - collision_count) / epsilon;

I don't quite understand what your two images are showing. Can you please elaborate?

taby said:
I don't quite understand what your two images are showing. Can you please elaborate?

Say we have 3 numbers: 2,3,9. And we want to know the maximum, but we're not allowed to do comparisions (becasue they are discontinuous).

One way is to sum up powers of the numbers, say the 10th: 1024 + 59049 + 3486784401 = 3486844474.
The average is 3486844474 / 3 = 1162281491
Then we convert back using pow (1162281491, 1/10) = 8.06

So we are pretty close to the expected value of 9.
Idk how this kind of math is called, but sadly it does not work if we want the minimum. We could try a negative exponent, and it's not totally bad, but it gives those ripples, bumps, or bands you see in my plot.

The problem is that the numbers are close to the origin of zero, causing oscillations from minor differences, i guess.
So i could try to use another origin.
Or i could use log and exp instead pow, and see how that goes…

taby said:
In terms of geometry,

I still miss them.

Did not look at the code in detail, but i assume you do this:

Create an emitter of gravity, e.g. a galaxy. You fit an ellipsoid to the mass distribution of the galaxy, so it would be very flat, almost like a disc. And you call that a dimension of 2.01.

Create a receiver for the gravity. A simple sphere.
To calculate how many gravitons are received, you generate rays from the emitter this way: Create uniform samples on the unit sphere, apply the non uniform scaling to the rays which turns the unit sphere to the emitter ellipsoid, and shoot the stretched rays.
Then you count how many rays hit the receiver sphere.

Probably i'm wrong, but that's the kind of 'geometrical explanation' i would need. ; )

Edit: Just found out the solution to my problem is to use the inverse of the numbers. \:D/

You’re right about things. Let’s focus on D = 3. Then I am getting a pseudorandom unit vector, and seeing if that vector intersects with the receiver sphere. I get the gradient of this intersection count and multiply it by r^3 G M, which is the math behind Newtonian gravitation.

I’m so glad to hear that you made a bit of progress with regard to your problem. Sometimes it just takes a teacher-student Q&A session to clear things up. 🙂

Here is a picture which shows the banding. It shows the vertices (red crosses), and also the lines that connect the vertices (blue lines). n = 100,000,000.

Thanks for all of your input. The paper will be called ‘Numerical Newtonian gravitation: the simplest model’. I won't publish until I figure out the banding curves, which is the last thing that I need to figure out.

I'm sorry for the delay in my response.

As for the central object being a ellipsoid, this is not quite what I had in mind.

In the case of the Milky Way, the central bulge is spherical, and the anisotropic gravitation kicks in when it comes to the disk. The further from the bulge radius, the more disk-like the Galaxy becomes. The ellipsoidal shape that I use in my code is only one possible shape… the emitter can still be spherical and emit gravitons non-spherically: as long as the pressure density is next to zero.

Sorry if this does not answer your question:

#include <iostream>
#include <cmath>
using namespace std;

int main(void)
{
	const double c = 299792458; // Speed of light in vacuum
	const double G = 6.674e-11; // Gravitational constant
	const double M = 1e41; // Galactic bulge mass

	const double start_distance = 6.1495e19; // Galactic bulge radius
	const double end_distance = 3.086e20; // Solar orbit radius (about 10 kiloparsecs)

	double v = sqrt(G * M / start_distance); // Speed with dark matter

	const size_t resolution = 1000;	
	const double step_size = (end_distance - start_distance) / (resolution - 1);

	for (double r = start_distance; r <= end_distance; r += step_size)
	{
		const double v_N = sqrt(G * M / r); // Speed without dark matter

		if (v < v_N)
			v = v_N;

		const double D = 3.0 - log(v / v_N) / log(c); // Dimension

		//const double G_a = v * v * r / M;
		//cout << r << " " << G_a << endl;

		cout << r << " " << D << endl;
	}

	return 0;
}

JoeJ said:
Edit: Just found out the solution to my problem is to use the inverse of the numbers. \:D/

Hehe, no. Exact same results. I've learned: pow(x,-4) = pow(1/x,4). I did not know before. :D

But got it working in the end, it seems. Weeks of SM research are over, and i only need to optimize implementation.
I don't trust it yet, though. Soon i may discover new forms of unacceptable artifacts, putting me again back to the start…

taby said:
The ellipsoidal shape that I use in my code is only one possible shape…

But if it's good enough for a proof, you could replace the ray sampling with something more accurate: Calculate the intersection of the ellipsoid with the cone bounding the receiver sphere, then divide the intersection area by the total area of the ellipsoid:

But if performance is no problem, it surely isn't worth the work.

taby said:
Sorry if this does not answer your question:

Well, i'm no physics guy, so ofc. not. ; )

But there is this branch (again): if (v < v_N) v = v_N;

If i see this, i immideately think this can't be physics. God does not branch. He floats. He never decides, he just follows the flow of his own. ;D

This is god, suffering from the anisotropic gravity from one of his galaxies.

OK, so I've got the code where D = 3.0 by definition. It's quite small:

https://github.com/sjhalayka/numerical_newtonian_gravity/blob/main/main.cpp

#include "main.h"

int main(int argc, char** argv)
{
	//cout << setprecision(10);

	const size_t n = 100000000; // Oscillator count

	cout << "Allocating memory for oscillators" << endl;
	unit_vectors.resize(n);

	for (size_t i = 0; i < n; i++)
	{
		unit_vectors[i] = RandomUnitVector();

		static const size_t output_mod = 10000;

		if (i % output_mod == 0)
			cout << "Getting pseudorandom locations: " << i << " of " << n << endl;
	}

	string filename = "3.0.txt";

	ofstream out_file(filename.c_str());
	out_file << setprecision(30);

	const real_type start_distance = 10;
	const real_type end_distance = 100;

	const size_t distance_res = 1000;

	const real_type distance_step_size = (end_distance - start_distance) / (distance_res - 1);

	for (real_type r = start_distance; r <= end_distance; r += distance_step_size)
	{
		const vector_3 receiver_pos(r, 0, 0);
		const real_type receiver_radius = 1.0;
		const real_type epsilon = 1.0;

		// https://en.wikipedia.org/wiki/Directional_derivative
		const vector_3 receiver_pos_plus(receiver_pos.x + epsilon, 0, 0);
		const long long signed int collision_count_plus = get_intersecting_line_count(receiver_pos_plus, receiver_radius);
		const long long signed int collision_count = get_intersecting_line_count(receiver_pos, receiver_radius);
		const vector_3 gradient(static_cast<real_type>(collision_count_plus - collision_count) / epsilon, 0, 0);
		const real_type gradient_strength = gradient.length() * pow(receiver_pos.x, 3.0);

		cout << " r: " << r << " " << gradient_strength << endl;

		out_file << r << " " << gradient_strength << endl;
	}

	out_file.close();

	return 0;
}

Well the simplest way would be to find a pseudorandom location on the circle formed by the receiver sphere’s cross section, and then trace that back to the emitter sphere’s centre. However, that’s not necessary. The hit ratio is just circle steradians /4pi, but we threw the baby (quantization) out with the bath water. 🙂

Thanks again so much for all of your insight and questions. I feel that I am getting closer to understanding the banding.

P.S. how does one go about calculating the solid angle?

This is what I'm trying:

#include "main.h"


int main(int argc, char** argv)
{
	//cout << setprecision(10);

	const size_t n = 100000000; // Oscillator count

	string filename = "3.0.txt";

	ofstream out_file(filename.c_str());
	out_file << setprecision(30);

	const real_type start_distance = 10;
	const real_type end_distance = 100;

	const size_t distance_res = 1000;

	const real_type distance_step_size = (end_distance - start_distance) / (distance_res - 1);

	for (real_type r = start_distance; r <= end_distance; r += distance_step_size)
	{
		const vector_3 receiver_pos(r, 0, 0);
		const real_type receiver_radius = 1.0;
		const real_type epsilon = 1.0;

		// https://en.wikipedia.org/wiki/Directional_derivative
		const vector_3 receiver_pos_plus(receiver_pos.x + epsilon, 0, 0);
		const long long signed int collision_count_plus = get_intersecting_line_count(n, receiver_pos_plus, receiver_radius);
		const long long signed int collision_count = get_intersecting_line_count(n, receiver_pos, receiver_radius);
		const vector_3 gradient(static_cast<real_type>(collision_count_plus - collision_count) / epsilon, 0, 0);
		const real_type gradient_strength = gradient.length() * pow(receiver_pos.length(), 3.0);

		cout << "r: " << r << " " << gradient_strength << endl;

		out_file << r << " " << gradient_strength << endl;
	}

	out_file.close();

	return 0;
}

where:

size_t get_intersecting_line_count(
	size_t n,
	const vector_3 sphere_location,
	const real_type sphere_radius)
{
	double spherical_cap_area = 2.0 * pi * sphere_radius * sphere_location.x;
	real_type sphere_area = 4.0 * pi * sphere_location.x * sphere_location.x;

	real_type ratio = spherical_cap_area / sphere_area;

	return static_cast<size_t>(n * ratio);
}

The result is:

Nice fractal pattern there. But, what I was expecting was a flat line.