def _density_par_approx_higherorder(self,Oparb,lowbindx,
                                        _return_array=False,
                                        gaussxpolyInt=None):
        """Contribution from non-linear spline terms"""
        spline_order= self._kick_interpdOpar_raw._eval_args[2]
        if spline_order == 1: return 0.
        # Form all Gaussian-like integrals necessary
        ll= (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])\
            /numpy.sqrt(2.*self._sortedSigOEig[2])
        ul= (Oparb[:-1]-self._kick_interpdOpar_poly.c[-1]-self._meandO)\
            /numpy.sqrt(2.*self._sortedSigOEig[2])
        if gaussxpolyInt is None:
            gaussxpolyInt=\
                self._densMoments_approx_higherorder_gaussxpolyInts(\
                ll,ul,spline_order+1)
        # Now multiply in the coefficients for each order
        powers= numpy.tile(numpy.arange(spline_order+1)[::-1],
                           (len(ul),1)).T
        gaussxpolyInt*= -0.5*(-numpy.sqrt(2.))**(powers+1)\
            *self._sortedSigOEig[2]**(0.5*(powers-1))
        powers= numpy.tile(numpy.arange(spline_order+1)[::-1][:-2],
                           (len(ul),1)).T
        for jj in range(spline_order+1):
            gaussxpolyInt[-jj-1]*= numpy.sum(\
                self._kick_interpdOpar_poly.c[:-2]
                *self._timpact**powers
                /(1.+self._kick_interpdOpar_poly.c[-2]*self._timpact)
                **(powers+1)
                *special.binom(powers,jj)
                *(Oparb[:-1]-self._kick_interpdOpar_poly.c[-1]-self._meandO)
                **(powers-jj),axis=0)
        if _return_array:
            return numpy.sum(gaussxpolyInt,axis=0)
        else:
            return numpy.sum(gaussxpolyInt[:,:lowbindx+1])
def _density_par_approx_higherorder ( self , Oparb , lowbindx , _return_array = False , gaussxpolyInt = None ) : spline_order = self . _kick_interpdOpar_raw . _eval_args [ 2 ] if spline_order == 1 : return 0. ll = ( 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 ] ) / numpy . sqrt ( 2. * self . _sortedSigOEig [ 2 ] ) ul = ( Oparb [ : - 1 ] - self . _kick_interpdOpar_poly . c [ - 1 ] - self . _meandO ) / numpy . sqrt ( 2. * self . _sortedSigOEig [ 2 ] ) if gaussxpolyInt is None : gaussxpolyInt = self . _densMoments_approx_higherorder_gaussxpolyInts ( ll , ul , spline_order + 1 ) powers = numpy . tile ( numpy . arange ( spline_order + 1 ) [ : : - 1 ] , ( len ( ul ) , 1 ) ) . T gaussxpolyInt *= - 0.5 * ( - numpy . sqrt ( 2. ) ) ** ( powers + 1 ) * self . _sortedSigOEig [ 2 ] ** ( 0.5 * ( powers - 1 ) ) powers = numpy . tile ( numpy . arange ( spline_order + 1 ) [ : : - 1 ] [ : - 2 ] , ( len ( ul ) , 1 ) ) . T for jj in range ( spline_order + 1 ) : gaussxpolyInt [ - jj - 1 ] *= numpy . sum ( self . _kick_interpdOpar_poly . c [ : - 2 ] * self . _timpact ** powers / ( 1. + self . _kick_interpdOpar_poly . c [ - 2 ] * self . _timpact ) ** ( powers + 1 ) * special . binom ( powers , jj ) * ( Oparb [ : - 1 ] - self . _kick_interpdOpar_poly . c [ - 1 ] - self . _meandO ) ** ( powers - jj ) , axis = 0 ) if _return_array : return numpy . sum ( gaussxpolyInt , axis = 0 ) else : return numpy . sum ( gaussxpolyInt [ : , : lowbindx + 1 ] )
Contribution from non-linear spline terms