#!/usr/bin/env python

import numpy as np
import pylab as py
from cosmology import Cosmology

# 2011 Oct 5, 25; Nov 15, 21; Dec 8, 12; 2012 Jan 8
datelabels=["2011 Oct 5", "Oct 25", "Nov 15", "Nov 21", "Dec 8", "Dec 12", "2012 Jan 8"]
obsmjd=[ 55839, 55859, 55880, 55886, 55903, 55907, 55934 ]

for o, s in zip(obsmjd, datelabels):
   py.axvline(o,linewidth=5, color='k')
   py.annotate(s, (o, 18.0), xycoords='data',
            xytext=(-20,0), textcoords='offset points', rotation=90,
            fontsize='x-large',
            arrowprops=None)

# py.ylim([0,1])
# Starting Sep 15, 2011
py.xlim([55819,55941])
py.xlabel('MJD (55819 = Sep 15, 2011)')

# Read in template and overplot using Oct 1 as the date of maximum.
file='JHK_template.txt'
(phase, jmag, jmagerr, hmag, hmagerr, kmag, kmagerr) = np.loadtxt(file, unpack=True)

# Magnitude of a Type Ia SNeIa
# Wood-Vasey08 (disagrees a bit with Krisciunas05)
absmagJ=-18.29 
absmagH=-18.08  
absmagK=-18.32

c = Cosmology()

file = 'oct25_z_phase.txt'
#If there is no redshift, I assume 0.04.
#If there is no phase, I assume phase = 15d on the discovery date
z0, phaseSN = np.loadtxt(file,usecols=(1,2), unpack=True)
name = np.loadtxt(file, unpack=True, usecols=(0,0), dtype=('S'))

outputfile = open('EstMag_EstTime.txt', "w")
for i in range(0,len(z0)):
   J = jmag[phaseSN[i]]+absmagJ+c.DistMod(z0[i])
   H = hmag[phaseSN[i]]+absmagH+c.DistMod(z0[i])
   K = kmag[phaseSN[i]]+absmagK+c.DistMod(z0[i])
   
   #Calculate integration time based on signal to noise and magnitude
   Nsj = 183000. #ADU/s
   Nsh = 195000. #ADU/s
   Nsk = 109000. #ADU/s
   
   Nbj =  8.00   #ADU/s-pixel
   Nbh =  25.0    #ADU/s-pixel
   Nbk =  70.0    #ADU/s-pixel
   
   g = 3.4     #e/ADU
   Npix = 78.5    #pix
   Asq =1.6*1.6 #arcsec^2
   rn = 20.0     #e
   pixscale = 0.0986 #arcsec
   
   s2n = 10.0
   
   C =s2n*s2n*(rn*rn*Npix*Asq)
   B = Nbj*g*Npix*Asq*s2n*s2n + 10**(-0.4*(J-10))*Nsj*g
   A = -(10**(-0.4*(J-10))*Nsj*g)**2.0
   timej = (-B - np.sqrt(B*B-4*A*C))/(2*A)/60.0
   
   C =s2n*s2n*(rn*rn*Npix*Asq)
   B = Nbh*g*Npix*Asq*s2n*s2n + 10**(-0.4*(H-10))*Nsh*g
   A = -(10**(-0.4*(H-10))*Nsh*g)**2.0
   timeh = (-B - np.sqrt(B*B-4*A*C))/(2*A)/60.0
   
   C =s2n*s2n*(rn*rn*Npix*Asq)
   B = Nbk*g*Npix*Asq*s2n*s2n + 10**(-0.4*(K-10))*Nsk*g
   A = -(10**(-0.4*(K-10))*Nsk*g)**2.0
   timek = (-B - np.sqrt(B*B-4*A*C))/(2*A)/60.0
   
   s = name[0][i] + " %5.4f" %z0[i] + " %5.0f" %phaseSN[i] +" %5.2f" %J + " %5.2f" %H + " %5.2f" %K + " %5.2f" %timej +" %5.2f" %timeh +" %5.2f" %timek +"\n"
   outputfile.write(s)

outputfile.close()