vlsr.py 3.79 KB
Newer Older
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
1
2
3
4
5
6
7
8
9
10
from astropy.time import Time
from astropy.coordinates import SkyCoord, EarthLocation
import astropy.units as u
from astropy.units import Quantity
import astropy.constants

import numpy as np

# Direction of motion of the Sun. Info from
# http://herschel.esac.esa.int/hcss-doc-15.0/load/hifi_um/html/hifi_ref_frame.html#hum_lsr_frame
11
12
psun = SkyCoord(ra="18:03:50.29", dec="+30:00:16.8", frame="icrs", unit=(u.hourangle, u.deg))
vsun = -20.0 * u.km / u.s
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
13

14
15

def vlsr(t, loc, psrc, verbose=False):
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
16
17
    """Compute the line of sight radial velocity

18
19
20
21
22
23
24
    Args:
        psrc: SkyCoord object or source
        loc: EarthLocation object of observer
        t: Time object

    Returns:
        Line of sight radial velocity (in km/s) as an astropy Quantity
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
25
26
    """
    # Radial velocity correction to solar system barycenter
27
    vsrc = psrc.radial_velocity_correction(obstime=t, location=loc)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
28
29

    # Projection of solar velocity towards the source
30
    vsun_proj = psrc.cartesian.dot(psun.cartesian) * vsun
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
31
32

    if verbose:
33
34
        print("Barycentric radial velocity: {0:+8.3f}".format(vsrc.to(u.km / u.s)))
        print("Projected solar velocity:    {0:+8.3f}".format(vsun_proj.to(u.km / u.s)))
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
35

36
37
38
39
40
    return vsun_proj - vsrc


def doppler_frequency(psrc, t, rest_frequency, loc, verbose=False):
    """Doppler corrected frequency, taking into account the line of sight radial velocity.
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
41
42

    Args:
43
44
        psrc (SkyCoord): sky location for correction
        t (float): time for correction
45
        rest_frequency (Union[Quantity, float]): rest frequency (frequency which was transmitted)
46
        loc (EarthLocation): observer location
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
47
48

    Returns:
49
        Observable frequency in Hz as astropy Quantity
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
50
    """
51
    v1 = vlsr(t, loc, psrc, verbose=verbose)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
52

53
    return (1 - v1 / astropy.constants.c) * rest_frequency
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
54

55
56

def test_vlsr():
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
57
58
59
    freq_hi = 1420.405751

    # Test 1: Source towards direction of solar motion, Earth moving prependicular
60
61
62
63
64
    loc = EarthLocation.from_geodetic(lat=52 * u.deg, lon=6 * u.deg, height=0 * u.m)
    psrc = SkyCoord(ra=18 * u.hourangle, dec=30.0 * u.deg, frame="icrs")
    t = Time("2018-06-21T12:00:00", scale="utc", format="isot")
    v = vlsr(t, loc, psrc, verbose=True)
    print("LSR velocity:                {0:+8.3f}".format(v.to(u.km / u.s)))
65
    f = doppler_frequency(psrc, t, freq_hi, loc)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
66
67
    print("Doppler HI: {}".format(f))

68
69
    cmd = "~harm/bin/vlsr %f %f 2000 %s %f %f 0" %\
          (psrc.ra.rad, psrc.dec.rad, t.unix, loc.lat.rad, loc.lon.rad)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
70
    print(cmd)
71
72
73
    cmd = "~harm/bin/doppler_bl %f %f 2000 %s %f %f 0 %f" %\
          (psrc.ra.rad, psrc.dec.rad, t.unix, loc.lat.rad, loc.lon.rad, freq_hi)
    print(cmd + "\n\n")
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
74
75

    # Test 2: Source perpendicular to direction of solar motion, Earth moving towards source
76
77
78
79
80
    loc = EarthLocation.from_geodetic(lat=52 * u.deg, lon=6 * u.deg, height=0 * u.m)
    psrc = SkyCoord(ra=6 * u.hourangle, dec=0.0 * u.deg, frame="icrs")
    t = Time("2018-06-21T12:00:00", scale="utc", format="isot")
    v = vlsr(t, loc, psrc, verbose=True)
    print("LSR velocity:                {0:+8.3f}".format(v.to(u.km / u.s)))
81
    f = doppler_frequency(psrc, t, freq_hi, loc)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
82
    print("Doppler HI: {}".format(f))
83
84
    cmd = "~harm/bin/vlsr %f %f 2000 %s %f %f 0" %\
          (psrc.ra.rad, psrc.dec.rad, t.unix, loc.lat.rad, loc.lon.rad)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
85
    print(cmd)
86
87
88
89
    cmd = "~harm/bin/doppler_bl %f %f 2000 %s %f %f 0 %f" %\
          (psrc.ra.rad, psrc.dec.rad, t.unix, loc.lat.rad, loc.lon.rad, freq_hi)
    print(cmd + "\n\n")

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
90
    # Test 3: Narrabri calculator
91
92
93
94
95
    loc = EarthLocation.from_geodetic(lat=-30.3128 * u.deg, lon=149.5502 * u.deg, height=0 * u.m)
    psrc = SkyCoord(ra=12 * u.hourangle, dec=0.0 * u.deg, frame="icrs")
    t = Time("2018-04-16T12:00:00", scale="utc", format="isot")
    v = vlsr(t, loc, psrc, verbose=True)
    print("LSR velocity:                {0:+8.3f}".format(v.to(u.km / u.s)))