def build_plot(self, usePCA):
        import pylab

        self.normaliser = MorphologyForRenderingOperator(self.morph, usePCA=usePCA)

        # Find the Point that is the furthest distance way from
        # the centre when the cell is centred and rotated:
        rotator = lambda s: self.normaliser(s.get_distal_npa3())

        rotated_section_dict = DictBuilderSectionVisitorHomo(morph=self.morph, functor=rotator)()

        # Add in the parents manually:
        dummy_scetion = self.morph._dummysection
        rotated_section_dict[dummy_scetion] = self.normaliser(dummy_scetion.get_distal_npa3())

        max_axis = max([numpy.linalg.norm(rotated_section_pt) for rotated_section_pt in rotated_section_dict.values()])
        plot_lims = (max_axis * -1.1, max_axis * 1.1)

        max_x = max([numpy.fabs(rotated_section_pt[0]) for rotated_section_pt in rotated_section_dict.values()])
        max_y = max([numpy.fabs(rotated_section_pt[1]) for rotated_section_pt in rotated_section_dict.values()])
        max_z = max([numpy.fabs(rotated_section_pt[2]) for rotated_section_pt in rotated_section_dict.values()])

        maxes = [max_x, max_y, max_z]

        # allMax = max(maxes)
        for i in self.plot_views:
            maxes[i] = maxes[i] + 0.2 * max([max_x, max_y, max_z])

        self.fig = pylab.figure(**self.fig_kwargs)  # figsize=(7, 7))
        self.fig.subplots_adjust(left=0.05, top=0.95, right=0.95, bottom=0.05, wspace=0.15, hspace=0.15)

        self.subplots = {}
        for i in self.plot_views:
            self.subplots[i] = self.build_draw_sub_plot(rotated_section_dict, self.fig, i, plot_lims)

    def setup_sampling(self, 
                       gradient, 
                       soln, 
                       linear_randomization,
                       quadratic_coef):

        self.accept_beta, self.total_beta = 0, 0

        random_direction = 2 * quadratic_coef * soln + linear_randomization 
        negative_subgrad = gradient + random_direction

        self.active_set = (soln != 0)
        self.initial_parameters = np.empty(1, self.dtype)
        self.initial_parameters['signs'] = np.sign(soln[self.active_set])
        abs_l1part = np.fabs(soln[self.active_set])
        l1norm_ = abs_l1part.sum()

        self.initial_parameters['simplex'] = (abs_l1part / l1norm_)[:-1]
        subgrad = -negative_subgrad[self.inactive_set]
        supnorm_ = np.fabs(negative_subgrad).max()

        if self.lagrange is not None: 
            self.initial_parameters['cube'] = subgrad / self.lagrange
            self.initial_parameters['scale'] = l1norm_
        else:
            if self._active_set.sum() != self.shape:
                self.initial_parameters['cube'] = subgrad / supnorm_
                self.initial_parameters['scale'] = supnorm_
        
        if self.lagrange is None:
            raise NotImplementedError("only lagrange form is implemented")

        return soln[self.active_set], subgrad

    def _calc_shift_ranges(self, x_shift, y_shift):
        '''Calculate shift indices for input and output arrays.
        '''
        LOGGER.debug("Calculating shift ranges.")
        # width of the portion to be moved
        width = self._img_shape[1] - int(np.fabs(x_shift))
        # height of the portion to be moved
        height = self._img_shape[0] - int(np.fabs(y_shift))

        # Calculate the corner indices of the area to be moved
        if x_shift < 0:
            n_x1, n_x2 = 0, width
            o_x1, o_x2 = -1*x_shift, -1*x_shift+width
        else:
            n_x1, n_x2 = x_shift, x_shift+width
            o_x1, o_x2 = 0, width
        if y_shift < 0:
            n_y1, n_y2 = 0, height
            o_y1, o_y2 = -1*y_shift, -1*y_shift+height
        else:
            n_y1, n_y2 = y_shift, y_shift+height
            o_y1, o_y2 = 0, height

        output_ranges = ((n_x1, n_x2), (n_y1, n_y2))
        input_ranges = ((o_x1, o_x2), (o_y1, o_y2))

        return (output_ranges[0], output_ranges[1],
                input_ranges[0], input_ranges[1])

def ransac(kernel, threshold):
    '''Robustly fit a model to data.

    >>> x = np.array([1., 2., 3.])
    >>> y = np.array([2., 4., 7.])
    >>> kernel = TestLinearKernel(x, y)
    >>> ransac(kernel, 0.1)
    (2.0, array([0, 1]), 0.10000000000000001)
    '''
    max_iterations = 1000
    best_error = float('inf')
    best_model = None
    best_inliers = []
    i = 0
    while i < max_iterations:
        try:
            samples = kernel.sampling()
        except AttributeError:
            samples = random.sample(range(kernel.num_samples()),
                                kernel.required_samples)
        models = kernel.fit(samples)
        for model in models:
            errors = kernel.evaluate(model)
            inliers = np.flatnonzero(np.fabs(errors) < threshold)
            error = np.fabs(errors).clip(0, threshold).sum()
            if len(inliers) and error < best_error:
                best_error = error
                best_model = model
                best_inliers = inliers
                max_iterations = min(max_iterations,
                    ransac_max_iterations(kernel, best_inliers, 0.01))
        i += 1
    return best_model, best_inliers, best_error

def trazo(ini, fin, ancho):
    '''Recibe la coordenada de inicio del trazo, el final y el ancho del
    trazo, devuelve un arreglo con las coordenadas de los puntos que lo
    conformman. Sólo funciona para trazos rectos.
    (recibe los puntos más izquierdos, o más  abajo del trazo)'''
    ancho = int(ancho)+1
    actual = ini
    if es_horizontal(ini, fin):
        paso = np.array([1, 0])
        ancho_1 = np.array([0, 1])
        rango = int(np.fabs(fin[0]-ini[0]))+1
    else:
        paso = np.array([0, 1])
        ancho_1 = np.array([1, 0])
        rango = int(np.fabs(fin[1]-ini[1]))

    trazo = np.zeros([(rango)*(ancho), 2])
    for a in range(ancho):
        for i in range(rango):
            trazo[i+(rango)*a] = np.array([actual])
            actual += paso
        actual = ini
        actual += (a+1)*ancho_1
    
    return trazo

def test_estimateDelta():

    #_____initialize some KKSPot_____
    Delta = 1.0
    pot = KuzminKutuzovStaeckelPotential(ac=20.,Delta=Delta,normalize=True)

    #_____initialize an orbit (twice)_____
    vxvv = [1.,0.1,1.1,0.01,0.1]
    o= Orbit(vxvv=vxvv)

    #_____integrate the orbit with C_____
    ts= numpy.linspace(0,101,100)
    o.integrate(ts,pot,method='leapfrog_c')


    #____estimate Focal length Delta_____
    #for each time step individually:
    deltas_estimate = numpy.zeros(len(ts))
    for ii in range(len(ts)):
        deltas_estimate[ii] = estimateDeltaStaeckel(pot,o.R(ts[ii]),o.z(ts[ii]))

    assert numpy.all(numpy.fabs(deltas_estimate - Delta) < 10.**-8), \
            'Focal length Delta estimated along the orbit is not constant.'

    #for all time steps together:
    delta_estimate = estimateDeltaStaeckel(pot,o.R(ts),o.z(ts))
    
    assert numpy.fabs(delta_estimate - Delta) < 10.**-8, \
            'Focal length Delta estimated from the orbit is not the same as the input focal length.'

    return None

    def __init__(self, num):
        global h, edge_times, edge_probs, state_density, states, Phdelta
        states, state_density, edge_probs, edge_times, Phdelta = [],[],[],[],[]
        
        h = (zmax-zmin)/float(num)
        self.tot_vert = num+1
        states = np.linspace(zmin, zmax, self.tot_vert)
        for s in states:
            c = process_var + h*np.fabs(self.drift(s))
            htime = h*h/c
            edge_times.append(htime)
        for i in range(self.tot_vert):
            s = states[i]
            c = process_var + h*np.fabs(self.drift(s))
            if (i != 0) and (i != self.tot_vert-1):
                edge_probs.append([(process_var/2 + h*np.fabs(self.drift(s)))/c, process_var/2/c])
            elif (i == self.tot_vert -1):
                edge_probs.append([1.0, 0.])
            elif (i == 0):
                edge_probs.append([0., 1.0])
        
        """ 
        # get filtering one_step transition probabilities using matrices
        self.delta = min(edge_times)*0.999
        P1 = np.zeros((self.tot_vert, self.tot_vert))
        P0 = np.zeros((self.tot_vert, self.tot_vert))
        for i in range(self.tot_vert):
            pt = edge_times[i]/(self.delta + edge_times[i])
            if( (i!=0) and (i!=self.tot_vert-1)):
                P1[i,i-1] = edge_probs[i][0]*pt
                P1[i,i+1] = edge_probs[i][1]*pt
                P0[i,i-1] = edge_probs[i][0]*(1-pt)
                P0[i,i+1] = edge_probs[i][1]*(1-pt)
            elif(i ==0):
                P1[i,i+1] = edge_probs[i][1]*pt
                P0[i,i+1] = edge_probs[i][1]*(1-pt)
            else:
                P1[i,i-1] = edge_probs[i][0]*pt
                P0[i,i-1] = edge_probs[i][0]*(1-pt)
        Phdelta = np.linalg.inv(eye(self.tot_vert) - P0)*P1
        """

        self.delta = min(edge_times)*0.999
        #self.delta = 0.0001
        #print 'min_htime: ', min(edge_times),' delta: ', self.delta
        if(min(edge_times) < self.delta):
            print "Add less nodes"
            sys.exit(0)
        # explicit method
        Phdelta = np.zeros((self.tot_vert, self.tot_vert))
        for i in range(self.tot_vert):
            ps = 1 - self.delta/edge_times[i]
            Phdelta[i,i] = ps 
            if( (i!=0) and (i!=self.tot_vert-1)):
                Phdelta[i,i+1] = edge_probs[i][1]*(1- ps)
                Phdelta[i,i-1] = edge_probs[i][0]*(1- ps)
            elif(i ==0):
                Phdelta[i,i+1] = edge_probs[i][1]*(1- ps)
            else:
                Phdelta[i,i-1] = edge_probs[i][0]*(1- ps)

def truncate_hist1( self, xmin, xmax ):
    buf   = get_buffer_hist1( self )
    sbuf  = get_err_buffer_hist1( self )
    edges, fixed = get_bin_edges_axis( self.GetXaxis(), type=True )

    e1 = numpy.fabs(edges[:-1]-xmin)<1.e-9
    e2 = numpy.fabs(edges[1:]-xmax)<1.e-9
    assert numpy.any( e1 ) and numpy.any( e2 ), 'Invalid new histogram limits'
    i1 = numpy.nonzero( e1 )[0][0]
    i2 = numpy.nonzero( e2 )[0][-1]+1

    if fixed:
        newhist = self.__class__( self.GetName(), self.GetTitle(), i2-i1, xmin, xmax )
    else:
        newhist = self.__class__( self.GetName(), self.GetTitle(), i2-i1, edges[i1:i2] )

    newbuf = get_buffer_hist1( newhist )
    if sbuf is None:
        newsbuf = None
    else:
        newhist.Sumw2()
        newsbuf = get_err_buffer_hist1( newhist )

    newbuf[:] = buf[i1:i2]
    if not sbuf is None:
        newsbuf[:] = sbuf[i1:i2]

    newhist.SetEntries( newhist.Integral() )

    return newhist

def geigen(Amat, Bmat, Cmat):
    """
    generalized eigenvalue problem of the form

    max tr L'AM / sqrt(tr L'BL tr M'CM) w.r.t. L and M

    :param Amat numpy ndarray of shape (M,N)
    :param Bmat numpy ndarray of shape (M,N)
    :param Bmat numpy ndarray of shape (M,N)

    :rtype: numpy ndarray
    :return values: eigenvalues
    :return Lmat: left eigenvectors
    :return Mmat: right eigenvectors

    """
    if Bmat.shape[0] != Bmat.shape[1]:
        print("BMAT is not square.\n")
        sys.exit(1)

    if Cmat.shape[0] != Cmat.shape[1]:
        print("CMAT is not square.\n")
        sys.exit(1)

    p = Bmat.shape[0]
    q = Cmat.shape[0]

    s = min(p, q)
    tmp = fabs(Bmat - Bmat.transpose())
    tmp1 = fabs(Bmat)
    if tmp.max() / tmp1.max() > 1e-10:
        print("BMAT not symmetric..\n")
        sys.exit(1)

    tmp = fabs(Cmat - Cmat.transpose())
    tmp1 = fabs(Cmat)
    if tmp.max() / tmp1.max() > 1e-10:
        print("CMAT not symmetric..\n")
        sys.exit(1)

    Bmat = (Bmat + Bmat.transpose()) / 2.
    Cmat = (Cmat + Cmat.transpose()) / 2.
    Bfac = cholesky(Bmat)
    Cfac = cholesky(Cmat)
    Bfacinv = inv(Bfac)
    Bfacinvt = Bfacinv.transpose()
    Cfacinv = inv(Cfac)
    Dmat = Bfacinvt.dot(Amat).dot(Cfacinv)
    if p >= q:
        u, d, v = svd(Dmat)
        values = d
        Lmat = Bfacinv.dot(u)
        Mmat = Cfacinv.dot(v.transpose())
    else:
        u, d, v = svd(Dmat.transpose())
        values = d
        Lmat = Bfacinv.dot(u)
        Mmat = Cfacinv.dot(v.transpose())

    return values, Lmat, Mmat

def replicate_line_with_threshold(g_expr,  g_expr_c,  M,  K,  threshold):
    '''Compare whether two lines are the same lines.'''
    coeffi_1 = np.zeros((M,  K))
    coeffi_2 = np.zeros((M,  K))
    coeff1_constant = g_expr[-1]
    coeff2_constant = g_expr_c[-1]
    #calculate the sum of all coefficients
    sum_coeff1 = coeff1_constant
    sum_coeff2 = coeff2_constant
    for i in xrange(M):
        for j in xrange(K):
            coeffi_1[i][j] = g_expr[i*K + j]
            coeffi_2[i][j] = g_expr_c[i*K + j]
            sum_coeff1 += coeffi_1[i][j]
            sum_coeff2 += coeffi_2[i][j]
    
    #check constant
    if np.fabs(sum_coeff1 * coeff2_constant - sum_coeff2 * coeff1_constant) > threshold:
        return 0
    #check all other coefficients
    for i in xrange(M):
        for j in xrange(K):
            if np.fabs(sum_coeff1 * coeffi_2[i][j] - sum_coeff2 * coeffi_1[i][j]) > threshold:
                #sum_coeff1 * coeffi_2[i][j] != sum_coeff2 * coeffi_1[i][j]
                return 0 
    #1 represents the two lines are the same equation
    return 1

def weighted_mean(_line):
    max_weight = 50
    
    # print _line.shape
    
    median_2d = bottleneck.nanmedian(_line, axis=1).reshape(_line.shape[0],1).repeat(_line.shape[1], axis=1)
    std = bottleneck.nanstd(_line, axis=1)
    std_2d = std.reshape(_line.shape[0],1).repeat(_line.shape[1], axis=1)
    
    weight_2d = numpy.fabs(std_2d / (_line - median_2d))
#    weight_2d[weight_2d > max_weight] = max_weight
    weight_2d[numpy.isinf(weight_2d)] = max_weight
    
    for i in range(3):
        avg = bottleneck.nansum(_line*weight_2d, axis=1)/bottleneck.nansum(weight_2d, axis=1)
        avg_2d = avg.reshape(_line.shape[0],1).repeat(_line.shape[1], axis=1)
        
        std = numpy.sqrt(bottleneck.nansum(((_line - avg_2d)**2 * weight_2d), axis=1)/bottleneck.nansum(weight_2d, axis=1))
        std_2d = std.reshape(_line.shape[0],1).repeat(_line.shape[1], axis=1)
        
        weight_2d = numpy.fabs(std_2d / (_line - avg_2d))
        #weight_2d[weight_2d > max_weight] = max_weight
        weight_2d[numpy.isinf(weight_2d)] = max_weight
    
    return bottleneck.nansum(_line*weight_2d, axis=1)/bottleneck.nansum(weight_2d, axis=1)

 def visiting(self, x, step, temperature):
     dim = x.size
     if step < dim:
         # Changing all coordinates with a new visting value
         visits = np.array([self.visit_fn(
             temperature) for _ in range(dim)])
         upper_sample = self.rs.random_sample()
         lower_sample = self.rs.random_sample()
         visits[visits > self.tail_limit] = self.tail_limit * upper_sample
         visits[visits < -self.tail_limit] = -self.tail_limit * lower_sample
         x_visit = visits + x
         a = x_visit - self.lower
         b = np.fmod(a, self.b_range) + self.b_range
         x_visit = np.fmod(b, self.b_range) + self.lower
         x_visit[np.fabs(
             x_visit - self.lower) < self.min_visit_bound] += 1.e-10
     else:
         # Changing only one coordinate at a time based on Markov chain step
         x_visit = np.copy(x)
         visit = self.visit_fn(temperature)
         if visit > self.tail_limit:
             visit = self.tail_limit * self.rs.random_sample()
         elif visit < -self.tail_limit:
             visit = -self.tail_limit * self.rs.random_sample()
         index = step - dim
         x_visit[index] = visit + x[index]
         a = x_visit[index] - self.lower[index]
         b = np.fmod(a, self.b_range[index]) + self.b_range[index]
         x_visit[index] = np.fmod(b, self.b_range[
             index]) + self.lower[index]
         if np.fabs(x_visit[index] - self.lower[
                 index]) < self.min_visit_bound:
             x_visit[index] += self.min_visit_bound
     return x_visit

    def test_get_distance(self):
        """
        Test GridMol.get_distance.
        """
        self.mol.add_atom((1, 2, 1), 1.6)
        self.mol.add_atom((1, 1, 1), 1.6)
        distances = self.mol.get_distance()

        # confirm that negative values are inside atoms
        # this tests distance correspondence with occupancy
        mask = self.mol.get_occupancy()
        assert np.all(distances[mask] <= 0)
        assert np.all(distances[~mask] > 0)

        # check for sane positive distances
        assert np.amax(distances) < max(self.mol.get_real_shape())

        # check that negative distances are not too large
        # min should be no larger than the largest atom radius plus the probe
        # radius (by definition)
        assert np.fabs(np.amin(distances)) <= (
            self.mol.atoms[0].radius + self.mol.probe_radius)

        # check that most distances are significantly less than max
        threshold = max(self.mol.get_real_shape()) / 2.
        assert np.count_nonzero(np.fabs(distances) < threshold) > (
            0.9 * distances.size)

def fresnelR(n1, n2, theta):
    temp1 = np.sqrt(1.0-((n1/n2)*np.sin(theta))**2)
    temp2 = np.cos(theta)
    R_s = np.fabs( (n1*temp2-n2*temp1)/(n1*temp2+n2*temp1) )**2
    R_p = np.fabs( (n1*temp1-n2*temp2)/(n1*temp1+n2*temp2) )**2
    R = (R_s + R_p)/2
    return R

    def remove_duplicate_candidates(self, verbosity=1):
        """Remove lower-significance 'duplicate' (i.e. same period)
            candidates from a list of candidates.  For the highest
            significance candidate, include a list of the DMs (and SNRs)
            of all the other detections.

            Inputs:
                verbosity: Verbosity level. (Default: 1)

            Ouputs:
                None
        """
        if verbosity >= 1:
            print "  Sorting the %d candidates by frequency..." % \
                        self.get_numcands()
        self.cands.sort(cmp_freq)
        if verbosity >= 1:
            print "  Searching for dupes..."
        ii = 0
        # Find any match
        while ii < self.get_numcands():
            jj = ii + 1
            if jj < self.get_numcands() and \
                        Num.fabs(self.cands[ii].r-self.cands[jj].r) < r_err:
                # Find others that match
                jj += 1
                while jj < self.get_numcands() and \
                        Num.fabs(self.cands[ii].r-self.cands[jj].r) < r_err:
                    jj += 1
                matches = self.cands[ii:jj]
                matches.sort(cmp_sigma)
                bestindex = self.cands.index(matches[0])
                #sigmas = [c.sigma for c in matches]
                #bestindex = Num.argmax(sigmas)+ii
                # flag the duplicates
                bestcand = self.cands[bestindex]
                # Add other matching cands as hit of highest-sigma cand
                for matchind in reversed(range(ii, jj)):
                    if matchind == bestindex:
                        # The current candidate is the highest-sigma cand
                        # Don't remove it
                        continue
                    match = self.cands[matchind]
                    bestcand.add_as_hit(match)
                    match.note = "This candidate is a duplicate of %s:%d" % \
                                (bestcand.filename, bestcand.candnum)
                    self.duplicate_cands.append(self.cands.pop(matchind))
                    if verbosity >= 2:
                        print "Removing %s:%d (index: %d)" % \
                                (match.filename, match.candnum, matchind)
                        print "    %s" % match.note
                # If the best candidate isn't at the same freq
                # as ii, then it's possible even more hits should
                # be added. So we don't increment the index
                # (note that the best cand has moved into position ii).
            else:
                ii += 1 # No candidates to be added as hits, move on
        if verbosity >= 1:
            print "Found %d candidates.\n" % self.get_numcands()
        self.cands.sort(cmp_sigma)

    def weights(self, z):
        """
        Hampel weighting function for the IRLS algorithm

        The psi function scaled by z

        Parameters
        ----------
        z : array-like
            1d array

        Returns
        -------
        weights : array
            weights(z) = 1                            for \|z\| <= a

            weights(z) = a/\|z\|                        for a < \|z\| <= b

            weights(z) = a*(c - \|z\|)/(\|z\|*(c-b))      for b < \|z\| <= c

            weights(z) = 0                            for \|z\| > c

        """
        z = np.asarray(z)
        a = self.a; b = self.b; c = self.c
        t1, t2, t3 = self._subset(z)
        v = (t1 +
            t2 * a/np.fabs(z) +
            t3 * a*(c-np.fabs(z))/(np.fabs(z)*(c-b)))
        v[np.where(np.isnan(v))]=1. # for some reason 0 returns a nan?
        return v

    def bartlettTest(self):
        M = -(self.n-1/2*(self.p+self.q+3))

        for i in range(len(self.s)):
            #有意水準を求める
            alf = self.sigAlf
            sig = sp.special.chdtri((self.p-i)*(self.q-i), alf)
            test = 1

            #ウィルクスのラムダ
            for j in range(len(self.s)-i):
                test = test*(1-self.s[len(self.s)-j-1])
            chi = M*math.log(test)

            #帰無仮説棄却/採用
            if chi > sig:
                print  "test["+str(i)+"]:"+str(chi) +" > sig("+str(alf)+"):"+str(sig)
                run = np.fabs(self.A[:,i:i+1])
                rvn = np.fabs(self.B[:,i:i+1])
                arg_u = np.argmax(run)
                arg_v = np.argmax(rvn)
                self.eigArray.append(str(i))
                val = [arg_u, arg_v]
                self.eigValArray.append(val)
                print "eigen:"+str(np.sqrt(self.s[i]))
                print "ru-max arg:"+str(arg_u)
                print "rv-max arg:"+str(arg_v)
            else:
                break

def plot(filename, x=0, y=1, x_abs=False, y_abs=False,
         xscale='linear', yscale='linear', xlabel='x', ylabel='y',
         scatter=False, size=12, title='', xmin=-np.inf, xmax=np.inf,
         ymin=-np.inf, ymax=np.inf):

    filtercols = (x, y)
    mins = (xmin, ymin)
    maxs = (xmax, ymax)
    data_filtered = __filter_outranged_x(filename, filtercols, mins, maxs)

    xs, ys = np.loadtxt(data_filtered, usecols=(x, y), unpack=True)

    if x_abs:
        xs = np.fabs(xs)
    if y_abs:
        ys = np.fabs(ys)

    f, ax = plt.subplots(1, 1)
    ax.set_xscale(xscale)
    ax.set_yscale(yscale)
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)

    if title:
        plt.title(title)

    if scatter:
        ax.scatter(xs, ys, s=size)
    else:
        ax.plot(xs, ys, linewidth=size)

    plt.axis('tight')
    plt.show()
    plt.close()

    def get_weights_cs2ll(self, dst, alpha, beta, panel, gids):
        cs_obj = self.cs_obj

        (a1,b1), (a2,b2), (a3,b3), (a4,b4) = \
                [cs_obj.alpha_betas[gid] for gid in gids]
        
        assert np.fabs(a1-a3)<1e-15
        assert np.fabs(a2-a4)<1e-15
        assert np.fabs(b1-b2)<1e-15
        assert np.fabs(b3-b4)<1e-15
        assert flge(a1,alpha,a2), "dst={}, a1={}, a2={}, alpha={}".format(dst,a1,a2,alpha)
        assert flge(b1,beta,b3), "dst={}, b1={}, b3={}, beta={}".format(dst,b1,b3,beta)

        panels = [cs_obj.gq_indices[gid,0] for gid in gids]
        for p in panels:
            if p != panel:
                print("(alpha,beta) ({},{})".foramt(alpha, beta))
                print("panel: {}, {}".format(panel, panels))
                print("dst: {}".format(dst))
                print("gids: {}".format(gids))
                sys.exit()

        # weights
        x, y = alpha, beta
        x1, x2 = a1, a2
        y1, y2 = b1, b3

        return self.get_bilinear_weights(dst, (x,y), (x1,y1), (x2,y2))

def test_orbitIntegrationC():

    #_____initialize some KKSPot_____
    Delta = 1.0
    pot = KuzminKutuzovStaeckelPotential(ac=20.,Delta=Delta,normalize=True)

    #_____initialize an orbit (twice)_____
    vxvv = [1.,0.1,1.1,0.,0.1]
    o_P= Orbit(vxvv=vxvv)
    o_C= Orbit(vxvv=vxvv)

    #_____integrate the orbit with python and C_____
    ts= numpy.linspace(0,100,101)
    o_P.integrate(ts,pot,method='leapfrog')  #python
    o_C.integrate(ts,pot,method='leapfrog_c')#C

    for ii in range(5):
        exp3= -1.7
        if   ii == 0: Python, CC, string, exp1, exp2 = o_P.R(ts) , o_C.R(ts) , 'R' , -5., -10.
        elif ii == 1: Python, CC, string, exp1, exp2 = o_P.z(ts) , o_C.z(ts) , 'z' , -3.25, -4.
        elif ii == 2: Python, CC, string, exp1, exp2 = o_P.vR(ts), o_C.vR(ts), 'vR', -3., -10.
        elif ii == 3: Python, CC, string, exp1, exp2, exp3 = o_P.vz(ts), o_C.vz(ts), 'vz', -3., -4., -1.3
        elif ii == 4: Python, CC, string, exp1, exp2 = o_P.vT(ts), o_C.vT(ts), 'vT', -5., -10.

        rel_diff = numpy.fabs((Python-CC)/CC) < 10.**exp1
        abs_diff = (numpy.fabs(Python-CC) < 10.**exp2) * (numpy.fabs(Python) < 10.**exp3)
        assert numpy.all(rel_diff+abs_diff), \
            'Orbit integration for '+string+' coordinate different in ' + \
            'C and Python implementation.'

    return None