def _determine_stream_spreadLB(self,simple=_USESIMPLE,
                                   ro=None,vo=None,
                                   R0=None,Zsun=None,vsun=None):
        """Determine the spread in the stream in observable coordinates"""
        if not hasattr(self,'_allErrCovs'):
            self._determine_stream_spread(simple=simple)
        if ro is None:
            ro= self._ro
        if vo is None:
            vo= self._vo
        if R0 is None:
            R0= self._R0
        if Zsun is None:
            Zsun= self._Zsun
        if vsun is None:
            vsun= self._vsun
        allErrCovsLB= numpy.empty_like(self._allErrCovs)
        obs= [R0,0.,Zsun]
        obs.extend(vsun)
        obskwargs= {}
        obskwargs['ro']= ro
        obskwargs['vo']= vo
        obskwargs['obs']= obs
        self._ErrCovsLBScale= [180.,90.,
                               self._progenitor.dist(**obskwargs),
                               numpy.fabs(self._progenitor.vlos(**obskwargs)),
                               numpy.sqrt(self._progenitor.pmll(**obskwargs)**2.
                                          +self._progenitor.pmbb(**obskwargs)**2.),
                               numpy.sqrt(self._progenitor.pmll(**obskwargs)**2.
                                          +self._progenitor.pmbb(**obskwargs)**2.)]
        allErrCovsEigvalLB= numpy.empty((len(self._thetasTrack),6))
        allErrCovsEigvecLB= numpy.empty_like(self._allErrCovs)
        eigDir= numpy.array([numpy.array([1.,0.,0.,0.,0.,0.]) for ii in range(6)])
        for ii in range(self._nTrackChunks):
            tjacXY= bovy_coords.galcenrect_to_XYZ_jac(*self._ObsTrackXY[ii])
            tjacLB= bovy_coords.lbd_to_XYZ_jac(*self._ObsTrackLB[ii],
                                               degree=True)
            tjacLB[:3,:]/= ro
            tjacLB[3:,:]/= vo
            for jj in range(6):
                tjacLB[:,jj]*= self._ErrCovsLBScale[jj]
            tjac= numpy.dot(numpy.linalg.inv(tjacLB),tjacXY)
            allErrCovsLB[ii]=\
                numpy.dot(tjac,numpy.dot(self._allErrCovsXY[ii],tjac.T))
            #Eigen decomposition for interpolation
            teig= numpy.linalg.eig(allErrCovsLB[ii])
            #Sort them to match them up later
            sortIndx= numpy.argsort(teig[0])
            allErrCovsEigvalLB[ii]= teig[0][sortIndx]
            #Make sure the eigenvectors point in the same direction
            for jj in range(6):
                if numpy.sum(eigDir[jj]*teig[1][:,sortIndx[jj]]) < 0.:
                    teig[1][:,sortIndx[jj]]*= -1.
                eigDir[jj]= teig[1][:,sortIndx[jj]]
            allErrCovsEigvecLB[ii]= teig[1][:,sortIndx]
        self._allErrCovsLBUnscaled= allErrCovsLB
        #Interpolate the allErrCovsLB covariance matrices along the interpolated track
        #Interpolate the eigenvalues
        interpAllErrCovsEigvalLB=\
            [interpolate.InterpolatedUnivariateSpline(self._thetasTrack,
                                                      allErrCovsEigvalLB[:,ii],
                                                      k=3) for ii in range(6)]
        #Now build the interpolated allErrCovsXY using slerp
        interpolatedAllErrCovsLB= numpy.empty((len(self._interpolatedThetasTrack),
                                               6,6))
        interpolatedEigval=\
            numpy.array([interpAllErrCovsEigvalLB[ii](self._interpolatedThetasTrack) for ii in range(6)]) #6,ninterp
        #Interpolate in chunks
        interpolatedEigvec= numpy.empty((len(self._interpolatedThetasTrack),
                                         6,6))
        for ii in range(self._nTrackChunks-1):
            slerpOmegas=\
                [numpy.arccos(numpy.sum(allErrCovsEigvecLB[ii,:,jj]*allErrCovsEigvecLB[ii+1,:,jj])) for jj in range(6)]
            slerpts= (self._interpolatedThetasTrack-self._thetasTrack[ii])/\
                (self._thetasTrack[ii+1]-self._thetasTrack[ii])
            slerpIndx= (slerpts >= 0.)*(slerpts <= 1.)
            for jj in range(6):
                for kk in range(6):
                    interpolatedEigvec[slerpIndx,kk,jj]=\
                    (numpy.sin((1-slerpts[slerpIndx])*slerpOmegas[jj])*allErrCovsEigvecLB[ii,kk,jj]
                     +numpy.sin(slerpts[slerpIndx]*slerpOmegas[jj])*allErrCovsEigvecLB[ii+1,kk,jj])/numpy.sin(slerpOmegas[jj])
        for ii in range(len(self._interpolatedThetasTrack)):
            interpolatedAllErrCovsLB[ii]=\
                numpy.dot(interpolatedEigvec[ii],
                          numpy.dot(numpy.diag(interpolatedEigval[:,ii]),
                                    interpolatedEigvec[ii].T))
        self._interpolatedAllErrCovsLBUnscaled= interpolatedAllErrCovsLB
        #Also calculate the (l,b,..) -> (X,Y,..) Jacobian at all of the interpolated and not interpolated points
        trackLogDetJacLB= numpy.empty_like(self._thetasTrack)
        interpolatedTrackLogDetJacLB=\
            numpy.empty_like(self._interpolatedThetasTrack)
        for ii in range(self._nTrackChunks):
            tjacLB= bovy_coords.lbd_to_XYZ_jac(*self._ObsTrackLB[ii],
                                                degree=True)
            trackLogDetJacLB[ii]= numpy.log(numpy.linalg.det(tjacLB))
        self._trackLogDetJacLB= trackLogDetJacLB
        for ii in range(len(self._interpolatedThetasTrack)):
            tjacLB=\
                bovy_coords.lbd_to_XYZ_jac(*self._interpolatedObsTrackLB[ii],
                                            degree=True)
            interpolatedTrackLogDetJacLB[ii]=\
                numpy.log(numpy.linalg.det(tjacLB))
        self._interpolatedTrackLogDetJacLB= interpolatedTrackLogDetJacLB
        return None
Determine the spread in the stream in observable coordinates