SRT Data Analysis 2

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