pro sig_noise,sn,debug=debug

;Purpose: to determine the signal/noise ratio

;Calling Sequence: sig_noise 

;plots the array of redshifts against the array of radii, and
;determines the signal/noise ratio in every bin

;lighter squares = higher signal/noise ratio

window,0,xsize=800,ysize=600

vdisp = 200.0
area_qso = 8000.0

numz = 12
numr = 10 

zmin = .01
zmax = .3

rmin = .5
rmax = 8.0

deltaZ = (zmax - zmin)/(numz - 1.0)

z_min = findgen(numz)*deltaZ+zmin
z_max = z_min + deltaZ

grid_zmin = replicate(1.0,numr)##z_min
grid_zmax = replicate(1.0,numr)##z_max

deltaR = (rmax - rmin)/(numr - 1)

r_min = findgen(numr)*deltaR+rmin
r_max = deltaR + r_min

grid_rmin = r_min##replicate(1.0,numz)
grid_rmax = r_max##replicate(1.0,numz)

n = n_pairs(area_qso,grid_rmin,grid_rmax,grid_zmin,grid_zmax)

noise = .3/sqrt(n)

grid_z = (grid_zmin + grid_zmax)/2.0
grid_r = (grid_rmin + grid_rmax)/2.0

ds = isothermal_sigma(grid_r,vdisp)
invcrit = inv_sig_crit(grid_z)
shear = ds*invcrit

SN = shear/noise
SNtot=sqrt(total(sn^2))
Err_vdisp=vdisp/SNtot

if keyword_set(debug) then stop

print,'Total Sig/Noise : ',sqrt(total(sn^2))
print,'for Vdisp = ',vdisp
print,'Err(vdisp) = ',Err_vdisp

!p.color = !white

tvim2,sn,title='Signal/Noise Ratio',xtitle='Redshift',ytitle='Radius',xscale=xscale,yscale=yscale

return 
end