I don’t know the sugya well but I think there are plenty that hold that zmanim are based off when sunset and sunrise are actually based on you’re height, and the height at which the sun disappears. meaning you’d have to calculate the zmanim for youre exact location (city isn’t enough) and calculate if there’s any terrain blocking the path between yourself and the sun which there almost always is closer to nightfall and sunrise and calculate back finding the time sunset actually is. there’s no such program out there that does this as far as I’m aware but I tried creating one but since there’s no such program i can’t test it’s accuracy. I believe its extremely inaccurate. just putting it out there if anyone is interested checking it out and maybe improving
python sloppy code
import datetime
import pytz
import numpy as np
import re
from skyfield.api import load, Topos
from skyfield import almanac
import srtm # You may need to run: pip install srtm.py
import rasterio
— 1. Coordinate Parsing Helper —
def parse_dms(dms_str):
“”“Parses a DMS coordinate string (e.g., “N 41° 5’ 59.8"”) into decimal degrees.”“”
parts = re.split(‘[^\d\w.]+’, dms_str)
deg = float(parts[1])
mnt = float(parts[2])
sec = float(parts[3])
decimal = deg + (mnt / 60.0) + (sec / 3600.0)
if parts[0].upper() in ('S', 'W'):
return -decimal
return decimal
— 2. Core Calculation Functions —
def get_terrain_adjusted_event(event_type, latitude, longitude, altitude, date_obj):
“”"
Calculates the terrain-adjusted time for sunrise or sunset.
Args:
event_type (str): Either 'sunrise' or 'sunset'.
latitude (float): Latitude of the observer.
longitude (float): Longitude of the observer.
altitude (float): Altitude of the observer in meters.
date_obj (datetime.date): The date for the calculation.
Returns:
datetime.datetime: The terrain-adjusted event time in UTC, or None.
"""
# --- Setup ---
ts = load.timescale()
eph = load('de421.bsp')
observer_location = Topos(latitude_degrees=latitude, longitude_degrees=longitude, elevation_m=altitude)
# --- Find Standard Event Time (Sea Level) ---
t0 = ts.utc(date_obj.year, date_obj.month, date_obj.day)
t1 = ts.utc(date_obj.year, date_obj.month, date_obj.day + 1)
f = almanac.sunrise_sunset(eph, observer_location)
times, events = almanac.find_discrete(t0, t1, f)
event_code = 1 if event_type == 'sunrise' else 0
is_event = (events == event_code)
if not np.any(is_event):
print(f"Standard {event_type} not found for this date.")
return None
standard_event_time = times[is_event][0]
# --- Get Azimuth and Calculate Horizon ---
sun = eph['sun']
earth = eph['earth']
astrometric = (earth + observer_location).at(standard_event_time).observe(sun)
_, az, _ = astrometric.apparent().altaz()
event_azimuth = az.degrees
print(f"Calculating horizon profile for {event_type} at azimuth {event_azimuth:.2f}°...")
horizon_angle = calculate_horizon_angle(latitude, longitude, altitude, event_azimuth)
print(f"Calculated horizon angle: {horizon_angle:.2f}°")
# --- Search for Terrain-Adjusted Time ---
search_start = standard_event_time.utc_datetime() - datetime.timedelta(hours=1)
search_end = standard_event_time.utc_datetime() + datetime.timedelta(hours=1)
time_step = datetime.timedelta(minutes=1)
current_time = search_start
# For sunrise, we look for the moment the sun *crosses above* the horizon angle.
# For sunset, we look for the moment the sun *dips below* the horizon angle.
# Find the sun's altitude at the start of our search
t_start = ts.from_datetime(current_time.replace(tzinfo=pytz.utc))
alt_start, _, _ = (earth + observer_location).at(t_start).observe(sun).apparent().altaz()
is_below_horizon_at_start = alt_start.degrees < horizon_angle
while current_time <= search_end:
t = ts.from_datetime(current_time.replace(tzinfo=pytz.utc))
alt, _, _ = (earth + observer_location).at(t).observe(sun).apparent().altaz()
is_currently_below = alt.degrees < horizon_angle
if event_type == 'sunrise' and is_below_horizon_at_start and not is_currently_below:
# We started below and are now above: this is sunrise
return current_time
elif event_type == 'sunset' and not is_below_horizon_at_start and is_currently_below:
# We started above and are now below: this is sunset
return current_time
current_time += time_step
print(f"Could not find terrain-adjusted {event_type} in the search window.")
return None
def calculate_horizon_angle(lat, lon, alt, azimuth, max_distance_km=20):
“”"
Calculates the horizon angle in a specific direction using the elevation library.
“”"
# It will download necessary data and cache it in a file called 'srtm_1_30m_41_75.tif' or similar.
import elevation
import rasterio
# Ensure the elevation data for the area is downloaded.
# The output is the path to the data file.
elevation.clip(bounds=[lon-0.5, lat-0.5, lon+0.5, lat+0.5], output='dem.tif')
dem_path = '/home/johnw/.cache/elevation/SRTM1/dem.tif' # Use your downloaded DEM file
min_distance_m = 1000
num_points = 100
angles = []
with rasterio.open(dem_path) as dem:
for i in range(1, num_points + 1):
distance = min_distance_m + ((i / num_points) * (max_distance_km * 1000 - min_distance_m))
# Calculate the coordinates of the point to check
d_lat = distance * np.cos(np.radians(azimuth)) / 111320
d_lon = distance * np.sin(np.radians(azimuth)) / (111320 * np.cos(np.radians(lat)))
point_lat = lat + d_lat
point_lon = lon + d_lon
try:
# Get elevation from the downloaded DEM file
terrain_elevation = next(dem.sample([(point_lon, point_lat)]))[0]
if terrain_elevation > -1000: # Filter out invalid data
height_diff = terrain_elevation - alt
angle = np.degrees(np.arctan2(height_diff, distance))
angles.append(angle)
except Exception:
# This can happen if we are at the edge of the available data, just ignore the point
pass
return max(angles) if angles else 0.0
if name == ‘main’:
import elevation # Import here to use for altitude lookup
print("--- Terrain-Adjusted Zmanim Calculator ---")
# --- Get Location Input ---
lat_str = input("Enter Latitude (e.g., 41.0999 or 'N 41 5 59'): ") or '41.087529'
lon_str = input("Enter Longitude (e.g., -74.0083 or 'W 74 0 30'): ") or '-74.060390'
# --- Parse Coordinates ---
try:
if '°' in lat_str or re.search('[a-zA-Z]', lat_str):
latitude = parse_dms(lat_str)
longitude = parse_dms(lon_str)
else:
latitude = float(lat_str)
longitude = float(lon_str)
except Exception as e:
print(f"Error: Could not understand coordinates. Please try again. Details: {e}")
exit()
print(f"\nParsed Location -> Lat: {latitude:.4f}, Lon: {longitude:.4f}")
# --- Automatically Get Altitude ---
try:
print("Fetching altitude for your location...")
dem_path = elevation.clip(bounds=[longitude-0.1, latitude-0.1, longitude+0.1, latitude+0.1], output='dem_altitude.tif')
with rasterio.open(dem_path) as dem:
auto_altitude = round(next(dem.sample([(longitude, latitude)]))[0])
print(f"Detected altitude: {auto_altitude} meters.")
manual_altitude = input(f"Press Enter to accept, or type a different altitude (in meters): ")
if manual_altitude.strip():
altitude = int(manual_altitude)
else:
altitude = auto_altitude
except Exception as e:
print(f"Could not automatically determine altitude. Error: {e}")
altitude = int(input("Please enter your altitude in meters: ")) or 111
# --- Get Date and Timezone ---
default_date = datetime.date.today().strftime('%Y-%m-%d')
date_str = input(f"Enter Date (YYYY-MM-DD) [default: {default_date}]: ") or default_date
# You can find your timezone here: https://en.wikipedia.org/wiki/List_of_tz_database_time_zones
local_timezone_str = input("Enter your Timezone (default: America/New_York): ") or "America/New_York"
print("-" * 40)
# --- Calculations ---
target_date = datetime.datetime.strptime(date_str, '%Y-%m-%d').date()
local_tz = pytz.timezone(local_timezone_str)
sunrise_utc = get_terrain_adjusted_event('sunrise', latitude, longitude, altitude, target_date)
print("-" * 40)
sunset_utc = get_terrain_adjusted_event('sunset', latitude, longitude, altitude, target_date)
print("-" * 40)
# --- Display Results ---
if sunrise_utc and sunset_utc:
sunrise_local = sunrise_utc.astimezone(local_tz)
sunset_local = sunset_utc.astimezone(local_tz)
print("RESULTS (Terrain-Adjusted):")
print(f" Sunrise: {sunrise_local.strftime('%I:%M:%S %p %Z')}")
print(f" Sunset: {sunset_local.strftime('%I:%M:%S %p %Z')}")
# Calculate Sha'os Zmaniyos (Halachic Hours)
day_duration = sunset_utc - sunrise_utc
shaah_zmanis = day_duration / 12
# Calculate Zmanim based on Gr"a
szks_gra = sunrise_utc + (3 * shaah_zmanis)
szt_gra = sunrise_utc + (4 * shaah_zmanis)
szks_gra_local = szks_gra.astimezone(local_tz)
szt_gra_local = szt_gra.astimezone(local_tz)
print(f"\n Shaah Zmanis (Gr\"a): {str(shaah_zmanis).split('.')[0]}")
print(f" S\"Z Krias Shema (Gr\"a): {szks_gra_local.strftime('%I:%M:%S %p %Z')}")
print(f" S\"Z Tfilah (Gr\"a): {szt_gra_local.strftime('%I:%M:%S %p %Z')}")
else:
print("\nCould not calculate all zmanim due to missing sunrise or sunset time.")
