Commit 0c8ccc6f authored by Tammo Jan Dijkema's avatar Tammo Jan Dijkema
Browse files

Add vlsr.py from Cees

parent 5f47f9f1
#!/usr/bin/env python3
from astropy.time import Time
from astropy.coordinates import SkyCoord, EarthLocation
import astropy.units as u
from astropy.units import Quantity
import astropy.constants
from typing import Union
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
# vlsr routine
def vlsr(t,loc,psrc,verbose=False):
"""Compute the line of sight radial velocity
psrc: SkyCoord object or source
loc: EarthLocation object of observer
t: Time object
"""
# 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(t: Time, loc: EarthLocation, psrc: SkyCoord,
rest_frequency: Union[Quantity, float], verbose=False):
"""
Compute the Doppler corrected frequency, taking into account the line of sight radial velocity.
Args:
psrc: sky location for correction
loc: observer location
t: time for correction
frequency: observed frequency in LSR
Returns:
Quantity: Observable frequency
"""
v1 = vlsr(t, loc, psrc)
beta = v1/astropy.constants.c
return np.sqrt((1 + beta)/(1 - beta)) * rest_frequency
if __name__=="__main__":
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(t, loc, psrc, freq_hi)
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(t, loc, psrc, freq_hi)
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)))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment