pro make_fits_image, evtfil, imgfil, chan, chantype, centad=centad, $
    epoch=epoch, centpix=centpix, pixsz=pixsz, bin=bin, det=det, $
    matchim=matchim, img=img, hdr=hdr, rbands=rbands, plist=plist, trim=trim, $
    hist=hist, fudge=fudge
;+
; NAME: 
;
;    MAKE_FITS_IMAGE
;
; PURPOSE:
;
;    Creates a fits image imgfil out of events in evtfil using
;    channels chan(0) - chan(1), default = all.  The output image is
;    binned so that (centpix(0), centpix(1)) lies at (centad(0), centad(1)) 
;    with pixel size pixsz or binned by bin.  Alternatively, the
;    output image can be made to be the same size and location as
;    matchim, useful for overlays.
;
; CATEGORY:
;
;    Imaging
;
; CALLING SEQUENCE:
;
;    MAKE_FITS_IMAGE, evtfil, imgfil, [chan, chantype, centad=centad, 
;                     epoch=epoch, centpix=centpix, pixsz=pixsz, bin=bin, 
;                     det=det, matchim=matchim, img=img, plist=plist,
;                     trim=trim, hist=hist, fudge=fudge]
; 
; INPUTS:
;
;    evtfil - FITS event file
;    bin - binning factor for output image (optional if astrometry
;          info is supplied -- see optional inputs)
;
; KEYWORD PARAMETERS:
;
;   matchim - FITS image file containing astrometry to be matched
;   centpix - center x,y of output image
;   centad - ra and dec of centpix in degrees
;   pixsz - pixel size of output image in arcsecs, can be substituted
;           for bin above
;   epoch - epoch of output image
;   chan - range of PHA or PI (default) channels to include in binned
;          image, default = minmax(events channels) (i.e., include all events)
;   chantype - = 0  for PI (default), = 1 for PHA
;   det - /det if you want the output image in detector coords
;   plist - use plist for photons and only use evtfil for header info
;   trim - /trim excludes dy values more than 3200 det. pixels (~ 50') from
;          center of image to avoid problems with expmap
;   hist - string array of history lines to be added to header
;   /fudge - add a 0.5 pixel offset before binning, should NOT be nec.
;            this is included to allow compat. for cases where FITS
;            convension was not followed correctly

;
; OUTPUTS:
;
;   Creates fits image imgfil
;
; OPTIONAL OUTPUTS:
;
;   img = 2-dimensional image array saved in imgfil
;
; COMMON BLOCKS:
;
;   None
;
; SIDE EFFECTS:
;
;   None currently known
;
; RESTRICTIONS:
;
;   Requires idlaul library (http://idlastro.gsfc.nasa.gov)
;
; PROCEDURE:
;
;   Bins x,y coordinates using either bin or pixsz, outputs apropriate
;   astromemtry in imgfil header
;
; EXAMPLE:
;
;   make_fits_image, 'ngc3998_pspc.fits', 'ngc3998_pspc.img', pixsz=30., /trim
;   Creates a PSPC image with 30" pixels
;
; MODIFICATION HISTORY:
;   Created Jan 1996, A. Ptak 
;-


if n_params(0) eq 0 then begin
    print,'make_fits_image, evtfil, imgfil, chan, chantype, centad=centad, epoch=epoch,
    print,'  centpix=centpix, pixsz=pixsz, bin=bin, det=det, matchim=matchim, img=img,'
    print,'  hdr=hdr, rbands=rbands, plist=plist, h=h, trim=trim'
;    print,'Creates a fits image imgfil out of events in evtfil'
;    print,'Only use channels in the range chan(0) -  chan(1), default = all
;    print,'chantype = 0 for pha, 1 for pi (default)'
;    print,'The output image is binned so that (centpix(0), centpix(1)) lies at 
;    print,'  (centad(0), centad(1)) with pixel size pixsz.'
;    print,'Output image can be matched to same size and location as matchim'
;    print,'  Warning... make_emap_hdr may need to be modified to make cdelt1 neg.'
;    print,'Or specify:'
;    print,'  centpix:  Center (x,y) of image (def. = center of pix. in evt. file)'
;    print,'  centad:   (ra,dec) of centpix (def. = crval1, crval2 in evt header)'
;    print,'  epoch:    epoch of (ra,dec), default = value in evt header or 2000.0'
;    print,'  pixsz:    size of pixels in "'
;    print,'  bin:      binning factor of pixels (pixsz take precedence)'
;    print,'            If neither pixsz nor bin is given, bin=1 is used'
;    print,'  /det:     use detector coords rather than sky coords'
;    print,'            [sizex,sizey] for a rectangular image'
;    print,'itype=itype: Use itype for image pixel type (i.e, first element of size(img)) 
;    print,'plist=plist: use plist for photon list and only use evtfil for header info'
;    print,'             rather than integer'
;    print,'             1 = byte'
;    print,'             2 = integer (default)'
;    print,'             3 = long integer (integer*4)'
;    print,'             4 = float (real*4)'
;    print,'             5 = double (real*8)' 
;    print,'Returns resultant image in img and header in hdr'
;    print,'If PSPC events, chan. range can be given by specifying rbands as given in'
;    print,'  Snowden, et al 1994, ApJ, 424, 714  (i.e., rbands=[1,2] for R1+R2)'
     print,"Type doc_library, 'make_fits_image' for more info"
    retall
endif

pl = n_elements(plist) gt 0
fm = '(f7.2)'
fm4 = '(f9.4)'

if n_elements(chantype) eq 0 then chantype = 1
channame = ['PHA','PI']
chantype = (chantype>0) < 1  ; range check 0 le chantype le 1 

typename = ['byte','integer','long','float','double']
if n_elements(itype) eq 0 then itype = 2
itype = (itype>1) < 5
print,'Output image will be '+typename(itype-1)

; Get primary header
h0 = headfits(evtfil)
obj = sxpar(h0,'object')
print, obj

; Move to EVENTS extension
ext = 1
h = ''
a = ''
while strmid(a,0,6) ne 'EVENTS' and strmid(a,0,6) ne 'STDEVT' do begin
   h = headfits(evtfil,ext=ext)
   a = sxpar(h,'EXTNAME')
   ext = ext+1
endwhile

ext = ext-1
opt = [-999., -999.]

; Find out what instrument this is and read approriate x,y coords
instr = strmid(strtrim(sxpar(h, 'instrume')),0,4)
if not pl then begin
  print,'Reading ' + evtfil
  tab = readfits(evtfil, ext=ext)
endif

case 1 of
   instr eq 'PSPC' or instr eq 'HRI': begin
     print,'ROSAT event list'
; Find out what keywords are used for astrometry
     if sxpar(h, 'cdelt2') ne 0 then begin
       print,'Assuming you are using US Rev0 data - good luck'
       if keyword_set(det) then begin
         if pl then begin
           x = plist.dx
           y = plist.dy
         endif else begin
           x = tbget(h, tab, 'dx')
           y = tbget(h, tab, 'dy')
         endelse
         cdelt = abs(sxpar(h, 'xs-inpxx'))
         if cdelt ne sxpar(h, 'xs-inpxy') then begin
           print,'Warning: abs(xs-inpxx) does not equal xs-inpxy in ' + evtfil
           print,'Using abs(xs-inpxx) for plate scale'
         endif
         crval = [sxpar(h, 'xs-rapt'), sxpar(h, 'xs-decpt')]
         crpix = [sxpar(h, 'xs-xpt'), sxpar(h, 'xs-ypt')]
         opt = [sxpar(h, 'xs-xdopt'), sxpar(h, 'xs-ydopt')]
         rot = sxpar(h, 'xs-rsrot')
       endif else begin ;  sky coords
         if pl then begin
           x = plist.x
           y = plist.y
         endif else begin
           x = tbget(h, tab, 'x')
           y = tbget(h, tab, 'y')
         endelse
         crval = [sxpar(h, 'crval1'), sxpar(h,'crval2')]
         cdelt = sxpar(h, 'cdelt2')
         if cdelt eq 0.0 then begin
            print,'Could not find astrometry info in '+evtfil
            stop
         endif
         if cdelt ne -sxpar(h,'cdelt1') then begin
           print,'cdelt1 does not equal -cdelt2 in '+evtfil
           print,'Using abs(cdelt2) for plate scale'
           cdelt = abs(cdelt)
         endif
         crpix = [sxpar(h, 'crpix1'), sxpar(h, 'crpix2')]
       endelse
     rot = sxpar(h, 'crota2')
     ctype1 = sxpar(h, 'ctype1')
     ctype2 = sxpar(h, 'ctype2')
     ra_nom = sxpar(h, 'XS-RAPT')
     dec_nom = sxpar(h, 'XS-DECPT')
     droll = sxpar(h, 'XS-DANG')
     texp = sxpar(h, 'XS-LIVTI')
   endif else begin
     if keyword_set(det) then begin
       if pl then begin
         x = plist.dx
         y = plist.dy
       endif else begin
         x = tbget(h, tab, 'dx')
         y = tbget(h, tab, 'dy')
       endelse
       crval = [sxpar(h, 'tcrvl6'), sxpar(h, 'tcrvl7')]
       cdelt = sxpar(h,'TCDLT7')
       if cdelt ne -sxpar(h,'tcdlt6') then begin
         print,'tcdlt7 does not equal -tcdlt6 in '+evtfil
         print,'Using abs(tcdlt7) for plate scale'
         cdelt = abs(cdelt)
       endif
       crpix = [sxpar(h,'TCRPX6'),sxpar(h,'TCRPX7')]
       rot = sxpar(h,'tcrot6')
     endif else begin ; sky coords
       if pl then begin
         x = plist.x
         y = plist.y
       endif else begin
         x = tbget(h, tab, 'x')
         y = tbget(h, tab, 'y')
       endelse
       crval = [sxpar(h, 'tcrvl1'), sxpar(h,'tcrvl2')]
       cdelt = sxpar(h, 'tcdlt2')
       if cdelt ne -sxpar(h,'tcdlt1') then begin
         print,'tcdlt2 does not equal -tcdlt1 in '+evtfil
         print,'Using abs(tcdlt2) for plate scale'
         cdelt = abs(cdelt)
       endif
       crpix = [sxpar(h,'TCRPX1'),sxpar(h,'TCRPX2')]
       rot = sxpar(h,'tcrot2')
     endelse
     ctype1 = sxpar(h, 'TCTYP1')
     ctype2 = sxpar(h, 'TCTYP2')
     ra_nom = sxpar(h, 'ra_nom')
     dec_nom = sxpar(h, 'dec_nom')
     droll = sxpar(h,'tcrot6')
     texp = sxpar(h, 'ontime')
   endelse
   end
   strmid(instr,0,3) eq 'SIS' or strmid(instr,0,3) eq 'GIS': begin
     print,'ASCA '+instr
     if strmid(instr,0,3) eq 'SIS' then begin
       detxnum = '8' & detynum = '9'
       skyxnum = '3' & skyynum = '4'
     endif else begin
       detxnum = '10' & detynum = '11'
       skyxnum = '2' & skyynum = '3'
     endelse
     if keyword_set(det) then begin
;       xnum = detxnum & ynum = detynum
       xnum = skyxnum & ynum = skyynum
       if pl then begin
         x = plist.detx & y = plist.dety
       endif else begin
         x = tbget(h, tab, 'detx') & y = tbget(h, tab, 'dety')
       endelse
     endif else begin
       xnum = skyxnum & ynum = skyynum
       if pl then begin
         x = plist.skyx & y = plist.skyy
       endif else begin
         x = tbget(h, tab, 'x') & y = tbget(h, tab, 'y')
       endelse
     endelse
     crval = [sxpar(h, 'tcrvl'+xnum), sxpar(h, 'tcrvl'+ynum)]
     crpix = [sxpar(h, 'tcrpx'+xnum), sxpar(h, 'tcrpx'+ynum)]
     cdelt = sxpar(h, 'tcdlt'+ynum)
     if sxpar(h, 'tcdlt'+xnum) ne -cdelt then begin
         print,'tcdlt'+ynum+' does not equal -tcdlt'+xnum+' in '+evtfil
         print,'This is expected (see ftool 
         print,'Using abs(tcdlt'+ynum+') for plate scale'
         cdelt = abs(cdelt)
     endif
     ctype1 = sxpar(h, 'TCTYP'+xnum)
     ctype2 = sxpar(h, 'TCTYP'+ynum)
     ra_nom = sxpar(h, 'ra_nom')
     dec_nom = sxpar(h, 'dec_nom')
     droll = sxpar(h,'tcrot'+detynum)
     texp = sxpar(h, 'ontime')
   end
   ELSE: begin
       print, 'Unknown instrument ' + instr
       return
   end
endcase

print, 'Astrometry in ' + evtfil + ':'
print, 'crval = (' + strn(crval(0)) + ', ' + strn(crval(1)) + ')'
print, 'crpix = (' + strn(crpix(0)) + ', ' + strn(crpix(1)) + ')'
print, 'cdelt = ' + strn(cdelt*3600, format='(f6.2)') + '"/pix'

if pl then begin
  if chantype then ph = plist.pi else ph = plist.pha
endif else ph = tbget(h, tab, channame(chantype))

if keyword_set(trim) then begin
  print,'Trimming photon list of high and low dy values'
  if pl then begin
    dx = plist.dx
    dy = plist.dy
   endif else begin
    dx = tbget(h, tab, 'dx')
    dy = tbget(h, tab, 'dy')
   endelse
   w = where(dy lt 7297  and dy gt 641)
   x = x(w)
   y = y(w)
   ph = ph(w)
endif

if n_elements(chan) eq 0 then chan = minmax(ph)
print, 'Using ' + channame(chantype) + ' values ' + strn(chan(0)) + ' - ' + $
  strn(chan(1))

; Get equinox
ep = sxpar(h,'epoch')
if ep eq 0 then ep = sxpar(h,'equinox')
if ep eq 0 then begin
   print,'Warning: No equinox/epoch found in'+evtfil+', assuming 2000.0'
   ep = 2000.0
endif

;print,'Initial pixel size = ' + strn(cdelt)*3600 + '"'

; If matching to another image, load info
if n_elements(matchim) ne 0 then begin
  print,'Getting astrometry from header of '+matchim
  mh = headfits(matchim)
  centpix = [sxpar(mh, 'crpix1'), sxpar(mh, 'crpix2')]
  epoch = sxpar(mh,'epoch')
  if epoch eq 0 then epoch = sxpar(mh,'equinox')  
  if epoch eq 0 then begin
    print,'Warning: No equinox/epoch found in' + matchim + ', assuming 2000.0'
    epoch = 2000.0
  endif
  centad = [sxpar(mh,'crval1'), sxpar(mh,'crval2')]
  print, '(' + strn(centpix(0), format=fm) + ', ' + strn(centpix(1), $
    format=fm) + ') is at (' + strn(centad(0)) + ', ' + strn(centad(1)) + $
    ', ' + strn(epoch, format=fm) + ') in ' + matchim
  if epoch ne ep then begin
    print,'Precessing crval0,1 in ' + matchim + ' from ' + strn(epoch) + $
      ' to ' + strn(ep)
    ra = centad(0) & dec = centad(1)
    precess, ra, dec, epoch, ep
    centad(0) = ra & centad(1) = dec
  endif
  pixsz = sxpar(mh, 'cdelt2')*3600
  if pixsz ne -sxpar(mh, 'cdelt1')*3600 then begin
    print,'cdelt1 does not equal -cdelt2 in '+matchim
    stop
  endif
  print,'Pixel size in ' + matchim + ' = ' + strn(pixsz) + '"/pix'
endif

if n_elements(bin) eq 0 then bin = 1
if n_elements(pixsz) eq 0 then pixsz = bin*cdelt*3600
print,'Output image pixel size = ' + strn(pixsz, format='(f5.2)') + '"'
bin = pixsz/cdelt/3600.0
print,'Output image binning = ' + strn(bin)
; The convension used in the M82 HRI products image (binned by 16)
; appears to be to add 0.5 after binning.
if n_elements(centpix) eq 0 then centpix = (crpix - 0.5)/bin+0.5
if n_elements(centad) eq 0 then centad = crval

; Compute offset in output image pixels to apply as a result of center
; ra,dec being different
; It is assumed that N is +y and E is -x
offx = (crval(0)-centad(0))*cos(crval(1)/!radeg)*3600.0/pixsz
offy = (centad(1) - crval(1))*3600.0/pixsz
print,'$(a,2(f8.2,x))','Offset applied (in output pixels) = ', offx, offy
print,'Binning original image pixels by ' + strn(bin)
if n_elements(matchim) ne 0 then im = make_array(sxpar(mh, 'naxis1'), $
  sxpar(mh,'naxis2'),type=itype) else $
  img = make_array(2*centpix(0) - 1, 2*centpix(1) - 1, type=itype)

print,'Output image will be centered at (' + strn(centpix(0), format=fm)$
  + ', ' + strn(centpix(1), format=fm) + ') = (' + strn(centad(0),format=fm4)$ 
  + ', ' + strn(centad(1),format=fm4) + ', ' + strn(fix(ep)) + ')'

s = size(img)
print, 'Extracting ' + strn(s(1)) + ' by ' + strn(s(2)) + ' ' + $
  typename(itype - 1) + ' image...'

w = where(ph ge chan(0) and ph le chan(1), ncnts)

back = string(replicate(8b,4))
step = 2000

; ix, iy = output image pixel of event i
ix = (x(w)-crpix(0))/bin + centpix(0)
iy = (y(w)-crpix(1))/bin + centpix(1)
zero = 1.0
if keyword_set(fudge) then zero = zero-0.5
xx = round(ix-offx-zero) ; Subtract 1 for IDL convensions, add 0.5 to conform
yy = round(iy-offy-zero) ; to M82 HRI product image (which may not be correct).

ww = where(xx ge 0 and xx lt s(1) and yy ge 0 and yy lt s(2),ncnts2)
for i=0l, ncnts2-1 do begin
  if i mod step eq 0 then print,"$(i3,a,$)",100.*i/ncnts2,'%'+back
  
  img(xx(ww(i)), yy(ww(i))) = img(xx(ww(i)), yy(ww(i))) + 1
endfor

if total(opt) ne -1998 then print, 'crpix = (' + strn(crpix(0)) + ', ' + $
  strn(crpix(1)) + ') is ' + strn(sqrt(total((opt-crpix)^2))*cdelt*60) + $
  "' off-axis."

print, "Writing fits file"

fxhmake, hdr, img
fxaddpar, hdr, 'content', 'IMAGE'
fxaddpar, hdr, 'FILIN001', evtfil
fxaddpar, hdr, 'origin', 'here'
fxaddpar, hdr, 'HDUCLASS', 'ogip'
fxaddpar, hdr, 'HDUCLAS1', 'IMAGE'
fxaddpar, hdr, 'DATE', systime(0)
fxaddpar, hdr, 'TELESCOP', sxpar(h, 'TELESCOP')
fxaddpar, hdr, 'INSTRUME', instr
fxaddpar, hdr, 'EQUINOX', ep
fxaddpar, hdr, 'OBJECT', obj
fxaddpar, hdr, 'RA_NOM', ra_nom
fxaddpar, hdr, 'dec_NOM', dec_nom
fxaddpar, hdr, 'RA_NOM', sxpar(h, 'XS-RAPT')
fxaddpar, hdr, 'DEC_NOM', sxpar(h, 'XS-DECPT')
fxaddpar, hdr, 'DROLLANG', sxpar(h, 'XS-DANG')
fxaddpar, hdr, 'EXPOSURE', sxpar(h, 'XS-LIVTI')
fxaddpar, hdr, 'DATE-OBS', sxpar(h, 'DATE-OBS')
fxaddpar, hdr, 'TIME-OBS', sxpar(h, 'TIME-OBS')
fxaddpar, hdr, 'DATE-END', sxpar(h, 'DATE-END')
fxaddpar, hdr, 'TIME-END', sxpar(h, 'TIME-END')
fxaddpar, hdr, 'IMGBIN', bin
fxaddpar, hdr, 'CRPIX1', centpix(0)
fxaddpar, hdr, 'DRPIX1', crpix(0)
fxaddpar, hdr, 'CRPIX2', centpix(1)
fxaddpar, hdr, 'DRPIX2', crpix(1)
fxaddpar, hdr, 'CRVAL1', centad(0)
fxaddpar, hdr, 'CRVAL2', centad(1)
fxaddpar, hdr, 'CDELT1', -pixsz/3600.
fxaddpar, hdr, 'CDELT2', pixsz/3600.
fxaddpar, hdr, 'DRDELT1', -cdelt
fxaddpar, hdr, 'DUNIT1', 'deg'
fxaddpar, hdr, 'DRDELT2', cdelt
fxaddpar, hdr, 'DUNIT2', 'deg'
fxaddpar, hdr, 'CROTA2', sxpar(h, 'CROTA2')
fxaddpar, hdr, 'CTYPE1', ctype1
fxaddpar, hdr, 'CTYPE2', ctype2
fxaddpar, hdr, 'TOTCTS', ncnts
fxaddpar, hdr, 'chanmin', chan(0)
fxaddpar, hdr, 'chanmax', chan(1)
fxaddpar, hdr, 'chantype', channame(chantype)
fxaddpar, hdr, 'CREATOR', 'IDL make_fits_image'
if n_elements(hist) gt 0 then for i=0,n_elements(hist)-1 do $
  fxaddpar, hdr, 'HISTORY', hist(i)

fxwrite, imgfil, hdr, img

return
end