;; sed.pro - plot eazy results from binary output files
; 
;Usage:
;
;    sed,44
;
; will plot the SED and p(z) distribution of the 45th object (zero index)
; in the photometric catalog provided to EAZY.
;
; Keywords: 
; 
; NOPZ - don't plot p(z) distribution
; ROOT - rootname of binary files ("MAIN_OUTPUT_FILE")
; TEMPFILTF = "CACHE_FILE"
; DIR = "OUTPUT_DIRECTORY"
; FNU - plot SEDs in F_nu rather than F_lambda
; STOPME - stop program after plotting to have access to the variables 
;          read from the binary files
; ZS  - z_spec value, will be overplotted on the p(z) plot for comparison
; CHITOP - convert CHI^2 to a probability
; PRINTSED - print broadband SED fluxes/errors to the screen
; PEAK - Force maximum value of the SED 
; 
; Also accepts any additional keywords passed to the 'plot' command
;
; NOTE: Requires the GSFC IDLASTRO library in your $IDL_PATH
;       (http://idlastro.gsfc.nasa.gov/)
;

pro sed_example
    
    ;; Given the default eazy parameters and setting BINARY_OUTPUTS=y
    ;; in the 'inputs' directory where you executed 'eazy'
    !p.multi=[0,2,1]
    sed,44,root='photz',dir='OUTPUT',/xlog,zs=1.012,xrange=[2.e3,3.e4]
    ;;; everything beyond the '44' is optional
    
    ;; Should produce a plot like  
    ;; http://www.astro.yale.edu/eazy/?sed

end
  
pro sed,idx,nopz=nopz,root=root,tempfiltf=tempfiltf,dir=dir,fnu=fnu,$
    stopme=stopme,zs=zs,chitop=chitop,$
    printsed=printsed,peak=peak,_extra=_extra

	if keyword_set(root) eq 0 then root='photz'
	if keyword_set(dir) eq 0 then dir='OUTPUT'
	if keyword_set(tempfiltf) eq 0 then begin
		spawn,'grep "might" '+dir+'/'+root+'.readbin.pro',tf
		tempfiltf=(strsplit(tf,"'",/extract))[1]
	endif
	
	;;; Default to F_lambda
	funit=1 
	if keyword_set(fnu) then funit=0
	
	
; s[0:3] = NFILT NTEMP NZ NOBJ 
; 16 6 161 6308
; Read with IDL:
; Template & catalog fluxes
	 openr,lun,tempfiltf,/swap_endian,/get_lun ; might not need 'swap_endian'
	 s = lonarr(4) & readu,lun,s
	 NFILT=s[0] & NTEMP = s[1] & NZ = s[2] & NOBJ = s[3]
	 tempfilt = dblarr(NFILT,NTEMP,NZ)
	 lc = dblarr(NFILT) ; central wavelengths
	 zgrid = dblarr(NZ)
	 fnu = dblarr(NFILT,NOBJ)
	 efnu = dblarr(NFILT,NOBJ)
	 readu,lun,tempfilt,lc,zgrid,fnu,efnu
	 close,/all
	; 
	; Coeffs (OBS_SED)
	 openr,lun,dir+'/'+root+'.coeff',/get_lun,/swap_endian
	 s = lonarr(4) & readu,lun,s
	 NFILT=s[0] & NTEMP = s[1] & NZ = s[2] & NOBJ = s[3]
	 coeffs = dblarr(NTEMP,NOBJ)
	 izbest = lonarr(NOBJ)
	 tnorm = dblarr(NTEMP)
	 readu,lun,coeffs,izbest,tnorm 
	 close,/all
	; 
	; Full templates (TEMP_SED)
	 openr,lun,dir+'/'+root+'.temp_sed',/get_lun,/swap_endian
	 s = lonarr(3) & readu,lun,s
	 NTEMP=s[0] & NTEMPL = s[1] & NZ = s[2]
	 templam = dblarr(NTEMPL)
	 temp_seds = dblarr(NTEMPL,NTEMP)
	 da = dblarr(NZ)
	 db = dblarr(NZ)
	 readu,lun,templam,temp_seds,da,db 
	 close,/all
	; 
	; POFZ
	 openr,lun,dir+'/'+root+'.pz',/get_lun,/swap_endian
	 s = lonarr(2) & readu,lun,s
	 NZ = s[0] & NOBJ = s[1]
	 pz = dblarr(NZ,NOBJ)
	 readu,lun,pz 
	 if keyword_set(chitop) then begin
	 	pz = exp(-0.5*pz)
	 	pz/=trapz(zgrid,pz)
	 endif
	 close,/all

	if funit eq 1 then begin
        fconvert = 1./(lc/5500.)^2      
        ytitle = 'F_lam'
	endif else begin
        fconvert = lc*0.+1
	    ytitle = 'F_nu'
	endelse 

    ;; Broadband SED, matrix multiply the filters of 
    ;; each template by the computed normalization coefficients
	OBS_SED = tempfilt[*,*,izbest[idx]] # coeffs[*,idx]

    ;; Full template SED
	TEMP_SED = temp_seds # coeffs[*,idx] ;;; F_lam
	
	if keyword_set(printsed) then forprint,lc,obs_sed,fnu[*,idx],efnu[*,idx]
	
	; Observed frame
	TEMPLAMZ = templam*(1+zgrid[izbest[idx]])

	; IGM absorption
	lim1 = where(templam lt 912)
	lim2 = where(templam ge 912 and templam lt 1026)
	lim3 = where(templam ge 1026 and templam lt 1216)
	
	TEMP_SED[lim1] *= 0
	TEMP_SED[lim2] *= 1-db[izbest[idx]]
	TEMP_SED[lim3] *= 1-da[izbest[idx]]

    if funit eq 0 then begin
            TEMP_SED *= (templamz/5500.)^2/(1+zgrid[izbest[idx]])^2 ;;; F_nu
    endif else begin
            TEMP_SED /= (1+zgrid[izbest[idx]])^2 ;;; F_lambda
    endelse
	
;	!p.multi=[0,2,2]
	plotsym,0,1,/fill
;	window,0,xs=600,ys=600
	loadct,5,/silent
	
	; plot SED+template
	
	ymax = 1.1*max((fnu[*,idx]+efnu[*,idx])*fconvert)
	!x.range=[2500,9.e4]
    xst = !x.style
    !x.style=1
    yst = !y.style
    !y.style=1
    
	!y.range=[-0.1*ymax < min((fnu[*,idx]-efnu[*,idx])*fconvert) < 0,ymax]
	
	if keyword_set(peak) then peak = peak/max(fnu[*,idx]*fconvert[*]) else peak = 1
	
	ploterror,lc,fnu[*,idx]*fconvert[*]*peak,efnu[*,idx]*fconvert*peak,psym=8,$
	    xtitle='wavelength',ytitle=ytitle,_extra=_extra ;,title=strint(idx)
	oplot,lc,(obs_sed*fconvert)*peak,psym=8,color=109.
	oplot,templamz,temp_sed*peak,color=49.
	oploterror,lc,fnu[*,idx]*fconvert*peak,efnu[*,idx]*fconvert*peak,psym=8
	
	; plot p(z)
	if (keyword_set(nopz) eq 0) then begin
		!x.range=0
		!y.range=0
		plot,zgrid,pz[*,idx],xtitle='z',ytitle='p(z)'	
		if keyword_set(zs) then begin
		    if n_elements(zs) gt 1 then zsi = zs[idx] else zsi = zs
		    oplot,zsi*[1,1],[0,1.e10],color=109.
        endif
        
	endif
 
    !x.style=xst
    !y.style=yst
    ;; Stop to have access to the read-in variables
 	if keyword_set(stopme) then stop
 	
 end