Thanks for the function! I tried using asin, and it draws a nice curve, but it's not quite the same as the curve generated by the brute force line intersection method.

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

I've also just tried the obvious, but it doesn't quite work either. Basically, whatever the curve is, multiplying it by r^3 does not form a straight line, like it should (?):

long long unsigned int get_intersecting_line_count(
	const vector<vector_3>& unit_vectors,
	const vector_3 sphere_location,
	const real_type sphere_radius)
{
	real_type big_area = 4 * pi * sphere_location.x * sphere_location.x;
	real_type small_area = pi * sphere_radius * sphere_radius;
	real_type ratio = small_area / big_area;
	
	return static_cast<long long unsigned int>(static_cast<real_type>(unit_vectors.size()) * ratio);
}

Edit: Aha, right… this only works for large r.

Edit: OK, I'm back on track… I had to divide the gradient strength by mass squared (or area).

real_type get_intersecting_line_count(
	const real_type n,
	const vector_3 sphere_location,
	const real_type sphere_radius)
{
	const real_type big_area =
		4 * pi
		* sphere_location.x * sphere_location.x;

	const real_type small_area =
		pi
		* sphere_radius * sphere_radius;

	const real_type ratio =
		small_area
		/ big_area;

	return n * ratio;
}

…

real_type get_intersecting_line_count(
	const real_type n,
	const vector_3 sphere_location,
	const real_type sphere_radius)
{
	const real_type big_area = 4 * sphere_location.x * sphere_location.x;
	const real_type small_area = sphere_radius * sphere_radius;
	const real_type ratio = small_area / big_area;
	
	return n * ratio;
}


int main(int argc, char** argv)
{
	// Field line count
	const real_type n = 1e70;

	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;

		const real_type epsilon = 0.1;

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

		// https://en.wikipedia.org/wiki/Directional_derivative
		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 real_type gradient = (collision_count_plus - collision_count) / epsilon;
		const real_type gradient_strength = -gradient / (4 * pi * receiver_radius * receiver_radius);

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

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

	out_file.close();

	return 0;
}

I've started on the paper. The latest version is always at: https://github.com/sjhalayka/papers/blob/main/newtonian%20gravitation/bezier_escape.pdf

Now I just have to do the English and the math parts, now that the code is done.

OK, I've got it figured out now, I think. Here is the code that computes the gradient in two different ways:

int main(int argc, char** argv)
{
	// Field line count
	// re: holographic principle:
	const real_type emitter_radius = 1.0;
	const real_type emitter_area = 4.0 * pi * emitter_radius * emitter_radius;
	const real_type n = (c3 * emitter_area) / (log(2.0) * 4.0 * G * hbar);
	const real_type emitter_mass = c2 * emitter_radius / (2*G);

	const 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 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 real_type gradient =
			(collision_count_plus - collision_count)
			/ epsilon;

		const real_type gradient_strength =
			-gradient
			/ (receiver_radius * receiver_radius);

		//const real_type gradient_strength_2 = 
		//	n / (2.0 * pow(receiver_pos.x, 3.0));

		// This should be clsoe to 1.0:
		// cout << gradient_strength_2 / gradient_strength << endl;

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

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

	out_file.close();

	return 0;
}

Thank you again JoeJ for all of your help, insight, and time.

I also found that G can be calculated using constants, for example:

const real_type G_2 = (c3 * pi) / (log(2.0) * hbar * 1.73502e+70);

Here 1.73502e+70 is the 't Hooft-Susskind constant… the number of field lines for a black hole emitter of unit Schwarzschild radius, as per the holographic principle.

I noticed that the following 5 variables plot the same curve in gnuplot, at large distance receiver_pos.x:

// These three are roughly equal in value:
real_type gradient_strength = -gradient / (receiver_radius * receiver_radius);
real_type gradient_strength_ = n / (2.0 * pow(receiver_pos.x, 3.0));
real_type gradient_strength__ = (c3 * emitter_area) / (log(2.0) * 8.0 * G * hbar * pow(receiver_pos.x, 3.0));

// These two are not roughly equal in value, but still form the same curve:
real_type gradient_strength___ = 1.0 / pow(receiver_pos.x, 3.0);
real_type gradient_strength____ = G * emitter_mass / pow(receiver_pos.x, 3.0);

What is it called when they form the same curve?

I figured it out. I can now calculate the Newtonian acceleration:

real_type get_intersecting_line_count(
	const real_type n,
	const vector_3 sphere_location,
	const real_type sphere_radius)
{
	const real_type big_area =
		4 * pi
		* sphere_location.x * sphere_location.x;

	const real_type small_area =
		pi
		* sphere_radius * sphere_radius;

	const real_type ratio =
		small_area
		/ big_area;

	return n * ratio;
}


int main(int argc, char** argv)
{
	const real_type emitter_radius = 1.0;
	
	const real_type emitter_area = 
		4.0 * pi * emitter_radius * emitter_radius;

	// Field line count
	// re: holographic principle:
	const real_type n = 
		(c3 * emitter_area) 
		/ (log(2.0) * 4.0 * G * hbar);
	
	const real_type emitter_mass = c2 * emitter_radius / (2.0 * G);

	// 1.73502e+70 is the 't Hooft-Susskind constant:
	// the number of field lines for a black hole of
	// unit Schwarzschild radius
	//
	//const real_type G_ = 
	//	(c3 * pi) 
	//	/ (log(2.0) * hbar * 1.73502e+70);

	const string filename = "newton.txt";
	ofstream out_file(filename.c_str());
	out_file << setprecision(30);

	const real_type start_distance = 10.0;
	const real_type end_distance = 100.0;
	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 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 real_type gradient =
			(collision_count_plus - collision_count)
			/ epsilon;

		real_type gradient_strength =
			-gradient
			/ (receiver_radius * receiver_radius);

		const real_type gradient_strength_ = 
			n / (2.0 * pow(receiver_pos.x, 3.0));

		const real_type newton_strength = 
			n * c * hbar * log(2.0)
			/ (pow(receiver_pos.x, 2.0) * emitter_mass * 4.0 * pi);

		cout << newton_strength / (G*emitter_mass/pow(receiver_pos.x, 2.0)) << endl;

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

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

	out_file.close();

	return 0;
}

In Newtonian gravity, the isotropic emitter ensures that space is curved.

This is not quite the same as what happens in general relativity \cite{misner}.

In general relativity, time is also curved.

One way to try to model the curvature of time is to produce graviton overlap by using pseudorandom direction vectors.

So, yes, up close, there is a difference between using pseudorandom normals, versus normals that are equal to the oscillator positions on the unit sphere. Of course, at large r (e.g. r = 4), the regular old Newtonian gravitation occurs.

Edit: Sorry for the empty post.