set xzeroaxis
set logscale y

set xlabel 'R_{mag}'
set ylabel 'S/N (T_{exp} = 180 s)'

piii = 3.14159265

# SBIG ST-8
gain = 2.3
RON  = 14.0
pixscale = 0.94
mzero = 19.75   # KVA zero point with ST-8     

rap  = 5.0/pixscale
rsky = 16.0/pixscale
wsky = 8.0/pixscale
sky  = 20.6 # Sky brightness mag/sq. arcsec
Texp = 180.0
ndarks = 4

# Number of pixels in the aperture
npix = piii*rap**2

# Number of pixels in the sky ring
nsky = piii*((rsky+wsky)**2-rsky**2)

# Sky (electrons/pix/s)
skyeps = gain*10**((mzero-sky)/2.5)*(pixscale*pixscale)

# Dark electrons/pix/s (ST-8 @-12 deg C)
darkeps = 1.559

# Signal
A(x) = gain*10**((mzero-x)/2.5)*Texp

# Background (sky + dark)
bgn = skyeps+(1.0+1.0/sqrt(ndarks))*darkeps

# Noise
B(x) = sqrt( A(x) + (npix + npix**2/nsky)*(bgn*Texp+RON**2) )

s(x) = A(x)/B(x)

plot [] [] 'sndata_edited.dat' u 1:(1/(10**($3/2.5)-1)) pt 7 title 'Obs. data',\
s(x) title 'Model' lw 2 lc rgb '#ff2222'

pause -1