import numpy as np

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

#Calculate signal to noise given magnitude and integration time

mj = 21.8
mh = 21.3
mk = 20.1

t = 60.0*60.0

signalj = 10**(-0.4*(mj-10))*Nsj*t*g
signalh = 10**(-0.4*(mh-10))*Nsh*t*g
signalk = 10**(-0.4*(mk-10))*Nsk*t*g

noisej = np.sqrt(signalj + Nbj*t*g*Npix*Asq + rn*rn*Npix*Asq)
noiseh = np.sqrt(signalh + Nbh*t*g*Npix*Asq + rn*rn*Npix*Asq)
noisek = np.sqrt(signalk + Nbk*t*g*Npix*Asq + rn*rn*Npix*Asq)

print signalj/noisej
print signalh/noiseh
print signalk/noisek

#Calculate integration time based on signal to noise and magnitude

mj = 20.8
mh = 20.3
mk = 20.1

s2n = 10.0

C =s2n*s2n*(rn*rn*Npix*Asq)
B = Nbj*g*Npix*Asq*s2n*s2n + 10**(-0.4*(mj-10))*Nsj*g
A = -(10**(-0.4*(mj-10))*Nsj*g)**2.0
timej = (-B - np.sqrt(B*B-4*A*C))/(2*A)

C =s2n*s2n*(rn*rn*Npix*Asq)
B = Nbh*g*Npix*Asq*s2n*s2n + 10**(-0.4*(mh-10))*Nsh*g
A = -(10**(-0.4*(mh-10))*Nsh*g)**2.0
timeh = (-B - np.sqrt(B*B-4*A*C))/(2*A)

C =s2n*s2n*(rn*rn*Npix*Asq)
B = Nbk*g*Npix*Asq*s2n*s2n + 10**(-0.4*(mk-10))*Nsk*g
A = -(10**(-0.4*(mk-10))*Nsk*g)**2.0
timek = (-B - np.sqrt(B*B-4*A*C))/(2*A)

print timej/60.0
print timeh/60.0
print timek/60.0

#calculate magnitude based on signal to noise and integration time

s2n = 10.0
t = 15.0*60.0

A = -1.
B = s2n*s2n
C = s2n*s2n*(Nbj*g*t*Npix*Asq + rn*rn*Npix*Asq)
signalj=(-B - np.sqrt(B*B-4*A*C))/(2*A)
mj = np.log10(signalj/(Nsj*t*g))/(-0.4) + 10.0

C = s2n*s2n*(Nbh*g*t*Npix*Asq + rn*rn*Npix*Asq)
signalh=(-B - np.sqrt(B*B-4*A*C))/(2*A)
mh = np.log10(signalh/(Nsh*t*g))/(-0.4) + 10.0

C = s2n*s2n*(Nbk*g*t*Npix*Asq + rn*rn*Npix*Asq)
signalk=(-B - np.sqrt(B*B-4*A*C))/(2*A)
mk = np.log10(signalk/(Nsk*t*g))/(-0.4) + 10.0

print mj
print mh
print mk