#PlotSNRwObservingFractionCurve.py
#
#PURPOSE: 
#  To plot the curve
#
#    [D*t/(sqrt(D*t+S*t)+P*t)]*sqrt(F)
#
#  The contamination fractions are calculated using the script 
#  OptimizeCollimatorOpeningAngle.py
#
#  The SNR curve is calculated using 
#  SNRvsCollimatorOpeningAngle.py
#
############################################################################
import os
execfile(os.path.expanduser('~')+'/Desktop/Software/Python/CZTCubeSat/CZTCubeSatSetup.py')

NumFOVs           = 30
FOV_ArrayWZero    = 3*arange(NumFOVs+1)

#Any source is contamination
ContaminationFrac_InfThresh = \
[  1.00000000e+00,   9.18704512e-01,   7.43910723e-01,   5.46967492e-01,
   3.68243571e-01,   2.36365842e-01,   1.42940320e-01,   7.84813197e-02,
   3.79184862e-02,   1.24939350e-02,   3.17806890e-03,   2.18340611e-04,
   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
   0.00000000e+00,   0.00000000e+00,   0.00000000e+00]
ContaminationFrac_InfThresh = array(ContaminationFrac_InfThresh, dtype=float32)

#Sources above 1 counts/sec are contamination
ContaminationFrac_1CntPerSec = \
[ 1.        ,  0.93840369,  0.79354682,  0.62299854,  0.46193595,  0.32823872,
  0.22777778,  0.15907327,  0.10992237,  0.07411451,  0.04507521,  0.02372635,
  0.01220281,  0.0056526 ,  0.00150412,  0.        ,  0.        ,  0.        ,  0.,
  0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.,
  0.        ,  0.        ,  0.        ,  0.        ,  0.        ]
ContaminationFrac_1CntPerSec = array(ContaminationFrac_1CntPerSec, dtype=float32)

#Sources above 2 counts/sec are contamination
ContaminationFrac_2CntPerSec = \
[  1.00000000e+00,   9.69068413e-01,   8.91727317e-01,   7.82241630e-01,
   6.53056769e-01,   5.25788452e-01,   4.10868510e-01,   3.14216400e-01,
   2.36147501e-01,   1.77632217e-01,   1.34182436e-01,   9.87627365e-02,
   7.21979622e-02,   4.93449782e-02,   3.10771470e-02,   1.88258127e-02,
   1.04075691e-02,   4.34255216e-03,   7.03541970e-04,   0.00000000e+00,
   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
   0.00000000e+00,   0.00000000e+00,   0.00000000e+00]
ContaminationFrac_2CntPerSec = array(ContaminationFrac_2CntPerSec, dtype=float32)

#Sources above 3 counts/sec are contamination
ContaminationFrac_3CntPerSec = \
[ 1.        ,  0.98440078,  0.94209122,  0.88231441,  0.80931587,  0.72680738,
  0.6352984 ,  0.54473557,  0.46278506,  0.39133916,  0.33212033,  0.27991266,
  0.23289665,  0.19078117,  0.15225619,  0.11899563,  0.08888889,  0.06162057,
  0.03806405,  0.0216885 ,  0.00853954,  0.00111596,  0.        ,  0.,          0.,
  0.        ,  0.        ,  0.        ,  0.        ,  0.        ,  0.,        ]
ContaminationFrac_3CntPerSec = array(ContaminationFrac_3CntPerSec, dtype=float32)

#Sources above 5 counts/sec are contamination
ContaminationFrac_5CntPerSec = \
[ 1.        ,  0.9918001 ,  0.96812227,  0.93088307,  0.88575934,  0.83495876,
  0.77714702,  0.71438622,  0.64679767,  0.57899078,  0.51246967,  0.4455118,
  0.38326055,  0.32421155,  0.26940805,  0.22593401,  0.19005337,  0.15900049,
  0.13044639,  0.10414847,  0.0839641 ,  0.06717613,  0.05245027,  0.03951965,
  0.02833576,  0.01909267,  0.01152353,  0.00572538,  0.00181951,  0.,          0.        ]
ContaminationFrac_5CntPerSec = array(ContaminationFrac_5CntPerSec, dtype=float32)

#SNR Curve from 
#SNR = \
#[  0.00000000e+00,   5.66146991e-05,   2.02771815e-04,   3.97621683e-04,
#   6.11180381e-04,   8.32173566e-04,   1.05247181e-03,   1.26921351e-03,
#   1.48555601e-03,   1.70092750e-03,   1.92657788e-03,   2.14497186e-03,
#   2.36012368e-03,   2.57377489e-03,   2.79387995e-03,   3.01675429e-03,
#   3.23904282e-03,   3.46453069e-03,   3.67814302e-03,   3.88745824e-03,
#   4.09953017e-03,   4.31163423e-03,   4.52686520e-03,   4.74274205e-03,
#   4.95871156e-03,   5.17416978e-03,   5.38811181e-03,   5.59977908e-03,
#   5.81076229e-03,   6.01998763e-03,   6.22723578e-03]
#SNR Curve from run SNRvsCollimatorOpeningAngle --SpinPeriod=15230700
SNR = \
[ 0.           ,   92.278862,       365.95202637,    816.20471191,
  1438.03320312,   2226.06274414,   3174.53515625,   4277.6171875,
  5529.22460938,   6922.97314453,   8452.79589844,   10111.125,       
  11891.06640625,  13785.78613281,  15788.97070312,  17893.11914062,  
  20090.55273438,  22374.40039062,  24735.39648438,  27166.77734375,  
  29662.01367188,  32213.38085938,  34813.9296875,   37455.875,       
  40131.765625,    42834.12890625,  45555.54296875,  48288.60546875,  
  51026.45703125,  53761.87109375,  56487.89453125]

SNR = array(SNR, dtype=float32)

#Form the product curves F^1
SNRxFCurve_InfThres = SNR*ContaminationFrac_InfThresh
SNRxFCurve_1CntPerS = SNR*ContaminationFrac_1CntPerSec
SNRxFCurve_2CntPerS = SNR*ContaminationFrac_2CntPerSec
SNRxFCurve_3CntPerS = SNR*ContaminationFrac_3CntPerSec
SNRxFCurve_5CntPerS = SNR*ContaminationFrac_5CntPerSec

#Form the product curves with sqrt(F1) F1=Fraction of time not contaminated
SNRxFCurve_InfThres_sqrt = SNR*(ContaminationFrac_InfThresh)**.5
SNRxFCurve_1CntPerS_sqrt = SNR*(ContaminationFrac_1CntPerSec)**.5
SNRxFCurve_2CntPerS_sqrt = SNR*(ContaminationFrac_2CntPerSec)**.5
SNRxFCurve_3CntPerS_sqrt = SNR*(ContaminationFrac_3CntPerSec)**.5
SNRxFCurve_5CntPerS_sqrt = SNR*(ContaminationFrac_5CntPerSec)**.5

#Form the product curves with the sqrt(F2) F2=Fraction of time not blocked by earth
#Form the product curves with sqrt(F1) F1=Fraction of time not contaminated
SNRxFCurve_InfThres_sqrt2 =  SNRxFCurve_InfThres_sqrt*sqrt(.5)
SNRxFCurve_1CntPerS_sqrt2 =  SNRxFCurve_1CntPerS_sqrt*sqrt(.5)
SNRxFCurve_2CntPerS_sqrt2 =  SNRxFCurve_2CntPerS_sqrt*sqrt(.5)
SNRxFCurve_3CntPerS_sqrt2 =  SNRxFCurve_3CntPerS_sqrt*sqrt(.5)
SNRxFCurve_5CntPerS_sqrt2 =  SNRxFCurve_5CntPerS_sqrt*sqrt(.5)

#Now
#Plot the product curves
mplot.figure(0)
mplot.clf()
mplot.plot(FOV_ArrayWZero, SNRxFCurve_InfThres, label='No Threshold (140)')
mplot.plot(FOV_ArrayWZero, SNRxFCurve_1CntPerS, label=r'3 Ergs/cm$^{2}$/sr/s/keV (98)')
mplot.plot(FOV_ArrayWZero, SNRxFCurve_2CntPerS, label=r'6 Ergs/cm$^{2}$/sr/s/keV (48)')
mplot.plot(FOV_ArrayWZero, SNRxFCurve_3CntPerS, label=r'9 Ergs/cm$^{2}$/sr/s/keV (24)')
mplot.plot(FOV_ArrayWZero, SNRxFCurve_5CntPerS, label=r'15 Ergs/cm$^{2}$/sr/s/keV (12)')
mplot.legend()
mplot.xlabel('Collimator Opening Half-Angle (Degrees)')
mplot.ylabel('Product of SNR and Fraction of Uncontaminated Observing Time')
mplot.savefig('./Figures/SignalToNoise_ContaminationFractionProduct.png')

#Plot the product curves with the sqrt
mplot.figure(1)
mplot.clf()
mplot.plot(FOV_ArrayWZero, SNRxFCurve_InfThres_sqrt, label='No Threshold (140)')
mplot.plot(FOV_ArrayWZero, SNRxFCurve_1CntPerS_sqrt, label=r'3 Ergs/cm$^{2}$/sr/s/keV (98)')
mplot.plot(FOV_ArrayWZero, SNRxFCurve_2CntPerS_sqrt, label=r'6 Ergs/cm$^{2}$/sr/s/keV (48)')
mplot.plot(FOV_ArrayWZero, SNRxFCurve_3CntPerS_sqrt, label=r'9 Ergs/cm$^{2}$/sr/s/keV (24)')
mplot.plot(FOV_ArrayWZero, SNRxFCurve_5CntPerS_sqrt, label=r'15 Ergs/cm$^{2}$/sr/s/keV (12)')
mplot.legend()
mplot.xlabel('Collimator Opening Half-Angle (Degrees)')
mplot.ylabel('Product of SNR and Fraction of Uncontaminated Observing Time')
mplot.savefig('./Figures/SignalToNoise_SqrtContaminationFractionProduct.png') 

#Plot the product curves for the F1*F2
mplot.figure(2)
mplot.clf()
mplot.plot(FOV_ArrayWZero, SNRxFCurve_InfThres_sqrt, label='No Threshold (140)')
mplot.plot(FOV_ArrayWZero, SNRxFCurve_1CntPerS_sqrt, label=r'3 Ergs/cm$^{2}$/sr/s/keV (98)')
mplot.legend()
mplot.xlabel('Collimator Opening Half-Angle (Degrees)')
mplot.ylabel('Product of SNR and Fraction of Uncontaminated Observing Time')
mplot.savefig('./Figures/SignalToNoise_SqrtContaminationFractionProduct.png')
 
