PRO snr,epdiamcm,ccddiamcm,starmag,ccdwl,ccdsnr,fovdeg $
       ,optau=optau,gratau=gratau,nmperpix=nmperpix $
       ,ccdqe=ccdqe,rnoise=rnoise,dark=dark,show=show,debug=debug

IF N_Params() LT 3 THEN BEGIN
    Print,'PRO snr,epdiamcm,ccddiamcm,starmag,ccdwl,ccdsnr,fovdeg $'
    Print,'       ,optau=optau,gratau=gratau,nmperpix=nmperpix $'
    Print,'       ,ccdqe=ccdqe,rnoise=rnoise,dark=dark,show=show,debug=debug'
    Return
    ENDIF

; Define default parameter values
IF N_Elements(optau) GE 1 THEN optau=optau(0) ELSE optau=0.8
IF N_Elements(gratau) GE 1 THEN gratau=gratau(0) ELSE gratau=0.5
IF N_Elements(nmperpix) GE 1 THEN nmperpix=nmperpix(0) ELSE nmperpix=5.
IF N_Elements(ccdqe) GE 1 THEN ccdqe=ccdqe(0) ELSE ccdqe=0.5
IF N_Elements(rnoise) GE 1 THEN rnoise=rnoise(0) ELSE rnoise=30.
IF N_Elements(dark) GE 1 THEN dark=dark(0) ELSE dark=280.
IF N_Elements(debug) GT 0 THEN debug=1 ELSE debug=0

; First work out appropriate spectrum for specified star magnitude (only
; G7-type for now)

; Load photo-visual response curve from Allen, p. 178
lambda=Findgen(31)*200+3000.    ; 3000 - 10000 A
vlambda=[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.01,0.03 $
        ,0.18,0.67,0.93,0.98,0.28,0.011,0.0,0.0,0.0,0.0 $
        ,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]

; Load star spectrum for G-type star
; Note: spectrum has units of photons/cm^2 s A, and wavelength scale is
; in A.
OpenR,lun,'\usask\893\crossect\wavelgth.txt',/Get_LUN
starwave=DblArr(1084)
ReadF,lun,starwave
Free_LUN,lun
OpenR,lun,'\usask\893\crossect\g1--2v.txt',/Get_LUN
photflux=DblArr(1084)
ReadF,lun,photflux
Free_LUN,lun

; Convert spectrum to units of erg/cm^2 s A
planck=6.626E-34        ; J-s
light=2.99792458E8      ; m/s
ph_energy=(planck*light/(starwave*1E-10))*1E7   ; 1E7 to get ergs, not joules
ergflux=photflux*ph_energy

; Calculate integral of PV response * flux
dwl=Median(starwave-Shift(starwave,1))
pvfit=(Spline(lambda,vlambda,starwave,10)>0.)<1.        ; high 'tension'
pvflux=Int_Tabulated(starwave,pvfit*ergflux)
mpvref=-2.5*Alog10(pvflux)-14.14
Print,'Reference spectrum photovisual magnitude is ',mpvref

; Convert ref. spectrum to desired photo-visual magnitude
convfac=10^(-0.4*(starmag-mpvref))
testflux=convfac*photflux

; Restrict spectrum to 450 - 750 nm region
inrange=Where((4500. LE starwave) AND (starwave LE 7500.),nflux)
starwave=starwave(inrange)
testflux=testflux(inrange)

; Now work out statistics for specified photon flux, entrance aperture,
; and CCD size
entr_aper=!Pi*epdiamcm*epdiamcm/4.
IF Keyword_Set(show) THEN Help,epdiamcm,entr_aper
ph_in=entr_aper*testflux
ph1ord=ph_in*gratau
ph2ccd=ph1ord*optau

; Now integrate spectrum over specified spectral bin size.
nsamples=Round(nmperpix*10./dwl)
nkeep=(N_Elements(ph2ccd)/nsamples)*nsamples
keep=Indgen(nkeep)
phperpix=Total(Reform(ph2ccd(keep),nsamples,nkeep/nsamples),1)*dwl
ccdwl=Total(Reform(starwave(keep),nsamples,nkeep/nsamples),1)/nsamples

; Now determine deposited charge and SNR
eperpix=phperpix*ccdqe
ccdsnr=eperpix/Sqrt(eperpix+dark+rnoise*rnoise)
IF Keyword_Set(show) THEN BEGIN
    Plot_IO,ccdwl,eperpix>10,PSym=10 $
           ,XTitle='Wavelength ('+Angstrom()+')' $
           ,YTitle='Electrons per pixel'
    OPlot,ccdwl,ccdsnr>10,PSym=10,Linestyle=1
    ENDIF

; Determine possible field of view.
fovdeg=2*!RaDeg*(0.5*ccddiamcm/epdiamcm)        ; full field of view

IF debug THEN Stop

Return
END