code: def _density_par_approx(self,dangle,tdisrupt,_return_array=False,
                              higherorder=False):
          
          # First construct the breakpoints for this dangle
          Oparb= (dangle-self._kick_interpdOpar_poly.x)/self._timpact
          # Find the lower limit of the integration in the pw-linear-kick approx.
          lowbindx,lowx= self.minOpar(dangle,tdisrupt,_return_raw=True)
          lowbindx= numpy.arange(len(Oparb)-1)[lowbindx]
          Oparb[lowbindx+1]= Oparb[lowbindx]-lowx
          # Now integrate between breakpoints
          out= (0.5/(1.+self._kick_interpdOpar_poly.c[-2]*self._timpact)\
                             *(special.erf(1./numpy.sqrt(2.*self._sortedSigOEig[2])\
                                               *(Oparb[:-1]-self._kick_interpdOpar_poly.c[-1]-self._meandO))\
                                   -special.erf(1./numpy.sqrt(2.*self._sortedSigOEig[2])*(numpy.roll(Oparb,-1)[:-1]-self._kick_interpdOpar_poly.c[-1]-self._meandO\
                                                                                               -self._kick_interpdOpar_poly.c[-2]*self._timpact*(Oparb-numpy.roll(Oparb,-1))[:-1]))))
          if _return_array:
              return out
          out= numpy.sum(out[:lowbindx+1])
          if higherorder:
              # Add higher-order contribution
              out+= self._density_par_approx_higherorder(Oparb,lowbindx)
          # Add integration to infinity
          out+= 0.5*(1.+special.erf((self._meandO-Oparb[0])\
                                    /numpy.sqrt(2.*self._sortedSigOEig[2])))
          return out
docs: Compute the density as a function of parallel angle using the 
          spline representation + approximations