-
Notifications
You must be signed in to change notification settings - Fork 16
/
astro.py
209 lines (186 loc) · 5.36 KB
/
astro.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
# Licensed under a 3-clause BSD style license - see LICENSE.rst
# -*- coding: utf-8 -*-
"""This module corresponds to the goddard/astro directory in idlutils.
"""
from __future__ import print_function
def airtovac(air):
"""Convert air wavelengths to wavelengths in vacuum.
Parameters
----------
air : array-like
Values of wavelength in air in Angstroms.
Returns
-------
array-like
Values of wavelength in vacuum in Angstroms.
Notes
-----
* Formula from `P. E. Ciddor, Applied Optics, 35, 1566 (1996)
<http://adsabs.harvard.edu/abs/1996ApOpt..35.1566C>`_.
* Values of wavelength below 2000 Å are not converted.
"""
from numpy import zeros
from astropy.units import Angstrom
try:
u = air.unit
except AttributeError:
u = None
try:
t = air.dtype
except AttributeError:
# Most likely, air is simply a float.
t = None
if t is None:
if air < 2000.0:
return air
vacuum = air
a = air
g = None
else:
try:
a = air.to(Angstrom).value
except AttributeError:
a = air
g = a < 2000.0
if g.all():
return air
vacuum = zeros(air.shape, dtype=t) + a
for k in range(2):
sigma2 = (1.0e4/vacuum)**2
fact = (1.0 + 5.792105e-2/(238.0185 - sigma2) +
1.67917e-3/(57.362 - sigma2))
vacuum = a * fact
if g is not None:
vacuum[g] = a[g]
if u is not None:
vacuum = (vacuum * Angstrom).to(u)
return vacuum
def gcirc(ra1, dec1, ra2, dec2, units=2):
"""Computes rigorous great circle arc distances.
Parameters
----------
ra1, dec1, ra2, dec2 : :class:`float` or array-like
RA and Dec of two points.
units : { 0, 1, 2 }, optional
* units = 0: everything is already in radians
* units = 1: RA in hours, dec in degrees, distance in arcsec.
* units = 2: RA, dec in degrees, distance in arcsec (default)
Returns
-------
:class:`float` or array-like
The angular distance. Units of the value returned depend on the
input value of `units`.
Notes
-----
The formula below is the one best suited to handling small angular
separations. See:
http://en.wikipedia.org/wiki/Great-circle_distance
"""
from numpy import arcsin, cos, deg2rad, rad2deg, sin, sqrt
if units == 0:
rarad1 = ra1
dcrad1 = dec1
rarad2 = ra2
dcrad2 = dec2
elif units == 1:
rarad1 = deg2rad(15.0*ra1)
dcrad1 = deg2rad(dec1)
rarad2 = deg2rad(15.0*ra2)
dcrad2 = deg2rad(dec2)
elif units == 2:
rarad1 = deg2rad(ra1)
dcrad1 = deg2rad(dec1)
rarad2 = deg2rad(ra2)
dcrad2 = deg2rad(dec2)
else:
raise ValueError('units must be 0, 1 or 2!')
deldec2 = (dcrad2-dcrad1)/2.0
delra2 = (rarad2-rarad1)/2.0
sindis = sqrt(sin(deldec2)*sin(deldec2) +
cos(dcrad1)*cos(dcrad2)*sin(delra2)*sin(delra2))
dis = 2.0*arcsin(sindis)
if units == 0:
return dis
else:
return rad2deg(dis)*3600.0
def get_juldate(seconds=None):
"""Returns the current Julian date.
Uses the MJD trick & adds the offset to get JD.
Parameters
----------
seconds : :class:`int` or :class:`float`, optional
Time in seconds since the UNIX epoch. This should only be used
for testing.
Returns
-------
:class:`float`
The Julian Day number as a floating point number.
Notes
-----
Do not use this function if high precision is required, or if you are
concerned about the distinction between UTC & TAI.
"""
from time import time
if seconds is None:
t = time()
else:
t = seconds
mjd = t/86400.0 + 40587.0
return mjd + 2400000.5
def get_juldate_main(): # pragma: no cover
"""Entry point for the get_juldate command-line script.
"""
jd = get_juldate()
print(jd)
return 0
def vactoair(vacuum):
"""Convert vacuum wavelengths to wavelengths in air.
Parameters
----------
vacuum : array-like
Values of wavelength in vacuum in Angstroms.
Returns
-------
array-like
Values of wavelength in air in Angstroms.
Notes
-----
* Formula from `P. E. Ciddor, Applied Optics, 35, 1566 (1996)
<http://adsabs.harvard.edu/abs/1996ApOpt..35.1566C>`_.
* Values of wavelength below 2000 Å are not converted.
"""
from numpy import zeros
from astropy.units import Angstrom
try:
u = vacuum.unit
except AttributeError:
u = None
try:
t = vacuum.dtype
except AttributeError:
# Most likely, vacuum is simply a float.
t = None
if t is None:
if vacuum < 2000.0:
return vacuum
air = vacuum
v = vacuum
g = None
else:
try:
v = vacuum.to(Angstrom).value
except AttributeError:
v = vacuum
g = v < 2000.0
if g.all():
return vacuum
air = zeros(vacuum.shape, dtype=t) + v
sigma2 = (1.0e4/v)**2
fact = (1.0 + 5.792105e-2/(238.0185 - sigma2) +
1.67917e-3/(57.362 - sigma2))
air = v / fact
if g is not None:
air[g] = v[g]
if u is not None:
air = (air * Angstrom).to(u)
return air