#SimulateCZTDetectorSpin.py                  Lance Simms, LLNL 2010
#
#PURPOSE:
#  To determine the optimum collimator opening angle for the Cosmic X-Ray Background Nanosat.
#  The program will go through a scan of the sky with a given FOV and determine what fraction of 
#  time it is contaminated by a source. 
#
#KEYWORDS:
#  DateTime: str
#    The date and time for which we want to calculate the observation.  This should be 
#    in the pyephem prefered format of "Year/Month/Day Hour:Minute:Second" in UT.
#  Lat: str
#    The latitude of the satellite when it makes the observation
#  Long: str
#    The longitude of the satellite when it makes the observation
#  FOV_Deg: float
#    The collimating angle of the instrument.  This will determine the threshold for 
#    angles at which a source is considered detected.
#  SpinPeriod: float
#    The number of seconds per full 360 degree rotation of the detector
#  NumAngles: int
#    The number of angles at which to record the hits in the detector
#  NumOrbits: int
#    The number of orbits to use for summing the angular information
#  CountRateThresh: int
#    The minimum count rate that we will use to consider as contamination 
#    (from the A4 catalog)
#
#EXAMPLES
#  run OptimizeCollimatorOpeningAngle.py 
#  run OptimizeCollimatorOpeningAngle.py --CountRateThresh=1
#  run OptimizeCollimatorOpeningAngle.py --CountRateThresh=3
#  run OptimizeCollimatorOpeningAngle.py --CountRateThresh=5
 
#
########
import os
execfile(os.path.expanduser('~')+'/Desktop/Software/Python/CZTCubeSat/CZTCubeSatSetup.py')

import ephem
from datetime import *

##PARSE KEYWORDS
keywords    = ['DateTime=', 'Lat=', 'Long=', 'FOV_Deg=', 'SpinPeriod=', \
               'NumAngles=', 'NumOrbits=', 'CountRateThresh=']

#DEFAULT VALUES FOR KEYWORDS
DateTime        = '2012/01/01 20:00:00'
Lat             = '37.00'
Long            = '-122.00'
FOV_Deg         = 20
SpinPeriod      = 60
NumAngles       = 360
NumOrbits       = 1
CountRateThresh = 0 

#NOW PARSE KEYWORDS
opts, extraparams = getopt.getopt(sys.argv[1:],'',keywords)
for o,p in opts:
  if o in ['--DateTime']:
    DateTime = p
  elif o in ['--Lat']:
    Lat = p
  elif o in ['--Long']:
    Long = p
  elif o in ['--FOV_Deg']:
    FOV_Deg = float(p)
  elif o in ['--SpinPeriod']:
    SpinPeriod = float(p)
  elif o in ['--NumAngles']:
    NumAngles = float(p)
  elif o in ['--NumOrbits']:
    NumOrbits = int(p)
  elif o in ['--CountRateThresh']:
    CountRateThresh = float(p)

#Read in the A4 Catalog first of all 
Names, RAs, DECs, CountRates, CountRateErrs, Fluxes, DataTabs, Observatories = \
  readcol('/Users/lances/Desktop/Software/Python/CZTCubeSat/A4_Catalog_NoHeader.txt', \
          fsep="|",twod=False)

#Convert the RAs and DECs to decimal
Cat_Size = size(RAs)
RAs_Deg  = zeros(Cat_Size, dtype=float64)
DECs_Deg = zeros(Cat_Size, dtype=float64)
RAs_Rad  = zeros(Cat_Size, dtype=float64)
DECs_Rad = zeros(Cat_Size, dtype=float64)

for i in arange(len(RAs)):
  RAs_Deg[i]  = hms2str_deg_space(RAs[i])
  DECs_Deg[i] = dms2str_deg_space(DECs[i])
  RAs_Rad[i]  = (pi/180)*RAs_Deg[i]
  DECs_Rad[i] = (pi/180)*DECs_Deg[i]

#Convert the Count Rates to our 1cm x 2cm = 2cm^2 CZT Detector
EffectiveAreaCZT     = 2.   #Area of CZT in cm^2
EffectiveAreaHEAO_A4 = 80.  #Effective area of A4 (see Levine et al. 1984) 
CountRatesCZT        = CountRates*(EffectiveAreaCZT/EffectiveAreaHEAO_A4)

#This is the DXRB Flux calculated from the script CalculateBackgroundFlux.py (from 10-100keV)
FluxFrom10_100keV = 0.65
CountRateDXRB     = FluxFrom10_100keV*(EffectiveAreaCZT*2*pi*(1-cos(FOV_Deg*pi/180)))

#Format the time appropriately
Year  = int(DateTime[0:4])
Month = int(DateTime[5:7])
Day   = int(DateTime[8:10])
Hours = int(DateTime[11:13])
Mins  = int(DateTime[14:16])
Secs  = int(DateTime[17:19])
DateTime = datetime(Year, Month, Day, Hours, Mins, Secs)

#Create the observer object and assign date and location
ObsTech = ephem.Observer()
ObsTech.long, ObsTech.lat = float(Lat), float(Long)
ObsTech.date = DateTime

#Compute the position of the Sun and the Moon
sun, moon = ephem.Sun(), ephem.Moon()
sun.compute(ObsTech)
moon.compute(ObsTech)

#Print out the info
print "Sun's location :"
print sun.ra, sun.dec
print "Phi=", float(sun.ra), "Theta=", float(sun.dec)

#Assign angles of the Sun for testing purposes
theta = float(sun.dec)
phi   = float(sun.ra)
#theta = 3*pi/6
#phi   = -1*pi/3# pi/2#

#Make an array of the CZT Pointings
AngularInterval    = 360./NumAngles
CZTPointingAngles  = AngularInterval*arange(NumAngles)
TotHitsPerPointing = zeros(NumAngles, dtype=float32)
TimePerPointing    = float(SpinPeriod)/NumAngles

#Determine the total number of angles we will test against
NumAngles = 0
for DEC in (arange(181)-90):

  #Get the declination in Radians
  ThisDEC_Rad = DEC*pi/180
  NumRAs = floor(359*cos(ThisDEC_Rad)+1)
  NumAngles+=NumRAs
  #Will equal 41220

#Array to fill up and index for it
NumFOVs            = 30
ContaminationFrac  = zeros(NumAngles, dtype=uint8)
ContaminationArray = zeros(NumFOVs+1)
FOV_Array          = 3*(arange(NumFOVs)+1)
AngleNum           = 0

#Change the collimator opening angle
for (FOV, FOV_Num) in zip(FOV_Array, arange(NumFOVs)+1):

  #Clear the AngleNum array 
  AngleNum = 0
  ContaminationFrac = zeros(NumAngles, dtype=uint8)

  #Print the FOV
  print "Field of View: " + str(FOV) + " Degrees"

  #Now go through and calculate these angles rotated to the sun pointed frame  
  for DEC in (arange(181)-90):

    #Get the declination in Radians
    ThisDEC_Rad = DEC*pi/180

    #Number of RAs
    NumRAs = floor(359*cos(ThisDEC_Rad)+1)
    RAArr  = arange(NumRAs)*360./NumRAs

    #Now find all of the objects within FOV_Deg degrees of pointing
    for RA in RAArr:
  
      #Get this right ascension in Radians    
      ThisRA_Rad = RA*pi/180
  
      #Now find all of the objects within FOV_Deg degrees of pointing
      for i in arange(Cat_Size):

        if (CountRates[i] > CountRateThresh) :    

          if ((ThisDEC_Rad-DECs_Rad[i])*180/pi < FOV) :
            #Determine separation between CZT Pointing and Object
            ThisAngSep  = Astro_Angular_Separation(ThisRA_Rad, ThisDEC_Rad, RAs_Rad[i], DECs_Rad[i]) 
            ThisAngSep  = (180/pi)*ThisAngSep
    
            #If object is in FOV, calculate flux (no cos dependence yet)
            if (ThisAngSep < FOV):
     
              #Now add to the total hits
              ContaminationFrac[AngleNum] = 1 
              break
 
      #Increment the angle number
      AngleNum+=1

  #Now add the Number of hits to the contamination array
  ContaminationArray[FOV_Num] = size(where(ContaminationFrac==0))

#Plot the ContaminationArray/Frac
#See PlotSNRwObservingFractionCurve.py for values obtained from this 
mplot.figure(0)
mplot.clf()
ContaminationArray[0]=41220
ContaminationFrac = ContaminationArray/41220.
FOV_ArrayWZero    = 3*arange(NumFOVs+1)
mplot.plot(FOV_ArrayWZero, ContaminationFrac)

SaveFig = 0
if SaveFig == 1:
  mplot.savefig('./Figures/ContaminationFracVsOpeningAngle_Theshold_'+str(CountRateThresh))
