;+
;
; NAME:
;	hxrs_response
;
; PURPOSE:
;       hxrs response function, count rate per channel as a function of
;		incident photon spectrum.
;
; CATEGORY:
;       HXRS, Spectral analysis, fitting
;
; CALLING SEQUENCE:
;
; 	hxrs_response, edges_out=edges_out, edges_in=edges_in, area=area, drm=drm, $
;	flux=flux, em_flux=em_flux, gain=gain, newflux=newflux, con_factors=con_factors, $
;	un_con=un_con, error=error
;
; CALLED BY:
; 	drm_4_spex
;
; CALLS TO:
;	none
;
; Inputs:
;   FLUX
;       - photons/cm2/s/keV defined at EM_FLUX,
;		  optional- only needed to obtain conversion factors
;	EM_FLUX
;		- energy in keV for FLUX
;		  optional- default is the energies of the loss matrix
;
;	GAIN
;		- HXRS gain state
;
; Inputs/Outputs:
;	EDGES_OUT(2,8)
;		- Energy edges of HXRS 8 channels in keV
;   EDGES_IN(2,500)
;		- Energy edges for HXRS input in keV
;   AREA
;       - HXRS area, 4.91 cm^2
;	DFILE
;		- Name of response matrix fits file to read (optional)
;	NEWFLUX(I)
;		- Photon flux defined midpoints of EDGES_IN,
;		- Interpolated to FLUX and EM_FLUX on midpoints of EDGES_IN
;		- if not input
;
; Outputs:
;	DRM(I,J)
;		- Probability of an input flux in units of photon/cm^2/s in
;		- bin J resulting in a rate of counts/cm2/keV/s in bin I
;	CON_FACTORS(I)
;		- Conversion factors.  Ratio of the (Integrated flux/bin width)
;		to NEWFLUX(I)
;	UN_CON(I)
;		- Uncertainties on the conversion factors.  An expression
;		of the uncertainty in the detector response
;	ERROR	- set if a parameter is entered out of range
;
; COMMON BLOCKS:
;	None.
;
; SIDE EFFECTS:
;	none
;
; RESTRICTIONS:
;	none
;
; PROCEDURE:
;
;  CALCULATE THE HXRS RESPONSE MATRIX AND/OR COUNT RATE.
;
;  Intended use:
;
;  This procedure calculates the pulse-height response.  The matrix is computed
;  from an energy-loss matrix obtained by running the program RESP_CALC, for a
;  set of 100 narrow line input functions.  The energy loss matrix is broadened
;  for the resolution of the detector to yield the pulse-height response matrix.
;  The pulse-height response matrix may be  thought as the detector pulse-height
;  spectrum in each of the 8 energy channels for an incident spectrum of a mono-
;  energetic beam at one of the 100 incident energies.
;
;  This procedure may be used in several ways:
;
;  To obtain a pulse-height response matrix:
;
;  HXRBS_RESPONSE, EDGES_OUT=EDGES_OUT, EDGES_IN=EDGES_IN, DRM=DRM
;
;  This gives the pulse-height response matrix, DRM, the edges of its energy
;  input photon bins, EDGES_IN, the edges of the 8 HXRS channels, EDGES_OUT.  Using OMATRIX and
;  an incident flux integrated over the E_MATRIX bins, one can compute the count
;  rate in the 8 HXRS channels, i.e.
;
;  COUNTS (counts/keV/cm2/sec) = DRM#(FLUX (Ph/cm2/s/keV) * (EDGES_IN(1:*)-EDGES_IN))
;
;  This calculation of the count rate is extremely quick after DRM is
;  obtained.
;
;  CONVERSION FACTOR calculation
;
;  CONVERSION FACTORS allow a simple transformation from the photon flux at a
;  given energy to the count rate in a given energy channel for a specified
;  model spectrum.  In general, they are defined as the ratio
;  of the count rate spectrum in counts/s/keV/cm^2 to the photon spectrum in
;  photons/s/keV/cm^2 provided an incident spectrum is given.  Specifically,
;  the CONVERSION FACTORS used in the DECONVOLVE section of DCPFIT are the ratio
;  of two quantities for a given channel I:
;
;  		Count rate(counts/sec) in channel I
;  		-----------------------------------
;  		Model Photon Flux at EMID(I) * Area * Width(I)
;
;  where EMID(I) is the mid-point energy in channel I, Area is the HXRS area
;  of 4.91 cm^2, and Width(I) is the channel width in keV.
;
;  This procedure will also return the CONVERSION FACTORS as follows:
;
;  On the first pass:
;
;  HXRS_RESPONSE, EDGES_OUT=EDGES_OUT, EDGES_IN=EDGES_IN, GAIN=GAIN
;
;  Then, compute the photon flux, FLUX, on the midpoints of EDGES_IN and
;  the photon flux, NEWFLUX, on the midpoints of EDGES_OUT outside of this procedure.
;  The units of the photon flux must be photons/s/keV/cm^2.
;
;  you can get the conversion factors by calling:
;
;  HXRS_RESPONSE, FLUX=FLUX, NEWFLUX=NEWFLUX, $
;  		 CON_FACTORS=CON_FACTORS, UN_CON=UN_CON
;
;  where the 8 CONVERSION FACTORS and their fractional uncertainties are given by
;  CON_FACTORS and UN_CON, respectively.
;
;  If the energy edges are already known, then these arguments FLUX, NEWFLUX,
;  CON_FACTORS, and UN_CON can be entered on the first pass throught the program.
;
;  As with the calculation of the response matrix , DRM, detailed above,
;  this calculation is quick after the first pass.
;
;
;  NON-STANDARD LIBRARY ROUTINES NEEDED:
;	EDGE_PRODUCTS
;	X_EOUT_DRM
;
; MODIFICATION HISTORY:
; Paul Bilodeau, NASA GSFC / RITSS, 31-May-2000 modified hxrbs_response for the HXRS instrument
;-

pro hxrs_response, edges_out=edges_out, edges_in=edges_in, area=area, drm=drm, dfile=dfile, $
	flux=flux, em_flux=em_flux, gain=gain, newflux=newflux, con_factors=con_factors, $
	un_con=un_con, error=error
;stop
error = 1 ;assume an error

if n_elements(area) eq 0 then area=4.91 ;cm^2

if n_elements(edges_out) eq 0 then edges_hxrs, edges_out

;NOMINAL FRACTIONAL UNCERTAINTY IN CONVERSION FACTORS
un_con = make_array(8,/float,value=0.05)

; Call resp_calc to get all quantities of interest
nflux=500
elo=10.
ehi=1000.
detector='NAI'
d=0.45
z=[13]
gmcm=.263
func='fwhm'
func_par=[1,00,.41,0]

; Get the drm defined on edges_in energies for both input and output sides
if n_elements(drm) eq 0 then begin
	resp_calc,detector,area,func,func_par,d,z,gmcm,nflux,elo,ehi,eloss_mat,pls_ht_mat,edges_in,tmp_drm

	; interpolate the drm onto the final bins
	drm = x_eout_drm( in_drm=tmp_drm, in_e_out=edges_in, out_e_out=edges_out, nsub=20)
endif else resp_calc,detector,area,func,func_par,d,z,gmcm,nflux,elo,ehi,eloss_mat,pls_ht_mat,edges_in

; midpoint energies of energy-loss matrix
edge_products, eloss_mat, mean=em_matrix, width=ww_matrix

;SET THE DEFAULT INPUT FLUX ENERGIES
;the input matrix energies are the defaults for the input flux
;if em_flux isn't defined it is assumed to be on inmatrix energy mid-points
checkvar,em_flux, em_matrix

;!!!!!!!!!!!!!!!!!!!!!!
;IF A PHOTON FLUX HAS BEEN INPUT, APPLY THE RESPONSE MATRIX AND COMPUTE
;CONVERSION FACTORS. IF NOT, THEN JUST RETURN THE RESPONSE MATRIX.

if n_elements(flux) ne 0 then begin

	;IF THE EM_FLUX ARE DIFFERENT FROM THE EM_MATRIX, INTERPOLATE FLUX TO EM_MATRIX
	wdiff = where(em_matrix - em_flux ne 0.0, ndiff)
	if total(abs(em_matrix-em_flux)) ne 0. then flux2 = interpol( flux, em_flux, em_matrix) else $
		flux2 = flux

	;FROM PHOTON/CM2/SEC/KEV TO PHOTON/CM2/SEC
	phot_per_bin = flux2 * ww_matrix

	;CALCULATE THE COUNT RATE IN THE 10 HXRS CHANNELS
	cnts_per_kevcm2 = drm#phot_per_bin

	;EXPRESS IN UNITS OF COUNTS/CM2/SEC/KEV TO GET CONVERSION FACTORS

	newflux = interpol(flux,em_flux,emedges)
	con_factors = cnts_per_kevcm2 / newflux
endif

error = 0	;success!
ERROR_RETURN:
return
end