##// END OF EJS Templates
Add skynoise with python
Add skynoise with python

File last commit:

r30:d43c87bf9288
r30:d43c87bf9288
Show More
plots.py
188 lines | 5.9 KiB | text/x-python | PythonLexer
Juan C. Espinoza
Add skynoise with python
r30
import os
import time
import numpy
import scipy
import base64
from io import BytesIO
import scipy.interpolate
from matplotlib.figure import Figure
from utils import TimeTools
def skyNoise(jd, ut=-5.0, longitude=-76.87, filename='/app/utils/galaxy.txt'):
"""
hydrapos returns RA and Dec provided by Bill Coles (Oct 2003).
Parameters
----------
jd = The julian date of the day (and time), scalar or vector.
dec_cut = A scalar giving the declination to get a cut of the skynoise over Jica-
marca. The default value is -11.95.
filename = A string to specify name the skynoise file. The default value is skynoi-
se_jro.dat
Return
------
maxra = The maximum right ascension to the declination used to get a cut.
ra = The right ascension.
Examples
--------
>> [maxra,ra] = skynoise_jro()
>> print maxra, ra
139.45 -12.0833333333
Modification history
--------------------
Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009.
"""
# Defining date to compute SkyNoise.
[year, month, dom, hour, mis, secs] = TimeTools.Julian(jd).change2time()
is_dom = (month==9) & (dom==21)
if is_dom:
tmp = jd
jd = TimeTools.Time(year,9,22).change2julian()
dom = 22
# Reading SkyNoise
g = os.getcwd()
f = open(filename)
lines = f.read()
f.close()
nlines = 99
lines = lines.split('\n')
data = numpy.zeros((2,nlines))*numpy.float32(0.)
for ii in numpy.arange(nlines):
line = numpy.array([lines[ii][0:6],lines[ii][6:]])
data[:,ii] = numpy.float32(line)
# Getting SkyNoise to the date desired.
otime = data[0,:]*60.0
opowr = data[1,:]
hour = numpy.array([0,23])
mins = numpy.array([0,59])
secs = numpy.array([0,59])
LTrange = TimeTools.Time(year,month,dom,hour,mins,secs).change2julday()
LTtime = LTrange[0] + numpy.arange(1440)*((LTrange[1] - LTrange[0])/(1440.-1))
lst = TimeTools.Julian(LTtime + (-3600.*ut/86400.)).change2lst()
ipowr = lst*0.0
# Interpolating using scipy (inside max and min "x")
otime = otime/3600.
val = numpy.where((lst>numpy.min(otime)) & (lst<numpy.max(otime))); val = val[0]
tck = scipy.interpolate.interp1d(otime,opowr)
ipowr[val] = tck(lst[val])
# Extrapolating above maximum time data (23.75).
uval = numpy.where(lst>numpy.max(otime))
if uval[0].size>0:
ii = numpy.min(uval[0])
m = (ipowr[ii-1] - ipowr[ii-2])/(lst[ii-1] - lst[ii-2])
b = ipowr[ii-1] - m*lst[ii-1]
ipowr[uval] = m*lst[uval] + b
if is_dom:
lst = numpy.roll(lst,4)
ipowr = numpy.roll(ipowr,4)
new_lst = numpy.int32(lst*3600.)
new_pow = ipowr
return ipowr, LTtime, lst
def skynoise_plot(year, month, day):
"""
getPlot method creates a skynoise map over Jicamarca for a desired day. Additionally
save a PNG file of this plot.
Examples
--------
>> SkyNoisePlot(skypower,skytime,skytime_lst).getPlot()
Modification history
--------------------
Created by Freddy R. Galindo, ROJ, 18 October 2009.
"""
# Working with the time before to plot the SkyNoise
julian = TimeTools.Time(year, month, day).change2julday()
power, tm, time_lst = skyNoise(julian)
secs = TimeTools.Julian(tm).change2secs()
secs_lst = time_lst*1.
if time_lst.argmin()>0:
secs_lst[time_lst.argmin():] = secs_lst[time_lst.argmin():] + 24.
secs_lst = secs_lst*3600.
label_secs = time.localtime(secs[power.argmax()] + time.timezone)
label_secs_lst = time.localtime(secs_lst[power.argmax()] + time.timezone)
# Creating labels for main x-labelticks (Local Time):
snow = TimeTools.Time(year, month, day, 0, 0, 0).change2secs()
today = secs - snow
xticks_dpos = []
xticks_dval = []
for ix in [0,120,240,360,480,600,720,840,960,1080,1200,1320,1439]:
xticks_dpos.append(today[ix])
time_tuple = time.localtime(today[ix] + time.timezone)
xticks_dval.append(time.strftime('%H:%M',time_tuple))
if ix==1439:xticks_dval[12] = ''
# Creating labels for secondary x-labelticks (Sidereal Time):
xticks_upos = [secs_lst[0],secs_lst[-1]]
xticks_uval = ['','']
for ix in [0,120,240,360,480,600,720,840,960,1080,1200,1320]:
index = numpy.argmin(numpy.abs(time_lst - (ix/60. + 1.)))
xticks_upos.append(secs_lst[index])
time_tuple = time.localtime(secs_lst[index] + time.timezone)
if time_tuple[4]==59:
tmp_list = list(time_tuple)
tmp_list[4] = 0
tmp_list[3] = tmp_list[3] + 1
time_tuple = tuple(tmp_list)
xticks_uval.append(time.strftime('%H:%M',time_tuple))
# Creating SkyNoise Map.
fig = Figure()
fig.subplots_adjust(top=0.80)
main_mesg = "SKY BRIGHTNESS AT 50Mhz - Date: "
main_mesg = main_mesg + time.strftime("%d-%b-%Y (%j)",label_secs)
main_mesg = main_mesg + "\nGalaxy Pass at " + time.strftime("%H:%M:%S LT",label_secs)
main_mesg = main_mesg + time.strftime(" (%H:%M:%S LST)",label_secs_lst)
fig.suptitle(main_mesg, fontsize=12)
axl = fig.add_subplot(1,1,1)
axl.plot(today,power,'b-')
axr = axl.twinx()
axl.set_xlabel('Local Time (Hours)')
axl.set_ylabel('Signal Strengh (mA)')
axr.set_ylabel('Temperature [K]x10^3')
axl.set_xlim(0,24)
axl.set_ylim(0,115)
axr.set_ylim(0,50)
axl.set_xticks(xticks_dpos)
axl.set_xticklabels(xticks_dval, fontsize=8)
axl.grid()
axt = axl.twiny()
axt.set_xlabel('Local Sidereal Time (Hours)')
axt.set_xlim(secs_lst[0],secs_lst[-1])
axt.set_xticks(xticks_upos)
axt.set_xticklabels(xticks_uval, fontsize=8)
buf = BytesIO()
fig.savefig(buf, format="png")
return buf