Commit e6bcc19f authored by Tammo Jan Dijkema's avatar Tammo Jan Dijkema
Browse files

Reformat, checked with pycodestyle

parent 312029bc
...@@ -15,18 +15,27 @@ class CamrasHpibDevice(object): ...@@ -15,18 +15,27 @@ class CamrasHpibDevice(object):
if not CamrasHpibDevice.command_thread: if not CamrasHpibDevice.command_thread:
serial_port = serial.Serial("/dev/ttyUSB0") serial_port = serial.Serial("/dev/ttyUSB0")
CamrasHpibDevice.command_thread = serial.threaded.ReaderThread(serial_port, hpib.GPIBProtocol) CamrasHpibDevice.command_thread = serial.threaded.ReaderThread(serial_port,
hpib.GPIBProtocol)
CamrasHpibDevice.command_thread.start() CamrasHpibDevice.command_thread.start()
CamrasHpibDevice.command_thread._connection_made.wait() CamrasHpibDevice.command_thread._connection_made.wait()
CamrasHpibDevice.command_thread.protocol.init_hpib() CamrasHpibDevice.command_thread.protocol.init_hpib()
@property @property
def frequency(self): def frequency(self):
"""Returns the frequency of the device as an astropy Quantity. """Frequency of the device as an astropy Quantity.
Throws a RuntimeError if no response"""
freq_str = CamrasHpibDevice.command_thread.protocol.command("freq?", address=self.hpib_address) Returns:
if not freq_str or len(freq_str)==0: Quantity: frequency (in Hz) of the device
raise RuntimeError("Camras device at address {} is not responding".format(self.hpib_address))
Raises:
RuntimeError: If the device does not respond
"""
freq_str = CamrasHpibDevice.command_thread.protocol.command("freq?",
address=self.hpib_address)
if not freq_str or len(freq_str) == 0:
raise RuntimeError("Camras device at address {} is not responding".format(
self.hpib_address))
return int(float(freq_str)) * u.Hz return int(float(freq_str)) * u.Hz
@frequency.setter @frequency.setter
...@@ -35,14 +44,20 @@ class CamrasHpibDevice(object): ...@@ -35,14 +44,20 @@ class CamrasHpibDevice(object):
Args: Args:
freq (Union[float, Quantity]): new frequency. If no unit present, assume Hz freq (Union[float, Quantity]): new frequency. If no unit present, assume Hz
Raises:
RuntimeError: If the device sticks to a different frequency (e.g. if the specified
frequency is outside the accepted range for the device)
""" """
if isinstance(freq, Quantity): if isinstance(freq, Quantity):
freq = freq.to(u.Hz).value freq = freq.to(u.Hz).value
CamrasHpibDevice.command_thread.protocol.command("freq {:d} Hz".format(int(freq)), address=self.hpib_address) CamrasHpibDevice.command_thread.protocol.command("freq {:d} Hz".format(int(freq)),
address=self.hpib_address)
new_freq = self.frequency new_freq = self.frequency
if new_freq.to(u.Hz).value != int(freq): if new_freq.to(u.Hz).value != int(freq):
raise RuntimeError("Setting frequency failed: tried to set to {}, it is now {}".format( raise RuntimeError("Setting frequency failed: tried to set to {}, it is now {}".format(
freq, new_freq)) freq, new_freq))
class Receiver(CamrasHpibDevice): class Receiver(CamrasHpibDevice):
"""Wrapper around HPIB commands for the Rohde & Schwartz receiver""" """Wrapper around HPIB commands for the Rohde & Schwartz receiver"""
...@@ -50,12 +65,14 @@ class Receiver(CamrasHpibDevice): ...@@ -50,12 +65,14 @@ class Receiver(CamrasHpibDevice):
def __init__(self, **kwargs): def __init__(self, **kwargs):
super(Receiver, self).__init__(1) super(Receiver, self).__init__(1)
class LocalOscillator(CamrasHpibDevice): class LocalOscillator(CamrasHpibDevice):
"""Wrapper around HPIB commands for the local oscillator""" """Wrapper around HPIB commands for the local oscillator"""
def __init__(self, **kwargs): def __init__(self, **kwargs):
super(LocalOscillator, self).__init__(28) super(LocalOscillator, self).__init__(28)
class ClockGenerator(CamrasHpibDevice): class ClockGenerator(CamrasHpibDevice):
"""Wrapper around HPIB commands for the clock generator (should be at 140MHz)""" """Wrapper around HPIB commands for the clock generator (should be at 140MHz)"""
......
from typing import Any, Union
import astropy.units as u import astropy.units as u
from astropy.coordinates import SkyCoord, EarthLocation from astropy.coordinates import EarthLocation
from astropy.time import Time from astropy.time import Time
import numpy as np
import argparse
import threading import threading
import subprocess import subprocess
from vlsr import doppler_frequency from vlsr import doppler_frequency
from telescope import telescope from camrasdevices import LocalOscillator
from camrasdevices import Receiver, LocalOscillator
freq_hi = 1420.405751 * u.MHz freq_hi = 1420.405751 * u.MHz
def doppler_harm(sky_coordinate, time, tracking_frequency, location, variant): def doppler_harm(sky_coordinate, time, tracking_frequency, location, variant):
ra = sky_coordinate.ra.to(u.rad).value ra = sky_coordinate.ra.to(u.rad).value
dec = sky_coordinate.dec.to(u.rad).value dec = sky_coordinate.dec.to(u.rad).value
t = time.unix t = time.unix
dt_lat = location.lat.to(u.rad).value dt_lat = location.lat.to(u.rad).value
dt_lon = location.lon.to(u.rad).value dt_lon = location.lon.to(u.rad).value
dt_alt = location.height.to(u.m).value dt_alt = location.height.to(u.m).value
argstring = "{ra} {dec} 2000 {time} {dt_lat} {dt_lon} {dt_alt}".format( argstring = "{ra} {dec} 2000 {time} {dt_lat} {dt_lon} {dt_alt}".format(
ra=ra, dec=dec, time=t, dt_lat=dt_lat, dt_lon=dt_lon, dt_alt=dt_alt) ra=ra, dec=dec, time=t, dt_lat=dt_lat, dt_lon=dt_lon, dt_alt=dt_alt)
doppler_executable = "/home/harm/bin/doppler_"+variant doppler_executable = "/home/harm/bin/doppler_" + variant
doppler_cmd = doppler_executable + " " + argstring + " " + str(tracking_frequency.to(u.MHz).value) doppler_cmd = doppler_executable + " " + argstring \
freq_doppler = float(subprocess.Popen(doppler_cmd, stdout=subprocess.PIPE, shell=True).stdout.read()) + " " + str(tracking_frequency.to(u.MHz).value)
freq_doppler = float(subprocess.Popen(doppler_cmd,
stdout=subprocess.PIPE,
shell=True
).stdout.read())
return freq_doppler * u.MHz return freq_doppler * u.MHz
def doppler_harm_bl(sky_coordinate, time, tracking_frequency, location): def doppler_harm_bl(sky_coordinate, time, tracking_frequency, location):
return doppler_harm(sky_coordinate, time, tracking_frequency, location, 'bl') return doppler_harm(sky_coordinate, time, tracking_frequency, location, 'bl')
def doppler_harm_mb(sky_coordinate, time, tracking_frequency, location): def doppler_harm_mb(sky_coordinate, time, tracking_frequency, location):
return doppler_harm(sky_coordinate, time, tracking_frequency, location, 'mb') return doppler_harm(sky_coordinate, time, tracking_frequency, location, 'mb')
def track_doppler(dt, lo, def track_doppler(dt, lo,
tracking_frequency=freq_hi, tracking_frequency=freq_hi,
doppler_function=doppler_frequency, doppler_function=doppler_frequency,
...@@ -44,7 +49,7 @@ def track_doppler(dt, lo, ...@@ -44,7 +49,7 @@ def track_doppler(dt, lo,
timeout=0.3): timeout=0.3):
"""Sets the Local Oscillator to a frequency to correct for Doppler shift. When exit event is """Sets the Local Oscillator to a frequency to correct for Doppler shift. When exit event is
set, return the LO frequency to 1GHz. set, return the LO frequency to 1GHz.
Args: Args:
dt (Telescope): a Telescope instance to read the current pointing from. dt (Telescope): a Telescope instance to read the current pointing from.
lo (LocalOscillator): the local oscillator of which the frequency will be set. lo (LocalOscillator): the local oscillator of which the frequency will be set.
...@@ -59,14 +64,14 @@ def track_doppler(dt, lo, ...@@ -59,14 +64,14 @@ def track_doppler(dt, lo,
dt_lon = 6.396169 * u.deg dt_lon = 6.396169 * u.deg
dt_alt = 25.0 * u.m dt_alt = 25.0 * u.m
dt_loc = EarthLocation.from_geodetic(lat = dt_lat, lon = dt_lon, height = dt_alt) dt_loc = EarthLocation.from_geodetic(lat=dt_lat, lon=dt_lon, height=dt_alt)
while not exit_event.is_set(): while not exit_event.is_set():
sky_coordinate = dt.radec sky_coordinate = dt.radec
if sky_coordinate: if sky_coordinate:
freq_doppler = doppler_function(sky_coordinate, Time.now(), tracking_frequency, dt_loc) freq_doppler = doppler_function(sky_coordinate, Time.now(), tracking_frequency, dt_loc)
dfreq = 1*u.GHz + tracking_frequency - freq_doppler dfreq = 1 * u.GHz + tracking_frequency - freq_doppler
lo.frequency = dfreq lo.frequency = dfreq
exit_event.wait(timeout=timeout) exit_event.wait(timeout=timeout)
lo.frequency = 1*u.GHz lo.frequency = 1 * u.GHz
exit_event.clear() exit_event.clear()
#!/usr/bin/env python3
from astropy.time import Time from astropy.time import Time
from astropy.coordinates import SkyCoord, EarthLocation from astropy.coordinates import SkyCoord, EarthLocation
import astropy.units as u import astropy.units as u
...@@ -9,32 +8,36 @@ import numpy as np ...@@ -9,32 +8,36 @@ import numpy as np
# Direction of motion of the Sun. Info from # 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 # 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)) 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 vsun = -20.0 * u.km / u.s
# vlsr routine
def vlsr(t,loc,psrc,verbose=False): def vlsr(t, loc, psrc, verbose=False):
"""Compute the line of sight radial velocity """Compute the line of sight radial velocity
psrc: SkyCoord object or source Args:
loc: EarthLocation object of observer psrc: SkyCoord object or source
t: Time object 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 # Radial velocity correction to solar system barycenter
vsrc = psrc.radial_velocity_correction(obstime = t, location = loc) vsrc = psrc.radial_velocity_correction(obstime=t, location=loc)
# Projection of solar velocity towards the source # Projection of solar velocity towards the source
vsun_proj = psrc.cartesian.dot(psun.cartesian)*vsun vsun_proj = psrc.cartesian.dot(psun.cartesian) * vsun
if verbose: if verbose:
print("Barycentric radial velocity: {0:+8.3f}".format(vsrc.to(u.km/u.s))) 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))) 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): return vsun_proj - vsrc
"""
Compute the Doppler corrected frequency, taking into account the line of sight radial velocity.
def doppler_frequency(psrc, t, rest_frequency, loc, verbose=False):
"""Doppler corrected frequency, taking into account the line of sight radial velocity.
Args: Args:
psrc (SkyCoord): sky location for correction psrc (SkyCoord): sky location for correction
...@@ -43,46 +46,51 @@ def doppler_frequency(psrc, t, rest_frequency, loc, verbose=False): ...@@ -43,46 +46,51 @@ def doppler_frequency(psrc, t, rest_frequency, loc, verbose=False):
loc (EarthLocation): observer location loc (EarthLocation): observer location
Returns: Returns:
Quantity: Observable frequency Observable frequency in Hz as astropy Quantity
""" """
v1 = vlsr(t, loc, psrc, verbose=verbose) v1 = vlsr(t, loc, psrc, verbose=verbose)
beta = v1/astropy.constants.c beta = v1 / astropy.constants.c
return np.sqrt((1 + beta)/(1 - beta)) * rest_frequency return np.sqrt((1 + beta) / (1 - beta)) * rest_frequency
if __name__=="__main__":
def test_vlsr():
freq_hi = 1420.405751 freq_hi = 1420.405751
# Test 1: Source towards direction of solar motion, Earth moving prependicular # 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) 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") psrc = SkyCoord(ra=18 * u.hourangle, dec=30.0 * u.deg, frame="icrs")
t = Time("2018-06-21T12:00:00",scale = "utc",format = "isot") t = Time("2018-06-21T12:00:00", scale="utc", format="isot")
v = vlsr(t,loc,psrc,verbose=True) v = vlsr(t, loc, psrc, verbose=True)
print("LSR velocity: {0:+8.3f}".format(v.to(u.km/u.s))) print("LSR velocity: {0:+8.3f}".format(v.to(u.km / u.s)))
f = doppler_frequency(t, loc, psrc, freq_hi) f = doppler_frequency(t, loc, psrc, freq_hi)
print("Doppler HI: {}".format(f)) 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) 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) 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) cmd = "~harm/bin/doppler_bl %f %f 2000 %s %f %f 0 %f" %\
print(cmd+"\n\n") (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 # 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) 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") psrc = SkyCoord(ra=6 * u.hourangle, dec=0.0 * u.deg, frame="icrs")
t = Time("2018-06-21T12:00:00",scale = "utc",format = "isot") t = Time("2018-06-21T12:00:00", scale="utc", format="isot")
v = vlsr(t,loc,psrc,verbose=True) v = vlsr(t, loc, psrc, verbose=True)
print("LSR velocity: {0:+8.3f}".format(v.to(u.km/u.s))) print("LSR velocity: {0:+8.3f}".format(v.to(u.km / u.s)))
f = doppler_frequency(t, loc, psrc, freq_hi) f = doppler_frequency(t, loc, psrc, freq_hi)
print("Doppler HI: {}".format(f)) 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) 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) 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) cmd = "~harm/bin/doppler_bl %f %f 2000 %s %f %f 0 %f" %\
print(cmd+"\n\n") (psrc.ra.rad, psrc.dec.rad, t.unix, loc.lat.rad, loc.lon.rad, freq_hi)
print(cmd + "\n\n")
# Test 3: Narrabri calculator # Test 3: Narrabri calculator
loc = EarthLocation.from_geodetic(lat = -30.3128*u.deg, lon = 149.5502*u.deg, height = 0*u.m) 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") psrc = SkyCoord(ra=12 * u.hourangle, dec=0.0 * u.deg, frame="icrs")
t = Time("2018-04-16T12:00:00",scale = "utc",format = "isot") t = Time("2018-04-16T12:00:00", scale="utc", format="isot")
v = vlsr(t,loc,psrc,verbose=True) v = vlsr(t, loc, psrc, verbose=True)
print("LSR velocity: {0:+8.3f}".format(v.to(u.km/u.s))) 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