Hengefinder: Finding When the Sun Aligns With Your Street

Next week in Manhattan, the sunset will align perfectly with the east-west streets of the city grid. It’s beautiful, and people know it. Twice a year, crowds gather to see the brief moment when the sun sits perfectly on the horizon, framed by skyscrapers on either side. It’s called Manhattanhenge, after Stonehenge. I wanted to know how astronomers figure out when Manhattanhenge happens. And if I could figure that out, why limit it to Manhattan?

A shot of crowds taking photos of Manhattanhenge at 42nd St in NYC.

This was one of my first projects at the Recurse Center: Hengefinder, a tool that lets you find a henge pretty much anywhere the sun sets.

Source code Hengefinder website Hengefinder mobile app (A follow up to the website, created by fellow Recurser John Pribyl)

The basic steps (and some terminology I would soon learn):

find the angle of a road (its bearing, relative to true north) find the angle of the sun at sunset each day (its azimuth) find the dates when those angles match.

It was appealing to me that this was mostly made up of many sub problems I could either think through thoroughly, or choose to mostly skip (by handing off to a library, a brute-force solution, an approximation, etc.). It’s like a bunch of closed boxes. I left plenty of those boxes closed — I didn’t build my own astronomical model, for instance. Other boxes I opened. And then I found many things I thought would be straightforward that turned out to be somewhat less trivial in practice as my assumptions broke; Like, in reality, you can’t treat roads like flat lines, latitude and longitude don’t behave like a Cartesian grid, and the word “sunset” is more ambiguous than I initially thought. I’ve enjoyed closing the gap between my assumptions and reality. So, rather than just talking about the project as a whole, I’m going to walk through a couple challenges I ran into in building this, and how I went about solving them. One challenge for each of those supposedly “simple” steps above. Challenge #1: Finding the road bearing (and rediscovering the Earth is not flat) The first challenge was calculating the bearing of a street (its angle relative to true north).

If you take the latitude and longitude of one address, and the latitude and longitude of another address down the street, you end up with the coordinates of two points on Earth. We can then use some trig to get the angle. My initial (incorrect) guess for getting the road bearing was to take the difference in latitude and the difference in longitude, then just get the angle with atan2 (the inverse tangent). The problem is, that would only work if the Earth were flat.

A schematic of the orientation of latitude and longitude. Degrees of longitude represent smaller distances as you move away from the equator.

Latitude and longitude lines aren’t arranged the same way on the Earth. Latitude lines (which run east-west) are evenly spaced: one degree of latitude is basically the same distance everywhere on Earth. Longitude lines (which run north-south) aren’t. They converge as you move toward the poles. For example, one degree of longitude spans roughly 52 miles (84 km) in New York City, but way less at higher latitudes (about 33 miles / 53 km in Anchorage). How do you manage that? TL;DR, if you just want the practical takeaway Scale longitude by cos(latitude), which puts longitude and latitude in the same “units.” Only then do atan2. ”Ok, but why?” (And if you don’t really care about the math, feel free to skip to the “corrected approach” section below.) When you do tan = opposite / adjacent, you’re assuming the dimensions of the triangle are in consistent units. They aren’t here. Like I said, degrees of longitude shrink as you move toward the poles, while latitude stays consistent. That means the horizontal and vertical sides of the triangle are in different units.

The same east-west distance maps on to different changes in longitude depending on latitude.

For a perfect “henge” moment, the sun needs to be sitting directly on the horizon. I used the Python library Astral to compute the timing, altitude, and azimuth (angle) of solar events at a given location. But Astral’s definition of a sunset is slightly different from what I want here. Astral uses the standard astronomical definition of sunset, i.e. when the sun has fully dipped below the horizon. For my purpose, that’s too late. I want to see when the sun is sitting right on top of the horizon. Astral’s version is still useful to me, however, since time-wise, it gets me in the ballpark on a given day. For a moment to be a “henge,” I want the sun’s disk to be fully visible and just touching the horizon. I have a target “altitude” I want the sun to be at, above the horizon, based on the apparent size of the sun.

The technical definition of a sunset is too late for what is needed for a henge.

Why this is a boundary search, not a value search So I’m trying to find the last moment when the sun’s altitude is still above my target threshold altitude. You could just do a linear search minute-by-minute in some time window before Astral’s sunset. It’d probably be fine, but each evaluation is an Astral API call. Plus, that’s less fun. This is one of the “boxes” I wanted to open. Because the sun’s altitude changes monotonically as sunset approaches, this indicates a binary search. The textbook approach to binary search is this: def basic_binary_search(arr, target): left = 0 right = len(arr) - 1 while left <= right: mid = (left + right) // 2 if arr[mid] == target: return mid elif arr[mid] < target: left = mid + 1 else: right = mid - 1 # Not found return -1 But that formulation makes several assumptions that don’t really apply here. For one, we’re not searching for an exact target match. Because with minute-level resolution, we may never get the exact moment.

Instead, in looking for a henge, I’m really evaluating a boolean each minute: “Is the sun still above the target altitude at this minute?” That gives a pattern like: True, True, True, True, False, False, False Rather than looking for a specific value, now we’re just looking for the last True before it flips.

A boundary (last-true) binary search So, I used a different approach to the binary search. Here’s the core loop, searching within a window of time before Astral’s sunset: while left < right: mid = (left + right + 1) // 2 # upper-biased midpoint for a "last true" search if altitude > target: # still valid, keep it left = mid # this minute is valid, could be the last one else: # sun at or below horizon right = mid - 1 # this minute is invalid, look earlier This turns it into a “last true” binary search, which lets me find the final minute when the sun is still above the target altitude. (This mirrors the general binary search template in this great LeetCode article, which argues this is a less error-prone approach overall to binary search.) Challenge #3: Finding when bearing = azimuth (a two-phase search) Once I knew how to get (a) a road’s bearing and (b) the sun’s azimuth at sunset, I just had to find when the two matched to get a henge date. I couldn’t do a binary search across the next 365 days, since the sun’s azimuth over the year is non-monotonic. I could just do the brute force approach and compute the sunset azimuth for every day of the year and check when (if ever) it matches the road bearing. That would have worked fine, it’s another box I could have left closed. But it felt like the wrong way to approach the problem. And again, it would make a lot of unnecessary API calls. Relationship between sunset azimuth and date If you plot the sun’s azimuth at a particular location, you get a smooth curve:

At a given location, the sun's sunset azimuth varies smoothly over the year and reverses direction around the solstices.

As a result, depending on the road’s bearing, there may be two alignment dates, one, or none at all.

Streets closer to north-south alignment, for example, are just plainly out of reach at most latitudes. You could model this function analytically and solve for the exact alignment dates. That’s how astronomers might approach it. But doing that correctly would mean committing to a pretty deep astronomical model, beyond the scope of what I wanted (the curve might look like a sine wave, but it isn’t quite). My goal wasn’t to do this perfectly, it was to build a tool that was easy to reason about and would reliably find me some henges. I ended up using a compromise as a two-phase search. Phase 1: Coarse search I start by sampling at coarse intervals. Let’s say, once every 30 days. If you sample in this coarse way, the goal is to identify a period in which we may have “missed” a henge, so that we can switch to a fine-grained search over that period. I do this by tracking not just the sun’s azimuth, but which side of the road bearing the sunset azimuth is on, and which direction it’s moving. Here’s why. When sampling coarsely, a “missed” henge can happen in one of two ways. The first way you might miss a henge in the coarse search is that the sun’s azimuth might move steadily toward the road bearing and cross it between samples, continuing in the same direction. That is effectively the Intermediate Value Theorem:

If a continuous quantity changes sign over an interval, it must have crossed zero inside that interval.

If that happens, a henge definitely occurred. The second way to miss a henge is if the sun’s azimuth reverses direction between samples. It’s possible a henge occurred during that reversal. It’s also possible that there was no henge, and that the road bearing is just “out of reach” of the sun’s maximum azimuth (for instance, there is no day in the year that the sunset azimuth in NYC would match a road with a bearing of 320˚). Either way, I’ve identified a finer window worth exploring.

Two ways coarse sampling can miss a henge, each shown as a pair of azimuth samples. The pink dots are on either side of the road bearing. A henge date is guaranteed to be in that window, which is flagged for fine-grained search.