vlsr.py

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


def vlsr(t, loc, psrc, verbose=False):
    """Compute the line of sight radial velocity

    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
    """
    # Radial velocity correction to solar system barycenter
    vsrc = psrc.radial_velocity_correction(obstime=t, location=loc)

    # Projection of solar velocity towards the source
    vsun_proj = psrc.cartesian.dot(psun.cartesian) * vsun

    if verbose:
        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)))

    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.

    Args:
        psrc (SkyCoord): sky location for correction
        t (float): time for correction
        rest_frequency (Union[Quantity, float]): rest frequency (frequency which was transmitted)
        loc (EarthLocation): observer location

    Returns:
        Observable frequency in Hz as astropy Quantity
    """
    v1 = vlsr(t, loc, psrc, verbose=verbose)

    return (1 - v1 / astropy.constants.c) * rest_frequency


def test_vlsr():
    freq_hi = 1420.405751

    # Test 1: Source towards direction of solar motion, Earth moving prependicular
    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)))
    f = doppler_frequency(psrc, t, freq_hi, loc)
    print("Doppler HI: {}".format(f))

    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)
    print(cmd)
    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")

    # Test 2: Source perpendicular to direction of solar motion, Earth moving towards source
    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)))
    f = doppler_frequency(psrc, t, freq_hi, loc)
    print("Doppler HI: {}".format(f))
    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)
    print(cmd)
    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")

    # Test 3: Narrabri calculator
    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)))