OK, here's the code. I'm not sure where to go from here:
int main(int argc, char** argv)
{
cout << setprecision(20) << endl;
srand(0);
double dimension = 2.5;
const size_t n = 10000000;
if (dimension <= 2)
dimension = 2.001;
else if (dimension > 3)
dimension = 3;
const double disk_like = 3 - dimension;
const double falloff_exponent = 2 - disk_like;
// Start with pseudorandom oscillator locations
for (size_t i = 0; i < n; i++)
threeD_oscillators.push_back(RandomUnitVector());
// Spread the oscillators out, so that they are distributed evenly across
// the surface of the ellipsoid emitter
for (size_t i = 0; i < n; i++)
{
vector_3 ring;
ring.x = threeD_oscillators[i].x;
ring.y = 0;
ring.z = threeD_oscillators[i].z;
vector_3 s = slerp(threeD_oscillators[i], ring, disk_like);
threeD_oscillators[i] = s;
}
// Get position on oblate ellipsoid
for (size_t i = 0; i < n; i++)
{
vector_3 vec = threeD_oscillators[i];
const vector_4 rv = RayEllipsoid(vector_3(0, 0, 0), vec, vector_3(1.0, 1.0 - disk_like, 1.0));
threeD_oscillators[i] = vector_3(rv.y, rv.z, rv.w);
}
// Get position and normal on prolate ellipsoid
for (size_t i = 0; i < n; i++)
{
vector_3 vec = threeD_oscillators[i];
const vector_4 rv = RayEllipsoid(vector_3(0, 0, 0), vec, vector_3(1.0 - disk_like, 1.0, 1.0 - disk_like));
vector_3 collision_point = vector_3(rv.y, rv.z, rv.w);
vector_3 normal = EllipsoidNormal(collision_point, vector_3(1.0 - disk_like, 1.0, 1.0 - disk_like));
normals.push_back(normal);
line_segment_3 ls;
ls.start = threeD_oscillators[i];
ls.end = threeD_oscillators[i] + normals[i] * 1e30;
threeD_line_segments.push_back(ls);
}
// Get intersecting lines
const double start_distance = 10;
const double end_distance = 10000;
const size_t resolution = 100000;
const double step_size = (end_distance - start_distance) / (resolution - 1);
for (double r = start_distance; r <= end_distance; r += step_size)
{
vector_3 receiver_pos(r, 0, 0);
size_t collision_count = get_intersecting_line_segments(receiver_pos, 1.0, dimension);
cout << collision_count * pow(receiver_pos.x, falloff_exponent) << endl;
}
return 0;
}
Basically, the value collision_count * pow(receiver_pos.x, falloff_exponent)
is constant where D = 3.0 and D = 2.001. For D = 2.5, it is not constant.