FUNCTION DUSTEM_CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=chi, fname=fname, help=help

;+
; NAME:
;   dustem_create_rfield
;  
; PURPOSE:
;   generates a file with the radiation field in format suitable for
;   DUSTEM (ISRF.DAT). Units erg/cm2/s/Hz (4*!pi*I_nu)
;   RFIELD = ISRF + WDIL*PI*BB(TEMP) '
;  
; CATEGORY:
;    DustEMWrap, Distributed, LowLevel, Helper
;  
; CALLING SEQUENCE:
;   mathis = DUSTEM_MATHIS_FIELD(x,unit=unit)
;
; INPUTS:
; TEMP -- blackbody temperature (can be an array). If 0 only ISRF.
;
;  OPTIONAL INPUT PARAMETERS:
;  X  --  wavelength -grid for the ISRF in microns. Default is 200 pts over 0.01-10^5 microns.
;           (including WCUT point). You want to include WCUT for accurate edges.
; ISRF -- if set to 0 no ISRF is added, 1 is Mathis (default), 2 is Habing
; WDIL --  blackbody dilution factor (can be an array)
; FNAME -- filename for ISRF.DAT
; WCUT --- for wave < wcut radiation field is 0 
;
; OUTPUTS:
;
; OPTIONAL OUTPUT PARAMETERS:
;     G0   -- factor flux(6-13.6eV) wrt. Mathis field
;     CHI  -- scaling factor at 100nm wrt. Mathis field
;
; ACCEPTED KEY-WORDS:
;    help = print this help
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURES AND SUBROUTINES USED:
;   
; EXAMPLES
;    tt = dustem_create_rfield([2d4,5d4],wdil=[1.d-14,1.d-16],x=x,fname='ISRF.DAT')
;  
; MODIFICATION HISTORY:
;    Written by LV 2009
;    Initially distributed with the Fortran, incorporated into DustEMWrap 2022
;    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_rfield'
     rfield=0
     goto, the_end
  end
  
  ;;    print,'------------------------------------------------------------------------------------------------------'
  ;;    print,'FUNCTION DUSTEM_CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=chi, fname=fname'
  ;;    print,'------------------------------------------------------------------------------------------------------'
  ;;    print,''
  ;;    print,' generates the radiation field for DUSTEM (ISRF.DAT)'
  ;;    print,' in erg/cm2/s/Hz (4*!pi*I_nu)'
  ;;    print,' RFIELD = ISRF + WDIL*PI*BB(TEMP) '
  ;;    print,' ISRF from Mathis et al. 1983, A&A 128, 212'
  ;;    print,' NB: if blackbody is isotropic set WDIL to 4 (to get 4*!pi factor)'
  ;;    print,''  
  ;;    print,' TEMP (I): blackbody temperature (can be an array). If 0 only ISRF.'
  ;;    print,' X    (I): X-grid for the ISRF in microns. Default is 200 pts over 0.01-10^5 microns.'
  ;;    print,'           (including WCUT point). You want to include WCUT for accurate edges.'
  ;;    print,' ISRF (I): if set to 0 no ISRF is added, 1 is Mathis (default), 2 is Habing'
  ;;    print,' WDIL (I): blackbody dilution factor (can be an array)'
  ;;    print,' FNAME(I): filename for ISRF.DAT'
  ;;    print,' WCUT (I): for wave < wcut radiation field is 0 '
  ;;    print,' G0   (O): factor flux(6-13.6eV) wrt. Mathis field '
  ;;    print,' CHI  (O): scaling factor at 100nm wrt. Mathis field'
  ;;    print,''
  ;;    print,'Example: tt = dustem_create_rfield([2d4,5d4],wdil=[1.d-14,1.d-16],x=x,fname=''ISRF.DAT'')'
  ;;    print,''
  ;;    print,' Created Aug. 2009, L. Verstraete, IAS'
  ;;    print,' Force the WCUT point, May 2011, LV'
  ;;    print,''
  ;;    print,'------------------------------------------------------------------------------------------------------'
  ;;    RETURN,0
  ;; endif

; inits 
  unit='W'
  ly_limit =  9.11267101235d-2
  xi = [ ly_limit, 2.07d-1 ] ; for G0 6-13.6 eV in microns
  x1 = 1.d-1 ; for CHI

  if n_elements(WCUT) EQ 0 then wcut=ly_limit  
  if n_elements(x) EQ 0 then begin
     nx = 199
     xb = [ 1.d-2, 1.d5 ] ; wave boundaries in microns
     dx = (alog(xb(1))-alog(xb(0))) / (nx-1)
     x = alog(xb(0)) + dindgen(nx)*dx
     x = exp(x)
     x = [x, wcut] ; add wcut to avoid coarse edge
     x = x(UNIQ(x,SORT(x)))
     nx = n_elements(x)
     if x(nx-1) NE xb(1) then x(nx-1)=xb(1) ; check rounding
  endif else begin
     print,'(W) DUSTEM_CREATE-RFIELD : using input wave x'
     x = double(x)
     nx = n_elements(x)
     ix = WHERE( ABS(x-wcut)/x LE 0.01, cx )
     if cx EQ 0 then begin
        print,'(W) DUSTEM_CREATE_RFIELD: your x-grid does not contain wcut.'
        print,'                   Should be included if radiation is 0 below wcut.'
     endif 
  endelse

  rfield = dblarr(nx)
; get Habing for normalization
  rhabing = DUSTEM_HABING_FIELD(x,unit=unit)
  rhabing =  1.d3*x*rhabing/3.d14 ; erg/cm2/s/Hz

; get isrf
  rmathis = DUSTEM_MATHIS_FIELD(x,unit=unit)
  rmathis =  1.d3*x*rmathis/3.d14 ; erg/cm2/s/Hz

  if n_elements(ISRF) EQ 0 then begin
     print,'(W) DUSTEM_CREATE_RFIELD : ISRF keyword not set, using default Mathis'
     ISRF=1
  endif
  
 if ISRF EQ 0 then begin 
     print,'(W) DUSTEM_CREATE_RFIELD : not adding diffuse ISM ISRF component'
 endif else if ISRF EQ 1 then begin 
     print,'(W) DUSTEM_CREATE_RFIELD : adding Mathis ISRF'
     rfield = rmathis
  endif else if ISRF EQ 2 then begin 
     print,'(W) DUSTEM_CREATE_RFIELD : adding Habing ISRF'
     rfield = rhabing
  endif

; get blackbody  
  ntemp = n_elements(TEMP)
  if ntemp GT 0 then begin 
     bb = DUSTEM_BLACKBODY( x, temp, unit=unit, wdil=wdil )
     bb = bb * 1.d3 * !pi * x^2 / 3.d14   ; erg/cm2/s/Hz
     print,'(W) DUSTEM_CREATE_RFIELD : adding BB with T= ',temp, format='(A38,10(1E10.4,1x))'
     print,'                dilution factor wdil= ',wdil, format='(A38,10(1E10.4,1x))'
  endif else bb=0.D0
  rfield = rfield + bb

; apply cut
  ix = WHERE( x LT WCUT, cx )
  if cx GT 0 then rfield(ix)=0.d0

; get G0
 ig = where( x GE xi(0) AND x LE xi(1), cg)
 if cg GT 0 then begin
    xx=x(ig)
    yy=rfield(ig)
    rr=rmathis(ig)
    s1 = TOTAL( (yy(1:cg-1)+yy(0:cg-2))*(xx(1:cg-1)-xx(0:cg-2))/2.D0 )
    s2 = TOTAL( (rr(1:cg-1)+rr(0:cg-2))*(xx(1:cg-1)-xx(0:cg-2))/2.D0 )
    g0 = s1/s2
 endif else g0 = 0.d0
 print,'(W) DUSTEM_CREATE_RFIELD : G0 =',g0

; get chi
 ig = where( abs(2.*(x-x1)/(x+x1)) LE 1.d-1, cg) 
 if cg GT 0 then begin 
   tmp = min( abs(x(ig)-x1), i_min )   
   ig = ig(i_min)
   ig = ig(0)
   chi = rfield(ig) / rmathis(ig)  
 endif else chi = 0
 print,'(W) DUSTEM_CREATE_RFIELD : chi =',chi

; write file
 if n_elements(FNAME) NE 0 then begin
   OPENW, iu, fname, /get_lun
   printf, iu, '# DUSTEM: exciting radiation field featuring '
   if ISRF EQ 1 then printf, iu, '# Mathis ISRF'
   if NTEMP GT 0 then begin 
      a =  '# Blackbody with     T='
      a1 = '# dilution factor wdil='
      for i = 0, ntemp - 1 do begin
         a = a   + string( format='(2x,1E11.4)', temp(i) )
         a1 = a1 + string( format='(2x,1E11.4)', wdil(i) )
      endfor
      printf, iu, a
      printf, iu, a1
   endif
   printf, iu, '# Nbr of points'
   printf, iu, '# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)'
   printf, iu, nx, format='(i4)'
   for i=0,nx-1 do begin 
    printf, iu, x(i), rfield(i), format='(2(1E13.6,2x))'
   endfor
   FREE_LUN, iu
   print,'(W) DUSTEM_CREATE_RFIELD: radiation field written in ', strtrim(fname,2)
 ENDIF

 the_end:
 RETURN, RFIELD
END