"""
Adapted from read_srt_fft.py
Tom Essinger-Hileman
November 2013
Scripts to read data from the Johns Hopkins Small Radio Telescope.
"""
import numpy as np
import pylab as plt
from scipy.ndimage.filters import gaussian_filter1d
import scipy.stats
import scipy.stats.stats
def read_srt_data(filename, data_range = None):
# Open file and read data.
f = open(filename, 'r')
data = f.readlines()
#Get velocity with respect ot the LSR from first line of file
pointingData = data[0].split()
vlsr = float(pointingData[15])
# Get frequencies for the spectra from the second line of the file.
# (Assume that this does not change as the data file progresses.)
freqData = data[1].split() # splits line up by spaces
start = float(freqData[1])
spacing = float(freqData[5])
Nfreqs = int(freqData[16])
freqs = start + spacing*np.arange(Nfreqs) # In MHz.
# Create an array to hold the spectrum and load from last line
# The last line is the average spectrum (the file has a running average)
spectrum = np.zeros(Nfreqs)
spectrum_strings = data[-1].split()
for j in range(Nfreqs):
spectrum[j] = float(spectrum_strings[j])
if data_range != None:
spectrum = spectrum[data_range[0]:data_range[1]]
freqs = freqs[data_range[0]:data_range[1]]
return spectrum, freqs, vlsr
def freq_to_vel(freqs):
"""
Takes an array of frequencies in and converts to redshifted 21-cm velocity.
Velocity defined as positive if moving away from earth
Returns in units of km/s.
"""
HIfreq = 1420.40575177 # 21-cm frequency, in MHz
c = 299792.458 # speed of light in km/s
return -1*c*(freqs - HIfreq)/(HIfreq)
plt.clf()
file_to_process = 'marriage_data/g90.rad' # replace with your file
spectrum, freqs, vlsr = read_srt_data(file_to_process)
vels_wrt_earth = freq_to_vel(freqs)
vels_wrt_lsr= vels_wrt_earth - vlsr
plt.plot(vels_wrt_earth, spectrum)
plt.xlabel("Velocity (km/s)")
plt.ylabel("Spectrum Value (Arbitrary)")
plt.savefig("spectrum.png")
#Remove spike
spectrum, freqs, vlsr = read_srt_data(file_to_process, data_range=[20,-1])
vels_wrt_earth = freq_to_vel(freqs)
plt.plot(vels_wrt_earth, spectrum)
plt.savefig("spectrum.png")
def read_and_get_velocity(filename, data_range=None):
spectrum, freqs, vlsr = read_srt_data(filename, data_range=data_range)
vels_wrt_earth = freq_to_vel(freqs)
return vels_wrt_earth, spectrum, vlsr
#Estimate background template from reference spectrum
vels_wrt_earth_ref, spectrum_ref, vlsr_ref = read_and_get_velocity( \
'essinger_data/reference.rad', data_range=[20,-1])
plt.plot(vels_wrt_earth_ref, spectrum_ref)
plt.savefig('spectrum.png')
smoothing_length = 3
spectrum_ref_smooth=gaussian_filter1d(spectrum_ref, smoothing_length)
plt.plot(vels_wrt_earth_ref, spectrum_ref_smooth)
plt.savefig('spectrum.png')
#Interpolate over HI in reference spectrum
indices = np.where( (vels_wrt_earth < 0)*(vels_wrt_earth > -70) )
vels_sub = vels_wrt_earth[indices]
ref_sub = spectrum_ref_smooth[indices]
slope = (ref_sub[-1]-ref_sub[0])/(vels_sub[-1]-vels_sub[0])
interp = np.ones(len(ref_sub))*ref_sub[0]+slope*(vels_sub-vels_sub[0])
spectrum_ref_smooth[indices] = interp
plt.plot(vels_wrt_earth_ref, spectrum_ref_smooth)
plt.savefig('spectrum.png')
# Fit background template and remove background
# Fit background to data outside of signal region
# (Determine background region with vel_wrt_earth,
# which is the same regardless of pointing)
# use linear model data = a*bg_template+b
def fit_background(data, ref, vels, vel_less_than=-200., \
vel_greater_than=50.):
"""
Returns background estimate and background subtracted data
"""
indices = np.where((vels > vel_greater_than)+(vels < vel_less_than))
data_sub = data[indices]
ref_sub = ref[indices]
slope,intercept,r_val,p_val,err = scipy.stats.linregress(ref_sub, \
data_sub)
bg = intercept + slope*ref
bgsub = data - bg
return bg, bgsub
bg, spectrum_bgsub = fit_background(spectrum, spectrum_ref_smooth, \
vels_wrt_earth, vel_less_than=-200, \
vel_greater_than=50.)
plt.plot(vels_wrt_earth_ref, bg)
plt.plot(vels_wrt_earth_ref, spectrum_bgsub)
plt.savefig('spectrum.png')
# Determine most blue-shifted frequency
# Simple method: choose threshold above which we call it HI
threshold=5
indices = np.where(spectrum_bgsub>threshold)
vels_HI= vels_wrt_earth[indices]
v_max = max(vels_HI)
print "Maximum velocity %.0f km/s" % (v_max-vlsr)
plt.clf()
plt.plot(vels_wrt_earth-vlsr,spectrum_bgsub)
plt.axhline(threshold,ls='dashed',color='black')
plt.xlabel("Velocity (km/s)")
plt.ylabel("Spectrum Value (Arbitrary)")
plt.xlim([-200,200])
plt.savefig("spectrum_final_90.png" )