pro pvpvcomp_sbsp,stks,nmult,liml,limc,lctr,qmask,vcross,ctrline,contin, $ pvpv,plpl,pqpq,pupu,pipi,plt,ptt,azmth ; routine to compute the signed, integrated Stokes V signal (pvpv) ; from a Stokes V spectrum, and the approximate location of the Stokes ; V zero-crossing wavelength, and the continuum intensity ; definition: pvpv = [sign of blue lobe of Stokes V profile]* ; [sum over lines{integral over spectrum(|V|)}] ; INPUTS: ; stks = spectral Stokes profiles ; nmult = number of lines to be considered in the sum to get ; pvpv ; liml = array [2,nmult] containing limits of integration ; for each line, in pixels ; limc = array [2] for continuum integration to get Ic ; lctr = integer array [nmult] for guess for center position of ; Stokes V ; qmask = mask from QS data for Q after rotating to preferred frame ; OUTPUTS: ; vcross = array [nmult] of Stokes V zero-crossing wavelength ; (pixels) ; ctrline = array [nmult] intensity minimum pixel position ; contin = continuum intensity Ic ; pvpv = ratio of integrated, signed Stokes V relative to ; continuum intensity ; plpl = ratio of wavelength-integrated preferred-frame ; Stokes Q relative to continuum intensity ; pqpq = ratio of wavelength-integrated Stokes Q to continuum ; pupu = ratio of wavelength-integrated Stokes U to continuum ; pipi = ratio of wavelength-integrated Stokes I to continuum ; plt = ratio of wavelength-integrated linear polarization ; sqrt(Q^2+U^2) to continuum ; ptt = ratio of wavelength-integrated total polarization ; sqrt(Q^2+U^2+V^2) to continuum ; azmth = field azimuth from ruff application to profiles ; History: ; 23-jan-2008- Bruce Lites - new definition of pvpv, avoids ; offsets near noise floor ; Compute Ic contin = mean(stks(limc(0):limc(1),0),/double) vcross = fltarr(nmult) nw = sizeof_sbsp(stks,1) vec1 = fltarr(nw) & vec2 = fltarr(nw) ww = findgen(nw) ; wavelength index ctrline = fltarr(nmult) ; line center intensity minimum pixel ; compute the (net linear polarization)^2 at each wavelength veclin = (stks(*,1)^2 + stks(*,2)^2) ; compute the (total polarization)^2 at each wavelength vectot = (stks(*,1)^2 + stks(*,2)^2 + stks(*,3)^2) ; initialize some wavelength integration quantities plpl = 0. pqpq = 0. pupu = 0. pipi = 0. plt = 0. ptt = 0. sumlin = 0. sumtot = 0. pvpv = 0. rnum = 0. bnum = 0. sgnb = 0. sgnr = 0. sumi = 0. sumq = 0. sumu = 0. ; begin loop over multiplet lines for im = 0,nmult-1 do begin iw1 = liml(0,im) & iw2 = liml(1,im) ; Locate the Stokes V zero-crossing wavelength closest to the line ; center guess ; first smooth the V-profile to minimize effects of noise smthv = smooth(stks(*,3),5) ; find sign of the smoothed Stokes V. Set zero values to positive temp = smthv whr = where(temp eq 0.,kount) if kount gt 0 then temp(whr) = 1. sgn = fix((temp)/abs(temp)) ; determine wavelength pixels for zero crossing irev = 0 ixp = lctr(im)+1 ixm = lctr(im)-1 while irev eq 0 do begin ; look in red direction, check for bounds of search if sgn(ixp) ne sgn(ixp-1) and ixp le iw2 then begin ;print,'found reversal for ixp ',ixp zc1 = ixp-1 & zc2 = ixp goto,kilroy endif else begin ixp = ixp+1 endelse ; look in blue direction, check for bounds of search if sgn(ixm) ne sgn(ixm+1) and ixm ge iw1 then begin ;print,'found reversal for ixm ',ixm zc1 = ixm & zc2 = ixm+1 goto,kilroy endif else begin ixm = ixm - 1 endelse ; check to see that pixels fall in the range for this multiplet, else ; just use guess for line center if ixm lt iw1 and ixp gt iw2 then begin irev = 1 vcross(im) = lctr(im) ; print,'did not find zero crossing for line',im goto,jump2 endif endwhile ; determine fractional pixel of zero-crossing by linear interpolation kilroy: sv1 = smthv(zc1) & sv2 = smthv(zc2) ;print,'sv1,sv2: ',sv1,sv2 vcross(im) = float(zc1) + sv1/(sv1-sv2) jump2: ; locate the positions of minimum line intensity ; fit parabola to bottom of line ; first locate the pixel line minimum smthi = smooth(stks(*,0),5) prfmin = min(smthi(iw1:iw2),imdx) imdx = imdx + iw1 ; fitting parabola of three points is equivalent to inverting a 3x3 matrix ; to get the 2nd order polynomial coefficients yy = a + b x + c x^2 if (smthi(imdx-1) eq smthi(imdx)) and (smthi(imdx) eq smthi(imdx+1)) $ then begin ctrln = float(imdx) goto, jump3 endif parabofit,ww(imdx-1:imdx+1),smthi(imdx-1:imdx+1),ctrln jump3: ctrline(im) = ctrln ; pick the line center position closest to the guess vcenter = vcross(im) if abs(ctrline(im)-lctr(im)) lt abs(vcross(im)-lctr(im)) then $ vcenter = ctrline(im) ; print,lctr(im),vcross(im),ctrline(im),vcenter ; CALCULATE pvpv ; initialize weight vectors for this multiplet vec1(*) = 0. & vec2(*) = 0. ; get weights for blue (vec1) and red (vec2) portions of integration ; range ;ctrof = fltarr(nw) ;ctrof(*) = 0. ;ll = iw1 ;lu = iw2 ;ctrof(ll:lu) = ww(ll:lu)-vcross(im) ; ;vvec1 = .5-ctrof ;stop ;whr = where(vvec1 ge 1.) ;vvec1(whr) = 1. ;stop ;whr = where(vvec1 le 0.) ;vvec1(whr) = 0. ; ; ;vvec2 = .5+ctrof ;whr = where(vvec2 ge 1.) ;vvec2(whr) = 1. ;whr = where(vvec2 le 0.) ;vvec2(whr) = 0. for iw = iw1,iw2 do begin ctroff = iw-vcenter vec1(iw) = max([0.,min([1.,.5-ctroff])]) vec2(iw) = max([0.,min([1.,.5+ctroff])]) endfor ;stop sumb = total(vec1*stks(*,3)) bnum = bnum + total(vec1) sumr = total(vec2*stks(*,3)) rnum = rnum + total(vec2) ; old definition of pvpv ;pvpv = pvpv+abs(sumb)+abs(sumr) ;sgnb = sgnb+sumb ;sgnr = sgnr+sumr ; new definition of pvpv, avoids offsets near noise floor pvpv = pvpv + sumb - sumr ; average the squares of the signals, then take the sqrt ; of the spectrally averaged squares sumlin = sumlin + total(vec1*veclin + vec2*veclin) sumtot = sumtot + total(vec1*vectot + vec2*vectot) sumq = sumq + total(vec1*stks(*,1) + vec2*stks(*,1)) sumu = sumu + total(vec1*stks(*,2) + vec2*stks(*,2)) sumi = sumi + total(vec1*stks(*,0) + vec2*stks(*,0)) ; end loop over multiplets endfor pvpv = pvpv/((bnum+rnum)*contin) plt = sqrt(sumlin/(bnum+rnum))/contin ptt = sqrt(sumtot/(bnum+rnum))/contin pqpq = sumq/((bnum+rnum)*contin) pupu = sumu/((bnum+rnum)*contin) pipi = sumi/((bnum+rnum)*contin) ; sign adjustment no longer needed for new pvpv definition ; set sign of Stokes V ;if(sgnr gt sgnb) then pvpv = -pvpv ; find effective azimuth with ruffazm dvdmin = 0.95 ; recast the Stokes profiles in a form for ruffazm: 3-dimensionsal rufstks = fltarr(nw,1,4) for istk = 0,3 do rufstks(*,0,istk) = stks(*,istk) ; call ruffazm, use 6302.5 ceni = lctr(1) limblu = liml(0,1) limred = liml(1,1) ruffazm_sbsp,rufstks,limblu,ceni,limred,dvdmin,azmout ;print,'azmout: ',azmout azmth = azmout(0) ; plpl computation ; rotate to the preferred frame srot = stks s_azm = sin(2.*azmth) & c_azm = cos(2.*azmth) srot(*,1) = stks(*,1)*c_azm + stks(*,2)*s_azm srot(*,2) = -stks(*,1)*s_azm + stks(*,2)*c_azm ; find the masked Q, U intensities qsig = srot(*,1)*qmask & usig = srot(*,2)*qmask plpl = total(qsig)/contin return end