moon_position_table_new3.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 import sys
3 import math
4 import ephem
5 import time
6 from datetime import datetime
7 from pytz import timezone
8 import pytz
9 from optparse import OptionParser
10 import subprocess
11 from ROOT import TTree, TFile, AddressOf, gROOT
12 
13 # Make a tree
14 f = TFile('semfe_tmp.root','RECREATE')
15 ntp_tr = TTree('ntp','My test tree')
16 
17 # Create a struct
18 gROOT.ProcessLine(\
19  "struct MyStruct{\
20  ULong64_t Night_Twilight;\
21  ULong64_t Twilight_Day;\
22  ULong64_t Day_Twilight;\
23  ULong64_t Twilight_Night;\
24  };")
25 
26 from ROOT import MyStruct
27 
28 # Create branches in the tree
29 s = MyStruct()
30 ntp_tr.Branch('Night_Twilight',AddressOf(s,'Night_Twilight'),'Night_Twilight/l')
31 ntp_tr.Branch('Twilight_Day',AddressOf(s,'Twilight_Day'),'Twilight_Day/l')
32 ntp_tr.Branch('Day_Twilight',AddressOf(s,'Day_Twilight'),'Day_Twilight/l')
33 ntp_tr.Branch('Twilight_Night',AddressOf(s,'Twilight_Night'),'Twilight_Night/l')
34 
35 version = "1.0"
36 usage = "usage: %prog [options]"
37 parser= OptionParser(usage=usage,version=version)
38 
39 parser.add_option("-r","--resolution", dest="res",default=60,
40  help="Resolution in seconds for calculations (default: 60)")
41 parser.add_option("-d","--date", dest="startdate",default='2015/10/31',
42  help='Date to start table from (example "2015/10/31")')
43 parser.add_option("-n","--numdays", dest="numdays",default='31',
44  help='Number of days from start to calculate')
45 
46 (options,args)=parser.parse_args()
47 
48 resolution = float(options.res)
49 
50 # Our timezone information
51 utc = pytz.utc
52 central = timezone('US/Central')
53 # mytime = central.localize(datetime.fromtimestamp(input_time))
54 
55 # We care about the sun as seen from NOvA (far det)
56 moon = ephem.Moon()
57 sun = ephem.Sun()
58 nova = ephem.Observer()
59 
60 # Set the latitude and longitude of Nova
61 #nova.lat = '48:22.715780'
62 #nova.lon = '-92:49.876930'
63 
64 # These are the coordinates of Orr
65 nova.lat = '48:04.0'
66 nova.lon = '-92:50.0'
67 
68 # Time Constants
69 minute = 60
70 second = 1
71 hour = 60*60
72 oneMinute = 1.0/(1440)
73 oneSecond = 1.0/(86400)
74 
75 # What day?
76 #start_day = '2015/10/31'
77 start_day = options.startdate
78 numdays = int(options.numdays)
79 
80 def SetupAnalysisRegion( AnalysisRegions ) :
81  AnalysisRegions['day_time'] = 0
82  AnalysisRegions['night_time'] = 0
83  AnalysisRegions['twilight_time'] = 0
84  AnalysisRegions['total_day_time'] = 0
85  AnalysisRegions['total_night_time'] = 0
86  AnalysisRegions['total_twilight_time'] = 0
87  AnalysisRegions['n2t_time_boundary'] = 0
88  AnalysisRegions['t2d_time_boundary'] = 0
89  AnalysisRegions['d2t_time_boundary'] = 0
90  AnalysisRegions['t2n_time_boundary'] = 0
91 
92 
93 AnalysisRegionsSun = {}
94 AnalysisRegionsMoon = {}
95 
96 SetupAnalysisRegion( AnalysisRegionsSun )
97 SetupAnalysisRegion( AnalysisRegionsMoon )
98 
99 
100 last_region = 999
101 
102 theta=10
103 
104 def AnalRegion(a, theta) :
105  dayBoundary = theta* math.pi / 180.0
106  nightBoundary = -1 *theta * math.pi / 180.0
107  type = 0
108 
109 
110  if a < nightBoundary :
111  type=1
112  if a > dayBoundary :
113  type=2
114  if (a < dayBoundary) and (a > nightBoundary) :
115  type=3
116 
117 # print '%f %f, %f %d' %(dayBoundary, nightBoundary, a, type)
118 
119  return type
120 
121 def FillAnal(a, region, time, resolution) :
122  global last_region
123  #
124  # Added a resolution step to the counters
125  stepSize = oneSecond*resolution
126  if region == 1 :
127  a['night_time']+= stepSize
128  a['total_night_time']+= stepSize
129  if region == 2 :
130  a['day_time']+= stepSize
131  a['total_day_time']+= stepSize
132  if region == 3 :
133  a['twilight_time']+= stepSize
134  a['total_twilight_time']+= stepSize
135 
136  if (region == 1) and (last_region == 3) :
137  a['t2n_time_boundary'] = time
138  if (region == 3) and (last_region == 1) :
139  a['n2t_time_boundary'] = time
140  if region == 2 and last_region == 3 :
141  a['t2d_time_boundary'] = time
142  if region == 3 and last_region == 2 :
143  a['d2t_time_boundary'] = time
144  last_region = region
145 
146 def tl(date):
147  return ephem.Date(date - 5 * ephem.hour)
148 
149 def computeDayBoundaryMoon(AnalysisRegions, day, theta, resolution):
150  global oneMinute
151  AnalysisRegions['night_time']=0
152  AnalysisRegions['day_time']=0
153  AnalysisRegions['twilight_time']=0
154 # for i in xrange(1,int(86400/resolution)):
155 # nova.date = day + (oneSecond*resolution * i)
156  for i in xrange(1,1440):
157  nova.date = day + oneMinute * i
158  #moon.compute(nova)
159  #region = AnalRegion(moon.alt, theta)
160  sun.compute(nova)
161  region = AnalRegion(sun.alt, theta)
162  FillAnal(AnalysisRegions, region, nova.date, resolution)
163 
164  # Print the transition times in UTC
165  print '%20s, %20s, %20s, %20s' % ((AnalysisRegions['n2t_time_boundary']),(AnalysisRegions['t2d_time_boundary']),
166  (AnalysisRegions['d2t_time_boundary']),(AnalysisRegions['t2n_time_boundary']))
167 
168  tmp_1 = str(AnalysisRegions['n2t_time_boundary']) + "UTC"
169  proc_1 = subprocess.Popen(["/grid/fermiapp/products/nova/externals/novadaq/v04_05_00/slf6.x86_64.e9.debug/bin/NovaTimeConvert",
170  tmp_1], stdout=subprocess.PIPE)
171  semfe_1 = proc_1.stdout.read()
172  proc_1.wait()
173 
174  tmp_2 = str(AnalysisRegions['t2d_time_boundary']) + "UTC"
175  proc_2 = subprocess.Popen(["/grid/fermiapp/products/nova/externals/novadaq/v04_05_00/slf6.x86_64.e9.debug/bin/NovaTimeConvert",
176  tmp_2], stdout=subprocess.PIPE)
177  semfe_2 = proc_2.stdout.read()
178  proc_2.wait()
179 
180  tmp_3 = str(AnalysisRegions['d2t_time_boundary']) + "UTC"
181  proc_3 = subprocess.Popen(["/grid/fermiapp/products/nova/externals/novadaq/v04_05_00/slf6.x86_64.e9.debug/bin/NovaTimeConvert",
182  tmp_3], stdout=subprocess.PIPE)
183  semfe_3 = proc_3.stdout.read()
184  proc_3.wait()
185 
186  tmp_4 = str(AnalysisRegions['t2n_time_boundary']) + "UTC"
187  proc_4 = subprocess.Popen(["/grid/fermiapp/products/nova/externals/novadaq/v04_05_00/slf6.x86_64.e9.debug/bin/NovaTimeConvert",
188  tmp_4], stdout=subprocess.PIPE)
189  semfe_4 = proc_4.stdout.read()
190  proc_4.wait()
191 
192  print semfe_1.split()[13], semfe_2.split()[13], semfe_3.split()[13], semfe_4.split()[13]
193 
194  s.Night_Twilight = long(semfe_1.split()[13])
195  s.Twilight_Day = long(semfe_2.split()[13])
196  s.Day_Twilight = long(semfe_3.split()[13])
197  s.Twilight_Night = long(semfe_4.split()[13])
198  ntp_tr.Fill()
199 
200 # theTime = datetime.fromtimestamp(time.time())
201 
202  unixtimes = []
203  del unixtimes[:]
204  for key in ['n2t_time_boundary','t2d_time_boundary','d2t_time_boundary','t2n_time_boundary']:
205  t = ephem.Date(AnalysisRegions[key]).datetime()
206  # unixtime = time.mktime(theTime.timetuple())
207  unixtime = 0
208  try:
209  unixtime = time.mktime(t.timetuple())
210  except:
211  print 'Value out of range'
212  unixtimes.append(unixtime)
213  # print '%20s, %20s, %20s, %20s' % (unixtimes[0], unixtimes[1], unixtimes[2], unixtimes[3] )
214  sys.stdout.flush()
215 
216 ########################################
217 # Loop over the days we care about
218 ########################################
219 # Compute the positions
220 for theta in range(10,11):
221  AnalysisRegionsSun['total_night_time']=0
222  AnalysisRegionsSun['total_day_time']=0
223  AnalysisRegionsSun['total_twilight_time']=0
224  AnalysisRegionsMoon['total_night_time']=0
225  AnalysisRegionsMoon['total_day_time']=0
226  AnalysisRegionsMoon['total_twilight_time']=0
227  print '#'*80
228  print '# Angle Boundary: %s degrees' %(theta)
229  print '# Object: Moon'
230  print '# %18s, %20s, %20s, %20s' %('Night/Twilight', 'Twilight/Day', 'Day/Twilight','Twilight/Night')
231  print '#'*80
232  for ndays in range(1,numdays):
233  day = ephem.Date(start_day)
234  day += ndays
235  computeDayBoundaryMoon(AnalysisRegionsMoon, day, theta, resolution)
236 
237 f.Write()
238 f.Close()
def computeDayBoundaryMoon(AnalysisRegions, day, theta, resolution)
def FillAnal(a, region, time, resolution)
def SetupAnalysisRegion(AnalysisRegions)