SHARE
TWEET

Untitled

a guest Jan 19th, 2020 82 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3.  
  4. from astroquery.vizier import Vizier
  5. from astroquery.simbad import Simbad
  6. from astropy.coordinates import SkyCoord
  7. import astropy.units as u
  8. from astropy.coordinates import AltAz
  9. from astropy.coordinates import EarthLocation
  10. from astropy.time import Time
  11. from astropy.utils.iers import conf
  12. conf.auto_max_age = None
  13.  
  14. %matplotlib inline
  15.  
  16.  
  17.  
  18. def unix_timestamp_from_iso_date(datetime):
  19.     return Time(datetime, format = 'isot').unix
  20.  
  21. def iso_date_from_unix_timestamp(datetime):
  22.     return Time(datetime, format = 'unix').isot
  23.  
  24. def radec_to_altaz(lat, lon, datetime, ra, dec):
  25.     """
  26.     Convert equatorial coordinates to local altitude and azimuth. The algorithm is approximate, not accounting
  27.     for the atmosphere, precession, nutation and other higher order effects. For XX and XXI century dates, the
  28.     expected error is 0.3 deg in alt and az * cos(alt). Maximum errors of 1.4 deg are possible for early XX or
  29.     late XXI century dates
  30.  
  31.     args:
  32.         lat       :    Latitude of the observatory in degrees (-90 <= lat <= 90)
  33.         lon       :    Longitude of the observatory in degrees
  34.         datetime  :    ISO date and time of observation (e.g. 2000-01-01T12:00:00Z for J2000)
  35.         ra        :    Right Ascension of the target in degrees, J2000
  36.         dec       :    Declination of the target in degrees (-90 <= dec <= 90), J2000
  37.  
  38.     returns:
  39.         alt       :    Altitude of the target in degrees
  40.         az        :    Azimuth of the target in degrees
  41.     """
  42.     # Days since J2000
  43.     unix_timestamp_requested = unix_timestamp_from_iso_date(datetime)
  44.     unix_timestamp_J2000 = unix_timestamp_from_iso_date('2000-01-01T12:00:00Z')
  45.     J2000 = (unix_timestamp_requested - unix_timestamp_J2000) / (60 * 60 * 24)
  46.  
  47.     # Check that the longitude complies with the format required by this approximation
  48.     lon = lon % 360
  49.     if lon > 180.0:
  50.         lon = lon - 360
  51.     # Siderial time and hour angle
  52.     LST = 100.46 + 0.985647 * J2000 + lon + 15 * (J2000 - 0.5 - np.floor(J2000)) * 24
  53.     LST = LST % 360
  54.     HA = LST - ra
  55.  
  56.     # Hour angle
  57.     dec, lat, HA = np.radians((dec, lat, HA))
  58.     alt = np.arcsin(np.sin(dec) * np.sin(lat) + np.cos(dec) * np.cos(lat) * np.cos(HA))
  59.     az = np.arccos((np.sin(dec) - np.sin(alt) * np.sin(lat)) / (np.cos(alt) * np.cos(lat)))
  60.     alt, az = np.degrees((alt, az))
  61.     if np.sin(HA) > 0:
  62.         az = 360 - az
  63.     return alt, az
  64.  
  65. def ascent_time(lat, lon, datetime_start, datetime_end, ra, dec, alt):
  66.     """
  67.     Determine the times when the target reaches a certain altitude within a window of dates. See the docstring
  68.     of radec_to_altaz() for a note on the accuracy
  69.  
  70.     args:
  71.         lat             :    Latitude of the observatory in degrees (-90 <= lat <= 90)
  72.         lon             :    Longitude of the observatory in degrees
  73.         datetime_start  :    Start date of the window (ISO)
  74.         datetime_end    :    End date of the window (ISO)
  75.         ra              :    Right Ascension of the target in degrees, J2000
  76.         dec             :    Declination of the target in degrees (-90 <= dec <= 90), J2000
  77.         alt             :    Desired altitude
  78.  
  79.     returns:
  80.         ascents         :    List of all times (as ISO dates) when the target is at the desired altitude.
  81.                              Empty list if the target never reaches the altitude within the given window
  82.     """
  83.     # Get the hour angle of the target
  84.     alt, dec, lat = np.radians((alt, dec, lat))
  85.     cos_HA = (np.sin(alt) - np.sin(dec) * np.sin(lat)) / (np.cos(dec) * np.cos(lat))
  86.     # The target never rises to the desired altitude
  87.     if np.abs(cos_HA) >= 1.0:
  88.         return []
  89.     # Now, we have the cosine of the HA. Of course, there will be infinitely many values of HA that
  90.     # have that cosine of two families: arccos(cos_HA) + n*360 and -arccos(cos_HA) + n*360
  91.  
  92.     # Days since J2000
  93.     unix_timestamp_end = unix_timestamp_from_iso_date(datetime_end)
  94.     unix_timestamp_start = unix_timestamp_from_iso_date(datetime_start)
  95.     unix_timestamp_J2000 = unix_timestamp_from_iso_date('2000-01-01T12:00:00Z')
  96.     J2000_end = (unix_timestamp_end - unix_timestamp_J2000) / (60 * 60 * 24)
  97.     J2000_start = (unix_timestamp_start - unix_timestamp_J2000) / (60 * 60 * 24)
  98.  
  99.     ascents = []     # Store found ascents here
  100.  
  101.     for phase in [1, -1]: # "phase" refers to the family of HA values considered. We search for ascents in both
  102.         # Check that the longitude complies with the format required by this approximation
  103.         lon = lon % 360
  104.         if lon > 180.0:
  105.             lon = lon - 360
  106.         # Get the HA value of the right family without n*360 for now
  107.         HA = np.degrees(np.arccos(cos_HA)) * phase
  108.         # Get the LST without n*360 as well
  109.         LST = HA + ra
  110.         # Try different values of n until LST matches the window
  111.         while (LST - 100.46 + 180 - lon - 360 * 1) / 0.985647 > J2000_start:
  112.             LST -= 360
  113.         while (LST - 100.46 + 180 - lon) / 0.985647 < J2000_start:
  114.             LST += 360
  115.         # Get the UT time of the ascent on the first day of the window as well as the rate
  116.         # at which it will drift per day
  117.         UT_start = (LST - 100.46 + 180 - lon - 0.985647 * np.floor(J2000_start)) / 360 / (1 + 0.985647 / 360)
  118.         UT_offset = -0.985647 / 360
  119.  
  120.         # Collect all ascents
  121.         i = 0
  122.         while True:
  123.             if i > 1000:
  124.                 raise ValueError()
  125.             ascent = np.floor(J2000_start) + i + UT_start + i * UT_offset
  126.             i += 1
  127.             if ascent < J2000_start:
  128.                 continue
  129.             if ascent > J2000_end:
  130.                 break
  131.             ascents += [iso_date_from_unix_timestamp(ascent * ((60 * 60 * 24)) + unix_timestamp_J2000)]
  132.     return ascents
  133.  
  134. def climax(lat, lon, datetime_start, datetime_end, ra, dec):
  135.     """
  136.     Determine the maximum altitude reached by the target and the time when it is reached within a window
  137.     of dates. See the docstring of radec_to_altaz() for a note on the accuracy
  138.  
  139.     args:
  140.         lat             :    Latitude of the observatory in degrees (-90 <= lat <= 90)
  141.         lon             :    Longitude of the observatory in degrees
  142.         datetime_start  :    Start date of the window (ISO)
  143.         datetime_end    :    End date of the window (ISO)
  144.         ra              :    Right Ascension of the target in degrees, J2000
  145.         dec             :    Declination of the target in degrees (-90 <= dec <= 90), J2000
  146.  
  147.     returns:
  148.         alt             :    Maximum altitude within the window in degrees
  149.         datetime        :    ISO date of when the altitude is reached. If reached multiple times, only the
  150.                              first instance is returned
  151.     """
  152.     # Find the minimum and maximum altitudes attainable by the target
  153.     extrema = [-lat - dec - 90, dec - lat + 90]
  154.     for i, extremum in enumerate(extrema):
  155.         if extremum > 90 or extremum < -90:
  156.             extrema[i] = (180 - extremum)
  157.         if extremum < -90:
  158.             extrema[i] = extrema[i] - 360
  159.  
  160.     # Get all times when the object is at its highest
  161.     ascents = ascent_time(lat, lon, datetime_start, datetime_end, ra, dec, np.max(extrema))
  162.     # If the object does reach its highest possible altitude, it will be the desired climax
  163.     if len(ascents) != 0:
  164.         return np.max(extrema), ascents[0]
  165.     else:
  166.     # If not, either the start or the end of the window will be...
  167.         start_alt = radec_to_altaz(lat, lon, datetime_start, ra, dec)[0]
  168.         end_alt = radec_to_altaz(lat, lon, datetime_end, ra, dec)[0]
  169.         if start_alt > end_alt:
  170.             return start_alt, datetime_start
  171.         else:
  172.             return end_alt, datetime_end
  173.  
  174. def get_sun(datetime):
  175.     """
  176.     Estimate J2000 equatorial coordinates of the Sun at any given time. Use the approximate algorithm
  177.     from stargazing.net/kepler/sun.html
  178.  
  179.     args:
  180.         datetime        :    Date of interest (ISO)
  181.  
  182.     returns:
  183.         ra              :    Right ascension of the Sun in degrees
  184.         dec             :    Declination of the Sun in degrees
  185.     """
  186.     # Days since J2000
  187.     unix_timestamp_requested = unix_timestamp_from_iso_date(datetime)
  188.     unix_timestamp_J2000 = unix_timestamp_from_iso_date('2000-01-01T12:00:00Z')
  189.     J2000 = (unix_timestamp_requested - unix_timestamp_J2000) / (60 * 60 * 24)
  190.  
  191.     # Just following the algorithm from now on...
  192.     L = (280.461 + 0.9856474 * J2000) % 360
  193.     g = np.radians(357.528 + 0.9856003 * J2000)
  194.     lmbda = np.radians(L + 1.915 * np.sin(g) + 0.020 * np.sin(2 * g))
  195.     epsilon = np.radians(23.439 - 0.0000004 * J2000)
  196.     Y = np.cos(epsilon) * np.sin(lmbda)
  197.     X = np.cos(lmbda)
  198.  
  199.     ra = np.degrees(np.arctan(Y / X))
  200.     if X < 0:
  201.         ra = ra + 180
  202.     elif Y < 0 and X > 0:
  203.         ra = ra + 360
  204.     dec = np.degrees(np.arcsin(np.sin(epsilon) * np.sin(lmbda)))
  205.     return ra, dec
  206.  
  207. def twilight(lat, lon, datetime, alt, niter = 5):
  208.     """
  209.     Iteratively determine sunrise/sunset/twilight times on a given date. ascent_time() can determine when
  210.     a given pair of equatorial coordinates reaches a given altitude. However, in order to know the coordinates
  211.     of the Sun when its transiting said altitude, one must know the time when it happens and feed it into
  212.     get_sun(). This recursive loop can be solved iteratively.
  213.  
  214.     Sunset/sunrise/twilight events are searched within 24 hours of the given time. On the 0th iteration, the
  215.     average coordinates of the Sun are assumed for the search period. Once the events are found, the coordinates
  216.     of the Sun are updated and the search is repeated. Starting with the 1st iteration, each event is searched
  217.     separately within 2 hours of the event time known from the previous iteration.
  218.  
  219.     See the docstring of radec_to_altaz() for a note on the accuracy. Ultimately, a 1 degree error in the RA/Dec
  220.     to AltAz conversion algorithm can result in a 4 minute error in sunset/sunrise times in the worst case
  221.     scenario
  222.  
  223.     args:
  224.         lat             :    Latitude of the observatory in degrees (-90 <= lat <= 90)
  225.         lon             :    Longitude of the observatory in degrees
  226.         datetime        :    Date of interest (ISO). Events will be searched within 24 hours of this date
  227.                              (12 hours in each direction)
  228.         alt             :    Altitude of the event. "0" for sunset/sunrise, "-6" for civil twlight etc...
  229.         n_iter          :    Number of iterations. Typically, takes only 3 or so to converge to seconds
  230.  
  231.     returns:
  232.         sunrise         :    Date of the first event (e.g. sunrise) in ISO
  233.         sunset          :    Date of the second event (e.g. sunset) in ISO
  234.     """
  235.     # Get ISO dates defining the initial search window
  236.     datetime_unix = unix_timestamp_from_iso_date(datetime)
  237.     datetime_start = iso_date_from_unix_timestamp(datetime_unix - 12 * 60 * 60)
  238.     datetime_end = iso_date_from_unix_timestamp(datetime_unix + 12 * 60 * 60)
  239.  
  240.     # Look for events. We must find 2!
  241.     ascents = ascent_time(lat, lon, datetime_start, datetime_end, *get_sun(datetime), alt)
  242.     if len(ascents) != 2:
  243.         raise ValueError('Expected two ascension events within the window, {} found'.format(len(ascents)))
  244.  
  245.     for i, ascent in enumerate(ascents):
  246.         # Iterate for each event
  247.         for j in range(niter):
  248.             # Define the new search window (narrower on all subsequent iterations)
  249.             ascent_unix = unix_timestamp_from_iso_date(ascent)
  250.             datetime_start = iso_date_from_unix_timestamp(ascent_unix - 60 * 60)
  251.             datetime_end = iso_date_from_unix_timestamp(ascent_unix + 60 * 60)
  252.             # Search again
  253.             updated_ascents = ascent_time(lat, lon, datetime_start, datetime_end, *get_sun(ascent), alt)
  254.             # Due to the narrow window, we must only find one even this time
  255.             if len(updated_ascents) != 1:
  256.                 raise ValueError('Ascension events too frequent (more than one in two hours)'.format(len(ascents)))
  257.             ascent = updated_ascents[0]
  258.         ascents[i] = ascent
  259.     return ascents
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Top