FUNCTION DUSTEM_BLACKBODY, x, temp, unit=unit, wdil=wdil, check=check, g0=g0, chi=chi, $
		     bb_to_hab=bb_to_hab, NORM=norm
  if n_params() EQ 0 then begin 
     print,'------------------------------------------------------------------------------------------'
     print,'FUNCTION DUSTEM_BLACKBODY, x, temp, unit=unit, wdil=wdil, NORM=norm, check=check, g0=g0, chi=chi'
     print,'------------------------------------------------------------------------------------------'
     print,' computes a (or a sum of) blackbody brightness(es) B(x,T) in W/m2/(x unit)/sr '
     print,' X	(I): blackbody abscissa in cm-1, GHz, microns or eV'
     print,' TEMP  (I): blackbody temperature (can be an array)'
     print,' WDIL  (I): dilution factor for each TEMP. [Default=1].'
     print,' UNIT  (I): blackbody x unit ''K'' is B_nu (W/m2/cm-1/sr),   x in cm-1 '
     print,'                             ''F'' is B_nu (W/m2/Hz/sr),     x in GHz [default]'
     print,'	                         ''W'' is B_lambda (W/m2/um/sr), x in microns.' 
     print,'	 	  	         ''E'' is B_E (W/m2/eV/sr), x in eV'
     print,' NORM  (I): overall normalizing factor, returns NORM*BB/MAX(BB)'
     print,' CHECK (O): check energy conservation for the BB (power is sigma*T^4)'
     print,' CHI   (O): scaling factor @ 1000 A in ISRF unit (1.6e-3 erg/cm2/s)'
     print,' G0    (O): scaling factor for integrated 5-13.6 eV flux in ISRF unit.'
     print,''
     print,' Created 2001, L Verstraete IAS'
     print,' More units and checks, Oct 2010 LV'
     print,'---------------------------------------------------------------------------------------'
     return,0
  endif
 nx = n_elements( X )
 nt = n_elements( TEMP )
 nw = n_elements( WDIL )
 if nw EQ 0 then begin 
   wdil = dblarr(nt) + 1. 
 endif else begin 
   if nw LT nt then begin 
     print,'DUSTEM_BLACKBODY (F): WDIL and TEMP have different sizes'
     return, 0
   endif 
 endelse
 
 if n_elements( UNIT ) EQ 0 then unit='F' $
 else begin 
  unit = strupcase(strcompress(unit,/rem))
  if unit NE 'K' AND unit NE 'F' AND unit NE 'W' AND unit NE 'E' then begin 
    print,'DUSTEM_BLACKBODY (W): not a valid BB unit --> make a B_nu'
    unit = 'F'
  endif
 endelse

 if n_elements( CHECK ) GT 1 then check = 0

 x = double(x) 
 temp = double(temp) 
 wdil = double(wdil)

 xpi = 3.1415926535897932385d
 clight = 2.9979246d08          ; m/s
 kb = 1.3806503d-23
 hp = 6.62606876d-34
 hck = hp*clight / kb
 ze  = 1.60217653d-19           ; electron charge
 hev = hp/ze                    ; h in eV unit

 xly = 12.4
 isrf0 = 1.6d-6         ; 4*pi*nu*I_nu (W/m2) in solar neighborhood

;
; loop on temperatures
; 
 bbt = dblarr(nx)
 acheck = 0.

 FOR j = 0, nt-1 do begin 

    xi = [6., 13.6]

    CASE UNIT OF

       'K': begin      
          be = exp( -(1d2*x) * hck / temp(j) )
          bb = be * (1d2*x)^3 / (1.d0 - be)
          bb = bb * 2d0 * hp * clight^2 
          bb = bb * 1d2         ; m-1 to cm-1
          xi = xi / hev/clight/1d2   
          x1 = xly / hev/clight/1d2 
       end

       'F': begin      
          be = exp( -(1d9*x) * hp/kb/ temp(j) )
          bb = be * (1d9*x)^3 / (1.d0 - be)
          bb = bb * 2d0 * hp / clight^2
          xi = xi / hev / 1d9
          x1 = xly / hev / 1d9
       end

       'W': begin
          be = exp( -hp*clight/kb/ (1d-6*x) / temp(j) )
          bb = be / (1d-6*x)^5 / (1.d0 - be)
          bb = bb * 2d0 * hp * clight^2 
          bb = 1d-6 * bb        ; m-1 to micron-1
          xi = hev*clight*1d6 / xi
          x1 = hev*clight*1d6 / xly
       end

       'E': begin
          be = exp( -(ze*x)/kb / temp(j) )
          bb = be * (ze*x)^3 / (1.d0 - be)
          bb = bb * 2d0 / hp^3 / clight^2 
          bb = bb * ze          ; J-1 to eV-1
          xi = xi
          x1 = xly
       end

    ENDCASE

    ck = 0.
    for jx = 1L,nx-1 do $
       ck = ck + ( bb(jx)+bb(jx-1) )*( x(jx)-x(jx-1) ) / 2.d0
    ck = xpi* ck / 5.6697d-8 / temp(j)^4  
    
    if n_elements( CHECK ) NE 0 then begin 
       if check EQ 1 then $
          print,'DUSTEM_BLACKBODY (W): integrated BB('+strtrim(temp(j),2)+$
                ' K) is '+ strtrim(ck,2)+ ' times Sigma*T^4'
    endif
    acheck = [ acheck, ck ]
    
    bbt = bbt + bb*wdil(j)
    
 ENDFOR
 
 check = acheck(1:*)

; get G0
 ig = where( x GE MIN(xi) AND x LE MAX(xi), cg)
 if cg GT 0 then begin
    xx=x(ig)
    yy=xpi*x(ig)*bbt(ig)
    s1 = TOTAL( (yy(1:cg-1)+yy(0:cg-2))*(xx(1:cg-1)-xx(0:cg-2))/2.D0 )
    g0 = s1/isrf0
 endif else g0 = 0.d0

; 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 = xpi * x(ig) * bbt(ig) / isrf0
 endif else chi = 0

 if n_elements( NORM ) NE 0 then bbt = norm * bbt/max(bbt)
                                                
the_end:
 RETURN, BBT
END                             ;FUNCTION BLACKBODY