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