;+
; NAME:
;	IR_dust_model
; PURPOSE:
;	Model the infrared emission from dust in an image-spectrum data cube.
;	Each line of sight is modeled independentaly (no scattering) by
;	fitting a simple 1D radiative transfer model of single temperature
;	dust emission and absorption to the spectrum at each image pixel.
;	The user can modify the model control parameters with X-widget GUI
;	that is automatically presented, or just accept the default values.
;	Assumes that IDL/XDR save file "Dust_Kappa.idl" containing the dust
;	absorption coefficients is in current directory.
;
; CALLING:
;	specfits = IR_dust_model( wavelens, images )
;
; INPUTS:
;	wavelens = array of wavelengths, in microns.
;
;	images = 3D array (stack of images), one image for each wavelength,
;		calibrated in Jansky's / arcsec^2.
;
; KEYWORDS:
;	BOX_AVERAGE = radius of moving box of pixels to average over.
;	/PLOT : plot the data & fit result at each pixel.
;
; OUTPUTS:
;	Function returns array of structures (one for each image pixel)
;	containing spectrum fit results: Temperature, mass of dust, etc...
;
; COMMON BLOCKS:
;	common IR_spectrum_fit, fitparms, wghts
;
; EXTERNAL CALLS:
;	pro X_Var_Edit
;	function IR_spectrum_fit
; PROCEDURE:
;	Call function IR_spectrum_fit in a double loop.
; HISTORY:
;	Written: Frank Varosi, HSTX @ NASA/GSFC, 1996.
;-

function IR_dust_model, wavelens, images, PLOT=plotit, BOX_AVERAGE=boxav

   common IR_spectrum_fit, fitparms, wghts

	sz = size( images )

	if sz(0) LE 0 then begin
		print," "
		message,"missing 2nd argument",/INFO
		print,"syntax:  sf = IR_dust_model( wavelens, images )"
		print,"syntax:  sf = IR_dust_model( wavelens, fluxes, /PLOT )"
		return,0
	   endif

	if sz(sz(0)) NE N_elements( wavelens ) then begin
		print," "
		message,"calling argument mismatch",/INFO
		print,"syntax:  sf = IR_dust_model( wavelens, images )"
		print,"syntax:  sf = IR_dust_model( wavelens, fluxes, /PLOT )"
		print," images or fluxes can be 1D, 2D, or 3D array,"
		print," and # elements in wavelens must match # elements"
		print," in 1st, 2nd, or 3rd dimension respectively."
		return,0
	   endif

	restore,"Dust_Kappa.idl",/VERB
	zero = IR_spectrum_fit( wavelens, 0, DUST=Dust_Kappa )
	X_Var_Edit, fitparms

	if sz(0) EQ 1 then return, IR_spectrum_fit( wavelens, images, P=plotit )

	sf = IR_spectrum_fit( wavelens, fltarr( N_elements( wavelens ) ) )
	nx = sz(1)
	ny = sz(2)
	sx = 0
	Lx = nx-1
	sy = 0
	Ly = ny-1

	if keyword_set( boxav ) then begin
		b = round( boxav > 1 )
		sx = sx + b
		Lx = Lx - b
		sy = sy + b
		Ly = Ly - b
	   endif

	if sz(0) EQ 2 then begin

		specfit = replicate( sf, nx )
		if keyword_set( boxav ) then npb = (2*b+1)

		for ix=sx,Lx do begin

		   if keyword_set( boxav ) then begin
			fluxes = total( images(ix-b:ix+b,*), 1 )/npb
		     endif else fluxes = reform( images(ix,*) )

		   specfit(ix) = IR_spectrum_fit( wavelens, fluxes, PL=plotit )
		   print,ix, specfit(ix).nfitp, specfit(ix).parms

		 endfor

		return, specfit
	   endif

	specfit = replicate( sf, nx, ny )
	if keyword_set( boxav ) then npb = (2*b+1)^2

	for iy=sy,Ly do begin
	   for ix=sx,Lx do begin

		if keyword_set( boxav ) then begin
		   fluxes = $
		      total( total( images(ix-b:ix+b,iy-b:iy+b,*), 1 ), 1 )/npb
		 endif else fluxes = reform( images(ix,iy,*) )

		specfit(ix,iy) = IR_spectrum_fit( wavelens, fluxes, PL=plotit )
		print,ix,iy, specfit(ix,iy).nfitp, specfit(ix,iy).parms

	     endfor
	  endfor

return, specfit
end