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