Hide keyboard shortcuts

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 

2 

3""" 

4Coordinate conversions for IceCube. 

5""" 

6 

7 

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) 

14 

15 

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. 

21 

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. 

31 

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. 

39 

40 Examples 

41 -------- 

42 Pull an archival neutrino and compare its official RA/dec to the RA dec 

43 conversion provided by this function. 

44 

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 

69 

70 

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). 

76 

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. 

85 

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. 

94 

95 Examples 

96 -------- 

97 Pull an archival neutrino and compare its official zenith/azimuth to the 

98 conversion provided by this function. 

99 

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