Function approximation

I am successfully running a numerical algorithm to determine the collision between gravitational field lines and an axis-aligned bounding box – https://github.com/sjhalayka/ellipsoid_emitter

I am trying to come up with an approximation. Below is one of many ways I've tried to attack the problem. How would you do it differently?

MyBig get_intersecting_line_count_real(
	const MyBig n,
	const MyBig box_location_x,
	const MyBig box_length,
	const MyBig D)
{
	MyBig running_total_small = 0;
	MyBig running_total_big = 0;

	MyBig running_ratio = 0;

	const MyBig big_area =
		4 * pi
		* box_location_x * box_location_x;

	const MyBig big_length =
		2 * pi
		* box_location_x;

	const MyBig small_area =
		box_length * box_length;

	const MyBig small_length =
		box_length;

	const MyBig disk_like = 3 - D;

	const MyBig start_distance = 0;
	const MyBig end_distance = 1;
	const size_t distance_res = 10000;

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

	for (size_t step_index = 0;
		step_index < distance_res;
		step_index++)
	{
		const MyBig rn =
			start_distance +
			step_index * distance_step_size;

		if (rn <= disk_like)
			running_ratio += 1;
	}

	const MyBig r =
		running_ratio 
		/ static_cast<MyBig>(distance_res);

	const MyBig big = lerp_func(big_area, big_length, r);
	const MyBig small = lerp_func(small_area, small_length, r);

	return n * (small/big);
}