def plot_graph_3d(G):
    from mayavi import mlab

    for (node0, node1) in G.edges():
        (x0, y0, z0) = G.node[node0]["xyz"]
        (x1, y1, z1) = G.node[node1]["xyz"]
        mlab.plot3d([x0, x1], [y0, y1], [z0, z1], color=(1, 0, 0), tube_radius=0.0025)

def plot_eclipse(a,b,c,color=(1,1,1)):
    zs=np.linspace(-c*0.95,c*0.95,10)
    for z0 in zs:
        theta=np.linspace(0,pi*2,num=30)
        r=(1-z0**2/c**2)/(np.cos(theta)**2/a**2+np.sin(theta)**2/b**2)
        r=r**0.5
        mlab.plot3d(r*np.cos(theta),r*np.sin(theta),np.ones_like(r)*z0,color=color,tube_radius=None)

    def test_stericValue_1(self):
        print("m0, m1, m2", self.cg1.steric_value(["m0", "m1", "m2"]))
        from_ = np.amin(self.cg1.coords._coordinates)
        to_ = np.amax(self.cg1.coords._coordinates)
        x, y, z = np.mgrid[from_:to_:4, from_:to_:4, from_:to_:4]
        from mayavi import mlab
        s = np.zeros_like(x)
        for i, j, k in np.ndindex(x.shape):
            s[i, j, k] = self.cg1.steric_value(
                np.array([x[i, j, k], y[i, j, k], z[i, j, k]]), "r**-3")
        #mlab.contour3d(x,y,z,s, contours= [0.5, 1, 2, 5], opacity=0.3)
        src = mlab.pipeline.scalar_field(x, y, z, s)
        mlab.pipeline.volume(src)
        #mlab.pipeline.iso_surface(src, contours=[0.1, ], opacity=0.3)
        #mlab.pipeline.iso_surface(src, contours=[0.5, ], opacity=0.7)
        #mlab.pipeline.iso_surface(src, contours=[1, ])

        colors = {"s": (0, 1, 0), "h": (0, 0, 1), "m": (1, 0, 0), "i": (
            1, 1, 0), "f": (0.5, 0.5, 0.5), "t": (0.5, 0.5, 0.5)}
        for d in self.cg1.defines:
            x = self.cg1.coords[d][0][0], self.cg1.coords[d][1][0]
            y = self.cg1.coords[d][0][1], self.cg1.coords[d][1][1]
            z = self.cg1.coords[d][0][2], self.cg1.coords[d][1][2]
            mlab.plot3d(x, y, z, tube_radius=2, color=colors[d[0]])
        mlab.show()
        assert False

def drawArrows(file1, file2, figure, bbox, index, descendant):
    """
    Draw an 'arrow' from the cell 'index' to 'descendant'
    'descendant' is assumed to be a 1D np.array with size 0, 1 or 2.
    """
    #load the center of mass position
    if descendant.size > 0 and descendant[0] != 0:
        com1 = file1["features"][str(index[0])]["com"][:]
        com2 = file2["features"][str(descendant[0])]["com"][:]
        
        #write the cell label as text
        mlab.text3d(com1[2]-bbox[2]+1,com1[1]-bbox[1]+1,com1[0]-bbox[0]+1, str(index), color=(1,1,1), figure=figure)
        
        #plot a point where the current cell is
        mlab.points3d([com1[2]-bbox[2]+1],[com1[1]-bbox[1]+1],[com1[0]-bbox[0]+1],color=(0,0,1), figure=figure)
                
        #plot a line to the descendant's center
        mlab.plot3d([com1[2]-bbox[2]+1,com2[2]-bbox[2]+1],
                    [com1[1]-bbox[1]+1,com2[1]-bbox[1]+1],
                    [com1[0]-bbox[0]+1,com2[0]-bbox[0]+1],
                    tube_radius=0.2, color=(1,0,0), figure=figure)
        
        #plot a second line, if there is a split
        if descendant.size == 2:
            com3 = file2["features"][str(descendant[1])]["com"][:]
            mlab.plot3d([com1[2]-bbox[2]+1,com3[2]-bbox[2]+1],
                        [com1[1]-bbox[1]+1,com3[1]-bbox[1]+1],
                        [com1[0]-bbox[0]+1,com3[0]-bbox[0]+1],
                        tube_radius=0.2, color=(1,0,0), figure=figure)

def plot_matrix(connectmat_file, centers_file, threshold_pct=5, weight_edges=False,
                node_scale_factor=2, edge_radius=.5, resolution=8, name_scale_factor=1,
                names_file=None, node_indiv_colors=[], highlight_nodes=[], fliplr=False):
    """
    Given a connectivity matrix and a (x,y,z) centers file for each region, plot the 3D network
    """
    matrix = core.file_reader(connectmat_file)
    nodes = core.file_reader(centers_file)
    if names_file:
        names = core.file_reader(names_file,1)
    num_nodes = len(nodes)
    edge_thresh_pct = threshold_pct / 100.0
    matrix_flat = np.array(matrix).flatten()
    edge_thresh = np.sort(matrix_flat)[len(matrix_flat)-int(len(matrix_flat)*edge_thresh_pct)]
    
    matrix = core.file_reader(connectmat_file)
    ma = np.array(matrix)

    thresh = scipy.stats.scoreatpercentile(ma.ravel(),100-threshold_pct)
    ma_thresh = ma*(ma > thresh)
    
    if highlight_nodes:
        nr = ma.shape[0]
        subset_mat = np.zeros((nr, nr))
        for i in highlight_nodes:
            subset_mat[i,:] = 1
            subset_mat[:,i] = 1
        ma_thresh = ma_thresh * subset_mat
        
    if fliplr:
        new_nodes = []
        for node in nodes:
            new_nodes.append([45-node[0],node[1],node[2]]) # HACK
        nodes = new_nodes
    
    mlab.figure(bgcolor=(1, 1, 1), size=(400, 400))
    for count,(x,y,z) in enumerate(nodes):
        if node_indiv_colors:
            mlab.points3d(x,y,z, color=colors[node_indiv_colors[count]], scale_factor=node_scale_factor, resolution=resolution)
        else:
            mlab.points3d(x,y,z, color=(0,1,0), scale_factor=node_scale_factor, resolution=resolution)
        if names_file:
            width = .025*name_scale_factor*len(names[count])
            print width
            print names[count]
            mlab.text(x, y,names[count], z=z,width=.025*len(names[count]),color=(0,0,0))
    for i in range(num_nodes-1):    
        x0,y0,z0 = nodes[i]
        for j in range(i+1, num_nodes):
            #if matrix[i][j] > edge_thresh:
            if ma_thresh[i][j] > edge_thresh:
                x1,y1,z1 = nodes[j]
                if weight_edges:
                    mlab.plot3d([x0,x1], [y0,y1], [z0,z1],
                            tube_radius=matrix[i][j]/matrix_flat.max(),
                            color=(1,1,1))
                else:
                    mlab.plot3d([x0,x1], [y0,y1], [z0,z1],
                            tube_radius=edge_radius,
                            color=(1,1,1))

def save_poly(fname, lines):
  fname += "_poly.txt"
  f = open(fname, 'w')
  print ' ', fname
  for line in lines:
    points = calc_points(line)
    #spoints = bspline(points, n=points.shape[0], degree=20)
    ##m = len(points)
    m = len(points)/2
    if m<4: continue
    kx = 3
    ##if(m>3): kx = 3
    ##else: kx = m-1
    wx = np.ones(len(points))
    wx[0] = wx[-1] = 100
    tck,u=si.splprep(np.transpose(points),w=wx,k=kx,s=10)
    ##m /= 2
    ##if(m<4) : m=4
    spoints = np.transpose([si.splev(np.linspace(0,1,m),tck)])
    f.write("%2d " % m) 
    for spoint in spoints:
      for vert in spoint:
        f.write("%0.2f " % vert)
    f.write('\n')
    mylab.plot3d(points[:,0], points[:,1], points[:,2], color=(1,0,0))   
    mylab.plot3d(spoints[:,0], spoints[:,1], spoints[:,2], color=(0,1,0))
  f.close()
  mylab.show()
  return

def show_lidar_with_boxes(pc_velo, objects, calib,
                          img_fov=False, img_width=None, img_height=None): 
    ''' Show all LiDAR points.
        Draw 3d box in LiDAR point cloud (in velo coord system) '''
    if 'mlab' not in sys.modules: import mayavi.mlab as mlab
    from viz_util import draw_lidar_simple, draw_lidar, draw_gt_boxes3d

    print(('All point num: ', pc_velo.shape[0]))
    fig = mlab.figure(figure=None, bgcolor=(0,0,0),
        fgcolor=None, engine=None, size=(1000, 500))
    if img_fov:
        pc_velo = get_lidar_in_image_fov(pc_velo, calib, 0, 0,
            img_width, img_height)
        print(('FOV point num: ', pc_velo.shape[0]))
    draw_lidar(pc_velo, fig=fig)

    for obj in objects:
        if obj.type=='DontCare':continue
        # Draw 3d bounding box
        box3d_pts_2d, box3d_pts_3d = utils.compute_box_3d(obj, calib.P) 
        box3d_pts_3d_velo = calib.project_rect_to_velo(box3d_pts_3d)
        # Draw heading arrow
        ori3d_pts_2d, ori3d_pts_3d = utils.compute_orientation_3d(obj, calib.P)
        ori3d_pts_3d_velo = calib.project_rect_to_velo(ori3d_pts_3d)
        x1,y1,z1 = ori3d_pts_3d_velo[0,:]
        x2,y2,z2 = ori3d_pts_3d_velo[1,:]
        draw_gt_boxes3d([box3d_pts_3d_velo], fig=fig)
        mlab.plot3d([x1, x2], [y1, y2], [z1,z2], color=(0.5,0.5,0.5),
            tube_radius=None, line_width=1, figure=fig)
    mlab.show(1)

 def plotLines(self, X, color=None, size=None, scalars=None):
     
     if self.autoRemove: self.removeLines()
     
     if isinstance(X, list):
         Ndim = 3
     else:
         Ndim = len(X.shape)
     
     if color==None: color=(1,0,0)
     
     #~ data.scene.disable_render = True
     #~ view = data.scene.mlab.view()
     #~ roll = data.scene.mlab.roll()
     #~ move = data.scene.mlab.move()
     
     if Ndim==2:
         if scalars==None:
             self.lines.append(mlab.plot3d(X[:,0], X[:,1], X[:,2], color=color, tube_radius=size))
         else:
             self.lines.append(mlab.plot3d(X[:,0], X[:,1], X[:,2], scalars, tube_radius=size))
     elif Ndim==3:
         for x in X:
             if scalars==None:
                 self.lines.append(mlab.plot3d(x[:,0], x[:,1], x[:,2], color=color, tube_radius=size))
             else:
                 self.lines.append(mlab.plot3d(x[:,0], x[:,1], x[:,2], scalars, tube_radius=size))

 def draw_pontoons(self, fig, truss, R, freeboard):
     from mayavi import mlab
     nE = truss.shape[0]
     c = (0.5,0,0)
     for k in xrange(nE):
         if np.any(truss[k,2,:] > freeboard): continue
         mlab.plot3d(truss[k,0,:], truss[k,1,:], truss[k,2,:], color=c, tube_radius=R, figure=fig)

def show3d():
	# cube
	a = list(itertools.product(([0.5,1]),repeat=3))
	controls = [np.array(list(b)) for b in a]

	# random
	controls = np.random.uniform(0.0,1.0,(5,3))
	norms = np.apply_along_axis(np.linalg.norm, 1, controls)
	controls = np.transpose(controls.T/norms)

	cone = get_cone(controls,verbose=1)

	print "Controls:",controls
	print "Cone:",cone

	for c in controls:
		plan3d(c,np.array([0]*3))
	
	mlab.points3d(controls[:,0],controls[:,1],controls[:,2])

	for c in cone:
		l = np.array([c*0.0] + [c*1.0])
		x = l[:,0]
		y = l[:,1]
		z = l[:,2]
		mlab.plot3d(x,y,z,line_width=1.0)

	mlab.show()

def save_spline(f, end_label, line_label, size, pnts, slabels):

  llist = [] # get list of line indicies
  for t in np.where(pnts['label']==end_label)[0]: # find the end points
    ijk = pnts['ijk'][t].tolist()  # get end point indices
    slabels['hit'][tuple(ijk)] = 0 # clear end point hits
  next_point(slabels, ijk, end_label, line_label, llist) # walk from the last end point
 
  points = np.zeros((len(llist[0]), len(llist))) # indicies -> point coordinates
  for i in range(len(llist)): 
    #points[:,i] = i2xyz(llist[i], size)
    points[0,i] = 8 * 0.069 * (llist[i][0] + 0.5)
    points[1,i] = 8 * 0.069 * (llist[i][1] + 0.5)
    points[2,i] = 0.798 * (llist[i][2] + 0.5)

  s = 15.0 # smoothness parameter
  k = 3 # b-spline order
  nest = -1 # estimate of number of knots needed (-1 = maximal)
  tckp,u = splprep(points, s=s, k=k, nest=-1) # find b-spline knot points
  xnew,ynew,znew = splev(np.linspace(0, 1, 30), tckp) # interpolate smooth points

  f.write(str(len(xnew))+" ")
  for i in range(len(xnew)):
    f.write("%2.4f %2.4f %2.4f "% (xnew[i], ynew[i], znew[i]))
  f.write('\n')

  import mayavi.mlab as mylab
  mylab.plot3d(points[0], points[1], points[2], color=(1,0,0))
  mylab.plot3d(xnew, ynew, znew, color=(0,1,0))
  mylab.show()
  return

    def draw_column(self, fig, centerline, freeboard, h_section, r_nodes, spacingVec=None, ckIn=None):
        from mayavi import mlab
        npts = 20

        nsection = h_section.size
        z_nodes = np.flipud( freeboard - np.r_[0.0, np.cumsum(np.flipud(h_section))] )

        th = np.linspace(0, 2*np.pi, npts)
        for k in xrange(nsection):
            rk = np.linspace(r_nodes[k], r_nodes[k+1], npts)
            z  = np.linspace(z_nodes[k], z_nodes[k+1], npts)
            R, TH = np.meshgrid(rk, th)
            Z, _  = np.meshgrid(z, th)
            X = R*np.cos(TH) + centerline[0]
            Y = R*np.sin(TH) + centerline[1]

            # Draw parameters
            if ckIn is None:
                ck = (0.6,)*3 if np.mod(k,2) == 0 else (0.4,)*3
            else:
                ck = ckIn
            #ax.plot_surface(X, Y, Z, alpha=0.5, color=ck)
            mlab.mesh(X, Y, Z, opacity=0.7, color=ck, figure=fig)

            if spacingVec is None: continue
            
            z = z_nodes[k] + spacingVec[k]
            while z < z_nodes[k+1]:
                rk = np.interp(z, z_nodes[k:], r_nodes[k:])
                #ax.plot(rk*np.cos(th), rk*np.sin(th), z*np.ones(th.shape), 'r', lw=0.25)
                mlab.plot3d(rk*np.cos(th) + centerline[0], rk*np.sin(th) + centerline[1], z*np.ones(th.shape), color=(0.5,0,0), figure=fig)
                z += spacingVec[k]
                
                '''

def plotDensity(mToB,coeffs,n=100,plotVals="real",plotBnd=True):
    """Plot boundary density with Mayavi
    
       INPUT:
       mToB     - Mesh-to-Basis map
       coeffs   - Coefficients of boundary density
       n        - Discretisation points per element (default 100)
       plotVals - "real", "imag" or "abs" for real/imaginary part or absolute values
                  (default "real")
       plotBnd  - Also plot boundary of shape (default True)
    """

    from mayavi import mlab

    transforms={"real":numpy.real,"imag":numpy.imag,"abs":numpy.abs}
    
    x,f=evalDensity(mToB,coeffs,n)    
    f=transforms[plotVals](f)
                    
    xmin=min(x[0])
    xmax=max(x[0])
    ymin=min(x[1])
    ymax=max(x[1])
    zmin=max([-1,min(f)])
    zmax=min([1,max(f)])
    extent=[xmin,xmax,ymin,ymax,zmin,zmax]
    ranges=[xmin,xmax,ymin,ymax,min(f),max(f)]
    
    fig=mlab.figure(bgcolor=(1,1,1))
    mlab.plot3d(x[0],x[1],f,extent=extent,color=(0,0,0),tube_radius=None,figure=fig)
    if plotBnd:
        mlab.plot3d(x[0],x[1],numpy.zeros(x[0].shape),extent=[xmin,xmax,ymin,ymax,0,0],tube_radius=None,figure=fig)
    mlab.axes(extent=extent,ranges=ranges,line_width=3.0,figure=fig)
    mlab.show()

def plot_tree(rrt):
    from mayavi import mlab
    tree = rrt.tree
    leafs = [i for i in tree.nodes() if len(tree.successors(i)) == 0]
    accounted = set()
    paths = []

    for leaf in leafs:
        this_node = leaf
        this_path = []
        paths.append(this_path)
        while True:
            this_path.append(this_node)
            accounted.add(this_node)
            p = tree.predecessors(this_node)
            if len(p) == 0:
                break
            assert len(p) == 1
            this_node = p[0]
            if this_node in accounted:
                this_path.append(this_node) #repeat branching node
                break

    for p in paths:
        s = np.array([rrt.tree.node[i]['state'] for i in p])
        c = np.array([rrt.tree.node[i]['cost'] for i in p])
        mlab.plot3d(s[:,0],s[:,1],s[:,2]/50,c,tube_radius=0.002,colormap='Spectral')

def show_matches(pc1, pc2, samples1, samples2, matches, tube_radius=0.0003, N=50):
    axis = 0
    diff = 1
    bbx = utils.compute_bounding_box_std(pc1)
    pc1[:, axis] -= (bbx[axis][1] - bbx[axis][0])*1.5
    samples1[:, axis] -= (bbx[axis][1] - bbx[axis][0])*1.5
    pc2[:, axis] += (bbx[axis][1] - bbx[axis][0])*1.5
    samples2[:, axis] += (bbx[axis][1] - bbx[axis][0])*1.5
    show_pc(pc1, pc2)
    assert(samples1.shape[0] == samples2.shape[0] == matches.shape[0])
    if matches.shape[0] > N:
        samples11 = samples1[matches[0:N,0].astype(dtype=np.int16)]
        samples22 = samples2[matches[0:N,1].astype(dtype=np.int16)]
    else:
        samples11 = samples1[matches[:,0].astype(dtype=np.int16)]
        samples22 = samples2[matches[:,1].astype(dtype=np.int16)]
    assert (samples11.shape == samples22.shape)
    line = np.zeros((2, 3))
    i = 0
    correct = 0
    wrong = 0
    for sample1, sample2 in zip(samples11, samples22):
        line[0] = sample1
        line[1] = sample2
        if matches[i,0] == matches[i, 1]:
            mlab.plot3d(line[:, 0], line[:, 1], line[:, 2], color=(0, 1, 0), tube_radius=tube_radius)
            correct += 1
        else:
            mlab.plot3d(line[:, 0], line[:, 1], line[:, 2], color=(1, 0, 0), tube_radius=tube_radius)
            wrong += 1
        i += 1
    print 'correct: ', float(correct) / (correct + wrong), '    worng: ', float(wrong) / (correct + wrong)
    mlab.show()
    return -1

def plotOrs(*ptss):

    mlab.figure(34); mlab.clf()

    r = 1
    phi, theta = mgrid[0:pi:101j, 0:2*pi:101j]
    
    x = r*sin(phi)*cos(theta)
    y = r*sin(phi)*sin(theta)
    z = r*cos(phi)

    
    mlab.mesh(x,y,z, colormap='gray',opacity=.2)
    
    
    
    colors = [(1,0,0),(0,1,0),(0,0,1)]
    print len(colors)
    print len(ptss)
    for (pts,col) in izip(ptss,colors):
        print col        
        ors = normr(pts)
        # Create a sphere
        
        x,y,z = ors.T
        
        mlab.points3d(x,y,z,color=col)
        mlab.plot3d(x,y,z,color=col,tube_radius=.025)

def plot_matrix_metric(connectmat_file,centers_file,threshold_pct,grp_metrics=None,node_metric='bc',
                       weight_edges=0,node_scale_factor=2,edge_radius=.5,resolution=8,name_scale_factor=1,names_file=0,
                       red_nodes=None):
    """
    Given a connectivity matrix and a (x,y,z) centers file for each region, plot the 3D network
    """
    matrix=core.file_reader(connectmat_file)
    nodes=core.file_reader(centers_file)
    if names_file:
        names = core.file_reader(names_file,1)
    num_nodes = len(nodes)
    edge_thresh_pct = threshold_pct / 100.0
    matrix_flat=np.array(matrix).flatten()
    edge_thresh=np.sort(matrix_flat)[len(matrix_flat)-int(len(matrix_flat)*edge_thresh_pct)]
    
    if grp_metrics: # regional metrics caclulated elsewhere, loaded in
        node_colors = {}     # define colors for each metric
        node_colors['s'] = (1,  0.733,  0)
        node_colors['cc'] = (0.53, 0.81, 0.98)
        node_colors['bc'] = (0.5,  1,  0)
        node_colors['eloc'] = (1,  0 ,  1)
        node_colors['ereg'] = (1,  1,  0)
        
        node_metrics={}
        metrics = np.array(core.file_reader(grp_metrics))
        cols = np.shape(metrics)[1]
        for i in range(cols):
            colmean = np.mean(metrics[:,i])
            colscale = 2 / colmean
            metrics[:,i] = metrics[:,i] * colscale # make node mean size 2
        node_metrics['s'] = metrics[:,0] # strength
        node_metrics['cc'] = metrics[:,1] # clustering coefficient
        node_metrics['bc'] = metrics[:,2] # betweenness centrality
        node_metrics['eloc'] = metrics[:,3] # local efficiency
        node_metrics['ereg'] = metrics[:,4] # regional efficiency
    
    mlab.figure(bgcolor=(1, 1, 1), size=(800, 800))
    for count,(x,y,z) in enumerate(nodes):
        if grp_metrics:
            mlab.points3d(x,y,z, color=node_colors[node_metric],
                          scale_factor=node_metrics[node_metric][count], resolution=resolution)
        else:
            mlab.points3d(x,y,z, color=(0,1,0), scale_factor=node_scale_factor, resolution=resolution)
        if names_file:
            mlab.text(x, y,names[count], z=z,width=.02*name_scale_factor*len(names[count]),color=(0,0,0))
    
    for i in range(num_nodes-1):    
        x0,y0,z0=nodes[i]
        for j in range(i+1, num_nodes):
            if matrix[i][j] > edge_thresh:
                x1,y1,z1=nodes[j]
                if weight_edges:
                    mlab.plot3d([x0,x1], [y0,y1], [z0,z1],
                            tube_radius=matrix[i][j]/matrix_flat.max(),
                            color=(1,1,1))
                else:
                    mlab.plot3d([x0,x1], [y0,y1], [z0,z1],
                            tube_radius=edge_radius,
                            color=(1,1,1))

    def plot_3D_spectrum(self, xmin=None, xmax=None, xN=None, ymin=None,
                         ymax=None, yN=None, trajectory=False, tube_radius=1e-2,
                         part='imag'):
        """Plot the Riemann sheet structure around the EP.

            Parameters:
            -----------
                xmin, xmax: float
                    Dimensions in x-direction.
                ymin, ymax: float
                    Dimensions in y-direction.
                xN, yN: int
                    Number of sampling points in x and y direction.
                trajectory: bool
                    Whether to include a projected trajectory of the eigenbasis
                    coefficients.
                part: str
                    Which function to apply to the eigenvalues before plotting.
                tube_radius: float
                    Trajectory tube thickness.
        """
        from mayavi import mlab

        X, Y, Z = self.sample_H(xmin, xmax, xN, ymin, ymax, yN)
        Z0, Z1 = [Z[..., n] for n in (0, 1)]

        def get_min_and_max(*args):
            data = np.concatenate(*args)
            return data.min(), data.max()

        surf_kwargs = dict(colormap='Spectral', mask=np.diff(Z0.real) > 0.015)

        mlab.figure(0)
        Z_min, Z_max = get_min_and_max([Z0.real, Z1.real])
        mlab.surf(X.real, Y.real, Z0.real, vmin=Z_min, vmax=Z_max, **surf_kwargs)
        mlab.surf(X.real, Y.real, Z1.real, vmin=Z_min, vmax=Z_max, **surf_kwargs)
        mlab.axes(zlabel="Re(E)")

        mlab.figure(1)
        Z_min, Z_max = get_min_and_max([Z0.imag, Z1.imag])
        mlab.mesh(X.real, Y.real, Z0.imag, vmin=Z_min, vmax=Z_max, **surf_kwargs)
        mlab.mesh(X.real, Y.real, Z1.imag, vmin=Z_min, vmax=Z_max, **surf_kwargs)
        mlab.axes(zlabel="Im(E)")

        if trajectory:
            x, y = self.get_cycle_parameters(self.t)
            _, c1, c2 = self.solve_ODE()

            for i, part in enumerate([np.real, np.imag]):
                e1, e2 = [part(self.eVals[:, n]) for n in (0, 1)]
                z = map_trajectory(c1, c2, e1, e2)
                mlab.figure(i)
                mlab.plot3d(x, y, z, tube_radius=tube_radius)
                mlab.points3d(x[0], y[0], z[0],
                              # color=line_color,
                              scale_factor=1e-1,
                              mode='sphere')

        mlab.show()

def draw_line(pt_1, pt_2, color=(1.0, 0.0, 0.0)):
    """Draw a line between two points"""
    mlab.plot3d(
        np.hstack([pt_1[0], pt_2[0]]),
        np.hstack([pt_1[1], pt_2[1]]),
        np.hstack([pt_1[2], pt_2[2]]),
        color=color
    )

    def line(self,n,p,c=None,l=1,o=1,w=0.1): 
        s = np.linspace(0,l)
        l = [p+e*n for e in s]

        xs = [e[0] for e in l]
        ys = [e[1] for e in l]
        zs = [e[2] for e in l]
        mlab.plot3d(xs,ys,zs,figure=self.fig,color=c,opacity=o)