Reading-in SRT Data

gistfile1.txt

# Script here

plot.py

import scipy
import scipy.signal
import numpy as np
import matplotlib

matplotlib.use('agg')

import matplotlib.pyplot as plt
#%matplotlib inline

import astropy
from astropy.coordinates import SkyCoord, ICRS, EarthLocation, AltAz, GCRS
import astropy.units as u
from astropy.time import Time

plt.style.use('ggplot')

#import vlsr

import otter
import otter.bootstrap as bs
import otter.plot as ml


import jplephem
from jplephem.spk import SPK
kernel = SPK.open('/home/daniel/data/de430.bsp')
import astropy
from astropy.coordinates import SkyCoord, ICRS, FK5, EarthLocation, AltAz
from astropy.time import Time
import astropy.units as u
import astropy.constants as con

import math
import numpy as np


@u.quantity_input(canon_velocity=(u.meter / u.second))
def v_sun(source, apex = "18h03m50.29s +30d00m16.8s", 
                  canon_velocity = 20.0 * 1000 * u.meter / u.second):
    """
    Calculate the velocity of the sun relative to a point in the sky.
    
    Parameters
    ----------
    source : `astropy.coordinate.SkyCoord` object
        The point on the sky which the Earth's projected velocity should be 
        calculated for.
        
    apex : str, optional
        The point in the sky which the Sun is travelling towards. This should 
        be given as a string in the format "ra dec". The default value is
        18h03m50.29s +30d00m16.8s.
        
    canon_velocity : float
        The velocity at which the sun is travelling towards `apex`. This 
        defaults to 20 000 m/s. The value must be specified using astropy 
        units.
        
    Returns
    -------
    v_sun : `astropy.quantity`
        The velocity of the sun, relative to the source, with correct units.
    
    Notes
    -----
    The Sun is moving towards a point in the sky close to Vega at a velocity
    of around 20 km/s. In order to calculate the velocity of the Earth towards 
    any other point in the sky we need to project the velocity vector onto
    the position vector of that point in the sky, using a dot product.
    """
    if isinstance(source.obstime.value, str):
        equinox = source.obstime
    else:
        equinox = source.obstime[0]
    c = SkyCoord(apex, equinox="J2000", frame=FK5)
    c = c.transform_to(FK5(equinox=equinox, representation='cartesian'))
    c = c.cartesian
    source = source.cartesian
    v = canon_velocity * c.xyz
    return canon_velocity*np.dot(c.xyz,source.xyz)

def v_obs(source, observer):
    """
    Calculate the velocity of the observatory due to its rotation about the 
    Earth's barycentre in the direction of a point on the sky.
    
    Parameters
    ----------
    source : `astropy.coordinate.SkyCoord` object
        The point on the sky which the Earth's projected velocity should be 
        calculated for.
    observer : `astropy.coordinates.EarthLocation` object
        The location on the Earth of the observatory.
    
    Returns
    -------
    v_obs : `astropy.quantity` object
        The velocity of the observer relative to the source.
    """
    pos = [observer.x.value, observer.y.value, observer.z.value] * u.meter
    vel = [0, 0,  (2*np.pi)/(24.*3600.)]  / u.second # 459
    source = source.transform_to(AltAz(location=observer))
    return np.dot(np.cross(pos, vel), source.cartesian.xyz)*u.meter/u.second

def v_earth(source):
    """
    Calculate the velocity of the Earth relative to some point on the sky.
    
    Parameters
    ----------
    source : `astropy.coordinate.SkyCoord`
        The point on the sky which the Earth's projected velocity should be 
        calculated for.
        
    Returns
    -------
    v_earth : `astropy.quantity` object
        The velocity of the Earth relative to the source.
    """
    time = source.obstime.jd
    position, velocity = kernel[0,3].compute_and_differentiate(time)
    position2, velocity2 = kernel[3,399].compute_and_differentiate(time)
    velocity = velocity - velocity2
    source = source.cartesian
    return np.dot(((velocity * 1000 * u.meter / u.day).to(u.meter / u.second)).T, source.xyz) 

def v_lsr(source, observer):
    """
    Calculate the velocity of the local standard of rest for an observatory in
    the direction of a source on the sky.
    
    Parameters
    ----------
    source : `astropy.coordinate.SkyCoord` object
        The point on the sky which the Earth's projected velocity should be 
        calculated for.
    observer : `astropy.coordinates.EarthLocation` object
        The location on the Earth of the observatory.
    
    Returns
    -------
    v_lsr : `astropy.quantity` object
        The velocity of the observer relative to the local standard of rest.
    """
    v =  v_sun(source) + (v_earth(source) + v_obs(source, observer))
    return v

def doppler_fraction(v):
    return 1 + v / con.c

def produce_v(start, duration, rate, ra, dec):
    time_range = np.linspace(0, duration, duration/float(rate))*rate*u.second
    times = start + time_range
    source = SkyCoord(ra*u.degree, dec*u.degree, frame=FK5, obstime=times)

def find_nearest(array,value):
    idx = (np.abs(array-value)).argmin()
    return array[idx], idx


report = otter.Otter('/home/daniel/www/data/radio/srt2017/index.html',
                     theme='./themes/daniel',
                     title = "SRT HI Observations",
                     subtitle = '2017 Honours Lab')


import glob
import re

# Form the calibration curve
calibration = scipy.fromfile('/home/daniel/data/radio/srt2017/calibration.cal', dtype=np.float32)
calibration = (np.reshape(calibration, [-1,256]))
cal = [np.abs(cal/np.max(cal)) for cal in calibration]
calc = np.mean(cal[6:-1], axis=0)

files = glob.glob('/home/daniel/data/radio/srt2017/*.dat')
a = re.compile("ra(.*)dec(.*)time(.*).dat")

counter = 0
for datafile in files:

    #try:

        #counter+= 1
        #if counter>5: break
        print datafile
        data = scipy.fromfile(datafile, dtype=np.float32)
        bandwidth = 4e6
        centre = 1420.4e6
        filepath = datafile.split('/')
        s_ra = float(a.search(filepath[-1]).group(1))
        s_dec = float(a.search(filepath[-1]).group(2))
        s_time = a.search(filepath[-1]).group(3)

        frequencies = np.linspace(centre - bandwidth/2, centre + bandwidth/2, 256)/1e6
        data = np.reshape(data, [-1,256])

        # Average and normalise the data
        flatten = np.mean(data[8500:9000, :], axis=0)
        data = data/calc
        # Correct for the variation in the observations using the calibration line
        # This is probably a bit naiive, but here's a shot at it
        data /= data.max()

        start_time = Time(s_time, format='isot', scale='utc')
        time_range = np.linspace(0, (np.shape(data)[0])*5, np.shape(data)[0])*u.second
        times = start_time + time_range
        acre_road = EarthLocation(lat=55.9024278*u.deg, lon=-4.307582*u.deg, height=61*u.m)

        start = SkyCoord(s_ra*u.deg, s_dec*u.deg, frame=ICRS)
        start = start.transform_to(AltAz(obstime=times,location=acre_road))

        positions = SkyCoord(start.az[0], start.alt[0], frame = AltAz(obstime=times,location=acre_road), unit='deg')
        positions = positions.transform_to(ICRS)

        # v = v_lsr(positions[::10], acre_road)

        # output = []
        # print "Outputting VLSR"
        # with open(datafile+'.vlsr', 'w') as f: 
        #     for i in range(len(times)):
        #         print>>f, times[::10][i].value, v[i].value

        decl = positions.dec[0]

        rollamount = find_nearest(positions.ra.value, 0)[1]
        newend = find_nearest(positions.ra.value[-1000:-1], positions.ra.value[0])[1]

        data = np.roll(data[:-(1000 - newend)], -rollamount, axis=0)
        positions_r = np.roll(positions.ra[:-(1000 - newend)], -rollamount)

        f, ax = plt.subplots(2,1, figsize=(7.5, 2.5), sharey=True)
        plt.subplots_adjust(hspace=0)
        #ax.set_title('Drift scan at $\delta=${}'.format(decl))
        ax[0].plot(positions_r, np.sum(np.fliplr(data), axis=1))
        ax[1].imshow(np.fliplr((data)[:,20:236].T), aspect=20, cmap='inferno',extent = (positions_r[-1].value, positions_r[0].value,frequencies.min(), frequencies.max(),  ), vmax = np.max(data))
        ax[1].set_ylabel('Frequency (MHz)')
        ax[1].set_xlabel('Right Ascension (deg)')
        f.savefig('d{}.pdf'.format(int(decl.value)))
        plt.close()
        row = report + bs.Row(1)
        row[0] += "#Declination = {}".format(decl)
        row[0] += ml.Figure(report, f, dpi=100)
        row[0] += "<a href='{}'>VLSR Values</a>".format(datafile.split('/')[-1]+'.vlsr')
        row[0] += "<p>Observation start time: {} </p>".format(s_time)
        row[0] += "<p>Declination: {} </p>".format(decl)

    #except: 
    #    pass

report.show()