;+
;
; PROJECT:  CHIANTI
;
;       CHIANTI is an atomic database package for the calculation of
;       continuum and emission line spectra from astrophysical plasmas. It is a 
;       collaborative project involving the Naval Research Laboratory
;       (Washington DC, USA), the Arcetri Observatory (Firenze, Italy), and the
;       Cambridge University (United Kingdom).
;
; NAME:
;	FB_RAD_LOSS
;
; PURPOSE:
;
;	Calculate the free-bound (radiative recombination) continuum
;       radiative energy losses losses.  The calculation follows the discussion
;       in Radiative Processes in Astrophysics by Rybicki and Lightman except that 
;       we have used observed energy levels in place of hydrogenic energy levels.
;
;
; CALLING SEQUENCE:
;
;      FB_RAD_LOSS,Temperature,LossRate
;
;
; INPUTS:
;
;
;
; OPTIONAL INPUTS:
;
;	None
;	
; KEYWORD PARAMETERS:
;
;	NO_SETUP:   If the procedure setup_elements has already been called then
;                  the keyword /no_setup should be set to avoid repeating this step
;
;       MIN_ABUND:  If set, calculates the continuum only from those elements which 
;                   have an abundance greater than min_abund.  Can speed up the 
;                   calculations.  For example:
;                   abundance (H)  = 1.
;                   abundance (He) = 0.085
;                   abundance (C)  = 3.3e-4
;                   abundance (Si) = 3.3e-5
;                   abundance (Fe) = 3.9e-5
;
;
; OUTPUTS:
;
;	Temperature:	temperature in degrees Kelvin, can be a 1 D array
;       LossRate:       radiative energy loss rate in erg s^-1 cm^3
;                       (radiative loss rate per emission measure (N_e N_H V)
;
;
;
; COMMON BLOCKS:
;
;	common elements,abund,abund_ref,ioneq,ioneq_t,ioneq_ref
;
;
;
; EXAMPLE:
;
;             > fb_rad_loss,t,rad
;             > fb_rad_loss,t,rad,min_abund=3.e-5
;             > fb_rad_loss,t,rad,/no_setup,min_abund=1.e-6
;
;
; MODIFICATION HISTORY:
; 	Written by:	Ken Dere
;	April 2000:     Version 3.0
;
;-
;
pro fb_rad_loss,t,rad_loss,no_setup=no_setup,min_abund=min_abund,elapsed=elapsed
;
;  to calculate the free-bound radiation loss as a function of temperature
;  units are ergs cm^-3 s-1 
;
common elements,abund,abund_ref,ioneq,ioneq_t,ioneq_ref
;
;
;
if n_params(0) lt 2 then begin
  print,'   IDL> ff_rad_loss,temperature,rad_loss  [,/no_setup, min_abund= ]'
  print,'     use /no_setup when called inside a self-consistent program'
  print,'     free-bound calculations are somewhat slow, setting higher min_abund values can help'
  print,'  Allen abundances for some elements'
  print,'     abund(H)  = 1.'
  print,'     abund(He) = 8.5e-2'
  print,'     abund(C)  = 3.3e-4'
  print,'     abund(Si) = 3.3e-5'
  print,'     abund(Fe) = 4.0e-5'
  return
endif
;
if keyword_set(elapsed) then time0=systime(1)
;
kb=1.38062d-16   ;  erg deg-1
h=6.6262d-27
c=2.997925d+10
ryd=2.17992d-11  ; erg
;
if not keyword_set(no_setup) then begin
   setup_elements
endif
   read_ip,!xuvtop+'/ip/chianti.ip',ionpot,ipref
;
zmax=max(where(abund gt 0.))+1
;
;  zmax=6  ;***************************
;
;  print,zmax,abund(zmax-1)
;
n_ioneq_t=n_elements(ioneq_t)
t=10.^ioneq_t
;
;
;
rad_loss=fltarr(n_ioneq_t)
;
;
;
for iz=1,zmax do begin       ;**************************************
   z2element,iz,element
   elstr=ch_strpad(element,10)
   this_abund=abund(iz-1)
   if keyword_set(min_abund) then begin
        test2 = this_abund gt min_abund
   endif else test2 = this_abund gt 0.
;
;  if test2 then print,iz,' element = ',elstr,'  this_abund = ',this_abund
;
  for ion=1,iz do begin   ; neutral to hydrogenic species
;
;  zion2spectroscopic,iz,ion,snote
;
  z=float(ion-1)
  dielectronic=0
  this_ioneq=ioneq(*,iz-1,ion-1+dielectronic+1)  ; +1 because recombining higher ion
;    print,' iz, ion = ',iz,ion,' max ioneq = ',max(this_ioneq)
;
;
  ip=ionpot(iz-1,ion-1)
;    print,'                   ip = ',ionpot(iz-1,ion-2)
;
test1=(max(this_ioneq) gt 0.)
test3=ip gt 0.
;
;print,' test1,2,3 = ',test1, test2, test3
;
if test1 and test2 and test3 then begin
;
;   print,' recombining to form ',snote
;   print,' ionization potential = ',ip,1.e+8/ip
;   print,' iz, ion = ',iz,ion
;
;   get data for the recombined ion
;
  zion2filename,iz,ion,filename
  ename=filename+'.elvlc'
;  print,' filename = ',filename
;
if ch_file_exist(ename) then begin
;  print,' filename = ',ename
   read_elvlc_direct,ename,l1,term,conf,ss,ll,spd,jj,mult,ecm,eryd,ecmth,erydth,eref
   mult=float(mult)
   necm=n_elements(ecm)
   conf=strtrim(conf,2)
   conf=strlowcase(conf)
   pqn=intarr(necm)
   for iecm=0,necm-1 do begin
      cf=conf(iecm)
      conf2n,cf,n
      pqn(iecm)=float(n)
   endfor
endif else begin
   ecm=[0.] & ecmth=[0.] & mult=[1.] & pqn=[1.]  ;  mostly these are for the neutrals
endelse
;
;  get data for the recombining ion, mostly need the statistical weight (mult)
;
;
  zion2filename,iz,ion+1,filenamer
  enamer=filenamer+'.elvlc'
;  print,' filename = ',filenamer
;
if ch_file_exist(enamer) then begin
;  print,' filename = ',enamer
  read_elvlc_direct,enamer,l1r,termr,confr,ssr,llr,spdr,jjr,multr, $
       ecmr,erydr,ecmthr,erydthr,erefr
  multr=float(multr)
endif else begin
   multr=[1.]  ;  mostly these are for the completely ionized species
endelse
;
;
zed=where(ecm eq 0.,nzed)
for i=0,nzed-1 do ecm(zed(i))=ecmth(zed(i))
necm=n_elements(ecm)
;
;    avoid summing over levels above the ionization potential
;
   sum=0.
   en=h*c*(ip-ecm)
   gooden=where(en gt 0.,ngooden)
;
;   for igen=0,ngooden-1 do begin
;      iecm=gooden(igen)
;                 print,pqn(iecm),mult(iecm),en(iecm),ecm(iecm)
;      sum=sum+mult(iecm)*en(iecm)^2/(pqn(iecm)*multr(0))
;   endfor

;
   iecm=gooden(indgen(ngooden))
   sum=total(mult(iecm)*en(iecm)^2/(pqn(iecm)*multr(0)))
;
goodt=where(this_ioneq gt 0.,ngt)
;
for igt=0,ngt-1 do begin
   it=goodt(igt)
;
;  print,'t = ',t(it), ' abund = ', this_abund,' this_ioneq = ',this_ioneq(it), ' ip =',ip
;
   alpha=0.4737*this_abund*this_ioneq(it)/sqrt(t(it))

;
          rad_loss(it)=rad_loss(it)+alpha*sum
;
endfor ;  temperature
;
;
endif  ;  various tests
;
endfor;  ion
;
endfor;  iz
;
;
;for it=0,n_ioneq_t-1 do begin
;   print,ioneq_t(it),t(it),rad_loss(it)
;endfor
;
;
if keyword_set(elapsed) then begin
    time1=systime(1)
    dt=time1-time0
    print,' elapsed seconds in synthetic = ',dt
endif
;
end