pro ruffazm_sbsp, stks, limblu, ceni, limred, dvdmin, azmout

;	subroutine:  ruffazam
;
;	purpose:  Find prefered azimuth angle by least-squares method of
;		  Auer Heasley House (1977,Solar Physics, 55, 47).
;
;	author:  paul@ncar  (revised 3/97)
;
;	input:
;	      stks	- observed Stokes spectra for one slit position
;
;		limblu  - lower limit of band pass, pixels (integer).
;		ceni	- approximate spectra line center, pixels (real).
;		limred  - upper limit of band pass, pixels (integer).
;		dvdmin  - minimum doppler width, pixels (real).
;	output:
;		azmout	- prefered azimuth, radians.
;
;-------------------------------------------------------------------------------
;pi = 3.14159 26535 89793 23846

;  find dimension of the input arrays
nw = sizeof_sbsp(stks,1) & nsl = sizeof_sbsp(stks,2)
azmout = fltarr(nsl)

;  Start loop over slit positions
for isl = 0,nsl-1 do begin

;  select Stokes Q, U profiles
	obsq = stks(*,isl,1) & obsu = stks(*,isl,2)

;  find guess for line center based on minimum of Stokes I
	obsi = smooth(stks(*,isl,0),7)
	prfi = obsi(limblu:limred)
	imin = min(prfi,minp)
	jcen = minp+limblu


	


;				    !Temporary arrays.
;	parameter( itmp = 2048 )
;	dimension vec1(itmp), vec2(itmp)

;	if( itmp .lt. np ) then
;		write(6,'(''a.ruffazm.f: itmp dimension to small'')')
;		stop 'a.ruffazm.f: itmp dimension to small'
;	endif
;				    !Compute arrays used in least squares.
;	do 600 w = limblu,limred
;		vec1(w) = obsq(w)*obsu(w)
;		vec2(w) = obsq(w)**2-obsu(w)**2
;600	continue

	vec1 = obsq*obsu
	vec2 = obsq*obsq - obsu*obsu

	snum = total(vec1(limblu:limred))
	sden = total(vec2(limblu:limred))

;	snum = 0.
;	do 100 w = limblu,limred
;		snum = snum+vec1(w)
;100	continue
;
;	sden = 0.
;	do 110 w = limblu,limred
;		sden = sden+vec2(w)
;110	continue
;				    !Angle that gives minimum or maximum
;				    !sum of squared u.
	dirazm = !pi/8.
	if(sden ne 0.) then dirazm = .25*atan(2.*snum/sden)

;				    !Resolve pi/4 ambiguity.
;				    !If second derivative of sum u*u is
;				    !negative (for maximum) add pi/4.
	if(.5*sden*cos(4.*dirazm)+snum*sin(4.*dirazm) lt 0.) then $
     	dirazm = dirazm+!pi/4.
;				    !Resolve pi/2 ambiquity.
;				    !Look for dip in middle of q-profile.
;				    !If second derivative of q-profile is
;				    !negative in line core add pi/2.

;				    !Rotated q-profile to azimuth.
	sin2azm = sin(2.*dirazm)
	cos2azm = cos(2.*dirazm)


;	do 200 w = 1,np
;		vec2(w) = obsq(w)*cos2azm+obsu(w)*sin2azm
;200	continue
	vec2 = obsq*cos2azm + obsu*sin2azm


;				    !First derivative of q-profile.
;	do 211 w = 2,np-1
;		vec1(w) = vec2(w+1)-vec2(w-1)
;211	continue
	for iw = 1,nw-2 do vec1(iw) = vec2(iw+1) - vec2(iw-1)

;				    !Second derivative of q-profile.
;	do 212 w = 3,np-2
;		vec2(w) = vec1(w+1)-vec1(w-1)
;212	continue
	for iw = 2,nw-3 do vec2(iw) = vec1(iw+1) - vec1(iw-1)



;				    !Sum of second derivative about center
;				    !to get trend.
	sum = 0.
;	jcen = nint(ceni)
	i2dvd = round(2.*dvdmin)
;	do 213 w = max(limblu,jcen-i2dvd,3),min(limred,jcen+i2dvd,np-2)
;		sum = sum+vec2(w)
;213	continue

	iw1 = max([limblu,jcen-i2dvd,2])
	iw2 = min([limred,jcen+i2dvd,nw-3])
	sum = total(vec2(iw1:iw2))
		

	if(sum lt 0.) then dirazm = dirazm+!pi/2.

;				    !Put azimuth in range -pi/2 to pi/2.
;	if(dirazm gt  pi) dirazm =  amod( dirazm,pi)
;	if(dirazm lt -pi) dirazm = -amod(-dirazm,pi)
	if(dirazm gt  !pi) then dirazm =  dirazm mod !pi
	if(dirazm lt -!pi) then dirazm =  -(-dirazm mod !pi)
	
	if(dirazm gt  !pi/2.) then dirazm = dirazm-!pi
	if(dirazm lt -!pi/2.) then dirazm = dirazm+!pi

	azmout(isl) = dirazm

;  end loop over slit length
endfor


	return
	end