#SimulateSkyCoverage.py                  Lance Simms, LLNL 2012
#
#PURPOSE:
#  To map out the portion of sky covered by the CXBN as it does its observatoins.  
#  The satellite points its "head" at the sun and spins around this axis.  The 
#  purpose of this program is to determine which regions of the sky it sees. 
#
#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'
#
########

import os
execfile(os.path.expanduser('~')+'/Desktop/Software/Python/CZTCubeSat/CZTCubeSatSetup.py')

import ephem
from datetime import *
import cgkit.cgtypes.mat3 as mat3
import cgkit.cgtypes.vec3 as vec3
from math import sin, cos, tan, acos, atan2, exp, pi

##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     = 18
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(os.path.expanduser('~')+'/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 = 5.55
CountRateDXRB     = FluxFrom10_100keV*(EffectiveAreaCZT*2*pi*(1-cos(FOV_Deg*pi/180)))
FOV_Rad           = FOV_Deg*pi/180

#Make a grid of points for the RA and DEC
NumGridPoints = 60
RA_Grid_x  = frange(NumGridPoints-1)*2*pi/NumGridPoints-pi
DEC_Grid_y = frange(NumGridPoints-1)*pi/NumGridPoints-pi/2
RA_Grid, DEC_Grid = meshgrid(RA_Grid_x, DEC_Grid_y) 
SkyCoverageGrid = zeros([NumGridPoints, NumGridPoints], dtype=int32)

#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_sun     = float(sun.dec)
phi_sun       = float(sun.ra)
theta_sun     = 0 #3*pi/6
theta_sun_rot = pi/2-theta
phi_sun_rot   = pi/6 #-1*pi/3# pi/2#

#Now find the rotation matrix that brings the +z axis to the sun pointing vector
sun_rot_mat = mat3(cos(theta_sun_rot)*cos(phi), -sin(theta)*cos(phi)

#Figure out the vector to the sun
x_sun   = cos(phi)*sin(pi/2-theta_sun_rot)
if (x_sun < 1e-10): x_sun = 0
y_sun   = sin(phi)*sin(pi/2-theta_sun_rot)
if (y_sun < 1e-10): y_sun = 0
z_sun   = cos(pi/2-theta_sun_rot)
if (z_sun < 1e-10): z_sun = 0
sun_vec = vec3(x_sun, y_sun, z_sun)

#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):
  
    #Convert to radians
    CZTPointingAngleRad = CZTPointingAngle*pi/180


    #Account for poor handling of cosine and sine
    if (abs(ThisPerpVec.x) < 1e-10) : ThisPerpVec.x = 0
    if (abs(ThisPerpVec.y) < 1e-10) : ThisPerpVec.y = 0
    if (abs(ThisPerpVec.z) < 1e-10) : ThisPerpVec.z = 0

    #Now get the angles corresponding to the vector in the original frame
    DEC_perp = acos(ThisPerpVec.z)-pi/2
    RA_perp  = atan2(ThisPerpVec.y,ThisPerpVec.x)
#    print DEC_perp, RA_perp
    if (RA_perp < 0) : RA_perp = 2*pi+RA_perp
    if RA_perp == 0:
      if ((AngleNum >= 90) and (AngleNum < 270)) : RA_perp=RA_perp+pi

    #Figure out what areas of the sky are in the FOV
    DECIndsInFOV = (where(abs(DEC_perp-DEC_Grid_y) < FOV_Rad))[0]
    for ThisDECInd in DECIndsInFOV:
      ThisDEC = DEC_Grid_y[ThisDECInd]
      for ThisRAInd in arange(NumGridPoints):
        ThisRA = RA_Grid_x[ThisRAInd] 
    
        #Determine separation between CZT Pointing and Object
        ThisAngSep  = Astro_Angular_Separation(RA_perp, DEC_perp, ThisRA, ThisDEC) 
        if (abs(ThisAngSep) < FOV_Rad): SkyCoverageGrid[ThisDECInd, ThisRAInd]+=1


#Uncomment for plot that represents 600,000 seconds of integration
HardCode = 1
if HardCode : SpinPeriod = 6
if HardCode : NumOrbits  = 100000
  
#Now plot all the angles
Fig=mplot.figure(0, figsize=[10,5])
mplot.subplot(111,projection="mollweide")
mplot.grid(True)

#Plot some points
mplot.plot(0,pi/2,'ro')
mplot.plot(-pi/2,0,'bo')
mplot.plot(0,-pi/2,'go')
mplot.plot(pi/2,0,'yo')

#Set up an array
mplot.imshow(SkyCoverageGrid, extent=(-pi,pi,-pi/2,pi/2), aspect=0.5)

SaveFig = 1
if (SaveFig == 1):
  mplot.savefig('./Figures/TestSkyCoverage.png')
