pro amb_shear,bx,by,bpx,bpy,hmax=hmax,bins=bins, median=median,$
    medsize=medsize, thd=thd 

;    new ambiguity solver
;    basic assumption : shear angle  follows a normal distribution
;
;    input : bx, by, bpx, bpy : bp is potential component
;
;    hmax= most probable shear (th-th_p)
;    if it is not given, this program estimates it.
;
;    bin= bin size for dertermining hmax(default= 5)
;    median :  smoothing keyword (default=7, for medsize)
;              inner product between bt and smoothed bt >0
;    
;   output : bx and by
;   
;   typical calling procedure
;     1) when hmax was already determined,
;        amb_shear, bx, by, bpx, bpy, hmax= xxx, median=median
;     2) when hmax was not determined
;          amb_shear,bx,by,bpx,bpy,median=median
;     	  
;
;
;   Y.-J. Moon      2003-April
;
;    present remark:  seems to be reasonable for localized sheared regions
;
;
 

    thp=atan(bpy, bpx)
    th=atan(by,bx)
    ra=!pi/180.
    thd=(th-thp)/ra

    s1=where(thd gt 180., count1)
    if count1 ge 1 then thd(s1)=thd(s1)-360
    s2=where(thd lt -180., count2)
    if count2 ge 1 then thd(s2)=thd(s2)+360

if not keyword_set(hmax) then begin
if not keyword_set(bins) then  bins=5
    h=histogram(thd,bin=bins,max=180,min=-180)
    hx=-180+findgen(n_elements(h))*bins
    sh=where(h eq max(h))
    hmax=hx(sh(0))+bins/2.
    print, 'hmax=', hmax
endif

    s3=where(thd lt hmax-90, count3)
    if count3 ge 1 then begin
      thd(s3)=thd(s3)+180.
      bx(s3)=-bx(s3)
      by(s3)=-by(s3)
    endif
    
    s4=where(thd gt hmax+90, count4)
    if count4 ge 1 then begin
     thd(s4)=thd(s4)-180.
     bx(s4)=-bx(s4)
     by(s4)=-by(s4)
    endif
if keyword_set(median) then begin
if not keyword_set(medsize) then medsize=7
   mbx=median(bx,medsize)
   mby=median(by,medsize)

   inner= bx*mbx+by*mby
   s5=where(inner lt 0) 
   bx(s5)=-bx(s5)
   by(s5)=-by(s5)
   thd(s5)=thd(s5)-180.
   s6=where(thd lt -180)
   thd(s6)=thd(s6)+360.
   
  
endif
 
end