pro imgrism,phafile,psname = psname,sigrange= sigr,grismspec=grismspec, $
   lower=lower,upper=upper,yrange=yrange,mkwave=mkwave,centroid=centroid, $
   srcwid=srcwid,png=png
;+
; NAME:
;     IMGRISM
; PURPOSE:
;     Image display of the .pha spectrum files output by UVOTIMGRISM
; EXPLANATION:
;      The 2-d image from the 3rd extension of the .pha file is displayed,
;      with the gross and background slits overlayed, and a co-aligned count 
;      rate plot appended.
; CALLING SEQUENCE:
;     IMGRISM, phafile, [ psname= ,sigrange=,/grismspec, /lower, /upper ]
; REQUIRED INPUTS:
;     phafile - scalar string giving name of the .pha file created by 
;         UVOTIMGRISM (or dogrism.pro).   Including the .pha extension is
;         optional
;     X - numeric scalar giving X centroid of the zero order 
;     Y - numeric scalar giving Y centroid of the zero order
; OPTIONAL KEYWORD INPUTS:
;     /CENTROID - Follow the trace of the spectrum by centroiding perpendicular
;            to the dispersion.    Implies the /grismspec option is also set.
;     /PNG - if set, then create a png output file.    The name will be the 
;            same as the input phafile but with a .png extension
;     PSNAME - indicate that the image is to written to a postscript file.
;           Can be either a scalar string giving the name of the postscript
;           file or /PSNAME which create a file with the same name as the
;           input phafile but with a .ps extension.
;     SIGRANGE - two element numeric vector specifying the the minimum and 
;           maximum display values in units of the image sigma.  Default
;           is [-1,7], i.e. minimum = avg-sigma : maximum = avg + 7*sigma
;     /GRISMSPEC - Indicates that the procedure grismspec.pro should be used
;           to extract the spectrum rather than the UVOTIMGRISM default.
;           The main difference is the use of a robust polynomial fit to the
;           background, and the option to use the /LOWER and /UPPER keywords.
;     /LOWER, /UPPER - Options to use only the lower (upper) background instead
;           of averaging.   Ignored if /GRISMSPEC is not set.
;    YRANGE - 2 element vector giving the Y range of the net counts plot.  
;            Useful when a spurious spike makes the default Y range too large.
; EXAMPLE:
;     IDL> imgrism,'sw00054300009ugv_3.pha',/psname 
;
;     This example creates an output postscript file sw00054300009ugv_3.ps
; NOTES:
;     Currently hardcoded to a first wavelength of 1700 A (U grism) or 2700 A
;     (V grism)
; PROCEDURES CALLED:
;     PLOTIMAGE (from http://astrog.physics.wisc.edu/~craigm/idl/idl.html)
; MODIFICATION HISTORY:
;      Written, W. Landsman   SSAI   Aug 2005
;
compile_opt idl2
On_error,2  

if N_params() LT 1 then begin
      print,'Syntax - imgrism, file, PSNAME= , Sigrange= ,/Grismspec'
      print,'             /Lower, /Upper, /Centroid, Srcwid='
      return
endif

fdecomp,phafile,disk,dir,fname,ext
if ext EQ '' then ext = 'pha'
file = disk + dir + fname + '.' + ext

fits_open,file[0],fcb
outfile=  keyword_set(psname) or keyword_set(png) ;Write to an output file

if outfile then begin
  sv_device = !D.Name
   fdecomp,file,disk,dir,fname
   
    if keyword_set(psname) then begin 
   if size(psname,/tname) NE 'STRING' then psname = fname + '.ps'
   set_plot,'ps'
   device,/color,/isolatin,file=psname
   sv_font = !p.font
   !p.font=0
   message,'Creating file ' + psname,/inf
   endif else begin 
      set_plot,'Z'
      message,'Creating file ' + fname + '.png',/inf
   endelse   
endif else erase  

 fits_read,fcb,tab,htab,exten=2     ;2nd extension contains spectrum
 fits_read,fcb,im,h,exten=3         ;3rd extension contains image
 fits_close,fcb

;Perform independent extraction?
 if keyword_set(grismspec) or keyword_set(centroid) then begin
     grismspec,file,tab,lower=lower,upper=upper,mkwave=mkwave,centroid=centroid,$
       srcwid=srcwid
    w = tab.w
    f = tab.net
    yy = tab.y
  
 endif else begin 
   w = tbget(htab,tab,'LAMBDA')
   f = tbget(htab,tab,'NET')
endelse    

   fdecomp,file,disk,dir,fname
   tabb = mrdfits(disk + dir + fname + '_back.pha',1,htabb)
   b = reverse(tabb.counts)
dim = size(im,/dimen)
wp = dim[0] - N_elements(w)
sky,im[wp:*,*],avg,var,/Sil
im = rebin(im,dim[0],dim[1]*2)
ratio = dim[1]*2.0/dim[0]
pos = [0.0,0.7,1.0,1.0]
if N_elements(sigr) NE 2 then sigr = [-1,7]
hh = sxpar(fcb.hmain,'HISTORY')
tmp = gettok(hh,' ')
keyword = gettok(hh,'=')

filter = strtrim(sxpar(h,'FILTER'),2)
if FILTER EQ 'UGRISM' then begin
   wtick = 1700 + indgen(40)*100
   wtickn = replicate(' ',40)
   wtickn[[3,13,23,33]] = ['2000','3000','4000','5000']
 endif else begin
    wtick = 2700 + indgen(36)*100
     wtickn = replicate(' ',36)
    wtickn[[3,13,23,33]] = ['3000','4000','5000','6000']
endelse

do_radec =0
keyword = strtrim(keyword,2) & hh = strtrim(hh,2)
g = where(keyword EQ 'sourcex')
sourcex = strcompress(hh[g[0]],/remove_all)
if sourcex EQ 0 then begin 
     do_radec=1
     g = where(keyword EQ 'ra')
     ra = strcompress(hh[g[0]],/remove_all)
     g = where(keyword EQ 'dec')
     dec = strcompress(hh[g[0]],/remove_all)
endif else begin      
    g = where(keyword EQ 'sourcey')
    sourcey = strcompress(hh[g[0]],/remove_all)
endelse    
g = where(keyword EQ 'infile')
finfile = strcompress(hh[g[0]],/remove_all)
fdecomp,finfile,disk,dir,infile,ext
infile = infile + '.' + ext
g = where(keyword EQ 'ang')
angle = strcompress(hh[g[0]],/remove_all)
if N_elements(srcwid) EQ 0 then begin 
  g = where(keyword EQ 'srcwid')
  srcwid = strcompress(hh[g[0]],/remove_all)
endif
g = where(keyword EQ 'bkgwid1')
bkgwid1 = strcompress(hh[g[0]],/remove_all)
g = where(keyword EQ 'bkgoff1')
bkgoff1 = strcompress(hh[g[0]],/remove_all)
object = sxpar(h,'OBJECT')
dateo = sxpar(h,'DATE-OBS')
expt = strtrim( string(sxpar(h,'LIVETIME'),f='(F6.1)'),2) 
detx = sxpar(htab,'ZERODETX')-1
dety = sxpar(htab,'ZERODETY')-1
detx = string(detx,'(F7.2)')
dety = string(dety,'(F7.2)')
if !D.name EQ 'X' then begin 
     win = !D.window > 0
     tit= object + ' ' +  filter + ' ' + expt + 's' + ' ' + infile
     window,win,title='IDL ' + strtrim(win,2) + ' ' + tit                     
endif     
plotimage,255b-bytscl(im,avg+sigr[0]*var,avg+sigr[1]*var), $
            position=pos,/keep,xsty=5,/noerase,imgyrange=[0,dim[1]-1], $
	    yrange=[0,dim[1]-1]
	    
axis,xaxis=1,xsty=1,xtit='Pixels'
pixcolor,12,'g'
cen = dim[1]/2
wid2 = fix(srcwid)/2.
if keyword_set(centroid) then begin 
   windex = wp + indgen(N_elements(w) )
   plots,windex, yy-wid2,/data,col=12,thick=4
   plots,windex, yy+wid2,/data,col=12,thick=4
   
endif else begin 
plots,[0,dim[0]],[cen-wid2,cen-wid2],/data,col=12,thick=4
plots,[0,dim[0]],[cen+wid2,cen+wid2],/data,col=12,thick=4
endelse
len = (pos[2] - pos[0])*N_elements(w)/float(dim[0])
bkgoff = fix(bkgoff1)
bkgwid = fix(bkgwid1)
if keyword_set(centroid) then begin 
  b1 = (yy - srcwid/2 -bkgoff-bkgwid) > 0
  b2 = b1 + bkgwid
   plots,windex, b1,/data,col=12,thick=4,/lines
   plots,windex, b2,/data,col=12,thick=4,/lines
endif else begin 
   
b1 = cen - srcwid/2 -bkgoff-bkgwid
b2 = b1 + bkgwid
plots,[0,dim[0]],[b1,b1],/data,col=12,/lines,thick=4
plots,[0,dim[0]],[b2,b2],/data,col=12,/lines,thick=4
endelse 


if keyword_set(centroid) then begin 
   b1 = yy + srcwid/2 + bkgoff
   b2 = (b1 + bkgwid) < dim[1]-1
   plots,windex, b1,/data,col=12,thick=4,/lines
   plots,windex, b2,/data,col=12,thick=4,/lines
endif else begin 

b1 = cen + srcwid/2 + bkgoff
b2 = b1 + bkgwid
plots,[0,dim[0]],[b1,b1],/data,col=12,/lines,thick=4
plots,[0,dim[0]],[b2,b2],/data,col=12,/lines,thick=4
endelse
xyouts,0,pos[1]-0.1,infile  +'!C' + object + '!C' + dateo + $
 '!C' + expt + 's',chars=0.9,/norm 
 
if do_radec then line1 = 'RA = ' + ra + '!CDec = ' + dec else $
                 line1 = 'SourceX = ' + sourcex + '!CSourceY = ' + sourcey
xyouts,0,pos[1]-0.25, line1  + $ 
  '!CDet X = ' + detx + '!CDet Y = ' + dety + $
  '!CAngle = ' + angle + '!CSrcwid = ' + strtrim(srcwid,2) + $
  '!CBkgoff = ' + bkgoff1 + '!CBkgwid = ' + bkgwid1 + '!C!C' + $
   strmid(systime(),4,20) , chars=0.9,/norm
pos[0] = pos[2] - len
xx = wp + indgen(N_elements(w))
bad =where(wtick GT max(w), Nbad)
if Nbad GT 0 then remove,bad,wtick,wtickn
xticks = N_elements(wtick)-1
quadterp,w,xx,wtick,xpos


if not keyword_set(yrange) then yrange = [0,max(f)]
plot,xx,f,xsty=5,yran=yrange,position=[pos[0],0.20,pos[2],pos[1]], $
  /noerase,ytit='Net Counts',chars=0.9
  
oplot,!X.crange,[0,0],/lines
plot,xx,b,/xsty,/ynoz,position=[pos[0],0.08,pos[2],0.20],/noerase,yticks=2, $
 ytit = 'Bkg',xtit= 'Wavelength (' + string(197b) + ')',xtickv=xpos, $
 xtickname=wtickn, xticks=xticks
 g = where(strlen(strtrim(wtickn,2)) GT 0 ,Ng)
 
 axis,xaxis=0,xtickv=xpos[g],xticks=ng-1,xticklen=-0.04, $
   xtickname=replicate(' ',ng)
   
if keyword_set(grismspec) then oplot,xx,tab.gross-f
if N_elements(psname) GT 0 then begin 
    device,/close
    !p.font = sv_font
    set_plot,sv_device
endif
if keyword_set(png) then begin 
    TVLCT, r, g, b, /GET
   im = tvrd() 
   Device, Z_Buffer=0
   Set_Plot, sv_Device
  
  Write_PNG, fname + '.png', im, r, g, b 
endif
return
end