;+
 ;  PROJECT:
 ;    SDAC
 ;  NAME:
 ;    MEWE_KEV
 ;  PURPOSE:
 ;    This function returns a Mewe thermal spectrum (line + continuum) for EM=1.e44 cm^-3
 ;    Differs principally from MEWE_SPEC. in allowing units of "keV" and flux in
 ;    Earth detector units, i.e. (ph/cm2/sec/keV).
 ;  CALLING SEQUENCE:
 ;    Flux = mewe_kev(Te6,wave)       ; ph s-1
 ;    Flux = mewe_kev(Te6,wave,/photon)   ; ph s-1
 ;    Flux = mewe_kev(Te6,wave,/photon,/kev,/earth)
 ;               ;wave in kev, ph cm-2 kev-1 s-1
 ;    Flux = mewe_kev(Te6,wave,/erg)     ; erg s-1
 ;    Flux = mewe_kev(Te6,wave,/cosmic)   ; Cosmic abundances
 ;
 ;  INPUTS:
 ;    Te6    = Electron Temperature in MK (may be a vector)
 ;    Wave   = Vector of wavelengths (Ang) in ascending order
 ;       or, if keyword EDGES set,
 ;       a 2xN array of the 2 wavelength edges of N bins
 ;       in any order.  Must be a cont
 ;  CALLS:
 ;    XLINFLX, CONFLX, ANYTIM, GET_RB0P, LOC_FILE, PATH_DIR, RD_ASCII
 ;
 ;  OUTPUTS:
 ;    Flux   = Fluxes in ph s-1 or erg s-1
 ;      Fluxes = fltarr(N_elements(Te6),N_elements(wave))
 ;
 ;  OPTIONAL INPUT KEYWORDS:
 ;    PHOTON = If set, calculation is made in ph s-1 (default)
 ;    ERG    = If set, calculation is made in erg s-1
 ;    COSMIC = If set, use cosmic abundances (solar is default)
 ;    EDGES      = If set, interpret Wave as a 2XN set of channel boundaries
 ;    KEV        = Units of wave are in keV, output units will be ph keV-1
 ;     If KEV is set, assumes EDGES is set so EDGES must be 2XN
 ;    EARTH_FLUX = calculate flux in units cm-2 s-1 (bin unit, ang or keV)-1
 ;    DATE       = optional date for calculation of earth flux, def='2-apr-92'
 ;    FILE_IN    = explicit file name for xlinflx.
 ;    NOLINE     = Line fluxes are not included.
 ;    NOCONT     = Continuous flux not included.
 ;    REL_ABUN   = A 2XN array, where the first index gives the atomic number
 ;     of the element and the second gives its relative abundance
 ;     to its nominal value given by ABUN.
 ;  OPTIONAL OUTPUT KEYWORDS:
 ;    ABUN       = Abundances used for calculation
 ;    ELEM       = Elements corresponding to the abundances
 ;    CONT_FLUX  = Continuum flux alone (without line contribution).
 ;    DWAVE      = Widths of wavelength bins (Ang)
 ;    WAVE_RANGE = Interval over which lines are extracted (Ang)
 ;  RESTRICTIONS:
 ;    N.B. If both edges aren't specified in WAVE, then the bins of WAVE must
 ;   be equally spaced.
 ;  METHOD:
 ;    Calls xlinflx and conflx.  Xlinflx reads the line data from either
 ;    $DIR_SXT_SENSITIVE/mewe_cosmic.genx or $DIR_SXT_SENSITIVE/mewe_solar.genx
 ;  COMMON BLOCKS:
 ;   MEWE_KEV holds the value of the linflux normalization factor.
 ;  MODIFICATION HISTORY:
 ;    29-oct-92, Written, J. R. Lemen, LPARL
 ;    25-feb-93, JRL, Added /cosmic option
 ;    31-mar-93, RAS, Specification of wavelength bin edges.
 ;    17-mar-95, ras, uses get_rb0p for apparent solar radius to compute distance
 ;     for flux at earth
 ;    Version 5, 22-may-96, ras, added noline keyword
 ;    Version 6, 7-aug-96, ras, checked version of linflx
 ;    Version 7, richard.schwartz@gsfc.nasa.gov, 17-aug-1997, used loc_file with path_dir
 ;   instead of chkarg
 ;    Version 8, richard.schwartz@gsfc.nasa.gov, 11-nov-1998, used xlinflx instead of linflx.
 ;    19-jul-02, richard.schwartz@gsfc.nasa.gov, added to documentation about
 ;     KEV keyword.
 ;   29-aug-03, richard.schwartz@gsfc.nasa.gov, protect against no line energies
 ;   17-sep-03, richard.schwartz@gsfc.nasa.gov, changed dist to thisdist
 ;   to protect against dist.pro function
 ;  14-apr-2004, ras, don't call xlinflx when energies are out of range.
 ;	10-mar-2006, ras, fixed bug in /NOCONT option
 ;	22-mar-2006, ras, fixed bug so it works with a  single interval
 ;		removed string(7b) bell rings, changed format on date default
 ;	internally to avoid repeated calls to anytim
 ;  26-Jun-2006, Kim.  Made compatible with <6.1 versions of IDL (size structure doesn't
 ;      contain structure_name before 6.1)

 ;
 ;-

 Function mewe_kev,Te6,wave,photon=photon,erg=erg,elem=elem,abun=abun, $
 cont_flux=cont_flux,cosmic=cosmic,dwave=dwave,edges=edges, file_in=file_in, $
 wave_range=wave_range,kev=kev, earth_flux=earth, date=date, $
 nocont=nocont, noline=noline, rel_abun=rel_abun


 common mewe_kev, lcon

 ; With no calling arguments, display the calling sequence:

 if n_params() eq 0 then begin
     return,'   Flux= mewe_kev(Te6,wave[,/erg]) ; for EM=1.e44 cm-3'
     endif
 if keyword_set(kev) and keyword_set(erg) then begin
     print,' **** Error in mewe_kev ****';,string(7b);,string(7b)
     print,'      You cannot specify both /keV and /erg'
     return,-1
     endif else if keyword_set(kev) then begin
     photon = 1
     erg = 0
     endif


 if keyword_set(photon) and keyword_set(erg) then begin
     print,' **** Error in mewe_kev ****';,string(7b);,string(7b)
     print,'      You cannot specify both /photon and /erg'
     return,-1
     endif else if keyword_set(erg) then Units=1 else Units=0

 if n_elements(wave) eq 0 then begin
     print,' **** Error in mewe_kev ****';,string(7b);,string(7b)
     print,'      Wave or Te6 is undefined'
     help,Te6,Wave
     return,-1
     endif

 if n_elements(cosmic) eq 0 then cosmic = 0  ; Solar is default
 szwave = size(wave)
 ;Wavelength Bins
 if not keyword_set(Edges)  then begin

     ;Edges aren't specified, only channel centers
     if (szwave(0) eq  2) and (szwave(1) eq 2) then begin
         print,' **** Error in mewe_kev ****';,string(7b);,string(7b)
         print,'      Wave expected to be defined fltarr(N) or fltarr(1,N).'
         help,Wave
         return,-1
         endif
     ; Compute the bin sizes:
     dwave = (wave(2:*)-wave) /2.
     dwave = [dwave(0),dwave,dwave(n_elements(dwave)-1)]    ; End bins
     dwave = abs( rebin(transpose(dwave),N_elements(Te6),N_elements(wave)))
     wave_range = [min(wave),max(wave)]
     wave_centers = wave

     endif else begin
     if (not (last_item(szwave) eq 2)) then if (szwave(0) ne  2) or (szwave(1) ne 2)   then begin
         print,' **** Error in mewe_kev ****';,string(7b);,string(7b)
         print,'      Wave not defined for edges, fltarr(2,N)'
         help,Wave
         return,-1
         endif
     dwave = reform( wave(1,*)-wave(0,*))
     dwave = abs(  $
     rebin(transpose(dwave),N_elements(Te6),N_elements(wave(0,*))))
     wave_range = [min(wave),max(wave)]
     wave_centers = reform( wave(1,*) + wave(0,*) )/2.
     endelse

 ;ras 13-apr-94
 funits =  1.      ;default units
 if keyword_set(earth) then begin
     date = fcheck(date,{anytim2ints,time:0L,day:4843}) ; '4-apr-92'
     date_test = 1
     if size(date,/tname) ne 'STRUCT' then date_test=0
     if date_test eq 1 then if tag_names(date,/structure_name) ne 'ANYTIM2INTS' then date_test=0
     if date_test eq 0 then date = anytim(date)
     radius = 6.9598e10
     arcsec = 180./!pi*3600.
     thisdist = radius/sin((get_rb0p(date))(0)/arcsec)
     funits = thisdist(0)^2*4*!pi
     endif

 ;Interpret wave quantities as keV
 const = 12.3985
 if keyword_set(kev) then begin
     wave_range = const / reverse( wave_range) ;convert this to angstroms
     funits = funits * dwave
     dsiz = (size(dwave))(1:2)
     ;wave_centers in keV, per angstrom
     dwave = n_elements(dwave) eq 1 ? const * dwave / wave_centers^2 : const * dwave / $
     (n_elements(wave_centers) eq 1 ? wave_centers[0]^2 : $
     temporary(rebin(reform(wave_centers^2,1,dsiz(1)),dsiz(0),dsiz(1))))
     wave_centers = const / wave_centers       ;convert to angstroms
     endif

 ; Compute the continuum flux (Convert from EM=1.e50 to 1.e44 cm-3):
 cont_flux = (conflx(Te6,Wave_centers,Units+2)/1.e6) * dwave / funits
 nocont = keyword_set(nocont) ? 1.0 : 0.0
 flux = cont_flux * funits * (1-nocont)  ; Return continuum flux

;Don't call line function when energy range is outside the range of line centers
 if not keyword_set(noline) and max(wave_range) gt 1.29  and min(te6) lt 1e3 then begin
     ; Compute the line flux:
     ; Older versions are in units of 1e50 cm-3, newer in units of 1.e44

     if not keyword_set( lcon ) then begin
         proc = loc_file('xlinflx.pro',path=path_dir(), count=pcount)
         proc = rd_ascii(proc)
         if (where(strpos(proc(0:19),'Compute the Mewe line spectrum for EM=1.e50 cm^-3') ne -1))(0) ne -1 $
         then lcon=1.e6 else lcon = 1.0
         endif

     xlinflx,Te6,wave_line,Lflux,ion_line, erg=units, $
     wave_range=wave_range, elem=elem,abun=abun, $
     cosmic=cosmic, file_in=file_in, rel_abun=rel_abun
     Lflux = Lflux / lcon        ;convert from EM=1.e50 to 1.e44 cm-3, MAYBE!


     ; Rebin into the user supplied wavelength bins:
     if n_elements(Te6) eq 1 then Flux = transpose(Flux)
     if wave_line[0] gt 0 then begin
         if keyword_set(kev) then begin
             bin4wave = find_ix( [(wave(0,*))(*),wave(1,n_elements(wave)/2-1)], $
             const/wave_line)
             for i=0,n_elements(wave_line)-1 do $
             flux(0, bin4wave(i)) = flux(*, bin4wave(i)) + lflux(*,i)
             endif else $
         for i=0,n_elements(wave_line)-1 do begin
             ;changed to vector addition across Te6 from a loop, ras 93/3/31
             a=min(abs(wave_line(i)-Wave_centers))
             Flux(0,!c) = Flux(*,!c) + Lflux(*,i)
             endfor
         endif
     endif

 return,reform(flux/funits)       ; Line + continuum
 end