Read SRT Data and Plot

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