function azam_smooth, azm, fld, psi, sdat, nchg, x0, y0, aax, aay, its=its_
;+
;
;	function:  azam_smooth
;
;	purpose:  return WHERE() magnetic field azimuth should be
;		  changed based on point to point continuity
;
;	author:  paul@ncar, 3/95	(minor mod's by rob@ncar)
;
;==============================================================================
;
;       Check number of parameters.
;
if n_params() eq 0 then begin
	print
	print, "usage:	ambig = azam_smooth( azm, fld, psi, sdat, nchg $"
	print, "		x0, y0, aax, aay )"
	print
	print, "	Return WHERE() magnetic field azimuth should be"
	print, "	changed based on point to point continuity."
	print
	print, "	Arguments"
	print, "		azm	- input 2D array of line of sight"
	print, "			  azimuth (degrees)"
	print, "		fld	- input 2D array of field strength"
	print, "			  (gauss)"
	print, "		psi	- input 2D array of line of sight"
	print, "			  inclination (degrees)"
	print, "		sdat	- input 2D array true for data"
	print, "			  locations"
	print, "		nchg	- size of returned array"
	print, "		x0 y0	- lower left position of sub array"
	print, "		aax aay	- dimensions of image"
	print, "			  in calling program"
	print
	print, "	Keywords"
	print, "		its	- number iterations (def: 20)."
	print
	return, -1
endif
;-
;							5/25/94
;
; There is a new option (smooth) available with the interactive mouse.
; To enable click on the lower right display button till it says (smooth).
; (smooth) is then run by moving the mouse with the right button pressed.
;
; (smooth) works on point to point continuity of the magnetic field.
; All azimuth changes with (smooth) improve continuity.  Lets say
; the mouse is working on (8x8) sub arrays.  With (smooth) only the 
; inner 6x6 are allowed to change.  By holding the mouse boundary fixed
; continuity outside the mouse boundary is not effected.
; Continuity of the inner 6x6 is force to the boundary.
; 
; The (wads) algorithm does not preserve continuity outside the
; mouse boundary.  For the bulk of work involving continuity
; you should use (wads).  (wads) can move the ambiguous boundaries
; boldly.
;
; For fine details use (smooth).  I do not recomend you use (smooth)
; till you are about done.  It will effect only ambiguous boundaries.
; I recomend (8x8) sub arrays when running (smooth).
;
;	 						5/25/94
; 
; To resolve 180 ambiguity the University of Hawaii has as one step what
; is called "acute angle resolution".   I tried this portion of their
; code on ASP data.  The input ASP data had already been disambiguated
; with 'azam'.  The output showed some improvements in point to point
; continuity.
; 
; Reading the UH code I find "acute angle resolution" is done with
; dot products of field vectors point to point.  I picked up on the dot
; product idea.
; 
; 'azam' has been revised to use field vector dot products when working
; on point to point continuity.  Options effected are (wads) and (smooth).
; Of these options (smooth) has a sequence most like the UH code.
; The codes are only similar; it would be a mistake to say we are
; using the UH method.
; 
; The dot product gives a good measure of continuity.  Only the field 
; transverse to the line of sight need to be considered; the line of
; sight field component is unambiguous.  Let us represent the transverse
; component of two adjacent field vectors as
; 
;	S = Ai+Bj
;	T = Di+Ej
; 
; The ambiguous vectors are
; 
;	U = -Ai-Bj = -S
;	V = -Di-Ej = -T
; 
; Note the relation
; 
;	S dot T = AD+BE = |S||T|cos(angle between S & T)
;		= -(S dot V)
;		= -(U dot T)
;		= U dot V
; 
; For continuity one would pick a dot product that is positive. 
; 
; The trouble is that a point has more than one neighbor.  The point in
; the middle of a 3x3 array has 8 neighbors.  Let's take the dot product
; of the center point with each of its 8 neighbors and sum them.
; If sum of the 8 dot products is negative, it is likely the negative
; (or ambiguous) vector would be better for the center point.
; This is how the (smooth) option in 'azam' works.
;-----------------------------------------------------------------------------

				    ;Get dimension.
if n_dims( azm, xdim, ydim ) ne 2 then  return, -1

				    ;Test for null conditions.
nchg = 0
if  total(sdat) eq 0          then return, -1
if  xdim lt 2  or  ydim lt 2  then return, -1

				    ;Azimuth in radians and transverse
				    ;field strength.
ttt = azm*(!pi/180.)
tot = fld*sin(psi*(!pi/180.))*sdat
				    ;Compute x y field components.
bx  = round(tot*cos(ttt))
by  = round(tot*sin(ttt))
				    ;Save initial x, y field components.
savbx = bx
savby = by
				    ;Do iterations.
if keyword_set(its_) then  its=its_-1  else  its=19
for iter=0,its do begin
				    ;Initialize dot product sums to zero.
	tot = lonarr(xdim,ydim)
				    ;Add left to right dot products.
	ttt = bx(0:xdim-2,*)*bx(1:xdim-1,*)+by(0:xdim-2,*)*by(1:xdim-1,*)
	tot(0:xdim-2,*) = tot(0:xdim-2,*)+ttt
	tot(1:xdim-1,*) = tot(1:xdim-1,*)+ttt

				    ;Add botton to top dot products.
	ttt = bx(*,0:ydim-2)*bx(*,1:ydim-1)+by(*,0:ydim-2)*by(*,1:ydim-1)
	tot(*,0:ydim-2) = tot(*,0:ydim-2)+ttt
	tot(*,1:ydim-1) = tot(*,1:ydim-1)+ttt

				    ;Add lower left to upper right products.
	ttt = bx(0:xdim-2,0:ydim-2)*bx(1:xdim-1,1:ydim-1) $
	     +by(0:xdim-2,0:ydim-2)*by(1:xdim-1,1:ydim-1)
	tot(0:xdim-2,0:ydim-2) = tot(0:xdim-2,0:ydim-2)+ttt
	tot(1:xdim-1,1:ydim-1) = tot(1:xdim-1,1:ydim-1)+ttt

				    ;Add upper left to lower right products.
	ttt = bx(0:xdim-2,1:ydim-1)*bx(1:xdim-1,0:ydim-2) $
	     +by(0:xdim-2,1:ydim-1)*by(1:xdim-1,0:ydim-2)
	tot(0:xdim-2,1:ydim-1) = tot(0:xdim-2,1:ydim-1)+ttt
	tot(1:xdim-1,0:ydim-2) = tot(1:xdim-1,0:ydim-2)+ttt

				    ;Find where dot product sums are negative.
				    ;Exit iterations if none found.
	whr = where( tot lt 0, nwhr )
	if nwhr eq 0 then  goto, break0

				    ;Loop where product sums were negative.
	ttt = tot
	iwhr = 0L
	for iwhr=0L,nwhr-1 do begin
				    ;Array index of negative product sum. 
		y = whr(iwhr)/xdim
		x = whr(iwhr)-y*xdim

				    ;Remove dot product if neighbor is more
				    ;negative.
		if	min( ttt(( 0 > (x-1)) : ((x+1) < (xdim-1)) $
	 		        ,( 0 > (y-1)) : ((y+1) < (ydim-1))) ) $
		lt	ttt(x,y) $
		then	tot(x,y) = 0
	end
				    ;Preserve boundaries.
	if x0 gt 0 then  tot(0,*) = 0
	if y0 gt 0 then  tot(*,0) = 0
	if x0+xdim lt aax then  tot(xdim-1,*) = 0
	if y0+ydim lt aay then  tot(*,ydim-1) = 0

				    ;Find where dot product sum is negative.
				    ;Exit iterations if none found.
	whr = where( tot lt 0, nwhr )
	if nwhr eq 0 then  goto, break0

				    ;Swap azimuth where dot product sum
				    ;is negative.
	bx(whr) = -bx(whr)
	by(whr) = -by(whr)
end
break0:
				    ;Set where there are changes.
return, where( bx ne savbx  or  by ne savby, nchg )

end