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 for iwhr=0,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