;WRITTEN BY YANG LIU, STANFORD UNIVERSITY 2009.
PRO tiltonly, mag, b_error=b_error, north = north, tiltangle
;
;
if not keyword_set(b_error) then b_error = 20.0 * 1.8
arcsec2km = 725.0 & pixel_size = 2.0 ; MDI magnetogram is 2.0 arcsec/pixel
coef = pixel_size * arcsec2km ; km/pixel for MDI full disk magnetograms
area = coef * coef * 1e10 ; cm^2

tmp = size ( mag ) & x_size = tmp (1) & y_size = tmp (2)
xx = fltarr(x_size, y_size) & yy = fltarr(x_size, y_size)
for i= 0, x_size -1 do begin
    for j= 0, y_size -1 do begin
        xx(i, j) = 1.0 * i & yy(i, j) = 1.0 * j
    endfor
endfor

a_pos = where(mag gt b_error, count)
a_neg = where(mag lt -b_error, count)
mag_pos = total(mag(a_pos)) 
x_pos = total(xx(a_pos) * mag(a_pos))/mag_pos & y_pos = total(yy(a_pos) * mag(a_pos))/mag_pos
mag_neg = total(mag(a_neg))
x_neg = total(xx(a_neg) * mag(a_neg))/mag_neg & y_neg = total(yy(a_neg) * mag(a_neg))/mag_neg

;loadct, 0
;tvim, bytscl(mag, max=500., min=-500.)
;colour
;plots, x_pos, y_pos, psym=1, color = 3, /data
;plots, x_neg, y_neg, psym=2, color = 2, /data

delta_x = ( x_neg - x_pos ) & delta_y = ( y_neg - y_pos )
distance = coef * sqrt ( delta_x ^2 + delta_y ^2 ) / 1000.0
angle = atan ( delta_y, delta_x ) / !DTOR
if keyword_set(north) then begin
   if angle ge 0.0 then angle = 180.0 - angle
   if angle lt 0.0 then angle = -180.0 - angle
endif

tiltangle=angle
print, angle
end