;+
;
; NAME:
;   F_MULTI_SPEC
;
;
; PURPOSE:
;   This function returns the differential photon spectrum from a model of
;   a double thermal bremsstrahlung plus broken power-law with a low-energy
;   cutoff.
;
;
; CATEGORY:
;   XRAY, SPECTRA
;
;
; CALLING SEQUENCE:  output = f_multi_spec( energy, parameters )
;
;
; CALLED BY: get_conversion
;
;
; CALLS:  f_bpow, f_vth, edge_products
;
;
; INPUTS:
;       energy
;   a
;       See Procedure for a description
; OPTIONAL INPUTS:
;   none
;
; OUTPUTS:
;       function returns photon spectrum in each bin in (/cm2/s/keV)
;
; OPTIONAL OUTPUTS:
;   none
;
; COMMON BLOCKS:
;   function_com
;
; SIDE EFFECTS:
;   none
;
; RESTRICTIONS:
;   This function is essentially identical to f_2therm_pow and both are
;   preserved for the moment, May 1996
;
; PROCEDURE:
;   The spectrum is built from the sum of two calls to f_vth() and one to f_bpow().
;
;   Uses f_vth for the optically thin bremsstrahlung radiation function,
;   the differential spectrum seen at Earth in units of photon/(cm2 s keV).
;
;INPUTS for f_vth
;   E- energy vector in keV, 2XN edges or N mean energies
;   A- free parameters, where the elements of A have these meanings:
;   a(0)= emission measure in units of 10^49 cm-3
;   a(1)= KT   plasma temperature in keV
;
;And the second one
;   a(2)= emission measure in units of 10^49 cm-3
;   a(3)= KT   plasma temperature in keV
;
;
;Uses broken power law function, F_BPOW
;   see F_BPOW for more on defaults.
;   a(4) - normalization of broken power-law at Epivot (usually 50 keV)
;   a(5) - negative power law index below break
;   a(6) - break energy
;   a(7) - negative power law index above break
;   a(8) - low energy cutoff for power-law, spectrum goes as E^(-1.5)
;   a(9) - spectral index below cutoff between -1.5 and -2.0
;   default parameters - a(8) and a(9) default to 10 keV and -1.5
;   can be overwritten by values in common function_com
;   none
;!!!!!!!!!!!!N.B. For developmental purposes this routine can be modified through
;   the use of two environment variables.
;   SPEX_ATTEN- This can be used to set fixed values for additional attenuation.
;   There should be two values separated by a space or comma which are
;   interpreted as the thickness in cm and the atomic number, Z.
;   Set SPEX_ATTEN to 0 to use apar(6:7) [see below] for the attenuation
;   variables.
;   SPEX_NEWMODEL- Currently, this only has meaning if it is set to F_CON_LIN.
;   In this case the two thermal components are the two components of
;   F_CON_LIN, a thermal continuum (conflx) and thermal line spectrum (linflx)
;   from MEWE_SPEC which can have independent emission measures apar(0) and apar(2)
;   and can have independent or identical temperatures apar(1) and apar(3).  For
;   linked temperatures apar(3) must be set less than or equal to 0.  In this
;   case apar(6) and apar(7) are used for the attenuation constants, thickness
;   and Z if apar(7) is between 1 and 94 and an integer. Also, the broken
;   power-law doesn't have a high energy break, but is fixed to the value of
;   apar(5)
;
; MODIFICATION HISTORY:
;   ras, 8-March-1995
;   Version 2, richard.schwartz@gsfc.nasa.gov, 14-aug-1997, added
;   private controls through environment variables
;   Version 3, richard.schwartz@gsfc.nasa.gov, 28-oct-1997,
;   fixed bug added by private controls, was potentially adding
;   attenuation where none was wanted.  Now, for attenuation,
;   fixed or variable SPEX_ATTEN must be defined to some value.
;   ras, 13-aug-03, call f_vth and f_bpow with energy to preserve 2xN input.
;   ras, 29-sep-03, made sure attenuate path is only used by design
;     previously apar_use(6:7) were being set to 100 and apar_use(5) because
;     this had been designed for low energy work.  Now it's general again.
;
;-


function f_multi_spec, energy, a

if (size(energy))(0) eq 2 then edge_products, energy, mean=em else em=energy

@function_com

;The next environment variable is an ad hoc method of inputting additional
;attenuation
atten_const = fltarr(2) ;  [d in cm., Z atomic number]
attenuate = getenv('SPEX_ATTEN')
newmodel  = getenv('SPEX_NEWMODEL')
add_attenuation = attenuate ne ''
if keyword_set(attenuate) then  reads, attenuate, atten_const

;FUNCTION_COM - Common and common source for information on fitting functions
;
;common function_com, f_model, Epivot, a_cutoff

checkvar,apar,fltarr(10)

npar = n_elements(a)
apar(0) = a(0: (npar-1)< (n_elements(apar)-1) )
if total(abs(apar(8:9))) eq 0.0 then apar(8)=a_cutoff

apar_use = apar

if atten_const(0) eq 0 and attenuate ne '' then begin
    atten_const = apar(6:7)
    if atten_const(1) ne fix(atten_const(1)) or atten_const(1) lt 1 or $
    atten_const(1) gt 94 or atten_const(0) lt 0.0 then atten_const(1)=0.0
;
; Use apar(6) and apar(7) for the attenuation constants - length in cm and material Z.
;

    apar_use(6:7) = [100., apar(5)] ; no break in power-law.  This is only used at low
              ; energy so there should be no compromise

endif

if strupcase(newmodel) eq 'F_CON_LIN' then begin
    if apar_use(1) eq 0 then apar_use(0)= [0.,1.]
    ans = f_bpow(energy,apar_use(4:9)) + f_con_lin(energy, apar(0:3))
stop
endif else begin


    if apar_use(3) eq 0 then apar_use(2)= [0.,1.]
    if apar_use(1) eq 0 then apar_use(0)= [0.,1.]
    ans = f_bpow( energy, apar_use(4:9) )
    if apar_use(0) ne 0 then ans = ans + f_vth(energy, apar_use(0:1))
    if apar_use(2) ne 0 then ans = ans + f_vth(energy, apar_use(2:3))

endelse

if atten_const(1) ne 0.0 and add_attenuation $
    then ans = exp(-((xsec(em, atten_const(1),'pe')*atten_const(0))<40.0)) * ans

return,ans
end