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