2014
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)
file_to_process = '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.clf()
plt.plot(vels_wrt_earth, spectrum)
plt.xlabel("Velocity (km/s)")
plt.ylabel("Spectrum Value (Arbitrary)")
plt.savefig("spectrum.png")