FUNCTION HABING_FIELD, x, unit=unit
; computes the Habing interstellar radiation field SED (ISRF)
; in W/m2 (4*!pi*nu*I_nu)
; (from Draine & Bertoldi 1996)
; X	(I): X-grid for the ISRF
; UNIT	(I): unit of the ISRF. Choices are 
;	     'W': wave in um [Default]
;	     'F': frequency in cm-1
;	     'E': energy in eV 

 if n_elements( UNIT ) EQ 0 then unit = 'W' $
 else unit = strupcase( strcompress( unit,/rem) )

 x = double( x )

 CASE unit of 

   'W': x3 = 1.d1 * x 

   'F': x3 = 1.d5 / x

   'E': x3 = 12.4 / x

 ENDCASE

 field = - x3^3*4.1667 + x3^2*12.5 - x3*4.3333 
 field = 1.d-1 * 3.d8 * 1.d-14 * field 

 i_neg = where( field LT 0.d0, c_neg )
 field( i_neg ) = 0.d0

 RETURN, field
END                             ;FUNCTION HABING_FIELD


FUNCTION MATHIS_FIELD, x, unit=unit
; computes the Mathis interstellar radiation field SED (ISRF)
; in W/m2 (4*!pi*nu*I_nu)
; from Mathis et al. 1983, A&A 128, 212
; X	(I): X-grid for the ISRF
; UNIT	(I): unit of the ISRF. Choices are 
;	     'W': wave in um [Default]
;	     'F': frequency in cm-1
;	     'E': energy in eV

 if n_elements( UNIT ) EQ 0 then unit = 'W' $
 else unit = strupcase( strcompress( unit,/rem) )

 ly_limit = 9.11267101235d-2

 x = double( x )

;
; visible part 
;
 wdil = 4.d0 * [ 4.d-13, 1.d-13, 1.d-14]	; Mathis definition
						; 4*!pi*Wdil*B_nu
 field = !pi * x * BLACKBODY( x, [ 3.d3, 4.d3, 7.5d3 ], wdil=wdil, unit=unit )

;
; UV part
;

; first convert to (lambda / 1e3 AA)
 CASE unit of 

   'W': x3 = 1.d1 * x 

   'F': x3 = 1.d5 / x

   'E': x3 = 12.4 / x

 ENDCASE
 
 il = where( x3 LE 2.46, cl ) 
 if cl GT 0 then field(il) = 0.D0

 il = where( x3 GE 1d1*ly_limit AND x3 LT 1.11, cl )
 if cl GT 0 then field(il) = 1.4817d-6 * x3(il)^(4.4172)
 il = where( x3 GE 1.11 AND x3 LT 1.34, cl ) 
 if cl GT 0 then field(il) = 2.0456d-6 * x3(il)
 il = where( x3 GE 1.34 AND x3 LT 2.46, cl )
 if cl GT 0 then field(il) = 3.3105d-6 * x3(il)^(-0.6678)
  

 RETURN, field
END                             ;FUNCTION MATHIS_FIELD


FUNCTION CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=chi, fname=fname

  if N_PARAMS() LT 1 then begin 
     print,'------------------------------------------------------------------------------------------------------'
     print,'FUNCTION 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 = 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) 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) 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 = HABING_FIELD(x,unit=unit)
  rhabing =  1.d3*x*rhabing/3.d14 ; erg/cm2/s/Hz

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

  if n_elements(ISRF) EQ 0 then ISRF=1
  if ISRF EQ 1 then begin 
     print,'(W) CREATE_RFIELD : adding Mathis ISRF'
     rfield = rmathis
  endif else if ISRF EQ 2 then begin 
     print,'(W) CREATE_RFIELD : adding Habing ISRF'
     rfield = rhabing
  endif

; get blackbody  
  ntemp = n_elements(TEMP)
  if ntemp GT 0 then begin 
     bb = BLACKBODY( x, temp, unit=unit, wdil=wdil )
     bb = bb * 1.d3 * !pi * x^2 / 3.d14   ; erg/cm2/s/Hz
     print,'(W) 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) 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) 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) CREATE_RFIELD: radiation field written in ', strtrim(fname,2)
 ENDIF

 RETURN, RFIELD
END