-
Notifications
You must be signed in to change notification settings - Fork 0
/
velo.py
62 lines (44 loc) · 1.98 KB
/
velo.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
# Filename: time.py
# Aim: to convert utc to mjd or lst at location of the FAST site
import numpy as np
from astropy import units as u
from astropy.time import Time
from astropy.io import fits
#from astropy.coordinates import *
from astropy.coordinates import SkyCoord
from .location import get_fast_location
# Schonrich et al.(2010) Barycentric velocity relative to the LSR, which is defined in
# Galactic (right-handed) cartesian velocity components (U,V,W) = (11.1,12.24,7.25) km/s
# v_lsr = v_bary+U*cos(l)cos(b)+V*sin(l)*cos(b)+W*sin(b) # l,b are the Galactic coordinates.
def getVframeLSR(ra,dec,mjd):
fast_loc = get_fast_location()
obstime = Time(mjd,format='mjd')
sc = SkyCoord(ra=ra*u.degree,dec=dec*u.degree,frame='icrs')
b=sc.galactic.b.rad
l=sc.galactic.l.rad
barycorr = sc.radial_velocity_correction(obstime=obstime,location=fast_loc)
vlsr=barycorr.value/1000.0+11.1*np.cos(l)*np.cos(b)+12.24*np.sin(l)*np.cos(b)+7.25*np.sin(b)
return vlsr # in km/s
def vframe2fits(hdul, ra, dec, mjd, update=True):
vframe = getVframeLSR(ra,dec,mjd)
if update:
if 'VFRAME' in hdul[1].columns.names:
hdul[1].data['VFRAME'] = vframe
else:
new_col = fits.Column(name='VFRAME',format='1D',array=vframe)
newhdu = fits.BinTableHDU.from_columns(hdul[1].columns + fits.ColDefs([new_col]))
hdul[1].data = newhdu.data
hdul.flush()
return vframe
def freq2velo(freq,line_rest_freq): # freq in MHz
rest_freq = line_rest_freq * u.MHz
frequency = freq * u.MHz
relativistic_equiv = u.doppler_relativistic(rest_freq)
velo = frequency.to(u.km/u.s,equivalencies=relativistic_equiv)
return velo.value # km/s
def velo2freq(velo,line_rest_freq): #velo in km/s; freq in MHz
rest_freq = line_rest_freq * u.MHz
velocity = velo * u.km/u.s
relativistic_equiv = u.doppler_relativistic(rest_freq)
freq = velocity.to(u.MHz,equivalencies=relativistic_equiv)
return freq.value # freq in MHz