Search code examples
pythonastropyastronomy

Using astropy to generate solar eclipse conditions based on my location


This is a question for the astronomy-minded folks on here.

I am an amateur astrophotographer looking to develop a personal script to aid my photography of next year's total solar eclipse. I am developing a Python script to automate my photography, so that I may enjoy the eclipse with my own eyes while my DSLR clicks away. Here's the script I've developed so far. The script uses digicamcontrol to control the camera.

Right now in the script I have just develop automation based on the partial phase of the eclipse (first contact, C1) and the timing of the eclipse in UTC (as well as my own PC). But a thought occurred to me: What if I can't connect to the internet and get the exact timing of the solar eclipse based on my location? I'd like to be able to generate those times. Is there a more efficient method to utilize astropy for this task? Thanks in advance.

import digiCamControlPython as dccp
import time
from datetime import datetime
from astropy.time import Time 

local_time = Time.now()

utc_time_now = local_time.utc

def PartialEclipse(start_time:str, end_time:str):
    camera = dccp.Camera()
    camera.setIso(100)
    camera.setShutterspeed("1/50")
    camera.setFolder(r"C:\Users\My_Name\Pictures\digiCamControl")
    
    # Set the target capture time in astropy time format
    partial_eclipse_start = Time(start_time, format='isot')
    partial_eclipse_end = Time(end_time, format='isot')

    # Wait until the capture time
    while utc_time_now < partial_eclipse_start:
        time.sleep(1)

    # Start capturing images
    while utc_time_now < partial_eclipse_end:
        camera.capture()
        time.sleep(30)  # Capture an image every 30 seconds
        
PartialEclipse("2024-04-08T17:12:13", "2024-04-08T18:29:24") #times of partial eclipse start and T-15s before totality

EDIT: In the event anyone ever looks at this question, I did make some progress on this.

import numpy as np
import astropy.units as u
from astropy.coordinates import solar_system_ephemeris, AltAz, EarthLocation, SkyCoord
from astropy.coordinates import get_body, get_moon, get_sun
from astropy.time import Time

myLocation = EarthLocation(lat=26*u.deg, lon=-80*u.deg, height=0*u.m)

# set the time step (how often to check for a solar eclipse in seconds)
time_step = 3600 # 1 hour

# set the number of days to check for a solar eclipse
num_days = 365

# set the start and end times to check for a solar eclipse
start_time = Time.now()
end_time = start_time + num_days * u.day

# initialize a list to store the times of a solar eclipse
eclipse_times = []

# loop over the desired time range, checking for a solar eclipse every time_step seconds
with solar_system_ephemeris.set('jpl'):
  for t in np.arange(start_time.unix, end_time.unix, time_step):
    time = Time(t, format='unix')
    moon = get_body('moon', time, myLocation)
    sun = get_body('sun', time, myLocation)
    sun_coord = SkyCoord(sun.ra, sun.dec, sun.distance, frame='icrs')
    moon_coord = SkyCoord(moon.ra, moon.dec, moon.distance, frame='icrs')
    
    # check if the angular separation between the moon and sun is close to zero
    angular_separation = moon_coord.separation(sun_coord)
    if angular_separation < 0.6 * u.deg: #elongation where the partial eclipse begins
      eclipse_times.append(time)

# print the times of the next solar eclipse
if len(eclipse_times) > 0:
  print("The next solar eclipse is at: ", eclipse_times[0].iso)
else:
  print("No solar eclipses found in the specified time range.")

Solution

  • I think your approach is really good! If you want to increase the accuracy of your start time prediction without using a lot more computational power, you can use scipy.optimize.root_scalar to refine the start time you found.

    In my solution below, I've defined a function called distance_contact() whose root represents the start of the eclipse. This function is zero if the Sun and Moon are barely touching, positive if they are separated, and negative if they are overlapping. Then I define a grid of times with a timestep of 1 hour similar to your code, and pass it into this function to search for eclipses. It then finds the first time where distance_contanct is negative and uses that time and the timestep before as the search space for scipy.optimize.root_scalar.

    Also, instead of using 0.6 * u.deg as the separation distance for an eclipse to occur, I've calculated the angular radius of the Sun and Moon for the time argument to distance_contact to make the prediction as accurate as possible.

    import numpy as np
    import scipy.optimize
    import astropy.units as u
    import astropy.time
    import astropy.constants
    import astropy.coordinates
    
    
    def distance_contact(
            location: astropy.coordinates.EarthLocation,
            time: astropy.time.Time,
            eclipse_type: str,
    ) -> u.Quantity:
    
        radius_sun = astropy.constants.R_sun
        radius_moon = 1737.4 * u.km
    
        coordinate_sun = astropy.coordinates.get_sun(time)
        coordinate_moon = astropy.coordinates.get_moon(time)
    
        frame_local = astropy.coordinates.AltAz(obstime=time, location=location)
    
        alt_az_sun = coordinate_sun.transform_to(frame_local)
        alt_az_moon = coordinate_moon.transform_to(frame_local)
    
        angular_radius_sun = np.arctan2(radius_sun, alt_az_sun.distance).to(u.deg)
        angular_radius_moon = np.arctan2(radius_moon, alt_az_moon.distance).to(u.deg)
    
        if eclipse_type == 'total':
            separation_max = angular_radius_moon - angular_radius_sun
        elif eclipse_type == 'partial':
            separation_max = angular_radius_moon + angular_radius_sun
        else:
            raise ValueError("Unknown eclipse type")
    
        return (alt_az_moon.separation(alt_az_sun).deg * u.deg) - separation_max
    
    
    def calc_time_start(
            location: astropy.coordinates.EarthLocation,
            time_search_start: astropy.time.Time,
            time_search_stop: astropy.time.Time,
            eclipse_type: str = 'partial'
    ) -> astropy.time.Time:
    
        astropy.coordinates.solar_system_ephemeris.set("de430")
    
        # If we're only looking for a partial eclipse, we can accept a coarser search grid
        if eclipse_type == "partial":
            step = 1 * u.hr
        elif eclipse_type == "total":
            step = 1 * u.min
        else:
            raise ValueError("Unknown eclipse type")
    
        # Define a grid of times to search for eclipses
        time = astropy.time.Time(np.arange(time_search_start, time_search_stop, step=step))
    
        # Find the times that are during an eclipse
        mask_eclipse = distance_contact(location=location, time=time, eclipse_type=eclipse_type) < 0
    
        # Find the index of the first time that an eclipse is occuring
        index_start = np.argmax(mask_eclipse)
    
        # Search around that time to find when the eclipse actually starts
        time_eclipse_start = scipy.optimize.root_scalar(
            f=lambda t: distance_contact(location, astropy.time.Time(t, format="unix"), eclipse_type=eclipse_type).value,
            bracket=[time[index_start - 1].unix, time[index_start].unix],
        ).root
        time_eclipse_start = astropy.time.Time(time_eclipse_start, format="unix")
    
        return time_eclipse_start
    
    
    def test_calc_time_start():
    
        location = astropy.coordinates.EarthLocation(lat=26 * u.deg, lon=-80 * u.deg, height=0 * u.m)
        eclipse_type = 'partial'
        time_start = calc_time_start(
            location=location,
            time_search_start=astropy.time.Time.now(),
            time_search_stop=astropy.time.Time.now() + 0.9 * u.yr,
            eclipse_type=eclipse_type,
        )
        print(time_start.isot)
    

    which outputs:

    2023-10-14T15:57:38.068