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

  if N_PARAMS() LT 1 or keyword_set(help) then begin 
     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

 RETURN, RFIELD
END