;+
; NAME:
;	Rad_Trans_1D
; PURPOSE:
;	Compute single temperature dust emission/absorption spectrum
;	with line of sight (LOS) absorption, in Janskys/arcsec^2.
;	If third argument is present, compute partial derivatives for fitting.
; CALLING:
;	IR_spec = Rad_Trans_1D( wvd, Params, Pder )
; INPUTS:
;	wvd = wavelength array, in microns.
;	Params = [ cLOS, T1, csrc1, T2, csrc2 ]
;	Pder = matrix of partial derivatives of parameters at wavelengths.
; KEYWORDS:
;	/FULL_SPECTRUM : use full resoltuion wavelength array in common block.
;	EXTINCT = array of LOS extinction factor vs. wavelength (returned).
;	/JANSKYS : return IR intensity in Jy/arcsec^2, otherwise
;		default is pure energy [nu.I(nu)]: W/m^2/arcsec^2.
; OUTPUTS:
;	Returns the IR intensity in W/m^2 or Janskys, per arcsec^2.
;
;	I(w)=  B(w,T1)*(1-exp[-TAU_src1])*exp[-TAU_LOS] 
;	     + B(w,T2)*(1-exp[-TAU_src2])*exp[-TAU_LOS]
; COMMON BLOCKS:
;	common IR_spectrum_comp, Fsil_LOS, Fsil_src1, Fsil_src2
;	common kappa_interp, wki, kigrf, kisil
;	common kappa_full, wkf, kgrf, ksil
;	common Rad_Trans_1D, cLOS
; EXTERNAL CALLS:
;	pro match
; PROCEDURE:
;
; HISTORY:
;	Written: Frank Varosi & Eli Dwek, NASA/GSFC, 1996.
;-

function Rad_Trans_1D, wvd, Params, Pder, FULL_SPECTRUM=fullspec, $
		EXTINCT=ext_LOS, JYS=janskys, T1_ONLY=T1_only, T2_ONLY=T2_only

   common IR_spectrum_comp, Fsil_LOS, Fsil_src1, Fsil_src2
   common kappa_interp, wki, kigrf, kisil
   common kappa_full, wkf, kgrf, ksil
   common RT_1D_fixed, cLOS, Tbb

 	Nwv = N_elements( wvd )
	Np = N_elements( Params ) 

	CASE Np OF
		1: begin
			csrc = Params(0)
			if N_elements( cLOS ) ne 1 then cLOS=0.
			if N_elements( Tbb ) ne 1 then Tbb = 200.
		     end
		2: begin
			Tbb = Params(0)
			csrc = Params(1)
			if N_elements( cLOS ) ne 1 then cLOS=0.
		     end
		3: begin
			cLOS = Params(0)
			Tbb = Params(1)
			csrc = Params(2)
		     end
		4: begin
			Tbb = Params(0)
			csrc = Params(1)
			Tbb2 = Params(2)
			csrc2 = Params(3)
			if N_elements( cLOS ) ne 1 then cLOS=0.
		     end
		5: begin
			cLOS = Params(0)
			Tbb = Params(1)
			csrc = Params(2)
			Tbb2 = Params(3)
			csrc2 = Params(4)
		     end
	  ENDCASE

	if keyword_set( fullspec ) then begin

		wvs = wkf
 		Nwv = N_elements(wvs)
		ksrc1 = Fsil_src1 * ksil + (1-Fsil_src1) * kgrf
		kLOS = Fsil_LOS * ksil + (1-Fsil_LOS) * kgrf

	 endif else begin

		match, wvd, wki, md, mki, COUNT=nmatch
		if( nmatch LT Nwv ) then begin
			message,"stored interpolate wavelengths do not match " $
				+"those for data, must restore Dust_Kappa.idl" $
				+ string(7b),/INFO
			stop
		   endif
		wvs = wvd
		ksrc1 = Fsil_src1 * kisil(mki) + (1-Fsil_src1) * kigrf(mki)
		kLOS = Fsil_LOS * kisil(mki) + (1-Fsil_LOS) * kigrf(mki)
	  endelse

if !DEBUG GT 2 then help,cLOS,Tbb,csrc
if !DEBUG GT 4 then stop

	arcsec2 = ( !DTOR/60/60 )^2
	wavnum = 1/wvs
janskys=1
	if keyword_set( janskys ) then begin
		hcj = 2 * !cv.h * (1e-7 * 1e26) * !cv.c * 1e-2
		bhc = arcsec2 * hcj * 1e18 * wavnum^3
	 endif else begin
		hcw = 2 * !cv.h * (1e-7 * !cv.c * 1e4) * !cv.c * 1e-2 *1e10
		bhc = arcsec2 * hcw * 1e18 * wavnum^4
	  endelse

	hcok = !cv.h * 1e4 * !cv.c / !cv.k
	ext_LOS = EXP( -cLOS*kLOS )

	if keyword_set( T2_only ) then goto,T2

;First temperature B.B.:

	BT = fltarr(Nwv)
	xT = wavnum * ( hcok / Tbb )
	w = where( xT LT 88., nex )

	if nex gt 0 then begin
		eT = EXP( xT(w) )
		BT(w) = bhc(w)/(eT-1)
	   endif

	BT = BT * ext_LOS
	ext_SRC = EXP( -csrc*ksrc1 )
	Flux = BT * (1 - ext_SRC)

	if N_params() GE 3 then begin

		Pder = fltarr( Nwv, Np )

		if (Np eq 1) then begin
			Pder(*,0) =  BT * ext_SRC * ksrc1
			return, Flux
		   endif

		exprat = replicate( 1.0, Nwv )
		if nex gt 0 then exprat(w) = eT/(eT-1)

		if (Np eq 2) or (Np eq 4) then begin
			Pder(*,0) =  Flux * (xT/Tbb) * exprat
			Pder(*,1) =  BT * ext_SRC * ksrc1
		 endif else begin
			Pder(*,0) = -Flux * kLOS
			Pder(*,1) =  Flux * (xT/Tbb) * exprat
			Pder(*,2) =  BT * ext_SRC * ksrc1
		  endelse
	   endif

	if (Np LE 3) or keyword_set( T1_only ) then return, Flux

T2:	;optional second temperature B.B:

	if keyword_set( fullspec ) then $
		ksrc2 = Fsil_src2 * ksil + (1-Fsil_src2) * kgrf $
	  else	ksrc2 = Fsil_src2 * kisil(mki) + (1-Fsil_src2) * kigrf(mki)

	BT = fltarr(Nwv)
	xT = wavnum * ( hcok / Tbb2 )
	w = where( xT LT 88., nex )

	if nex gt 0 then begin
		eT = EXP( xT(w) )
		BT(w) = bhc(w)/(eT-1)
	   endif

	ext_SRC = EXP( -csrc2*ksrc2 )
	BT = BT * ext_LOS
	Flux2 = BT * (1 - ext_SRC)
	if N_elements( Flux ) LE 0 then Flux = 0
	Flux = Flux + Flux2 

	if N_params() GE 3 then begin
		exprat = replicate( 1.0, Nwv )
		if nex gt 0 then exprat(w) = eT/(eT-1)
		if (Np eq 4) then begin
			Pder(*,2) =  Flux2 * (xT/Tbb2) * exprat
			Pder(*,3) =  BT * ext_SRC * ksrc2
		 endif else begin
			Pder(*,0) = -Flux * kLOS
			Pder(*,3) =  Flux2 * (xT/Tbb2) * exprat
			Pder(*,4) =  BT * ext_SRC * ksrc2
		  endelse
	   endif

return, Flux
end