In [1]:
from matplotlib import pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = 12,8
In [2]:
import os
import sys
import numpy as np
from astropy.time import Time
import astropy.coordinates
from myPyLib.Mollweide import moll
#from matplotlib import pyplot as plt
import util
import schdutil
import Moon
In [3]:
tel, yr, mn, dy = 'bok', 2016, 9, 3
moon_dis_limit = 50.0
moon_dis_limit, airmass_limit, ha_limit=50.0,1.75,3.0
obs_begin, obs_end = None, None
# airmass lower limit, set this to avoid zenith
airmass_lbound = 1.005 # about 84 deg
In [4]:
# load site and telescope basic data
site = schdutil.load_basic(tel)
# load fields and plan
plans = schdutil.load_expplan(tel)
fields = schdutil.load_field(tel)
plancode = plans.keys()
plancode.sort()
nplan = len(plancode)
In [5]:
# night parameters
# mjd of 18:00 of site timezone, as code of tonight
mjd18 = int(schdutil.mjd(yr, mn, dy, 18, 0, 0, site.tz)) % 10000
# mjd of local midnight, as calculate center
mjd24 = schdutil.mjd(yr, mn, dy, 24 - site.lon / 15.0, 0, 0, 0)
tmjd24 = Time(mjd24, format="mjd")
# night length (hour)
night_len = schdutil.night_len(mjd24, site.lat)
dark_len = night_len - 2.5 # assume twilight 1.25 hours
# local sidereal time of midnight
lst24 = schdutil.fmst (yr, mn, dy, site.lon)
# timezone correction: between local time and timezone standard time
tzcorr = site.tz - site.lon / 15.0
# observation start and end time, in timezone time
sunset_time, sunrise_time = 24 + tzcorr - night_len / 2, 24 + tzcorr + night_len / 2
# if obs time is given, use given time
if obs_begin is None :
obs_begin = 24 + tzcorr - dark_len / 2
else :
if obs_begin < 12.0 : obs_begin += 24
if obs_end is None :
obs_end = 24 + tzcorr + dark_len / 2
else :
if obs_end < 12.0 : obs_end += 24
# moon position at midnight, as mean coord to calculate moon-object distance
mpos = astropy.coordinates.get_moon(tmjd24)
mphase = Moon.moon_phase(tmjd24)
# sun position at midnight
spos = astropy.coordinates.get_sun(tmjd24)
In [6]:
# find all obsed file, and mark them
obsedlist = schdutil.ls_files("{tel}/obsed/*/obsed.J*.lst".format(tel=tel))
schdutil.load_obsed(fields, obsedlist, [], plancode, None)
afields = np.array(fields.values())
ara = np.array([f.ra for f in afields])
ade = np.array([f.de for f in afields])
In [7]:
# mark fields near moon and sun
moon_dis = Moon.distance(mpos.ra.deg, mpos.dec.deg, ara, ade)
for f in afields[np.where(moon_dis < moon_dis_limit)]:
f.tag += 10
sun_dis = Moon.distance(spos.ra.deg, spos.dec.deg, ara, ade)
for f in afields[np.where(sun_dis < 75)]:
f.tag += 20
atag = np.array([f.tag for f in afields])
In [8]:
# keep only unfinished fields, and must away from moon
newfield = afields[np.where(atag <= 1)]
# count histogram
n_tag = len(afields)
n_tag_01 = sum(atag <= 1)
n_tag_2 = sum(atag % 10 == 2)
n_tag_10 = sum((atag % 10 <= 1) & (atag >= 10))
In [9]:
# blocks and unique blocks
newfieldblock = np.array([f.bk for f in newfield])
newblockset = set(newfieldblock)
n_block = len(newblockset)
In [12]:
newfieldblock
Out[12]:
In [13]:
# block parameter
newblock = {}
for b in newblockset :
f_in_b = newfield[np.where(newfieldblock == b)]
newblock[b] = schdutil.block_info(b, f_in_b)
In [10]:
clock_now = obs_begin
lst_now = lambda : (lst24 + clock_now) % 24.0 # lst of start, use lst_now() to call this
# define a lambda rank function
rank = lambda aa : aa.argsort().argsort()
In [ ]:
block_sn = 0 # block sn, count blocks, and also for output
time_skip = 0 # skipped time, sum when no good block available
span_skip = 1.0 / 60.0 # how long skipped each loop
In [71]:
rep_fmt = "{sn:02}: #{bn} ({ra:9.5f} {de:+9.5f}) {airm:4.1f} @ {clock:5}".format
sum_fmt = "{clock:5}: {sn:02} #{bn} ({ra:9.5f} {de:+9.5f}) {airm:4.1f}\n".format
chk_fmt = "{ord:04} {bn} ({ra:9.5f} {de:+9.5f}) {airm:4.2f} {key:>5.1f}\n".format
scr_fmt = (site.fmt + "\n").format
In [61]:
bra = np.array([b.ra for b in newblock.values()])
bde = np.array([b.de for b in newblock.values()])
bname = np.array(newblock.keys())
In [62]:
# calculate airmass for all available block
ha = np.abs(lst_now() - bra) / 15.0
airm = schdutil.airmass(site.lat, lst_now(), bra, bde)
In [15]:
plt.hist(airm, range=(0,2))
Out[15]:
In [63]:
ix_1 = np.where((airm < airmass_limit) & (airm > airmass_lbound) &
((ha < ha_limit) | (ha > 24.0 - ha_limit)))
In [64]:
minde = bde[ix_1].min()
ix_2 = np.where((airm < airmass_limit) & (airm > airmass_lbound) &
((ha < ha_limit) | (ha > 24.0 - ha_limit)) &
(bde < minde + 4.0))
len(ix_2[0])
Out[64]:
In [65]:
# make key for each block, key = airmass rank + ra rank + de rank
airm_2, bra_2, bde_2, bname_2 = airm[ix_2], bra[ix_2], bde[ix_2], bname[ix_2]
rank_2 = rank(airm_2) + rank(bra_2)/2.0 + rank(bde_2)/2.0
rank_2
Out[65]:
In [66]:
plt.plot(bra[ix_1]/15.0,bde[ix_1],".")
plt.xlim(0,4)
Out[66]:
In [67]:
lst_now(), clock_now
Out[67]:
In [68]:
plt.hist(ha)
Out[68]:
In [69]:
# choose the best block, and generate scripts, and generate check file
sort_2 = rank_2.argsort() # sort for generate check file
ix_best = rank_2.argmin() # the best block
bname_best = bname_2[ix_best]
block_best = newblock[bname_best]
# inc sn
block_sn += 1
In [72]:
print("#{ord:>3} {bn:>7} ({ra:^9} {de:^9}) {airm:>4} {key:>5}\n".format(
ord="ORD",bn="B_NAME", ra="RA", de="DEC", airm="AIRM",key="KEY")).strip()
for i in range(len(ix_2[0])) :
s = sort_2[i]
b = newblock[bname_2[s]]
print(chk_fmt(ord=i+1, bn=b.bname, ra=b.ra, de=b.de, airm=airm_2[s], key=rank_2[s])).strip()
In [73]:
block_time = 0 # time cost for this block, in seconds
#with open(scr_fn, "w") as scr_f :
if True:
# script: plan loop, field loop, do only factor < 1
for p in plancode :
for f in block_best.fields :
factor_work = 1.0 - f.factor[p]
nrepeat = int(np.ceil(factor_work / plans[p].factor))
for i in range(nrepeat) :
print(scr_fmt(e=schdutil.exposure_info.make(plans[p], f))).strip()
block_time += plans[p].expt + site.inter
for f in block_best.fields :
f.tag = 9 # mark as used in this night
In [59]:
block_time
Out[59]:
In [74]:
# remove used block from newblock
del newblock[bname_best]
# clock forward
clock_now += block_time / 3600.0
In [ ]: