def generate_plots_3d(self):
        self.ax = mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(800, 600))
        self.clf = mlab.clf()

        minS, maxS = maxint, 0
        contour_plots = []
        for cond in self.conductors.itervalues():

            minS, maxS, face_data = self.generate_plot_data_for_faces_3d(cond, minS, maxS)
            for (x, y, z, s) in face_data:
                if isinstance(cond, conductor_type_3d['Unstructured']):
                    pts = mlab.points3d(x, y, z, s, scale_mode='none', scale_factor=0.002)
                    mesh = mlab.pipeline.delaunay3d(pts)
                    contour_plots.append(mlab.pipeline.surface(mesh, colormap='viridis'))
                else:
                    if np.min(s) < 0.0:
                        contour_plots.append(mlab.mesh(x, y, z, color=(0, 0, 0), colormap='viridis'))
                    else:
                        contour_plots.append(mlab.mesh(x, y, z, scalars=s, colormap='viridis'))

        for cp in contour_plots:
            cp.module_manager.scalar_lut_manager.trait_set(default_data_range=[minS * 0.95, maxS * 1.05])

        mlab.draw()
        mlab.colorbar(object=contour_plots[0], orientation='vertical')
        mlab.show()

def plotvfonsph3D(theta_rad, phi_rad, E_th, E_ph, freq=0.0,
                     vcoord='sph', projection='equirectangular'):
    PLOT3DTYPE = "quiver"
    (x, y, z) = sph2crtISO(theta_rad, phi_rad)
    from mayavi import mlab
    
    mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(400, 300))
    mlab.clf()
    if PLOT3DTYPE == "MESH_RADIAL" :
        r_Et = numpy.abs(E_th)
        r_Etmx = numpy.amax(r_Et)
        mlab.mesh(r_Et*(x)-1*r_Etmx, r_Et*y, r_Et*z, scalars=r_Et)
        r_Ep = numpy.abs(E_ph)
        r_Epmx = numpy.amax(r_Ep)
        mlab.mesh(r_Ep*(x)+1*r_Epmx , r_Ep*y, r_Ep*z, scalars=r_Ep)
    elif PLOT3DTYPE == "quiver":
        ##Implement quiver plot
        s2cmat = getSph2CartTransfMatT(numpy.array([x,y,z]))
        E_r = numpy.zeros(E_th.shape)
        E_fldsph = numpy.rollaxis(numpy.array([E_r, E_ph, E_th]), 0, 3)[...,numpy.newaxis]
        E_fldcrt = numpy.rollaxis(numpy.matmul(s2cmat, E_fldsph).squeeze(), 2, 0)
        #print E_fldcrt.shape
        mlab.quiver3d(x+1.5, y, z,
                      numpy.real(E_fldcrt[0]),
                      numpy.real(E_fldcrt[1]),
                      numpy.real(E_fldcrt[2]))
        mlab.quiver3d(x-1.5, y, z,
                      numpy.imag(E_fldcrt[0]),
                      numpy.imag(E_fldcrt[1]),
                      numpy.imag(E_fldcrt[2]))              
    mlab.show()

def plot_both_mayavi(f, xbounds=(-1, 1), ybounds=(-1, 1), res=401):
    """ Plot the real and imaginary parts of
    the function 'f', given the bounds and resolution. """
    X, Y, vals = get_vals(f, xbounds, ybounds, res)
    ml.mesh(X, Y, vals.real)
    ml.mesh(X, Y, vals.imag)
    ml.show()

def PlotHorizon3d(tss):
    """
    Plot a list of horizons.

    Parameters
    ----------

    tss : list of trappedsurface
        All the trapped surfaces to visualize.
    """
    from mayavi import mlab
    cmaps = ['bone', 'jet', 'hot', 'cool', 'spring', 'summer', 'winter']
    assert len(cmaps) > len(tss)
    extents = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    for ts, cm in zip(tss, cmaps):
        mlab.mesh(ts.X, ts.Y, ts.Z, colormap=cm, opacity=0.4)
        extents[0] = min(extents[0], np.min(ts.X))
        extents[1] = max(extents[1], np.max(ts.X))
        extents[2] = min(extents[2], np.min(ts.Y))
        extents[3] = max(extents[3], np.max(ts.Y))
        extents[4] = min(extents[4], np.min(ts.Z))
        extents[5] = max(extents[5], np.max(ts.Z))
    mlab.axes(extent=extents)
    mlab.outline(extent=extents)
    mlab.show()

 def plot_mayavi(self):
     """Use mayavi to plot a phenotype phase plane in 3D.
     The resulting figure will be quick to interact with in real time,
     but might be difficult to save as a vector figure.
     returns: mlab figure object"""
     from mayavi import mlab
     figure = mlab.figure(bgcolor=(1, 1, 1), fgcolor=(0, 0, 0))
     figure.name = "Phenotype Phase Plane"
     max = 10.0
     xmax = self.reaction1_fluxes.max()
     ymax = self.reaction2_fluxes.max()
     zmax = self.growth_rates.max()
     xgrid, ygrid = meshgrid(self.reaction1_fluxes, self.reaction2_fluxes)
     xgrid = xgrid.transpose()
     ygrid = ygrid.transpose()
     xscale = max / xmax
     yscale = max / ymax
     zscale = max / zmax
     mlab.surf(xgrid * xscale, ygrid * yscale, self.growth_rates * zscale,
               representation="wireframe", color=(0, 0, 0), figure=figure)
     mlab.mesh(xgrid * xscale, ygrid * yscale, self.growth_rates * zscale,
               scalars=self.shadow_prices1 + self.shadow_prices2,
               resolution=1, representation="surface", opacity=0.75,
               figure=figure)
     # draw axes
     mlab.outline(extent=(0, max, 0, max, 0, max))
     mlab.axes(opacity=0, ranges=[0, xmax, 0, ymax, 0, zmax])
     mlab.xlabel(self.reaction1_name)
     mlab.ylabel(self.reaction2_name)
     mlab.zlabel("Growth rates")
     return figure

def plot_sphere_func(f, grid='Clenshaw-Curtis', theta=None, phi=None, colormap='jet', fignum=0):

    # Note: all grids except Clenshaw-Curtis have holes at the poles

    import matplotlib
    matplotlib.use('WxAgg')
    matplotlib.interactive(True)
    from mayavi import mlab

    if grid == 'Driscoll-Healy':
        b = f.shape[0] / 2
    elif grid == 'Clenshaw-Curtis':
        b = (f.shape[0] - 2) / 2
    elif grid == 'SOFT':
        b = f.shape[0] / 2
    elif grid == 'Gauss-Legendre':
        b = (f.shape[0] - 2) / 2

    if theta is None or phi is None:
        theta, phi = meshgrid(b=b, convention=grid)

    phi = np.r_[phi, phi[0, :][None, :]]
    theta = np.r_[theta, theta[0, :][None, :]]
    f = np.r_[f, f[0, :][None, :]]

    x = np.sin(theta) * np.cos(phi)
    y = np.sin(theta) * np.sin(phi)
    z = np.cos(theta)

    mlab.figure(fignum, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(600, 400))
    mlab.clf()
    mlab.mesh(x, y, z, scalars=f, colormap=colormap)

    # mlab.view(90, 70, 6.2, (-1.3, -2.9, 0.25))
    mlab.show()

def zoncaview(m):
    """
    m is a healpix sky map, such as provided by WMAP or Planck.
    """

    nside = hp.npix2nside(len(m))
    vmin = -1e3; vmax = 1e3

    # Set up some grids:
    xsize = ysize = 1000
    theta = np.linspace(np.pi, 0, ysize)
    phi   = np.linspace(-np.pi, np.pi, xsize)
    longitude = np.radians(np.linspace(-180, 180, xsize))
    latitude = np.radians(np.linspace(-90, 90, ysize))

    # Project the map to a rectangular matrix xsize x ysize:
    PHI, THETA = np.meshgrid(phi, theta)
    grid_pix = hp.ang2pix(nside, THETA, PHI)
    grid_map = m[grid_pix]

    # Create a sphere:
    r = 0.3
    x = r*np.sin(THETA)*np.cos(PHI)
    y = r*np.sin(THETA)*np.sin(PHI)
    z = r*np.cos(THETA)

    # The figure:
    mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(400, 300))
    mlab.clf()

    mlab.mesh(x, y, z, scalars=grid_map, colormap="jet", vmin=vmin, vmax=vmax)

    mlab.draw()

    return

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 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 addDepthMap(figure,
                input_atoms,
                vmin,
                vmax,
                scale_factors=None,
                colormap='hot',
                resolution=32,
                ):
    from mayavi import mlab
    if scale_factors is None:
        scale_factors=np.ones(len(input_atoms))

    numbers = input_atoms.get_atomic_numbers()
    collections = set(zip(numbers, scale_factors))

    for number, scale_factor in collections:
        take1 = numbers == number
        take2 = scale_factors == scale_factor
        take = np.logical_and(take1, take2)
        atoms = input_atoms[take]
        points = atoms.positions
        radius = my_radii[number]/2.0
        radius *= scale_factor

        points = points[points[:, 2] > vmin - radius]

        for point  in points:
            x, y, z = sphere(point, radius, resolution=resolution)
            mlab.mesh(x,y,z,
                      scalars=z,
                      vmin=vmin,
                      vmax=vmax,
                      colormap=colormap,
                      figure=figure,
                      )

 def display(self,bgcolor=None,showAxes=False):
     mlab.figure(bgcolor=bgcolor)
     for x,y,z,op,col in zip(self.xMesh,self.yMesh,self.zMesh,self.meshOpacity,self.meshColor):
         mlab.mesh(x,y,z,opacity=op,color=col)
     if showAxes:
         mlab.axes()
     mlab.show()

def plotView(rays,pts=None, **kwargs):
      
    if not pts is None:
        x = scipy.zeros((len(rays)+1,pts))
        y = scipy.zeros(x.shape)
        z = scipy.zeros(x.shape)
        for i in rays:
            i.norm.s = scipy.linspace(i.norm.s[0],i.norm.s[-1],pts)
    else:
        x = scipy.zeros((len(rays)+1,len(rays[0].norm.s)))
        y = scipy.zeros(x.shape)
        z = scipy.zeros(x.shape)
    

    for i in xrange(len(rays)):
        if rays[i]._origin.flag:
            rays[i] = rays[i].c()
        x[i] = rays[i].x()[0]
        y[i] = rays[i].x()[1]
        z[i] = rays[i].x()[2]

    x[-1] = rays[0].x()[0]
    y[-1] = rays[0].x()[1]
    z[-1] = rays[0].x()[2]

    mlab.mesh(x,y,z,**kwargs)

    def draw_ballast(self, fig, centerline, freeboard, h_section, r_nodes, h_perm, h_water):
        from mayavi import mlab
        npts = 40
        th = np.linspace(0, 2*np.pi, npts)
        z_nodes = np.flipud( freeboard - np.r_[0.0, np.cumsum(np.flipud(h_section))] )

        # Permanent ballast
        z_perm = z_nodes[0] + np.linspace(0, h_perm, npts)
        r_perm = np.interp(z_perm, z_nodes, r_nodes)
        R, TH = np.meshgrid(r_perm, th)
        Z, _  = np.meshgrid(z_perm, th)
        X = R*np.cos(TH) + centerline[0]
        Y = R*np.sin(TH) + centerline[1]
        ck = np.array([122, 85, 33]) / 255.0
        ck = tuple(ck.tolist())
        mlab.mesh(X, Y, Z, color=ck, figure=fig)

        # Water ballast
        z_water = z_perm[-1] + np.linspace(0, h_water, npts)
        r_water = np.interp(z_water, z_nodes, r_nodes)
        R, TH = np.meshgrid(r_water, th)
        Z, _  = np.meshgrid(z_water, th)
        X = R*np.cos(TH) + centerline[0]
        Y = R*np.sin(TH) + centerline[1]
        ck = (0.0, 0.1, 0.8) # Dark blue
        mlab.mesh(X, Y, Z, color=ck, figure=fig)

def draw_encoders(encoders, angle=None, colors=None):
    """Add encoders to the scene.

    Can either supply a collection of encoders, or an integer giving the
    number of encoders to randomly generate.

    encoders: int or collection
    angle: angle giving size of the cap drawn to represent an encoder
    colors: collection of normalized rgb colors to paint the encoders
    """

    if isinstance(encoders, int):
        encoders = dists.UniformHypersphere(3, surface=True).sample(encoders)

    if colors is None:
        colors = [black for e in encoders]

    if not angle:
        angle = 0.01 * np.pi

    for enc, color in zip(encoders, colors):
        front = color
        back = black

        cap = threed.make_cap(r=1.01*radius, cap_angle=angle, direction=enc)
        mlab.mesh(*cap, color=front, opacity=1.0)

        cap = threed.make_cap(r=1.007*radius, cap_angle=angle, direction=enc)
        mlab.mesh(*cap, color=back, opacity=1.0)

def drawstick(p1, p2, radius=1, clr = (1,0,0), alpha = 1, samplenum = 6, ):
	if isinstance(radius,tuple):
		clr = radius[1]
		alpha = radius[2]
		samplenum = radius[3]
		radius = radius[0]

	pi = np.pi
	cos = np.cos
	sin = np.sin
	theta = np.linspace(0,2*pi,samplenum)
	normv = (p1-p2) / np.linalg.norm(p1-p2)
	u = np.empty((3,1))
	if normv[0]!= 0:
		u[1] = u[2] = 1
		u[0] = -(normv[1]+normv[2])/normv[0]
	elif normv[1]!= 0:
		u[0] = u[2] = 1
		u[1] = -(normv[0]+normv[2])/normv[1]
	else:# normv[2]!= 0
		u[0] = u[1] = 1
		u[2] = -(normv[1]+normv[0])/normv[2]
	u = u / np.linalg.norm(u)
	circle = radius * np.multiply(cos(theta),u) + radius * np.multiply( sin(theta), np.cross(normv.T, u.T).T )
	p1arr = np.tile(p1,circle.shape[1])
	p2arr = np.tile(p2,circle.shape[1])
	c1 = circle + p1arr
	c2 = circle + p2arr
	x = np.vstack((p1arr[0],c1[0],c2[0],p2arr[0]))
	y = np.vstack((p1arr[1],c1[1],c2[1],p2arr[1]))
	z = np.vstack((p1arr[2],c1[2],c2[2],p2arr[2]))

	mlab.mesh(x,y,z, color=clr, opacity=alpha)

    def plot(self, Ps):
        from mayavi import mlab

        mlab.figure(mlab.figure(), bgcolor=(1, 1, 1))
        for s in range(len(Ps)):
            mlab.mesh(Ps[s][:, :, 0], Ps[s][:, :, 1], Ps[s][:, :, 2], color=(65 / 256, 105 / 256, 225 / 256))
        mlab.show()

    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 makefinegrid():
    coor=make_coor(20j)

    mlab.figure()
    mlab.mesh(coor.x, coor.y, coor.z)
    mlab.points3d(coor.x, coor.y, coor.z, 
                      scale_factor=0.1)
    savefig('finegrid.png')

    def plane(self,n,p,c=None):
        d = -np.sum(p*n)
        x = np.linspace(-2,2)
        y = np.linspace(-2,2)
        [xx,yy]=np.meshgrid(x,y);
        zz = (-n[0]*xx - n[1]*yy - d)/n[2]

        mlab.mesh(xx,yy,zz,opacity=0.15,figure=self.fig)

def plot_alpha_sphere():
  x,y,z = make_sphere(50, 1.0)

  mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(400, 300))
  mlab.clf()

  s = np.zeros(x.shape)
  mlab.mesh(x, y, z, scalars=s, colormap='jet', opacity=0.1)
  mlab.show(stop=True)