;+
; NAME:
;	NGAO_BKG
;
; PURPOSE:
;	Model thermal background and point source sensitivity of the
;	Keck NGAO system.
;
; EXPLANATION:
;	For a given optical configuration and temperatures, computes the
;	transmission and emissivity spectra at any point in an optical system.
;	Includes both on-axis and off-axis (thermal) components.  Accesses
;	a library of transmission and reflectivity data at:
;		/Users/abouchez/ngao/bkg_data
;
; NOTES:
;	All surface coatings should be defined in terms of their reflectivity
;	(coatings used in transmission will use trans. = 1 - refl.).
;
; OPTIONAL INPUT KEYWORDS:
;	design = String defining the optical design to be modeled.  Currently
;		supported options are:
;			'K2AO', 'SPLIT-NIRI', 'SPLIT-DIFU',
;			'NGAO-AM2','NGAO-OAPIR','P3K'
;
;	band = String defining wavelength band to be modeled.  Supported
;		bands are listed in the file:  'ngao_bkg_filters.dat'
;
;	res = Scalar defining spectral resolution at which to compute
;		background and transmission spectra.
;
;	offaxis = Scalar defining solid angle of off-axis thermal radiation
;		impinging on the science detector (in steradians).
;
;	sampwvl = Sampling wavelength, in m.  If defined, will output
;		the ratio of background to sky background at that wavelength.
;
;	aotemp = Override the default AO temperature, in K.
;		Note that you may also define surface temperatures individually
;		in the "optics" definition file (eg. '*_optics_k2ao.dat').
;
;	dirt = Added emissivity to simulate imprefect or dirty coatings.
;		If scalar, added to all optics.  If 2-element vector, first
;		is added to telescope, second to AO system optics.  Default
;		for K2AO system is dirt=0.008
;
;	/verbose Set this flag for maximum textual output.
;
;	/save	Save PNG versions of plots to 'pdir' directory, defined below.
;
; OUTPUT:
;	dat	Structure containing emission and transmission spectra
;		after each surface in the optical system.  For a case with
;		464 elements in the spectral sampling (band='K', res=2e3):
;
;			WVL             FLOAT     Array[464]
;			DWVL            FLOAT       1.10000e-09
;			FILT            STRUCT    -> <Anonymous> Array[1]
;			FTRAN           FLOAT     Array[464]
;			OBS             STRUCT    -> <Anonymous> Array[1]
;			TEL             STRUCT    -> <Anonymous> Array[1]
;			AO              STRUCT    -> <Anonymous> Array[1]
;			OPT             STRUCT    -> <Anonymous> Array[15]
;			SKY_EMIS        FLOAT     Array[464]
;			SKY_TRAN        FLOAT     Array[464]
;			EMIS            FLOAT     Array[464, 15]
;			OEMIS           FLOAT     Array[464, 15]
;			TRAN            FLOAT     Array[464, 15]
;			TOT_EMIS        FLOAT     Array[464]
;			TOT_TRAN        FLOAT     Array[464]
;
; CALLING SEQUENCE:
;	NGAO_BKG [, design=, band=, res=, offaxis=]
;
; PROCEDURES USED:
;	RESAMPLE
;
; EXAMPLE
;	Plot the background and transmission spectra and band-averaged
;	values at the science focal plane for the point design NGAO system.
;		NGAO_BKG, band='H'
;
; MODIFICATION HISTORY:
;	Original version writen in Jan. 2007; A. Bouchez, Caltech Optical Obs.
;	Added PNG plot output; Jul. 2007, AHB,
;	Released Version 1.0 on NGAO TWiki; Aug. 2007, AHB,
;	Added ability to define temps. of individual optics; 10/08 AHB.
;-
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
PRO	NGAO_BKG, dat, design=design, band=band, res=res, $
  aotemp=aotemp, offaxis=offaxis, sampwvl=sampwvl, dirt=dirt, $
  verbose=verbose, save=save

  if not KEYWORD_SET(design) then design = 'K2AO'
  if not KEYWORD_SET(band) then band = 'K'
  if not KEYWORD_SET(res) then res = 2000.
  if not KEYWORD_SET(dirt) then dirt = 0.
  if not KEYWORD_SET(offaxis) then offaxis = 0.		; ster.
  if not KEYWORD_SET(verbose) then verbose = 0
  if not KEYWORD_SET(save) then save = 0
  inst_tran = 0.90			; instrument transm. to detector
  tim_reads = 1.4			; detector readout time (s)
  max_reads = 64			; max number of non-destructive reads
  ho_wfe = 170e-9			; hight order wavefront error (m)
  r0 = 0.18				; at 0.5um

  DEVICE,decompose=0
  TEK_COLOR

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; DEFINE CONSTANTS, BANDPASSES, OPTICAL MATERIALS, TELESCOPE AND AO SYSTEM

;;; Constants
  ddir = '/Users/abouchez/ngao/bkg_data/'
  pdir = '/Users/abouchez/ngao/bkg_model/'
  c = 2.99792458d8		; m s-1
  h = 6.6260755d-34		; J s
  k = 1.380658d-23		; J k-1

;;; Atmosphere, telescope, and AO parameters
  case STRUPCASE(STRMID(design,0,1)) of
    'P':begin
      if not KEYWORD_SET(aotemp) then aotemp = 288.0
      obs = {id:'Palomar', za:0.0, mphase:0.0}
      tel = {id:'Hale', diam:5.10, t:288.0, area:17.8, pscl:381e-6}
      ao =  {id:design, t:aotemp}
    end
    'T':begin						; Needs Checking!
      if not KEYWORD_SET(aotemp) then aotemp = 243.0
      obs = {id:'TMT', za:0.0, mphase:0.0}
      tel = {id:'TMT', diam:30.0, t:275.6, area:650., pscl:4e-3}
      ao =  {id:design, t:aotemp}
    end
    else:begin
      if not KEYWORD_SET(aotemp) then aotemp = 278.0
      obs = {id:'Mauna Kea', za:0.0, mphase:0.0}
      tel = {id:'Keck 2', diam:10.95, t:275.6, area:79.0, pscl:1.375e-3}
      ao =  {id:design, t:aotemp}
    end
  endcase

;;; Read filter data
  fname = 'ngao_bkg_filters.dat'
  filt_struc = {id:'', cwvl:0., fwhm:0., pflx:0., sky:0., ext:0., $
      file:'', skip:0, wunit:0., funit:0.}
  fmt = '(A3,X,3(E9.3,X),F5.2,X,F4.2,X,A24,X,I2.2,X,2(E7.1,X))'
  filt = READSTRUC(ddir+fname, filt_struc, format=fmt, skip=2)

;;; Read optical materials data
  fname = 'ngao_bkg_materials.dat'
  mat_struc = {id:'',  vis:0., nir:0., file:'', skip:0, wunit:0., funit:0.}
  fmt = '(A12,X,2(E8.2,X),A24,X,I2.2,X,2(E7.1,X))'
  mat = READSTRUC(ddir+fname, mat_struc, format=fmt, skip=2)

;;; Read atmospheric emission and absorption spectra
  fname = 'ngao_bkg_atmspec.dat'
  atm_struc = {id:'', am:0., mphase:0., h2o:0., min:0., max:0., $
      file:'', skip:0L, wunit:0., funit:0.}
  fmt = '(A12,X,2(F4.2,X),3(E8.2,X),A25,X,I2.2,X,2(E7.1,X))'
  atm = READSTRUC(ddir+fname, atm_struc, format=fmt, skip=2)

;;; Read optical design data (based on 'design' keyword)
  fname = 'ngao_bkg_optics_' + STRLOWCASE(design) + '.dat'
  opt_struc = {id:'', mat:'', thk:0., t:''}
  fmt = '(A12,X,A12,X,E9.2,X,A12,X)'
  opt = READSTRUC(ddir+fname, opt_struc, format=fmt, skip=2)
  if SIZE(opt,/type) ne 8 then RETURN
  MESSAGE, /info, 'Optical design: ' + STRUPCASE(design)

;;; Read detector parameters
  fname = 'ngao_bkg_detectors.dat'
  det_struc = {id:'', minwvl:0., maxwvl:0., rsig:0., nreads:0L, $
    dcur:0., qe:0., wdepth:0.}
  fmt = '(A12,X,2(E8.2,X),F4.1,X,I1.1,X,F5.3,X,F4.2,X,E8.2)'
  det = READSTRUC(ddir+fname, det_struc, format=fmt, skip=2)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; LOAD FILTER DATA AND SET UP OUTPUT DATA STRUCTURE

;;; Identify filter of interest and set up data structure
  wf = WHERE(STRTRIM(filt.id,2) eq STRUPCASE(band))
  if wf[0] ne -1 then begin
    wf = wf[0]
    dwvl = filt[wf].cwvl / res
    wvl_mm = filt[wf].cwvl + filt[wf].fwhm * [-0.75, 0.75]
    nwvl = (1.5 * filt[wf].fwhm / dwvl) + 1
    wvl = filt[wf].cwvl - 0.75 * filt[wf].fwhm + FINDGEN(nwvl) * dwvl
  endif else begin
    MESSAGE, /info, "Band '" + band + "' not defined."
    RETURN
  endelse

;;; Set up output data structure
  nopt = N_ELEMENTS(opt)		; number of optical elements
  dat = { wvl:		wvl, $			; m
          dwvl:		dwvl, $			; m
          filt:		filt[wf], $		; filter data
          ftran:	FLTARR(nwvl), $		; filter transmission
          obs:		obs, $			; observatory data
          tel:		tel, $			; telescope data
          ao:		ao, $			; AO system data
          opt:		opt, $			; optical design data
          sky_emis:	FLTARR(nwvl), $		; ph s-1 m-3 arcsec-2
          sky_tran:	FLTARR(nwvl), $		; fractional
          emis:		FLTARR(nwvl,nopt), $	; ph s-1 m-3 arcsec-2
          oemis:	FLTARR(nwvl,nopt), $	; " (optics only)
          aemis:	FLTARR(nwvl,nopt), $	; " (AO optics only)
          tran:		FLTARR(nwvl,nopt), $	; fractional
          otran:	FLTARR(nwvl,nopt), $	; fractional (optics only)
          tot_emis:	FLTARR(nwvl), $		; includes off-axis light
          tot_tran:	FLTARR(nwvl)}		; includes filter

;;; Read filter curve data
  fname = STRTRIM(filt[wf].file, 2)
  if STRLEN(fname) gt 0 then begin
    MESSAGE, /info, 'Reading filter data: ' + fname
    READCOL, ddir+fname, c0, c1, skip=filt[wf].skip, /silent
    ww = WHERE((c0*filt[wf].wunit ge MIN(wvl)) and $
               (c0*filt[wf].wunit le MAX(wvl)))
    if ww[0] ne -1 then begin
      mtxt = STRING(dwvl*1e9,f='(F4.1)') + ' nm'
      if MIN(ABS(DERIV(c0[ww]*filt[wf].wunit))) lt dwvl then begin
        if verbose then MESSAGE, /info, '  Resampling to ' + mtxt
        dat.ftran = $
          RESAMPLE(c1[ww]*filt[wf].funit, c0[ww]*filt[wf].wunit, wvl) > 0
      endif else begin
        if verbose then MESSAGE, /info, '  Interpolating to ' + mtxt
        dat.ftran = $
          INTERPOL(c1[ww]*filt[wf].funit, c0[ww]*filt[wf].wunit, wvl) > 0
      endelse
    endif else $
      MESSAGE, /info, '  No valid points in spectrum!'
  endif

;;; Compute equivalent rectangular filter function
  if total(dat.ftran) eq 0 then begin
    MESSAGE, /info, 'Using rectangular filter approximation.'
    ww = WHERE((wvl ge filt[wf].cwvl-0.5*filt[wf].fwhm) and $
               (wvl lt filt[wf].cwvl+0.5*filt[wf].fwhm))
    dat.ftran[ww] = 1.0
  endif

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; LOAD SKY EMISSION AND TRANSMISSION SPECTRA

;;; Resample sky emission spectrum to output wavelengths
  wa = WHERE((STRPOS(atm.id,'emis') ne -1) and $
             (atm.min le filt[wf].cwvl-filt[wf].fwhm/2) and $
             (atm.max ge filt[wf].cwvl+filt[wf].fwhm/2))
  if wa[0] ne -1 then begin
    wa = wa[0]
    fname = STRTRIM(atm[wa].file, 2)
    if STRLEN(fname) gt 0 then begin
      MESSAGE, /info, 'Reading atm. emission spectrum: ' + fname
      READCOL, ddir+fname, c0, c1, skip=atm[wa].skip, /silent
      ww = WHERE((c0*atm[wa].wunit ge MIN(wvl)) and $
                 (c0*atm[wa].wunit le MAX(wvl)))
      if ww[0] ne -1 then begin
        if verbose then MESSAGE, /info, $
          '  Resampling to ' + STRING(dwvl*1e9,f='(F4.1)') + ' nm'
        emis = (c1[ww] * atm[wa].funit / atm[wa].am) / COS(obs.za)
        dat.sky_emis = RESAMPLE(emis, c0[ww]*atm[wa].wunit, wvl)
      endif
    endif
  endif

;;; ... or approximate from band-averaged sky brightness.
  if TOTAL(dat.sky_emis) eq 0 then begin
    MESSAGE, /info, 'Using band-averaged atm. emission.'
    dat.sky_emis = dat.dwvl * filt[wf].pflx * 10^(-0.4*filt[wf].sky)
  endif

;;; Resample sky transmission spectrum to output wavelengths
  wa = WHERE((STRPOS(atm.id,'tran') ne -1) and $
             (atm.min le filt[wf].cwvl-filt[wf].fwhm/2) and $
             (atm.max ge filt[wf].cwvl+filt[wf].fwhm/2))
  if wa[0] ne -1 then begin
    wa = wa[0]
    fname = STRTRIM(atm[wa].file, 2)
    if STRLEN(fname) gt 0 then begin
      MESSAGE, /info, 'Reading atm. transmission spectrum: ' + fname
      READCOL, ddir+fname, c0, c1, skip=atm[wa].skip, /silent
      if atm[wa].funit gt 0 then $
        tran = (c1 * atm[wa].funit / atm[wa].am) / COS(obs.za) else $
        tran = (10^(-0.4*c1) / atm[wa].am) / COS(obs.za)
      if verbose then MESSAGE, /info, $
        '  Resampling to ' + STRING(dwvl*1e9,f='(F4.1)') + ' nm'
      if STRUPCASE(filt[wf].cwvl) lt 1e-6 then $
        dat.sky_tran = INTERPOL(tran, c0*atm[wa].wunit, dat.wvl) else $
        dat.sky_tran = RESAMPLE(tran, c0*atm[wa].wunit, dat.wvl)
    endif
  endif

;;; ... or approximate from band-averaged extinction.
  if TOTAL(dat.sky_tran) eq 0 then begin
    MESSAGE, /info, 'Using band-averaged atm. transmission.'
    dat.sky_tran = 10^(-0.4*filt[wf].ext)
  endif

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; ADD TRANSMISSION AND EMISSIVITY EFFECTS OF EACH SURFACE

  nopt = (SIZE(opt))[1]
  for i=0,nopt-1 do begin

;;; Read optical transmission/reflectivity data for surface
    mtxt = 'Surface ' + STRING(i+1,f='(I2.2)') + ':' + opt[i].id
    wm = WHERE(STRTRIM(mat.id,2) eq STRTRIM(STRUPCASE(opt[i].mat),2))
    if wm[0] ne -1 then begin
      fname = STRTRIM(mat[wm].file, 2)
      if STRLEN(fname) ne 0 then begin
        mtxt = mtxt + '; Using material spectrum ' + fname
        if verbose then MESSAGE, /info, mtxt
        READCOL, ddir+fname, c0, c1, skip=mat[wm].skip, /silent
        tmp = INTERPOL(c1 * mat[wm].funit, c0 * mat[wm].wunit, dat.wvl)
      endif else begin

;;; If not available, use Vis/IR average values
        if filt[wf].cwvl lt 1e-6 then $
          tmp = mat[wm].vis else $
          tmp = mat[wm].nir
        mtxt = mtxt + '; Using avg. trans. =' + STRING(tmp,f='(F6.3)')
        if verbose then MESSAGE, /info, mtxt
      endelse
    endif else begin
      mtxt = mtxt + '; Material undefined.  Ignoring.'
      if verbose then MESSAGE, /info, mtxt
      tmp = 1.0
    endelse

;;; Determine transmission spectrum for surface
    if opt[i].thk lt 0 then tran = tmp				; reflective
    if opt[i].thk eq 0 then tran = 1 - tmp			; transmissive
    if opt[i].thk gt 0 then tran = 1 - (tmp * opt[i].thk)	; bulk

;;;; ADJUSTED EMISSIVITY TO MATCH NIRC2
    if (opt[i].thk le 0) then begin
      if (i le 2) then begin
        if verbose then $
          MESSAGE, /info, 'Adding e=' + STRING(dirt[0],f='(F5.3)')
        tran = tran - dirt[0]
      endif else begin
        if verbose then $
          MESSAGE, /info, 'Adding e=' + $
            STRING(dirt[N_ELEMENTS(dirt)-1],f='(F5.3)')
        tran = tran - dirt[N_ELEMENTS(dirt)-1]
      endelse
    endif

;;; Compute emission spectrum
    if not ISNUMBER(opt[i].t) then begin
      case STRTRIM(opt[i].t, 2) of
        'AO':temp = ao.t
        'TEL':temp = tel.t
        else:begin
          MESSAGE, /info, 'Temperature zone '+opt[i].t+' undefined!'
          temp = 0.0
        end
      endcase
    endif else $
      temp = FLOAT(opt[i].t)

    e = 1 - tran
    fdens = ((2 * c^2 * h) / dat.wvl^5) * $		; W  m-3 sr-1
      (1 / (EXP((h * c) / (dat.wvl * k * temp)) - 1))
    pflx = e * fdens * dat.wvl / (h * c)		; ph s-1 m-3 sr-1
    emis = pflx / 206265.^2				; ph s-1 m-3 arcsec-2

;;; Attenuate incoming beam and add thermal emission

    if (i eq 0) then begin
      dat.tran[*,i] = dat.sky_tran * tran
      dat.otran[*,i] = tran
      dat.emis[*,i] = dat.sky_emis + emis
      dat.oemis[*,i] = emis
    endif else begin
      dat.tran[*,i] = dat.tran[*,i-1] * tran
      dat.otran[*,i] = dat.otran[*,i-1] * tran
      dat.emis[*,i] = dat.emis[*,i-1] * tran + emis
      dat.oemis[*,i] = dat.oemis[*,i-1] * tran + emis
      dat.aemis[*,i] = dat.aemis[*,i-1] * tran + emis
    endelse
    if (i eq 3) then begin
      dat.aemis[*,0:2] = 0
      dat.aemis[*,i] = emis
    endif
  endfor

;;; Add off-axis thermal emission and filter transmission

  fdens = ((2 * c^2 * h) / dat.wvl^5) * $		; W  m-3 sr-1
    (1 / (EXP((h * c) / (dat.wvl * k * ao.t)) - 1))
  pflx0 = 1.0 * fdens * dat.wvl / (h * c)		; ph s-1 m-3 sr-1
  pflx1 = pflx0 * tel.pscl^2 * offaxis / tel.area	; ph s-1 m-3 arcsec-2

  dat.tot_tran = dat.tran[*,nopt-1]
  dat.tot_emis = (dat.emis[*,nopt-1] + pflx1)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; COMPUTE BAND-AVERAGE EMISSION AND TRANSMISSION

  if filt[wf].pflx ne 0 then begin			; for single bands...
    sky_emis = dat.sky_emis * dat.ftran * $
      dat.filt.fwhm / TOTAL(dat.ftran * dat.dwvl)
    sky_mag = -2.5 * ALOG10( TOTAL(sky_emis * dat.dwvl) / filt[wf].pflx )
    MESSAGE, /info, STRUPCASE(band) + ' sky magnitude:       ' + $
      STRING(sky_mag,f='(F5.2)')

    avg_tran = TOTAL(dat.otran[*,nopt-1] * dat.ftran)/TOTAL(dat.ftran)
    avg_tran_tel = TOTAL(dat.otran[*,2] * dat.ftran)/TOTAL(dat.ftran)
    avg_tran_ao = avg_tran / avg_tran_tel

    bkg_emis = dat.tot_emis * dat.ftran * $
      dat.filt.fwhm / (TOTAL(dat.ftran * dat.dwvl) * avg_tran)
    bkg_mag = -2.5 * ALOG10( TOTAL(bkg_emis * dat.dwvl) / filt[wf].pflx )
    MESSAGE, /info, STRUPCASE(band) + ' background after AO: ' + $
      STRING(bkg_mag,f='(F5.2)')
    MESSAGE, /info, 'Avg. TEL trans.:        ' + STRING(avg_tran_tel,f='(F5.3)')
    MESSAGE, /info, 'Avg. AO trans.:         ' + STRING(avg_tran_ao,f='(F5.3)')
    MESSAGE, /info, 'Avg. Total trans.:      ' + STRING(avg_tran,f='(F5.3)')
  endif

  if KEYWORD_SET(sampwvl) then begin
    wsamp = WHERE(ABS(dat.wvl - sampwvl)/dat.dwvl lt 1.0)
    if wsamp[0] ne -1 then begin
      wsamp = wsamp[0]
      bkgr = dat.tot_emis[wsamp] / dat.sky_emis[wsamp]
      MESSAGE, /info, 'BKG/SKY at ' + STRING(sampwvl*1e9,f='(I0.0)') + $
        ' nm:   ' + STRING(bkgr,f='(F6.2)')
      bkgt = dat.tot_emis[wsamp] / dat.emis[wsamp,2]
      MESSAGE, /info, 'BKG/TEL at ' + STRING(sampwvl*1e9,f='(I0.0)') + $
        ' nm:   ' + STRING(bkgt,f='(F6.2)')
      MESSAGE, /info, 'I(BKG) ' + STRING(sampwvl*1e9,f='(I0.0)') + $
        ' nm:       ' + STRING(dat.tot_emis[wsamp]*1e-9,f='(F6.2)') + ' phot'
    endif else $
      MESSAGE, /info, 'Sample wavelength outside of band!'
  endif

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; COMPUTE POINT SOURCE SENSITIVITY IN BAND.  INCOMPLETE!!!

  if filt[wf].pflx ne 0 then begin			; for single bands...
    wd = WHERE((det.minwvl lt filt[wf].cwvl) and $
               (det.maxwvl gt filt[wf].cwvl))
    if wd[0] ne -1 then begin
      det = det[wd[0]]

;;; Determine max integration time for background, at 2*Nyquist pixel size

      pscl = 206265 * filt[wf].cwvl / (2 * tel.diam)		; arcsec/pix
;      MESSAGE, /info, 'Platescale:           ' + $
;        STRING(pscl*1e3,f='(F5.1," mas/pix")')
      bkg_emis = dat.tot_emis * dat.ftran * $		; ph s-1 m-2 "-2 nm-1
        dat.filt.fwhm / TOTAL(dat.ftran * dat.dwvl)
      bkg_flx = TOTAL(bkg_emis * dat.dwvl) * $		; ph s-1 pix-1
        tel.area * pscl^2 * inst_tran
      maxint = 0.80 * det.wdepth / bkg_flx			; s (80% full)
      nfr = ROUND(3600. / maxint)				; no. frames
;      MESSAGE, /info, 'Integration:           ' + $
;        STRING([nfr,maxint],f='(I0.0," x ",F5.1," s")')

;;; Compute light distribution in halo/core

      strehl = EXP(-1*(2*!pi*ho_wfe/filt[wf].cwvl)^2)
;      MESSAGE, /info, 'Strehl:                 ' + $
;        STRING(strehl,f='(F5.3)')
      r0_wvl = r0 * (filt[wf].cwvl/5.0e-7)^(6/5.)	; m
      hlo_fwhm = 206265 * filt[wf].cwvl / r0_wvl	; "
      hlo_pix = (!pi * (hlo_fwhm/2)^2) / pscl^2		; pix
      hlo_rat = (1-strehl) / hlo_pix			; pix-1 (star fraction)

;;; Noise sources per pixel, per integration

      mag = 20+FINDGEN(30)/2.
      star_flx = tel.area * filt[wf].pflx * 10^(-0.4*mag)	; ph s-1
      
      if (det.nreads gt 1) then $
        nreads = FIX(maxint/tim_reads) < max_reads else $
        nreads = 1
      rd_sig = det.rsig / SQRT(nreads)
      bkg_sig = SQRT((det.dcur + bkg_flx + star_flx*hlo_rat) * maxint)

;;; Computer SNR in one frame

      tot_sig = SQRT(rd_sig^2 + bkg_sig^2) * SQRT(4.)
      snr_frm = (star_flx * strehl * maxint) / tot_sig
      snr_tot = snr_frm * SQRT(nfr)
;plot,mag,snr_tot,/ylog,yr=[1,10]
;retall

    endif else $
      MESSAGE, /info, 'No relevant detector found'
  endif

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; PLOT RESULTS

  wsz = [600, 400]
  if MAX(STRUPCASE(filt[wf].id) eq ['VIS','NIR','L-M']) eq 1 then $
    wsz = [1000,400]

;;; PLOT BACKGROUND FLUX

  if save then begin
    wr = 2
    WINDOW, /free, xs=wsz[0]*wr, ys=wsz[1]*wr, /pixmap
  endif else begin
    wr = 1
    WINDOW, 12, xs=wsz[0]*wr, ys=wsz[1]*wr
  endelse

  TV, REPLICATE(1B, wsz[0]*wr, wsz[1]*wr)
  if MAX(STRUPCASE(filt[wf].id) eq ['VIS','NIR','L-M']) eq 1 then begin
    pos = [0.08, 0.12, 0.98, 0.98]
    xyx = 0.11
  endif else begin
    pos = [0.10, 0.12, 0.98, 0.98]
    xyx = 0.13
  endelse
  xr = (filt[wf].cwvl + filt[wf].fwhm * [-0.5, 0.5]) * 1e9
  w = WHERE((dat.wvl gt xr[0]*1e-9) and (dat.wvl lt xr[1]*1e-9))
  yr = [0.6,1.2]*MINMAX(dat.sky_emis[w]*1e-9)
  xtxt = 'Wavelength (nm)'
  ytxt = 'photons s!E-1!N m!E-2!N arcsec!E-2!N nm!E-1!N'
  PLOT, dat.wvl*1e9, dat.sky_emis*1e-9, pos=pos, $
    /xsty, xr=xr, xtit=xtxt, $
    /ylog, /ysty, yr=yr, ytit=ytxt, $
    color=0, chars=1.5*wr, xthick=1.5*wr, ythick=1.5*wr, $
    thick=1.0*wr, charthick=1.5*wr, /noerase
  OPLOT, dat.wvl*1e9, $
    dat.oemis[*,2]*(dat.otran[*,nopt-1]/dat.otran[*,2])*1e-9, thick=1.5*wr, co=3
  OPLOT, dat.wvl*1e9, dat.aemis[*,nopt-1]*1e-9, thick=1.5*wr, co=4
  OPLOT, dat.wvl*1e9, dat.tot_emis*1e-9, thick=1.5*wr, co=2

;; Code for plotting ratio of sky+telescope / sky.  Use with design=telonly.
;  yr = [0.6,1.2]*MINMAX(dat.tot_emis/dat.sky_emis)
;  ytxt = 'Sky + Telescope / Sky'
;  PLOT, dat.wvl*1e9, dat.tot_emis/dat.sky_emis, pos=pos, $
;    /xsty, xr=xr, xtit=xtxt, $
;    /ylog, /ysty, yr=yr, ytit=ytxt, $
;    color=0, chars=1.5*wr, xthick=1.5*wr, ythick=1.5*wr, $
;    thick=1.5*wr, charthick=1.5*wr, /noerase

  if KEYWORD_SET(sampwvl) then begin
    if wsamp[0] ne -1 then begin
      OPLOT, sampwvl*1e9*[1,1], $
        [dat.sky_emis[wsamp], dat.tot_emis[wsamp]]*1e-9, $
        thick=2.0*wr, co=4
    endif
  endif

  XYOUTS, xyx, 0.94, 'Black: Atmosphere', /norm, $
    align=0, chars=1.5*wr, charth=1.5*wr, color=0
  XYOUTS, xyx, 0.90, 'Green: Telescope', /norm, $
    align=0, chars=1.5*wr, charth=1.5*wr, color=3
  XYOUTS, xyx, 0.86, 'Blue:  AO system', /norm, $
    align=0, chars=1.5*wr, charth=1.5*wr, color=4
  XYOUTS, xyx, 0.82, 'Red:   Total', /norm, $
    align=0, chars=1.5*wr, charth=1.5*wr, color=2

;;; Save to PNG file

  if save then begin
    im = REBIN(TVRD(/true), 3, wsz[0], wsz[1])
    WDELETE
    WINDOW, 12, xs=wsz[0], ys=wsz[1]
    TV, REBIN(im, 1, wsz[0], wsz[1])

    sname = 'bkg_' + STRLOWCASE(design) + '_' + $
      STRLOWCASE(band) + '_' + STRING(res,f='(I4.4)') + $ 
      '_' + STRING(ao.t,f='(I3.3)') + '.png'
    MESSAGE, /info, 'Writing file: ' + sname 
    WRITE_PNG, pdir+sname, im
  endif

;;; PLOT SYSTEM TRANSMISSION

  if save then $
    WINDOW, /free, xs=wsz[0]*wr, ys=wsz[1]*wr, /pixmap else $
    WINDOW, 13, xs=wsz[0]*wr, ys=wsz[1]*wr

  TV, REPLICATE(1B, wsz[0]*wr, wsz[1]*wr)
  if MAX(STRUPCASE(filt[wf].id) eq ['VIS','NIR','L-M']) eq 1 then begin
    pos = [0.08, 0.12, 0.98, 0.98]
    xyx = 0.71
  endif else begin
    pos = [0.10, 0.12, 0.98, 0.98]
    xyx = 0.63
  endelse
  xr = (filt[wf].cwvl + filt[wf].fwhm * [-0.5, 0.5]) * 1e9
  yr = [0,1]
  xtxt = 'Wavelength (nm)'
  ytxt = 'transmission'
  PLOT, dat.wvl*1e9, dat.sky_tran, pos=pos, $
    /xsty, xr=xr, xtit=xtxt, $
    /ysty, yr=yr, ytit=ytxt, $
    color=0, chars=1.5*wr, xthick=1.5*wr, ythick=1.5*wr, $
    thick=1.5*wr, charthick=1.5*wr, /noerase
  OPLOT, dat.wvl*1e9, dat.tran[*,2] / dat.sky_tran, thick=1.5*wr, co=3
  OPLOT, dat.wvl*1e9, dat.tot_tran / dat.tran[*,2], thick=1.5*wr, co=4
  OPLOT, dat.wvl*1e9, dat.tot_tran, thick=1.5*wr, co=2

  XYOUTS, xyx, 0.27, 'Black: Atmosphere', /norm, $
    align=0, chars=1.5*wr, charth=1.5*wr, color=0
  XYOUTS, xyx, 0.23, 'Green: Telescope', /norm, $
    align=0, chars=1.5*wr, charth=1.5*wr, color=3
  XYOUTS, xyx, 0.19, 'Blue:  AO system', /norm, $
    align=0, chars=1.5*wr, charth=1.5*wr, color=4
  XYOUTS, xyx, 0.15, 'Red:   Total', /norm, $
    align=0, chars=1.5*wr, charth=1.5*wr, color=2

;;; Save to PNG file

  if save then begin
    im = REBIN(TVRD(/true), 3, wsz[0], wsz[1])
    WDELETE
    WINDOW, 13, xs=wsz[0], ys=wsz[1]
    TV, rebin(im, 1, wsz[0], wsz[1])

    sname = 'tra_' + STRLOWCASE(design) + '_' + $
      STRLOWCASE(band) + '_' + STRING(res,f='(I4.4)') + $ 
      '_' + STRING(ao.t,f='(I3.3)') + '.png'
    MESSAGE, /info, 'Writing file: ' + sname 
    WRITE_PNG, pdir+sname, im
  endif

END