#SNRvsCollimatorOpeningAngle.py                  Lance Simms, LLNL 2010
#
#PURPOSE:
#  To simulate a full 360 revolution of the CZT detector aboard the 
#  DXRB nanosat.  The satellite points at the sun, and the CZT detector
#  spins in a plane perpendicular to the sun pointing axis.
#
#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
#
#EXAMPLES
#  run SimulateCZTDetectorSpin.py --SpinPeriod=360 --NumAngles=360 --FOV_Deg=20 --DateTime='2011/09/01 18:00:00'
#
########

execfile('/Users/lances/Desktop/Software/Python/CZTCubeSat/CZTCubeSatSetup.py')

import ephem
from datetime import *

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

#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

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

#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

#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

#Go through a number of orbits
for OrbitNum in arange(NumOrbits):

  #Print out the orbit number
  print 'Orbit Number : ' + str(OrbitNum)

  #Now go through and calculate these angles rotated to the sun pointed frame  
  for (AngleNum,CZTPointingAngle) in zip(arange(NumAngles), CZTPointingAngles):
  
    #Determine the angle in the unrotated plane
    x=cos(CZTPointingAngle*pi/180)
    y=sin(CZTPointingAngle*pi/180)
  
    #Rotate about the Declination axis with an angle of ra
    x_prime  = x*math.cos(phi)-y*math.sin(phi)
    y_prime  = x*math.sin(phi)+y*math.cos(phi)
  
    #Rotate about the Right Ascension=0 axis with an angle of dec
    #y-axis points to 0 RA, so rotate about x-axis
    x_dprime = x_prime
    y_dprime = y_prime*math.cos(theta)
    z_dprime = y_prime*math.sin(theta)
  
    #Now get the angles corresponding to the vector in the original frame
    DEC_perp = math.acos(z_dprime)-pi/2 
    RA_perp  = math.atan(y_dprime/x_dprime)
  
    #Determine number of hits from background
    ThisNumHits = 0
    NumHitsFromDXRB = CountRateDXRB*TimePerPointing
  
    #If number of hits is more than 1, use shot noise, else use probabilistic treatment
    if(NumHitsFromDXRB >= 1):
      ThisNumHits += numpy.random.poisson(NumHitsFromDXRB)
    else: 
      #Draw a random sample.  If the sample is less than HitsInPointing, count hit 
      RandomSample = numpy.random.random()
      if (RandomSample < NumHitsFromDXRB):
         ThisNumHits += 1
  
    #Now add to the total hits
    TotHitsPerPointing[AngleNum] += ThisNumHits
    if (AngleNum < 10):
       print 'Angle : '+str(AngleNum)+' Num ' + str(TotHitsPerPointing[AngleNum]) 
    #Now find all of the objects within FOV_Deg degrees of pointing
    for i in arange(Cat_Size):
  
      #Determine separation between CZT Pointing and Object
      ThisAngSep  = Astro_Angular_Separation(RA_perp, DEC_perp, RAs_Rad[i], DECs_Rad[i]) 
      ThisAngSep  = (180/pi)*ThisAngSep
  
      #If object is in FOV, calculate flux (no cos dependence yet)
      if (ThisAngSep < FOV_Deg):
  
        #See how many hits we expect
        HitsInPointing = CountRatesCZT[i]*TimePerPointing
        #If more than 1, use shot noise.  If less than one, use probabilistic treatment for hit
        if (HitsInPointing >= 1):
          ThisNumHits += numpy.random.poisson(HitsInPointing)
        else:
          #Draw a random sample.  If the sample is less than HitsInPointing, count hit 
          RandomSample = numpy.random.random()
          if (RandomSample < HitsInPointing): 
            ThisNumHits += 1
  
        #Now add to the total hits
        TotHitsPerPointing[AngleNum] += ThisNumHits 
        #print DEC_perp, RA_perp, RAs[i], DECs[i], HitsInPointing, ThisNumHits
  
  #Now plot all the angles
Fig=mplot.figure(0)
mplot.clf()
Fig.add_axes([0.100,0.100,0.800,0.6750])
mplot.plot(TotHitsPerPointing, label=str(SpinPeriod))
mplot.xlabel('Spin Phase (Degrees)')
mplot.ylabel('Counts')
mplot.figtext(0.05,0.94,r'Date      : '+str(DateTime))
#mplot.figtext(0.28,0.94, DateTime)
mplot.figtext(0.05,0.89,r'Sun RA   : '+str(sun.ra))
mplot.figtext(0.05,0.84,r'Sun DEC : '+str(sun.dec))
mplot.figtext(0.50,0.94,r'Spin Period : '+str(SpinPeriod)+' seconds')
mplot.figtext(0.50,0.89,r'Num Orbits  : '+str(NumOrbits))

SaveFig = 1
if (SaveFig == 1):
  mplot.savefig('./Figures/CZTSpinSimulation_Period_'+str(SpinPeriod)+'_NumOrbits_'+str(NumOrbits)+\
                '_'+str(DateTime.year)+'_'+str(DateTime.month)+'_'+str(DateTime.day)+'.png')
