pro bapp_sbsp,polnetin,pvin,tbbb,tpolnet,lbbb,lpv,polnetsm,pvsim

;  routine to convert the output linear, circ. polarization measures from 
;  the stksimages.pro output, polnet, to apparent longitudinal,
;  transverse flux density

;  INPUTS:	polnetin= the measure of linear polarization of Stokes Q
;			  rotated into the preferred frame, and applied
;			  the mask function
;		pvin	= the measure of line-integrated circular polarization
;		tbbb	= set of input field strengths for ME Stokes profile
;			  computations, used in calibration
;		tpolnet = qmask linear polarization computation (see
;			  stksimages_sbsp.pro), for ME profiles
;		lbbb	= set of input field strengths for ME Stokes profile
;			  computations, used in calibration
;		lpv	= wavelength-integrated Stokes V signal, ME profiles

;  OUTPUT:	polnetsm= the equivalent transverse apparent flux density,
;			  Mx/cm-2
;		pvsim	= the equivalent longitudinal apparent flux density,
;			  Mx/cm-2

;  LINEAR POLARIZATION: TRANSVERSE FIELDS
;  Current version of lincalib90 works with the Stokes Q integrated signal
;  using a Quiet sun mask derived from deep-mode spectra (20070227).
;  Thus, no longer uses net linear polarization.
;  only go to 1500 G since curve turns over after 1700 G
nsyn = 20
lbbb = tbbb(0:nsyn-1)
lpolcal = tpolnet(0:nsyn-1)
;  determine the polarization signal for 100 Mx cm-2
whr100 = where(lbbb le 100.,countbbb)
lpol100 = lpolcal(countbbb)

;  for the range 0-100 Mx cm-2, use sqrt function fit
;  fit the observed curve with a square root function to get the
;  amplitude of the square-root relationship in Mx cm-2, sqrtamp
sqrtamp = 3000. ; initial guess
weights = fltarr(countbbb+1) & weights(*) = 1.
result=curvefit(lpolcal(0:countbbb),lbbb(0:countbbb),weights,sqrtamp, $
        function_name='sqrtfn_sbsp',/noderivative)
;  above 100 Mx cm-2, use spline fit


;  get the effective transverse field polnetsm
polnetin = abs(polnetin)
;  gentle median filter the net linear polarization to remove few spurious points
;  perform the median filtering along the slit by transposing the polnetin
;  array, then transposing back.  Increase the threshold to 5-sigma
; medfilt_sbsp,polnetin,3,3,polnetsm
if (n_elements(polnetin) gt 7) then medfilt_sbsp,transpose(polnetin),3,5,polnetsm
polnetsm = transpose(polnetsm)
;  calculate the transverse field from the calibration.  Avoid negative values
;  by taking the absolute value of the wavelength-integrated Q-polarization
;  calculate the effective transverse fields
;  get the linear polarization calibration
plnet = abs(polnetsm)
;  convert linear polarization to Mx cm-2, use spline interpolation
;  of points
srtindx = sort(plnet)
linsort = plnet(srtindx)
splinelin = spline(lpolcal,lbbb,linsort)

;  for fields le 100 Mx cm-2, use sqrt function
whr100 = where(linsort le lpol100,count100)
if count100 gt 0 then splinelin(whr100) = sqrtamp*sqrt(linsort(whr100))

;  for polarizations gt maximum [lpolcal(nsyn-1)], set polnetsm
;  to 1500 Mx cm-2.  This analysis does not give meaningful results
;  above that polarization.
whrlmax = where(linsort gt lpolcal(nsyn-1),countlmax)
if countlmax gt 0 then splinelin(whrlmax) = lbbb(nsyn-1)

;  restore the Btapp image
linsort(srtindx) = splinelin
polnetsm(*,*) = linsort


;  CIRCULAR POLARIZATION: LONGITUDINAL FIELDS
;  now do the same for circular polarization for the 45-degree simulations 
pvv = abs(lpv)
;bb45 = lbbb/sqrt(2.)
bb45 = tbbb/sqrt(2.)
npoly = 3
;pvcoef=poly_fit(pvv,bb45,npoly,yfit=yfit)
pvsim = pvin & pvsim(*,*) = 0.
pvtmp = abs(pvin)

;  find maximum of this spline fit, and do not spline fit
;  points above that maximum -- will use a linear relationship
;  extrapolated from the final two points at highest field
;  strengths of the bb45-pvv relationship
maxpol = max(pvv)

;  slope and intercept for points with pv gt maxpol for 
;  extrapolation to high field strengths
npvv = n_elements(pvv)
slop = (bb45(npvv-1) - bb45(npvv-2))/(pvv(npvv-1) - pvv(npvv-2))
xinte = bb45(npvv-1) - slop*pvv(npvv-1)

;  determine the polarization signal for 100 Mx cm-2
whr100v = where(bb45 le 100.,countbbv)
vpol100 = pvv(countbbv)
calcoef = bb45(countbbv)/vpol100
;  convert integrated Stokes V to Mx cm-2, use spline interpolation
;  of points
srtindx = sort(pvtmp)
circsort = pvtmp(srtindx)
splinec = spline(pvv,bb45,circsort)

;  for fields le 100 Mx cm-2, use linear relation
whr100 = where(circsort le vpol100,count100)
if count100 gt 0 then splinec(whr100) = calcoef*circsort(whr100)

;  for pixels where observed pv > max(pvv) use linear extrapolation
whrmax = where(circsort gt maxpol,countmxx)
if countmxx gt 0 then splinec(whrmax) = xinte + slop*circsort(whrmax)

;  recover the Blapp image
circsort(srtindx) = splinec
pvsim(*,*) = circsort

;  recover the signs of the original integrated polarization
whrneg = where(pvin lt 0.)
pvsim(whrneg) = -pvsim(whrneg)


return
end