OK, so I have it drawing a straight line. I'm closer:

#include "main.h"

real_type get_intersecting_line_count(
	size_t n,
	const vector_3 sphere_location,
	const real_type sphere_radius)
{
	real_type 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<real_type>(n) * ratio);
}

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 = 0.01;

		// https://en.wikipedia.org/wiki/Directional_derivative
		const vector_3 receiver_pos_plus(receiver_pos.x + epsilon, 0, 0);
		const real_type collision_count_plus = get_intersecting_line_count(n, receiver_pos_plus, receiver_radius);
		const real_type collision_count = get_intersecting_line_count(n, receiver_pos, receiver_radius);
		const vector_3 gradient((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;
}

I also did away with the sphere-ray intersection test. I use a dot product comparison now. It's not much faster, but it is simpler to follow.

long long unsigned int get_intersecting_line_count(
	const vector_3 sphere_location,
	const real_type sphere_radius)
{
	long long unsigned int count = 0;

	for (size_t i = 0; i < unit_vectors.size(); i++)
	{
		vector_3 veca = unit_vectors[i];

		vector_3 vecb(sphere_location.x, sphere_radius, 0);
		vecb.normalize();

		vector_3 vecc(sphere_location.x, 0, 0);
		vecc.normalize();

		real_type veca_dot = veca.dot(vecc);
		real_type vecb_dot = vecb.dot(vecc);

		if (veca_dot >= vecb_dot)
			count++;
	}

	return count;
}

I also added in repulsion, to space the unit vectors apart along the 2-sphere. It's glacially slow. I'm also not sure if it's needed where n is much larger than 1.

const size_t repulsion_rounds = 10;

for (size_t i = 0; i < repulsion_rounds; i++)
{
	cout << "Repulsing round " << i << endl;

	vector<vector_3> unit_vectors_backup = unit_vectors;

	for (size_t j = 0; j < n - 1; j++)
	{
		vector_3 accel;

		for (size_t k = j + 1; k < n; k++)
		{
			//if (j == k)
			//	continue;

			vector_3 grav_dir = unit_vectors_backup[k] - unit_vectors_backup[j];
				
			real_type distance = grav_dir.length();
			grav_dir.normalize();

			// Repulsion
			accel += -grav_dir / distance;
		}

		unit_vectors[j] += accel;
	}

	for (size_t j = 0; j < n; j++)
		unit_vectors[j].normalize();
}

OK, so the intersection count reduces with distance. The gradient of the intersection count also reduces with distance. It's important to remember that we're dealing with integers. Below is a plot of the gradient for the numerical simulation:

Below is a plot of 1.0/r^3, which is the same plot as above:

So, the multiplication by the distanced cubed, to get gradient strength, is no longer required. Along with the distance cubed, also goes the banding.

long long unsigned int get_intersecting_line_count(
	const vector<vector_3> &unit_vectors,
	const vector_3 sphere_location,
	const real_type sphere_radius)
{
	long long unsigned int count = 0;

	vector_3 cross_section_edge_dir(sphere_location.x, sphere_radius, 0);
	cross_section_edge_dir.normalize();

	vector_3 receiver_dir(sphere_location.x, 0, 0);
	receiver_dir.normalize();

	const real_type min_dot = cross_section_edge_dir.dot(receiver_dir);

	for (size_t i = 0; i < unit_vectors.size(); i++)
		if (unit_vectors[i].dot(receiver_dir) >= min_dot)
			count++;

	return count;
}


int main(int argc, char** argv)
{
	const size_t n = 10000000; // Oscillator count


	cout << "Allocating memory for oscillators" << endl;
	vector<vector_3> unit_vectors(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: " << static_cast<float>(i)/n << endl;
	}


	//const size_t repulsion_rounds = 10;

	//for (size_t i = 0; i < repulsion_rounds; i++)
	//{
	//	cout << "Repulsing round " << i << endl;

	//	vector<vector_3> unit_vectors_backup = unit_vectors;

	//	for (size_t j = 0; j < n - 1; j++)
	//	{
	//		vector_3 accel;

	//		for (size_t k = j + 1; k < n; k++)
	//		{
	//			//if (j == k)
	//			//	continue;

	//			vector_3 grav_dir = unit_vectors_backup[k] - unit_vectors_backup[j];
	//			
	//			real_type distance = grav_dir.length();
	//			grav_dir.normalize();

	//			// Repulsion
	//			accel += -grav_dir / distance;
	//		}

	//		unit_vectors[j] += accel;
	//	}

	//	for (size_t j = 0; j < n; j++)
	//		unit_vectors[j].normalize();
	//}


	string filename = "newton.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 (size_t step_index = 0; step_index < distance_res; step_index++)
	{
		const real_type r = start_distance + step_index * distance_step_size;

		const vector_3 receiver_pos(r, 0, 0);
		const real_type receiver_radius = 1.0;

		const real_type epsilon = 1.0;

		vector_3 receiver_pos_plus = receiver_pos;
		receiver_pos_plus.x += epsilon;

		// https://en.wikipedia.org/wiki/Directional_derivative
		const long long signed int collision_count_plus = get_intersecting_line_count(unit_vectors, receiver_pos_plus, receiver_radius);
		const long long signed int collision_count = get_intersecting_line_count(unit_vectors, receiver_pos, receiver_radius);
		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: " << gradient_strength << endl;

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

	out_file.close();

	return 0;
}

taby said:
Along with the distance cubed, also goes the banding.

So the bands are likely the result of precision loss from raising to some power (or other transcendental functions you use)?
That's what i've thought from the start, knowing how quickly the numbers become jaggy when doing so.
But then i wonder: Why does the high precision library not help against the banding?

Edit: Maybe the library does not implement some functions, but internally converts to double and back to use standard functions?

Btw, it sounds like you have just verified that you can calculate the ‘force’ of gravity between two bodies using projected area of the solid angle? Meaning, we could calculate gravity using the same methods we typically use to solve for radiosity?

If so, there seems an obvious mistake: Mass is not defined from volume alone, so the projection stuff would fail if the volume / mass ratio (density) of your objects changes, but their size remains the same?

As far as I can tell, it is not a precision loss issue, per se. Rather, it is because of the integers, and the bands are curved because this integer gradient strength decreases in a way that accelerates as distance increases.

What I'm doing is making field lines. There is only one situation where mass equals radius, and that's for a black hole, but that's kind of out of the realm of basic brute force Newtonian gravitation using field lines.

Please let me think about your valid objection. Thanks again for your input and questions!

JoeJ said:

Btw, it sounds like you have just verified that you can calculate the ‘force’ of gravity between two bodies using projected area of the solid angle? Meaning, we could calculate gravity using the same methods we typically use to solve for radiosity?

Well, you can definitely speed up the calculations using a cone like I did in a previous post:

real_type get_intersecting_line_count(
	const vector<vector_3>& unit_vectors,
	const vector_3 sphere_location,
	const real_type sphere_radius)
{
	real_type solid_angle = 2.0 * pi * (1 - sqrt(sphere_location.x * sphere_location.x - sphere_radius * sphere_radius) / sphere_location.x);

	//real_type rho = asin(sphere_radius / sphere_location.x);
	//real_type solid_angle = 2 * pi * (1 - cos(rho));

	//real_type area = 4 * pi * sphere_location.x * sphere_location.x;
	//real_type x = area / (sphere_radius * sphere_radius);
	//real_type solid_angle = x;
	
	real_type ratio = solid_angle / (4.0 * pi);
	return (static_cast<real_type>(unit_vectors.size()) * ratio);
}

Am I calculating the solid angle properly?

Yes, it's sort of like the Solar irradiance problem, except that this has to do with the gradient.

taby said:
Am I calculating the solid angle properly?

I did this recently to set the fov of my light so it bounds the bou8nding sphere of the scene:

lightFOV = asin(sceneRadius / lightDist);

Not sure if this gives the half or the full angle, though.