function dustem_create_ionfrac, key=key, val=val,help=help

;+
; NAME:
;   dustem_create_ionfrac
;  
; PURPOSE:
;   Calculates the neutral/positive/negative PAH fractions for a given
;   temperature, density and grain size distribution and writes it to
;   the appropriate .DAT file (according to !run_ionfrac). Deprecated
;  
; CATEGORY:
;    DustEMWrap, Distributed, MidLevel, Deprecated
;  
; CALLING SEQUENCE:
;   
; INPUTS:
;  key --  (1) gas temperature and (2) n_elec
;  val -- values corresponding to parameters passed by key
;  other physical parameters taken from the dust model as specified in
;  the grain structure, i.e. (*!dustem_params).grains
;  
; OPTIONAL INPUT PARAMETERS:
;
; OUTPUTS:
;
; OPTIONAL OUTPUT PARAMETERS:
;
; ACCEPTED KEY-WORDS:
;    help = print this help
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURES AND SUBROUTINES USED:
;   
; EXAMPLES
;  
; MODIFICATION HISTORY:
;    Written by NF 2007
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
;-

  
  if N_PARAMS() LT 1 or keyword_set(help) then begin 
     doc_library,'dustem_create_ionfrac'
     goto, the_end
  end

  

;PURPOSE :
;--------
; pour une Température et une densité donnée
; pour une liste de taille donnée (en nb de carbone ou angströms)
; donne les fractions de PAH neutre/cation/anion pour chaque taille
;
;INPUTS :
;-------
; gastemp = gas temperature
; sizes = PAHs size in Angstroms
; xe = electronic density relative to the hydrogen
; n_h = hydrogen density in cm-3
; n_elec = electronic density in cm-3
; file_isrf = ISRF file

; !run_ionfrac tells where to write the ionization fraction
; 1 : in MIX_xxx.DAT
; 2 : in SIZE_xxx.DAT

;Default values for the parameters
  gastemp = 100.d0
  n_elec = 3.d-4
;Check if these values have to be modified
  IF keyword_set(key) THEN BEGIN 
     a=where(key EQ 1,count1)
     b=where(key EQ 2,count2)
     IF count1 NE 0 then gastemp=(val(a))(0)
     IF count2 NE 0 then n_elec=(val(b))(0)
  ENDIF

;=== GET SIZES ===
;Check version of the model ...
;  IF getenv('DUSTEM_WHICH') NE 'VERSTRAETE' THEN BEGIN  
  IF !dustem_which NE 'VERSTRAETE' THEN BEGIN  
     message, "YOU SHOULD BE USING THE NEW VERSION OF THE MODEL AND THE INPUT FILES"
; ... and identify PAH0 and PAH1
  ENDIF ELSE BEGIN
    pah0 = (where(strmid((*!dustem_params).grains.type,0,4) eq 'PAH0'))(0)
    pah1 = (where(strmid((*!dustem_params).grains.type,0,4) eq 'PAH1'))(0)
  ENDELSE
;Get info on PAH0 and PAH1 sizes
  test_size = stregex( (*!dustem_params).keywords,'SIZE',/bool,/fold)
;   ... either from the SIZE_xxx.DAT files
  IF test_size THEN BEGIN
     pah0_nsize = (*(*!dustem_params).size(pah0)).nsize
     pah1_nsize = (*(*!dustem_params).size(pah1)).nsize
     pah0_sizes = (*(*!dustem_params).size(pah0)).size_st.a
     pah1_sizes = (*(*!dustem_params).size(pah1)).size_st.a
     if pah0_sizes ne pah1_sizes then message, "PAH0 and PAH1 should have the same size distribution. Please correct."
     sizes = pah0_sizes
  ENDIF ELSE BEGIN
;   ... or from the GRAIN.DAT
     pah0_nsize = (((*!dustem_params).grains.ndiscr)(pah0))(0)
     pah1_nsize = (((*!dustem_params).grains.ndiscr)(pah1))(0)
     alpha0 = (((*!dustem_params).grains.alpha)(pah0))(0)
     tmin0 = (((*!dustem_params).grains.tmin)(pah0))(0)
     tmax0 = (((*!dustem_params).grains.tmax)(pah0))(0)
     alpha1 = (((*!dustem_params).grains.alpha)(pah1))(0)
     tmin1 = (((*!dustem_params).grains.tmin)(pah1))(0)
     tmax1 = (((*!dustem_params).grains.tmax)(pah1))(0)
     if alpha0 ne alpha1 or tmin0 ne tmin1 or tmax0 ne tmax1 then message, "PAH0 and PAH1 should have the same size distribution. Please correct."
     sizes = fltarr(pah0_nsize)
     for i=0L,pah0_nsize-1 do begin
        sizes(i) = tmin0 * exp(i * alog(tmax0/tmin0) / (pah0_nsize-1)) * 1.d8
     endfor
  ENDELSE

;Compute other parameters
  Nc = (sizes/0.9)^2.           ;number of carbon atoms
  dim_n = pah0_nsize
  

;OUTPUT :
;--------
; fn = fraction of neutral PAH
; fc = fraction of cation PAH
; fa = fraction of anion PAH
  
;-----
; ISRF
;-----
  lam_isrf = (*!dustem_params).isrf.lambisrf
  isrf = (*!dustem_params).isrf.isrf
;convert into microns-1 and sort
  x = 1./lam_isrf
  isrf=isrf[sort(x)]
  x=x[sort(x)]
  dim_l = (size(x))(1)
  
;--------------------------------------
; absorption cross sections in UV (refer to L.Verstraete's HDR)
;--------------------------------------
  sig_uv_pah =  fltarr(dim_l)   ;in Mbarn/C
  
;below 5.9
  on = where(x lt 5.9)
  sig_uv_pah(on) = 6.6*(4.5/x(on)*(3.5*x(on))^2/((x(on)^2-4.5^2)^2+(3.5*x(on))^2))
;between 5.9 and 10.84
  on = where(x ge 5.9 and x lt 10.84)
  sig_uv_pah(on) = 6.6*(4.5/x(on)*(3.5*x(on))^2/((x(on)^2-4.5^2)^2+(3.5*x(on))^2)) + 0.85*(0.5392*(x(on)-5.9)^2+0.0564*(x(on)-5.9)^2)
;between 10.84 and 21.36
  on = where(x ge 10.84 and x lt 21.36)
  sig_uv_pah(on) = 18.6*(13.4*x(on))^2/((x(on)^2-12.8^2)^2+(13.4*x(on))^2)
;above 21.36
  on = where(x ge 21.36)
  sig_uv_pah(on) = 9.2*(21.36/x(on))^1.7
  
  C = fltarr(dim_n,dim_l)
  Ec = 8.2/sqrt(Nc)*1./sqrt(1.-sqrt(6./Nc)+2./Nc) ;in eV
  xc = Ec*1.6/(6.626d-1*3.)                       ;cutoff frequency in microns-1
  for j=0, dim_l-1 do begin
     for i=0, dim_n-1 do begin
        C(i,j) = 1./!pi*atan(1.d5*(x(j)/xc(i)-1.)^3) + .5 ;cutoff function
     endfor
  endfor
  
  sig_uv_pah_tab = fltarr(dim_n,dim_l)
  for j=0, dim_l-1 do begin
     for i=0, dim_n-1 do begin
        sig_uv_pah_tab(i,j) = sig_uv_pah(j) * C(i,j) * Nc(i) * 1.e6 ;in barn (good comparison with L.Verstraete's HDR)
     endfor
  endfor
  
;---------------------
;ionization parameters
;---------------------
  Icoronene = 7.3                         ; eV
  Icoronene = Icoronene*1.6/(6.626d-1*3.) ; microns-1
  Ipah_ev = fltarr(dim_n)
  Ipah_ev = 4.8+10.9/sizes         ; eV
  Ipah = Ipah_ev*1.6/(6.626d-1*3.) ; microns-1
  
  Yion = fltarr(dim_n,dim_l)
  for j=0, dim_l-1 do begin
     for i=0, dim_n-1 do begin
        Yion(i,j) = 0.8*exp(-0.00128*((12.-Icoronene)/(12.-Ipah(i))*(x(j)-12.))^4.) ; energies in microns-1
     endfor
  endfor
  
  E_isrf = (6.626d-7*3.)*x*1.d6/1.6 ; lam_isrf in eV, sort in ascending order
  n_ph = isrf*1./6.626d-27/E_isrf   ; photon density in s-1.cm-2.eV-1 sort in ascending order of E_isrf
  
  sig_ion = fltarr(dim_n,dim_l)
  for j=0, dim_l-1 do begin
     for i=0, dim_n-1 do begin
        sig_ion(i,j) = Yion(i,j) * sig_uv_pah_tab(i,j) * 1.d-24 ; cm2
     endfor
  endfor

  f = fltarr(dim_n,dim_l)
  for j=0, dim_l-1 do begin
     for i=0, dim_n-1 do begin
        f(i,j) = (double(sig_ion(i,j)*n_ph(j)))
     endfor
  endfor
  
;-----------------------------------
;definition of reaction efficiencies
;-----------------------------------
;recombination in cm3/s
  krec = 1.e-6*sqrt(300./gastemp) ; cm3/s
  krec_s = krec*n_elec            ; s-1
  
;electronic attachement
  ke = fltarr(dim_n)
  ke = 1.e-7*(Nc)^(3./4.)       ; cm3/s
  ke_s = ke*n_elec              ; s-1
  
;ionization
  kion = fltarr(dim_n)
  for i=0, dim_n-1 do begin
     on = where(E_isrf ge Ipah_ev(i) and E_isrf lt 13.6, count)
     kion(i) = int_tabulated(E_isrf(on), (f(i,*))(on)) ;BORNES !!!!
  endfor
  
;electronic photodetachement
  kp = fltarr(dim_n)
  Eae = 3                       ;energie d'affinite electronique en eV
  on_eae = where(E_isrf ge 3 and E_isrf le 13.6)
  E_isrf2 = E_isrf(on_eae)      ;energie isrf de 3 a 13.6 eV
  n_ph2 = (n_ph)(on_eae)        ;densite de photon entre 3 et 13.6 eV
  
  Nph = int_tabulated(E_isrf2, n_ph2) ;nombre de photon d'energie comprise entre 3 et 13.6 eV
;kp = sig_uv_pah*Yion*Nc*Nph     ;taux de photodetachement electronique
  sigion = 1.6d-18
  kp = sigion * Nc * Nph
  
  
  k_str = {krec:0., ke:fltarr(dim_n), kion:fltarr(dim_n), kp:fltarr(dim_n)}
  k_str.krec = krec_s & k_str.ke = ke_s & k_str.kion = kion & k_str.kp = kp
  
;--------------------------------------------
;fractions of PAH neutral, cations and anions
;--------------------------------------------
  fn = (kp*krec*n_elec)/(kp*krec*n_elec+kp*kion+krec*ke*n_elec*n_elec)
  fc = (kp*kion)/(kp*krec*n_elec+kp*kion+krec*ke*n_elec*n_elec)
  fa = (krec*ke*n_elec*n_elec)/(kp*krec*n_elec+kp*kion+krec*ke*n_elec*n_elec)
  
  f0 = fn
  f1 = fa + fc

; WRITE RESULTS INTO PROPER "FILES"
;   in MIX_xxx.DAT
  if abs(!run_ionfrac) eq 1. then begin
     (*(*!dustem_params).mix(pah0)).fmix = f0
     (*(*!dustem_params).mix(pah1)).fmix = f1
;     dustem_write_mix, getenv('DUSTEM_DAT'), (*!dustem_params).mix
     dustem_write_mix, !dustem_dat, (*!dustem_params).mix
  endif
;   in SIZE_xxx.DAT
  if abs(!run_ionfrac) eq 2. then begin
     (*(*!dustem_params).size(pah0)).size_st.f_mix = f0
     (*(*!dustem_params).size(pah1)).size_st.f_mix = f1
;     dustem_write_size_lv, getenv('DUSTEM_DAT'), (*!dustem_params).size
     dustem_write_size_lv, !dustem_dat, (*!dustem_params).size
  endif

the_end:

  return,0
  
end