Hard cutoff in Atmospheric Refraction near the horizon

[Dmytro-Malets]

Hi. I am using Skyfield to calculate precise celestial object positions for a visualization. I have observed a significant, non-physical jump (discontinuity) in the calculated apparent altitude precisely when the object crosses the geometric altitude threshold of $-1.0^\circ$.

This issue is visible as a sudden, sharp drop in the object's apparent position, which creates an artifact in any system that requires smooth, continuous altitude tracking.

This behavior seems to stem from the hard cutoff implemented in the core refraction logic, where the full correction d is applied only within a strict range:

def refraction(alt_degrees, temperature_C, pressure_mbar):
    """Given an observed altitude, return how much the image is refracted.

    Zero refraction is returned both for objects very near the zenith,
    as well as for objects more than one degree below the horizon.

    """
    r = 0.016667 / tan((alt_degrees + 7.31 / (alt_degrees + 4.4)) * DEG2RAD)
    d = r * (0.28 * pressure_mbar / (temperature_C + 273.0))
    return where((-1.0 <= alt_degrees) & (alt_degrees <= 89.9), d, 0.0)


def refract(alt_degrees, temperature_C, pressure_mbar):
    """Given an unrefracted `alt` determine where it will appear in the sky."""
    alt = alt_degrees
    while True:
        alt1 = alt
        alt = alt_degrees + refraction(alt, temperature_C, pressure_mbar)
        converged = abs(alt - alt1).max() <= 3.0e-5
        if converged:
            break
    return alt

When alt_degrees is $<-1.0^\circ$, the refraction correction is instantly set to $0.0^\circ$, causing the jump.

Why was the specific value of $-1.0^\circ$ chosen as the strict lower limit for applying the refraction correction, rather than a deeper limit like $-3^\circ$ or $-5^\circ$? Does this choice stem from a known limitation of the underlying Bennett formula at that specific point? Given that the refraction effect does not vanish instantaneously in the physical world, why was a smooth transition (e.g., tapering, blending, or damping) not implemented to bridge the gap between the full correction at $-1.0^\circ$ and zero correction further below the horizon?

For applications requiring smooth altitude progression, this discontinuity is problematic. I have already tested a solution using a cosine damping function that gradually fades the refraction correction to zero between $-3.0^\circ$ and $-5.0^\circ$ of geometric altitude.

$>89.9^\circ$: No refraction
$-3.0^\circ$ to $89.9^\circ$: Full refraction (Bennett's formula).
$-3.0^\circ$ to $-5.0^\circ$: Smooth cosine damping from 100% to 0%.
$<-5.0^\circ$: No refraction.

This method entirely removes the jump. If you are interested in this approach, I can share the implementation details.

[okkine]

I'm intrigued. Are you completely bypassing Skyfields's refraction calculations, and implementing your own?

I'm building a sensor that updates the elevation at every 0.5 degrees, by using find_descrete to find when it crosses the next 0.5 degree threshold, and then calling altaz() at the returned time. Works perfectly, except at elevations between about -1, and -0.4. I suspect the problem is related to Skyfield's hard cutoff for refraction at -1.

I definitely like how refraction was handled better in pyephem.

Edit: I believe I know what you're saying, and I'm working my own plan to accomplish the same! Thanks for the idea!

[brandon-rhodes]

Why was the specific value of − 1.0 ∘ chosen as the strict lower limit…

When Skyfield was first written, its goal was to give exactly the positions returned by the science-grade USNO library NOVAS which, as it happens, cuts off refraction at –1.0° because by that point the object is below the horizon and isn't visible.

And so Skyfield did the same, because it meant that the two libraries could be compared, for testing purposes, to very high precision. Which was useful because NOVAS was written by actual professionals in astrometry, whereas I wrote Skyfield merely as an amateur astronomer who happened to love Python. 🙂

But the sudden jump is indeed annoying. Maybe it could be kept, but be hidden behind a compatibility flag (like usno=True) that would be False by default but that test code could use when comparing Skyfield results to NOVAS?

In your own code, I assume that –3 and –5 are arbitrary design choices, just as –1 was for the NOVAS authors?

[okkine]

I can't speak for Dmytro, but that sounds like a good solution for me. I've built my own workaround, but having built in smoothing would be ideal. Even if it's disabled by default, but available to those who need it. Whatever you had in Pyephem worked well for my needs.

[Dmytro-Malets]

In your own code, I assume that –3 and –5 are arbitrary design choices, just as –1 was for the NOVAS authors?

Yes, the -3.0° and -5.0° thresholds are empirical design choices. I derived them based on visual observations from software like Stellarium and data consistency from TimeAndDate, specifically aiming for visual smoothness in animations.

I think your suggestion of a compatibility flag (e.g., usno=True) is a great solution. I have already implemented this locally to verify the concept. I updated refraction() to use cosine damping by default (usno=False) and added the usno=True flag to the tests in test_against_novas.py.

I am ready to submit a Pull Request with these changes. Please let me know if you would like me to proceed!

Below is a comparison of the calculated Apparent Altitude of the Sun during sunset on Dec 11, 2025, from Maputo, demonstrating the difference in behavior:

Local Time	Geometric Alt	usno=True	usno=False
18:35:00	-0.8109	-0.1932	-0.1932
18:35:30	-0.9116	-0.2745	-0.2745
18:36:00	-1.0122	-1.0122	-0.3552
18:36:30	-1.1129	-1.1129	-0.4354
18:37:00	-1.2134	-1.2134	-0.5152
18:37:30	-1.3139	-1.3139	-0.5946
18:38:00	-1.4144	-1.4144	-0.6737
18:38:30	-1.5148	-1.5148	-0.7525
18:39:00	-1.6152	-1.6152	-0.8313
18:39:30	-1.7155	-1.7155	-0.9100
18:40:00	-1.8158	-1.8158	-0.9890
18:40:30	-1.9160	-1.9160	-1.0684
18:41:00	-2.0162	-2.0162	-1.1485
18:41:30	-2.1163	-2.1163	-1.2296
18:42:00	-2.2164	-2.2164	-1.3121
18:42:30	-2.3164	-2.3164	-1.3966
18:43:00	-2.4164	-2.4164	-1.4836
18:43:30	-2.5163	-2.5163	-1.5741
18:44:00	-2.6162	-2.6162	-1.6693
18:44:30	-2.7160	-2.7160	-1.7708
18:45:00	-2.8158	-2.8158	-1.8812
18:45:30	-2.9155	-2.9155	-2.0041
18:46:00	-3.0151	-3.0151	-2.1452
18:46:30	-3.1147	-3.1147	-2.3130
18:47:00	-3.2143	-3.2143	-2.5153
18:47:30	-3.3138	-3.3138	-2.7427
18:48:00	-3.4133	-3.4133	-2.9664
18:48:30	-3.5127	-3.5127	-3.1812
18:49:00	-3.6120	-3.6120	-3.3842
18:49:30	-3.7113	-3.7113	-3.5596
18:50:00	-3.8105	-3.8105	-3.7118
18:50:30	-3.9097	-3.9097	-3.8473
18:51:00	-4.0089	-4.0089	-3.9714
18:51:30	-4.1079	-4.1079	-4.0793
18:52:00	-4.2070	-4.2070	-4.1835
18:52:30	-4.3059	-4.3059	-4.2875
18:53:00	-4.4048	-4.4048	-4.3910
18:53:30	-4.5037	-4.5037	-4.4939
18:54:00	-4.6025	-4.6025	-4.5961
18:54:30	-4.7013	-4.7013	-4.6976
18:55:00	-4.8000	-4.8000	-4.7983
18:55:30	-4.8986	-4.8986	-4.8982
18:56:00	-4.9972	-4.9972	-4.9972
18:56:30	-5.0957	-5.0957	-5.0957
18:57:00	-5.1942	-5.1942	-5.1942