;+
; NAME:
;      DT_DUST
;
;
; PURPOSE:
;      The DT_DUST procedure reads the dust brightness per pixel 
;      and the estimated dust temperature per pixel from the 
;      IRAS/DIRBE 100 micron maps for the given values of 
;      Right Ascension and Declination, then works out the 
;      extrapolation of the dust brightness to any frequency.
;      The dust emission is modeled as a grey body with emissivity
;      index alpha. 
;      The temperature of the grey body is read for every pixel 
;      from the maps.
;      The result is converted into equivalent CMB temperature 
;      fluctuation: this is calculated by interpreting the dust 
;      emission as a linear perturbation term to a 2.276 K blackbody 
;      emission and calculating the corresponding temperature fluctuation.
;      Printed on screen are the mean, variance and standard 
;      deviation calculated on the entire scan.
;
;
; CALLING SEQUENCE:
;      DT_DUST,ra,dec,nu,alpha,value
;
; 
; INPUTS:
;      ra:  Right Ascension (hours)
;      dec: declination (degrees) 
;      nu:  frequency of interest (Hz)
;      alpha: spectral index of the grey body
;
; OUTPUTS:
;      value: a vector containing the values of the 
;             equivalent CMB temperature fluctuation 
;             per sample due to dust, in K, for the given scan
;
; EXAMPLE:
;      Calculate the dust contamination per sample for the set of
;      coordinates 'ra', 'dec', at frequency 150 GHz, for a spectral index
;      1.6, and store the result in the vector 'dust':
;
;      IDL> DT_DUST,ra,dec,150.e9,1.6,dust
;
; PROCEDURES USED:
;      DUST_GETVAL: reads dust intensity and temperature from the maps
;      EULER: converts between coordinates systems
;      PRECESS: precesses equinoxes
;
; MODIFICATION HISTORY:
;	Created, Amedeo Balbi, April 1998 
;       Bug fixed December 1998 (input ra in PRECESS were in hours instead
;       of degrees)
;
;-
pro DT_DUST,ra,dec,nu,alpha,value

; converting from RA-dec to galactic coordinates

np=n_elements(ra)
galb=fltarr(np) 
gall=fltarr(np)
ra0=ra*15.
dec0=dec
PRECESS,ra0,dec0,2000,1950
EULER,ra0,dec0,gall,galb,1

; definition of constants: 
h=6.62618e-27                   ; Planck's constant (erg/s)
c=2.99792458e10                 ; speed of light (cm/s)
k=1.38066e-16                   ; Boltzmann constant (erg/K)
pi=3.1415926535                 ; pi
lambda0=100./1.e4               ; IRAS wavelength in cm
nu0=c/lambda0                   ; IRAS frequency in Hz

; looking at the dust brightness from the maps 
value=DUST_GETVAL(gall,galb,map='I100',/noloop)

; looking at the corresponding dust temperature from the maps
t0=DUST_GETVAL(gall,galb,map='T',/noloop)

; working out extrapolation
index=3.0+alpha
value=value*(nu/nu0)^index * (exp(h*nu0/k/t0)-1.0)/(exp(h*nu/k/t0)-1.0)
value=value/1.e17               ; 1 MJy/Sr = 1.e-17 erg/s/cm2/Hz/Sr

; converting to equivalent CMB temperature fluctuation
tcmb=2.726
factor=h*nu/k/tcmb
x=exp(factor)
factor=2.*factor*h*nu^3/c^2/tcmb
dbdt=factor*x/(x-1.)
dbdt=dbdt/(x-1.)
tt=value/dbdt
value=tt

mean=total(value)/n_elements(value)
qdev=(value-mean)^2
var=total(qdev)/(n_elements(qdev)-1.)
sdev=sqrt(var)
print,'mean: ',mean,' K'
print,'var: ',var,' K^2'
print,'sdev: ',sdev,' K'

return
end