• API
• FAQ
• Tools
• Archive
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
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.
Top