;+
; NAME:
;	analyze_fluxset
; PURPOSE:
;	Set fluxes at 2 points in stack of mosaic images,
;	for comparing (mosaic math) different wavelengths,
;	or averaging mosaics of same wavelength.
; CALLING:
;	analyze_fluxset
; INPUT & OUTPUT:
;	all thru common math_mosaics.
; COMMON BLOCKS:
;	common math_mosaics, math_List, math_images, math_imscaled
; EXTERNAL CALLS:
;	pro match_mos_apply
;	pro box_draw2
;	pro box_erase2
; HISTORY:
;	Written: Frank Varosi NASA/GSFC 1990.
;	F.V. 1991, use match_mos_apply to factor and offset the mosaic images.
;	F.V. 1999, option to use TOTAL or AVERAGE fluxes in SOURCE/SKY box.
;-

pro analyze_fluxset

   common math_mosaics, math_List, math_images, math_imscaled

	Nmath = N_struct( math_List )
	if (Nmath LT 1) then return
	old_386i = (!version.arch EQ '386i')

	if (Nmath LT 2) AND (old_386i) then begin
		print," need at least 2 mosaics ...
		return
	   endif

	border_images, math_List, INUM=indgen( Nmath )
	Magf = math_List(0).Magf
	input = ""
	format = "(A30,G10.3)"
	point = "SOURCE"
	BELL = string( 7b )
	LF = string( 10b )
	menuf = [ "Select method using:", "AVERAGE flux in box","TOTAL flux"]
	ftype = next_word( menuf( wmenu( menuf, INIT=1, TIT=0 ) ) )
	init = 3
	menub = [ strtrim( 2*indgen(27)+3, 2 ), "ALL pixels" ]
PICK:
	menu = [ "Select " + point + " Box Size", menub ]
	sel = wmenu( menu, INIT=1, TIT=0 )
	if (sel LE 0) then goto,PICK

	if (next_word( menu(sel) ) EQ "ALL") then begin
		boxsel = max( math_List.size_image(1:2) )
	  endif else boxsel = fix( menu(sel) )

	printw, ["select POINT to set "+point+" fluxes",$
		 "or RIGHT button to abort"], /ERASE

	imos = select_image( math_List, x,y, /ALL, /CURSET )
	printw,[" "," "],/ERASE

	if (imos(0) LT 0) then return

	if (N_elements( imos ) LT 2) AND (old_386i) then begin
		print,LF+" point must be in at least 2 images"+BELL
		print," try again..."
		wait,1
		goto,PICK
	   endif

	boxhs = boxsel/2.
	box = boxhs*Magf
	box_draw2, POS=[x,y], RADIUS=box

	if (point EQ "SOURCE") then BEGIN
			imosave = imos
			xsrc = x
			ysrc = y
			boxsrc = boxsel-1
			boxhsrc = boxhs
			point = "SKY"
			init=7
			print," select POINT for SKY fluxes"
			goto,PICK
	  endif else if (point EQ "SKY") then BEGIN
			w = where( imos NE imosave, Ndiff )
			if (Ndiff GT 0) then begin
				print,LF+" SKY box must be in same mosaics"+BELL
				box_erase2
				goto,PICK
			   endif
			xsky = x
			ysky = y
			boxsky = boxsel-1
			boxhsky = boxhs
	   endif else begin
		message,"error: point="+point,/INFO
		return
	    endelse

	if N_elements( imos ) GT 1 then $
		imos = imos( sort( math_List(imos).Level ) )
	names = strmid( math_List(imos).name, 0,30 )

	mf = float( Magf )
	xsrc = xsrc/mf - math_List(imos).Locx
	ysrc = ysrc/mf - math_List(imos).Locy
	xsky = xsky/mf - math_List(imos).Locx
	ysky = ysky/mf - math_List(imos).Locy

	Npick = N_elements( imos )
	srcvals = fltarr( Npick )
	skyvals = fltarr( Npick )

; Note:	array <math_sizes> points to the Locations
;	of the selected mosaic image data in <math_images>,
;	and must be used sequentially,
;	since <imos> subscript mapping has already been applied.

	math_sizes = math_List(imos).size_image		;Locations of data!
	xsrcmin = round_off( xsrc - boxhsrc ) > 0
	xsrcmax = ( xsrcmin + boxsrc ) < (math_sizes(1,*)-1)
	ysrcmin = round_off( ysrc - boxhsrc ) > 0
	ysrcmax = ( ysrcmin + boxsrc ) < (math_sizes(2,*)-1)
	xskymin = round_off( xsky - boxhsky ) > 0
	xskymax = ( xskymin + boxsky ) < (math_sizes(1,*)-1)
	yskymin = round_off( ysky - boxhsky ) > 0
	yskymax = ( yskymin + boxsky ) < (math_sizes(2,*)-1)

	for im=0,Npick-1 do begin
		mosaic = image_extract( im, math_images, math_sizes )
		msrc =  mosaic( xsrcmin(im):xsrcmax(im),$
				ysrcmin(im):ysrcmax(im)  )
		msky =  mosaic( xskymin(im):xskymax(im),$
				yskymin(im):yskymax(im)  )
		if ftype eq "TOTAL" then begin
			srcvals(im) = total( msrc )
			skyvals(im) = total( msky )
		 endif else begin
			srcvals(im) = total( msrc )/N_elements( msrc )
			skyvals(im) = total( msky )/N_elements( msky )
		  endelse
	  endfor

	if !DEBUG then stop
	point = "SOURCE"
	boxhs = boxhsrc
	print, ftype + " fluxes in ",point," box, size =",byte(2*boxhs+1)
	for i=0,Npick-1 do print,names(i),srcvals(i),FORM=format

	read," enter new SOURCE fluxes in same order: ",input
	srcnew = float( get_words( input ) )
	nval = N_elements( srcnew )

	if (input EQ "") OR (nval LT 1) then begin
		box_erase2
		print," retry ..."
		goto,PICK
	   endif

	while (nval LT Npick) do begin
		read,input
		srcnew = [ srcnew, float( get_words( input ) )]
		nval = N_elements( srcnew )
	  endwhile

	point = "SKY"
	boxhs = boxhsky
	print, ftype + " fluxes in ",point," box, size =",byte(2*boxhs+1)
	for i=0,Npick-1 do print,names(i),skyvals(i),FORM=format

	read," enter new SKY fluxes in same order: ",input
	skynew = float( get_words( input ) )
	nval = N_elements( skynew )

	if (input EQ "") OR (nval LT 1) then begin
		box_erase2
		print," retry ..."
		goto,PICK
	   endif

	while (nval LT Npick) do begin
		read,input
		skynew = [ skynew, float( get_words( input ) )]
		nval = N_elements( skynew )
	  endwhile

	Factor = (srcnew-skynew)/(srcvals-skyvals)
	Offset = skynew - Factor*skyvals
	print," Factors :",Factor
	print," Offsets :",Offset
	if !DEBUG then stop
	box_erase2

	match_mos_apply, math_List, math_images, math_imscaled,  $
						imos, Factor, offset, /APPLY
end