Coverage for /opt/anaconda/lib/python3.7/site-packages/llama/files/i3/utils.py : 100%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# (c) Stefan Countryman 2019
3"""
4Coordinate conversions for IceCube.
5"""
8def icecube_coords():
9 """Get the ``astropy.coordinates.EarthLocation`` of the IceCube
10 detector."""
11 import astropy.units as u
12 from astropy.coordinates import EarthLocation
13 return EarthLocation(lat=-89.9944*u.deg, lon=-62.6081*u.deg, height=0*u.m)
16def zen_az2ra_dec(mjd, zen_deg, az_deg):
17 """Take the zenith and azimuth of the neutrino (these are earth-based
18 coordinates measured from the ICECUBE site at the south pole) and convert
19 them into the right ascension/declination equatorial coordinates used by
20 LIGO and others.
22 Parameters
23 ----------
24 mjd : float or array-like
25 The Modified Julian Date at which to make the conversion.
26 zen_deg : float or array-like
27 The zenith angle of the point on the sky as measured from IceCube's
28 location in degrees.
29 az_deg : float or array-like
30 Same, but for the azimuthal angle.
32 Returns
33 -------
34 ra_deg : float or array-like
35 The right-ascension in degrees of the input point or points. Will be a
36 scalar or vector value based on the input.
37 dec_deg : float or array-like
38 Same, but for the declination.
40 Examples
41 --------
42 Pull an archival neutrino and compare its official RA/dec to the RA dec
43 conversion provided by this function.
45 >>> import numpy as np
46 >>> from llama.files.i3.json import get_archive
47 >>> ns = get_archive(58392, 58393)
48 >>> mjds = np.array([n['mjd'] for n in ns])
49 >>> ras = np.array([n['ra'] for n in ns])
50 >>> decs = np.array([n['dec'] for n in ns])
51 >>> azs = np.array([n['azimuth'] for n in ns])
52 >>> zens = np.array([n['zenith'] for n in ns])
53 >>> cras, cdecs = zen_az2ra_dec(mjds, zens, azs)
54 >>> max(abs(cras - ras)) < 0.001
55 True
56 >>> max(abs(cdecs - decs)) < 0.001
57 True
58 """
59 import astropy.units as u
60 from astropy.time import Time
61 from astropy.coordinates import SkyCoord, AltAz, ICRS
62 i3coord = icecube_coords()
63 times = Time(mjd, format='mjd', scale='ut1')
64 alt = (90 - zen_deg)*u.degree
65 ap_az = (90-az_deg)*u.degree - i3coord.lon
66 radec = AltAz(alt=alt, az=ap_az, obstime=times, location=i3coord,
67 pressure=0).transform_to(ICRS)
68 return radec.ra.deg, radec.dec.deg
71def ra_dec2zen_az(mjd, ra_deg, dec_deg):
72 """Take the right ascension/declination equatorial coordinates (used by
73 LIGO and others) of the neutrino and convert them into the zenith and
74 azimuth measured from IceCube (the earth-based coordinates measured from
75 the ICECUBE site at the south pole).
77 Parameters
78 ----------
79 mjd : float or array-like
80 The Modified Julian Date at which to make the conversion.
81 ra_deg : float or array-like
82 The right-ascension in degrees of the point or points on the sky.
83 dec_deg : float or array-like
84 Same, but for the declination.
86 Returns
87 -------
88 zen_deg : float or array-like
89 The zenith angle of the point on the sky as measured from IceCube's
90 location in degrees. Will be a scalar or vector value based on the
91 input.
92 az_deg : float or array-like
93 Same, but for the azimuthal angle.
95 Examples
96 --------
97 Pull an archival neutrino and compare its official zenith/azimuth to the
98 conversion provided by this function.
100 >>> import numpy as np
101 >>> from llama.files.i3.json import get_archive
102 >>> ns = get_archive(58392, 58393)
103 >>> mjds = np.array([n['mjd'] for n in ns])
104 >>> ras = np.array([n['ra'] for n in ns])
105 >>> decs = np.array([n['dec'] for n in ns])
106 >>> azs = np.array([n['azimuth'] for n in ns])
107 >>> zens = np.array([n['zenith'] for n in ns])
108 >>> czens, cazs = ra_dec2zen_az(mjds, ras, decs)
109 >>> max(abs(czens - zens)) < 0.001
110 True
111 >>> max(abs(cazs - azs)) < 0.001
112 True
113 """
114 # https://astropy4cambridge.readthedocs.io/en/latest/_static/Astropy%20-%20Celestial%20Coordinates.html
115 # https://code.icecube.wisc.edu/projects/icecube/browser/IceCube/meta-projects/combo/trunk/astro/resources/examples/compare_with_astropy.py
116 import astropy.units as u
117 from astropy.time import Time
118 from astropy.coordinates import SkyCoord, AltAz
119 coords = SkyCoord(ra=ra_deg*u.degree, dec=dec_deg*u.degree, frame='icrs')
120 times = Time(mjd, format='mjd', scale='ut1')
121 i3coord = icecube_coords()
122 altaz = coords.transform_to(AltAz(obstime=times, location=i3coord,
123 pressure=0))
124 zen_deg = (90*u.degree - altaz.alt).degree
125 az_deg = (90*u.degree-altaz.az-i3coord.lon).wrap_at(360*u.degree).degree
126 return zen_deg, az_deg