Ring Shadows


I needed to draw the shadow of a planet on its rings. To render the rings, each point on a disk around the planet is drawn separately, and the inputs are the angle of the line from the center to the point and the length of that line. That turns out to be a polar coordinate system, where each point is parametrized by some length \(r\) and angle \(\theta\).

Now the shadows. Since the light source (Sun) is pretty far away (like a trillion meters or something) we can consider its rays parallel. That means that in an ideal world the planet's "shadow" (the portion of space that light doesn't reach) is a cylinder the same radius as the planet, and its caps perpendicular to the light rays. Since our rings are almost parallel to the planet's orbit, the shadow intersects the ring disk in a sort of rectangular shape.

\(R_p\) is the radius of the planet (in red) relative to the rings (so that the outer disk radius is \(1\))

Notice how the angles of the two lines from the origin to the shadow's edge have different angles? That might feel complicated at first, since straight lines in polar coordinates are kind of weird! But let's simplify things and take a single circular "slice" of the disk, at radius \(r\).

\(\alpha\) is the angle of the Sun, \(\beta\) and \(\beta'\) are the angles of the shadow's edges.

Now, since \(r\) is fixed, we only need to check if \(\theta\) falls into our shadow. That boils down to figuring out if \(\theta\) is between the angles \(\beta\) and \(\beta'\). After a bit of trigonometry, we can conclude that

\(\beta,\beta' = \alpha \pm \arcsin\left(\frac{R_p}{r}\right)\)

Ok, looks like we're done! Yay! Here's the shader code:

bool is_in_shadow(
  float radius,   // = r
  float angle,    // = θ
  float sun_angle // = α
) {
  float d = acos(relative_planet_radius / radius);
  float l = angle - d; // = β or β'
  float r = angle + d; // = β' or β
  return l <= angle && angle >= r;
}

If we actually run this, it works! But let's fast-forward a couple hundred years (Saturn has a big orbit) just to be sure that it doesn't break.

I'm too lazy sorry :( Saturn from above, but only half of the shadow is visible.

Oh...

Let's try to figure out what's happening: the shadow is getting cut off at a certain angle and then it appears again. Oh! (after like an hour) Looks like we forgot that \(\alpha\) and \(\theta\) are in the range \([-\pi, \pi]\), since that's what atan2 returns1, but the \(\beta\)s can be outside of it! Let's try to come up with an example that breaks our check. Say \(\theta\) is \(-3\), and \((\beta, \beta')=(\pi - 0.5, \pi + 0.5)\). On a unit circle it would look like this:

Notice how \(\theta\) looks like its inside of the blue region, but numerically, the angles don't overlap: \(\pi-0.5< -3 < \pi+0.5\) is false. To actually check if the angle is inside, we need to add \(\pm 2\pi\) to it, to bring it in the same period as the edges. Interestingly, we don't even need to check for this situation explicitly (ifs are costly on the GPU)! We can compare \(\theta\), \(\theta+2\pi\), and \(\theta-2\pi\) with the \(\beta\)s regargless of what the values are, and if any of the results is true, the point is between the angles. This is how it would look in code:

// returns true if a <= b <= c.
bool is_ordered(float a, float b, float c) {
  return a <= b && b <= c;
}

bool is_in_shadow(
  float radius,   // = r
  float angle,    // = θ
  float sun_angle // = α
) {
  float d = acos(relative_planet_radius / radius);
  float l = angle - d; // = β or β'
  float r = angle + d; // = β' or β
  return
    is_ordered(l, angle - 2.0 * PI, r) ||
    is_ordered(l, angle, r) ||
    is_ordered(l, angle + 2.0 * PI, r);
}

Here's the pretty image:

I'm too lazy sorry :( Saturn from the side, camera facing the sun.

1. I'm using atan2 to get the angles, but they could just as well have been constrained to \([0, 2\pi]\), the problem would still exist, just at different angles.