;+
;
; 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
;       Cambridge University (United Kingdom).
;
; NAME:
;	FREEFREE
;
; PURPOSE:
;
;	Calculate the free-free continuum from an hot, low density plasma
;
;       Uses the free-free gaunt factor calculations of Sutherland, 1998, 
;       MNRAS, 300, 321.
;
;       Note that Sutherland's Eq.(15) has units of erg/cm^3/s. Comparing 
;       with Rybicki & Lightman's Eq.5.14(a) (in their book 'Radiative 
;       Processes in Astrophysics'), suggests that Sutherland's units 
;       should be erg/cm^3/s/sr/Hz. We are assuming the latter to be 
;       correct in this routine.
;
;       When using the DEM_INT optional input, FREEFREE expects the 
;       differential emission measure to have been derived from a product 
;       of N_e*N_H rather than N_e*N_e. This can be important when dealing 
;       with a regime (typically T < 10^4.5) where H and He are not fully 
;       ionised.
;
; CALLING SEQUENCE:
;
;       FREEFREE,temperature, wavelength, intensity
;
;
; INPUTS:
;
;	T	   Temperature in degrees Kelvin, can be a 1 D array
;       WVL        Wavelength in Angstroms
;
;
; 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.  
;                  This can speed up the calculations. ;
;
;	DEM_INT    The intensity array is multiplied by a DEM number for 
;                  each temperature. DEM_INT needs to be of the same size 
;                  as TEMPERATURE. It is needed for the synthetic spectrum 
;                  routines.
;
;       PHOTONS    The output spectrum is given in photon units rather 
;                  than ergs.
;
;       SUMT       The default is to output the intensity array as an array 
;                  of size (nwvl x nT). Setting this keyword performs a sum 
;                  over the temperatures to yield a vector of same size as 
;                  the input wavelengths, thus producing the complete 
;                  free-free spectrum.
;	
; OUTPUTS:
;
;	RAD        Free-free continuum intensity in units
;                  10^-40 erg cm^3 s^-1 str^-1 angstrom^-1 per unit 
;                  emission measure int(n^2 dh) in cm^-5. If T is given as 
;                  a 1-D array, then RAD will be output as a 2-D array, 
;                  with one element for each temperature and wavelength 
;                  (but also see SUMT).
;
;
; PROGRAMMING NOTES
;
;     The gaunt factors from Sutherland (MNRAS 300, 321, 1998) are a 
;     function of uu and gg (see his Eq. 14). uu is a function of both 
;     wavelength and T, while gg is a function of T only.
;
;     The gaunt factor (gff) is tabulated for values of uu and gg at 
;     fixed intervals in log(uu) and log(gg). The log(uu) values go from 
;     -4 to 4 in 0.2 steps; the log(gg) values go from -4 to 4 in 0.1 steps.
;
;     A particular (input) temperature and wavelength give rise to values 
;     uu0 and gg0, the logs of which lie between -4 and 4. To derive the 
;     corresponding gff0 value, I use the IDL routine BILINEAR. 
;
;     BILINEAR requires, not uu and gg values as input, but indices. 
;     E.g., the indices corresponding to the tabulated values of uu are 
;     0 (=-4.0), 1 (=-3.8), 2 (=-3.6), etc. Thus, if log(uu0)=-3.76, then 
;     i_uu0=1.20 is the index of uu0.
;
;     In order to make significant time-savings, I give BILINEAR all of 
;     the wavelengths and temperatures in the same call for a particular 
;     ion. To do this, I make my i_uu and i_gg values 2-D arrays of size 
;     (nwvl x nT), and BILINEAR then produces a (nwvl x nT) array 
;     containing the gff values.
;
;     A problem occurred if nT=1, as BILINEAR will turn the input uu and 
;     gg vectors into 2-D arrays of size (nwvl x nwvl). If there are a 
;     large number of wavelengths, this uses a lot of memory. To solve 
;     this I make a 2 element temperature vector whose values are 
;     identical, and then change this back to a 1 element vector after 
;     BILINEAR has been called. See the parts of the code where I use 
;     'tst1'.
;
;
; COMMON BLOCKS:
;
;	common elements,abund,abund_ref,ioneq,ioneq_t,ioneq_ref
;
; CALLS
;
;       READ_IP, READ_GFFGU, SETUP_ELEMENTS
;
; EXAMPLES:
;
;       IDL> freefree,1.e+6,wvl,int
;       IDL> freefree,1.e+6,wvl,int,min_abund=3.e-5
;       IDL> freefree,1.e+6,wvl,int,/no_setup,min_abund=1.e-6
;
;       IDL> wvl=findgen(5001)/10. + 50.
;       IDL> temp=10.^(findgen(41)/10. +4.)
;       IDL> freefree,temp,wvl,int
;
;
; MODIFICATION HISTORY:
; 	Written by:	Ken Dere
;	March 1999:     Version 2.0
;       September 1999:  Version 3.0
;
;       Ver.3.1, 11-Aug-00, Peter Young
;           Improved call to bilinear, allowing routine to solve for 
;           all temperatures in one go. This makes the routine quicker, 
;           and also lowers the memory usage of the routine when dealing 
;           with many wavelengths.
;
;       Ver.3.2, 16-Aug-00, Peter Young
;           Corrected expression for 'gg', replacing ip with z^2.
;
;       Ver.3.3, 16-Oct-00, Peter Young
;           Now deals with dem_int correctly
;
;       Ver.3.4, 10-Oct-01, Ken Dere
;           Corrected for labelling errors in Sutherland's gffgu.dat file
;           No longer reads ionization potential file
;
;       Ver.3.5  17-Dec-2001 Giulio Del Zanna
; 
;           Corrected for an error, a constant left out during the previous 
;           correction.  
; 
;-


pro freefree,t,wvl,rad,no_setup=no_setup,min_abund=min_abund, $
       dem_int=dem_int, photons=photons, sumT=sumT

COMMON elements,abund,abund_ref,ioneq,ioneq_t,ioneq_ref
;
;
if n_params(0) lt 3 then begin
  print,'   IDL> freefree,temperature,wavelength,intensity  [,/no_setup, min_abund=3.e-5]'
  print,'     use /no_setup when freefree is called inside a self-consistent program'
  print,'     continuum calculations are fairly 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

t1 = systime(1)

kb=1.38062d-16   ; erg/K
h=6.6262d-27     ; erg s
c=2.997925d+10   ; cm/s
ryd=2.17992d-11  ; erg
factor=5.44436d-39  ; this is f_k of Eq.15 of Sutherland
rescale=1.d+40
;
if not keyword_set(no_setup) then setup_elements
;
;
zmax=max(where(abund gt 0.))+1
;
;
n_ioneq_t=n_elements(ioneq_t)
dlnt=alog(10.^(ioneq_t(1)-ioneq_t(0)))
;
;  read free-free gaunt factors
;
read_gffgu,g2,u,gff
;

t = double(t)
wvl = double(wvl)

;
; If t and wvl only have 1 element, then need to convert them to arrays
;
IF n_elements(wvl) EQ 1 THEN wvl = [wvl]
IF n_elements(t) EQ 1 THEN t = [t]

nwvl=n_elements(wvl)
nt=n_elements(t)

ident_t = make_array(nt,val=1.,/double)
ident_wvl = make_array(nwvl,val=1.,/double)

;
; the isothermal here refers to the keyword used in the routine ch_synthetic
; when it is set, no DEM distribution is read in, and the so the 'demfactor' 
; is not used (see later)
;
IF n_elements(dem_int) EQ 0 THEN BEGIN
   dem_int = make_array(nt,val=1.) 
   isothermal = 1
ENDIF ELSE BEGIN
   dem_int = double(dem_int)
   isothermal = 0
ENDELSE


;
IF keyword_set(sumt) THEN rad = dblarr(nwvl) ELSE rad = dblarr(nwvl,nt)
;
;  if the # of wavelengths is large, this procedure can be time consuming
;
nwvls=nwvl
wvls=wvl


for iz=1,zmax do begin
   z2element,iz,element
   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.
;
;
  for ion=2,iz+1 do begin
    z=float(ion-1)
    this_ioneq=ioneq(*,iz-1,ion-1)
;
    goodt=where(this_ioneq gt 0.,ngoodt)
    if ngoodt ge 2 then begin
       mingoodt=min(ioneq_t(goodt))
       maxgoodt=max(ioneq_t(goodt))
    endif else begin
       mingoodt=-1. & maxgoodt=-2.
    endelse
;
    ind = WHERE((ALOG10(t) GE mingoodt) AND (ALOG10(t) LE maxgoodt),nind)
    IF ind[0] EQ -1 THEN BEGIN
       ioneq1 = 0.
    ENDIF ELSE BEGIN
       ioneq1 = SPLINE(ioneq_t,this_ioneq,ALOG10(t[ind]))
    ENDELSE

;
    test1=(ioneq1[0] NE 0.)
;
    if test1 and test2 then begin

      tst1 = 0
      temp = t[ind]
      IF n_elements(temp) EQ 1 THEN BEGIN
        temp = [temp,temp]
        tst1 = 1
      ENDIF

     ;
     ; uu is the u=hv/kT of Sutherland
     ;
      uu = h*c*1.d8/kb*(1/wvl)#(1/temp)
      i_uu = (ALOG10(uu)+4.d)*10.d

     ;
     ; gg is the gamma^2=Z^2 Ry/kT of Sutherland
     ;

      gg = z^2*ryd/kb* ident_wvl #(1/temp)

;      gg = z^2*ryd/109737.3d/kb* ident_wvl #(1/temp)
;      gg = ip*ryd/109737.3d/kb* ident_wvl #(1/temp)

      i_gg = (ALOG10(gg)+4.d)*5.d


; big is an array of size (nwvl x nt)
;
      big = bilinear(gff,i_gg,i_uu)

      newfactor=this_abund*ioneq1*factor*c*1.d+8*rescale

      IF isothermal EQ 0 THEN demfactor = temp*0.1/alog10(exp(1)) $
                         ELSE demfactor = make_array(nt,val=1.,/double)

      IF tst1 EQ 1 THEN BEGIN     ; just one temperature
         newbig = big*exp(-uu)
         newbig = newbig[*,0]
         newbig = newbig*z^2*(1/wvl^2)
         t_factor = newfactor*1/SQRT(temp[0])*dem_int[ind]*demfactor
      ENDIF ELSE BEGIN
         newbig = big*exp(-uu)*z^2*( (1/wvl^2)#ident_t)
         t_factor = newfactor*1/SQRT(temp)*dem_int[ind]*demfactor
      ENDELSE

      IF keyword_set(sumt) THEN BEGIN
         rad = rad+t_factor##newbig
      ENDIF ELSE BEGIN
         rad[*,ind] = rad[*,ind]+(ident_wvl#t_factor)*newbig
      ENDELSE


    ENDIF         ; various tests

  ENDFOR          ; ion

ENDFOR            ; iz


IF keyword_set(photons) THEN rad = rad*wvl/1.d8/h/c

t2 = systime(1)
print,format='("The call to freefree took ",f8.1," seconds")',t2-t1

END