PRO hxrs_drm2fits, edges_in, drm, edges_out=edges_out, file=file, origin=origin

IF NOT(keyword_set(filename)) THEN filename = 'hxrs_resp.rmf'

;Create a file and header
fxhmake, head_1, /extend, /initialize ;create a header

;Add keywords to the primary header
if not(keyword_set(origin)) then origin='CU-Boulder'
fxaddpar, head_1, 'ORIGIN', origin
sttim = anytim(!stime, /ecs)
fxaddpar, head_1, 'DATE', strmid(sttim, 0, 10)

;write the file
fxwrite, file, head_1

;Add the extension, number of rows = number of energies

; Get spex variables of interest
if not(keyword_set(edges_out)) then begin
    edges_hxrs,edges_out, width=delta_light
endif else begin
    delta_light=edges_out[1,*]-edges_out[0,*]
endelse

d2f_drm=drm
s = size(drm)
nchans = s[1]
dchans=fix(nchans)
nrows = s[2]

; XSpec wants a drm with counts/cm^2 per photons/cm^2.
; Spex uses a drm with counts/cm^2/keV per photons/cm^2.
for i=0,nchans-1 do d2f_drm[i,*]=drm[i,*]*delta_light[i]

fxbhmake, head_rmf, nrows, 'HXRS SPECTRAL RESPONSE RMF EXTENSION'

;Add keywords to the header
fxaddpar, head_rmf, 'EXTNAME', 'SPECRESP MATRIX'
fxaddpar, head_rmf, 'TELESCOP', 'MTI'
fxaddpar, head_rmf, 'INSTRUME', 'HXRS'
fxaddpar, head_rmf, 'FILTER', 'None'
fxaddpar, head_rmf, 'CHANTYPE', 'PHA'
fxaddpar, head_rmf, 'DETCHANS', dchans
fxaddpar, head_rmf, 'HDUCLASS', 'OGIP'
fxaddpar, head_rmf, 'HDUCLAS1', 'RESPONSE'
fxaddpar, head_rmf, 'HDUCLAS2', 'RSP_MATRIX'
fxaddpar, head_rmf, 'HDUCLAS3', 'DETECTOR'
fxaddpar, head_rmf, 'HDUVERS', '1.3.0'
fxaddpar, head_rmf, 'TLMIN4', 0

;Ok, start adding columns
fxbaddcol, 1, head_rmf, edges_in[0,0], 'ENERG_LO', tunit = 'keV'
fxbaddcol, 2, head_rmf, edges_in[1,0], 'ENERG_HI', tunit = 'keV'
fxbaddcol, 3, head_rmf, 1, 'N_GRP'
fxbaddcol, 4, head_rmf, 0, 'F_CHAN'
fxbaddcol, 5, head_rmf, dchans, 'N_CHAN'
fxbaddcol, 6, head_rmf, d2f_drm[*,0], 'MATRIX'

;Create the extension
fxbcreate, unit0, file, head_rmf

;Ok, now fill the file up
FOR j = 0, nrows-1 DO BEGIN
    rowno = j+1
    fxbwrite, unit0, edges_in[0,j], 1, rowno
    fxbwrite, unit0, edges_in[1,j], 2, rowno
    fxbwrite, unit0, 1, 3, rowno
    fxbwrite, unit0, 0, 4, rowno
    fxbwrite, unit0, dchans, 5, rowno
    fxbwrite, unit0, d2f_drm[*, j], 6, rowno
ENDFOR

;You're done, finish the extension
fxbfinish, unit0

;Write the EBOUNDS extension
;Add the extension, number of rows = number of energies
dchans=fix(nchans)
fxbhmake, head_e, dchans, ' SPECTRAL RESPONSE EBOUNDS EXTENSION'

;Add keywords to the header
fxaddpar, head_e, 'EXTNAME' , 'EBOUNDS'
fxaddpar, head_e, 'TELESCOP', 'MTI'
fxaddpar, head_e, 'INSTRUME', 'HXRS'
fxaddpar, head_e, 'FILTER'  , 'None'
fxaddpar, head_e, 'CHANTYPE', 'PHA'
fxaddpar, head_e, 'DETCHANS', dchans
fxaddpar, head_e, 'HDUCLASS', 'OGIP'
fxaddpar, head_e, 'HDUCLAS1', 'RESPONSE'
fxaddpar, head_e, 'HDUCLAS2', 'EBOUNDS'
fxaddpar, head_e, 'HDUVERS' , '1.2.0'

;add columns
fxbaddcol, 1, head_e, 1, 'CHANNEL'
fxbaddcol, 2, head_e, edges_out[0,0], 'E_MIN', tunit = 'keV'
fxbaddcol, 3, head_e, edges_out[1,0], 'E_MAX', tunit = 'keV'

;Create the extension
fxbcreate, unit0, file, head_e

;add the data, nchans rows
FOR j = 0, nchans-1 DO BEGIN
    rowno = j+1
    fxbwrite, unit0, j, 1, rowno
    fxbwrite, unit0, edges_out[0,j], 2, rowno
    fxbwrite, unit0, edges_out[1,j], 3, rowno
ENDFOR

fxbfinish, unit0
RETURN


EN