def plot(func, dataset, x1s, x1e, x2s, x2e, delta, levels, title=None, track=None, f_val=None):
    # t, a, b, c = dataset

    k1 = np.arange(x1s, x1e, delta)
    k2 = np.arange(x2s, x2e, delta)
    K1, K2 = np.meshgrid(k1, k2)
    FX = np.zeros(K1.shape)
    row, col = K1.shape
    for i in xrange(row):
        for j in xrange(col):
            FX[i, j] = func(array([K1[i, j], K2[i, j]]), *dataset)

    fig = plt.figure(title)
    subfig1 = fig.add_subplot(1, 2, 1)
    surf1 = plt.contour(K1, K2, FX, levels=levels, stride=0.001)

    if track:
        track = array(track)
        plt.plot(track[:, 0], track[:, 1])
    plt.xlabel("k1")
    plt.ylabel("k2")

    subfig2 = fig.add_subplot(1, 2, 2, projection='3d')
    surf2 = subfig2.plot_wireframe(K1, K2, FX, rstride=10, cstride=10, color='y')
    plt.contour(K1, K2, FX, stride=1, levels=levels)
    if track != None and f_val != None:
        f_val = array(f_val)
        subfig2.scatter(track[:, 0], track[:, 1], f_val)
    plt.show()

def plot_solution(problem_data, solution):
    from numpy import linspace
    from matplotlib.pylab import contour, colorbar, contourf, xlabel, ylabel, title, show
    from matplotlib.mlab import griddata

    print(" * Preparing for plotting...")
    NN = problem_data["NN"]

    # Extract node coordinates seperately for plotting
    x, y = [0] * NN, [0] * NN
    for i, node in enumerate(problem_data["nodes"]):
        x[i] = node[0]
        y[i] = node[1]

    # Refine the contour plot mesh for a "smoother" image, by generating a
    # 200*200 grid
    xi = linspace(min(x), max(x), 200)
    yi = linspace(min(y), max(y), 200)

    # Approximate the mid values from neighbors
    zi = griddata(x, y, solution, xi, yi)

    print(" * Plotting...")
    # Plot the contour lines with black
    contour(xi, yi, zi, 15, linewidths=0.5, colors='k')
    # Plot the filled contour plot
    plot = contourf(xi, yi, zi, 15, antialiased=True)

    colorbar(plot, format="%.3f").set_label("T")
    xlabel('X')
    ylabel('Y')
    title("Contour plot of T values for {0}".format(problem_data["title"]))

    show()

def plot_svc(X, y, mysvc, bounds=None, grid=50):
    if bounds is None:
        xmin = np.min(X[:, 0], 0)
        xmax = np.max(X[:, 0], 0)
        ymin = np.min(X[:, 1], 0)
        ymax = np.max(X[:, 1], 0)
    else:
        xmin, ymin = bounds[0], bounds[0]
        xmax, ymax = bounds[1], bounds[1]
    aspect_ratio = (xmax - xmin) / (ymax - ymin)
    xgrid, ygrid = np.meshgrid(np.linspace(xmin, xmax, grid),
                              np.linspace(ymin, ymax, grid))
    plt.gca(aspect=aspect_ratio)
    plt.xlim(xmin, xmax)
    plt.ylim(ymin, ymax)
    plt.xticks([])
    plt.yticks([])
    plt.hold(True)
    plt.plot(X[y == 1, 0], X[y == 1, 1], 'bo')
    plt.plot(X[y == -1, 0], X[y == -1, 1], 'ro')
    
    box_xy = np.append(xgrid.reshape(xgrid.size, 1), ygrid.reshape(ygrid.size, 1), 1)
    if mysvc is not None:
        scores = mysvc.decision_function(box_xy)
    else:
        print 'You must have a valid SVC object.'
        return None;
    
    CS=plt.contourf(xgrid, ygrid, scores.reshape(xgrid.shape), alpha=0.5, cmap='jet_r')
    plt.contour(xgrid, ygrid, scores.reshape(xgrid.shape), levels=[0], colors='k', linestyles='solid', linewidths=1.5)
    plt.contour(xgrid, ygrid, scores.reshape(xgrid.shape), levels=[-1,1], colors='k', linestyles='dashed', linewidths=1)
    plt.plot(mysvc.support_vectors_[:,0], mysvc.support_vectors_[:,1], 'ko', markerfacecolor='none', markersize=10)
    CB = plt.colorbar(CS)

 def plot(x,y,field,filename,c=200):
     plt.figure()
     # define grid.
     xi = np.linspace(min(x),max(x),100)
     yi = np.linspace(min(y),max(y),100)
     # grid the data.
     si_lin = griddata((x, y), field, (xi[None,:], yi[:,None]), method='linear')
     si_cub = griddata((x, y), field, (xi[None,:], yi[:,None]), method='linear')
     print np.min(field)
     print np.max(field)
     plt.subplot(211)
     # contour the gridded data, plotting dots at the randomly spaced data points.
     CS = plt.contour(xi,yi,si_lin,c,linewidths=0.5,colors='k')
     CS = plt.contourf(xi,yi,si_lin,c,cmap=plt.cm.jet)
     plt.colorbar() # draw colorbar
     # plot data points.
     #    plt.scatter(x,y,marker='o',c='b',s=5)
     plt.xlim(min(x),max(x))
     plt.ylim(min(y),max(y))
     plt.title('Lineaarinen interpolointi')
     #plt.tight_layout()
     plt.subplot(212)
     # contour the gridded data, plotting dots at the randomly spaced data points.
     CS = plt.contour(xi,yi,si_cub,c,linewidths=0.5,colors='k')
     CS = plt.contourf(xi,yi,si_cub,c,cmap=plt.cm.jet)
     plt.colorbar() # draw colorbar
     # plot data points.
     #    plt.scatter(x,y,marker='o',c='b',s=5)
     plt.xlim(min(x),max(x))
     plt.ylim(min(y),max(y))
     plt.title('Kuubinen interpolointi')
     plt.savefig(filename)

    def _plot_nullclines(self, resolution):
        """
        Plot nullclines.

        Arguments
            resolution
                Resolution of plot
        """
        x_mesh, y_mesh, ode_x, ode_y = self._get_ode_values(resolution)

        plt.contour(
            x_mesh, y_mesh, ode_x,
            levels=[0], linewidths=2, colors='black')
        plt.contour(
            x_mesh, y_mesh, ode_y,
            levels=[0], linewidths=2, colors='black',
            linestyles='dashed')

        lblx = mlines.Line2D(
            [], [],
            color='black',
            marker='', markersize=15,
            label=r'$\dot\varphi_0=0$')
        lbly = mlines.Line2D(
            [], [],
            color='black', linestyle='dashed',
            marker='', markersize=15,
            label=r'$\dot\varphi_1=0$')
        plt.legend(handles=[lblx, lbly], loc='best')

def plot_hyperplane(X, Y, model, K, plot_id, d = 500):
    I0 = np.where(Y==-1)[0]
    I1 = np.where(Y==1)[0]

    plt.subplot(plot_id)

    plt.plot(X[I1, 0], X[I1, 1], 'og')
    plt.plot(X[I0, 0], X[I0, 1], 'xb')

    min_val = np.min(X, 0)
    max_val = np.max(X, 0)

    clf = model()
    clf.train(X, Y, K)

    x0_plot = np.linspace(min_val[0, 0], max_val[0, 0], d)
    x1_plot = np.linspace(min_val[0, 1], max_val[0, 1], d)

    [x0, x1] = plt.meshgrid(x0_plot, x1_plot);

    Y_all = np.matrix(np.zeros([d, d]))

    for i in range(d):
        X_all = np.matrix(np.zeros([d, 2]))
        X_all[:, 0] = np.matrix(x0[:, i]).T
        X_all[:, 1] = np.matrix(x1[:, i]).T
        Y_all[:, i] = clf.predict(X_all)

    plt.contour(np.array(x0), np.array(x1), np.array(Y_all), levels = [0.0], colors = 'red')

def zsview(im, cmap=pl.cm.gray, figsize=(8,5), contours=False, ccolor='r'):
    z1, z2 = zscale(im)
    pl.figure(figsize=figsize)
    pl.imshow(im, vmin=z1, vmax=z2, origin='lower', cmap=cmap, interpolation='none')
    if contours:
        pl.contour(im, levels=[z2], origin='lower', colors=ccolor)
    pl.tight_layout()

def plot_cost_function(data):
  x = data[:, 0]
  y = data[:, 1]
  theta_0_mesh, theta_1_mesh = np.meshgrid(np.arange(-10, 10, 0.01), np.arange(-1, 4, 0.01))
  cost_all_points = np.zeros(theta_0_mesh.shape)
  cost_all_points = compute_costs_thetas(theta_0_mesh, theta_1_mesh, x, y)
  plt.contour(theta_0_mesh, theta_1_mesh, cost_all_points, 100)    
  plt.show()

 def sanity_steppar2(self):
   import numpy as np
   import matplotlib.pylab as plt
   from PyAstronomy import funcFit as fuf
   
   # Set up a Gaussian model
   # and create some "data"
   x = np.linspace(0,2,100)
   gf = fuf.GaussFit1d()
   gf["A"] = 0.87
   gf["mu"] = 1.0
   gf["sig"] = 0.2
   y = gf.evaluate(x)
   y += np.random.normal(0.0, 0.1, len(x))
   
   # Thaw parameters, which are to be fitted ...
   gf.thaw(["A", "mu", "sig"])
   # ... and "disturb" starting values a little.
   gf["A"] = gf["A"] + np.random.normal(0.0, 0.1)
   gf["mu"] = gf["mu"] + np.random.normal(0.0, 0.1)
   gf["sig"] = gf["sig"] + np.random.normal(0.0, 0.03)
   # Find the best fit solution
   gf.fit(x, y, yerr=np.ones(len(x))*0.1)
   
   # Step the amplitude (area of the Gaussian) and the
   # center ("mu") of the Gaussian through the given
   # ranges.
   sp = gf.steppar(["A", "mu"], ranges={"A":[0.8, 0.95, 20], \
                   "mu":[0.96,1.05,15]})
   
   # Get the values for `A`, `mu`, and chi-square
   # from the output of steppar.
   As = map(lambda x:x[0], sp)
   mus = map(lambda x:x[1], sp)
   chis = map(lambda x:x[2], sp)
   
   # Create a chi-square array using the
   # indices contained in the output.
   z = np.zeros((20, 15))
   for s in sp:
     z[s[3]] = s[2]
   
   # Find minimum chi-square and define levels
   # for 68%, 90%, and 99% confidence intervals.
   cm = min(chis)
   levels = [cm+2.3, cm+4.61, cm+9.21]
   
   # Plot the contours to explore the confidence
   # interval and correlation.
   plt.xlabel("mu")
   plt.ylabel("A")
   plt.contour(np.sort(np.unique(mus)), np.sort(np.unique(As)), z, \
               levels=levels)
   # Plot the input value
   plt.plot([1.0], [0.87], 'k+', markersize=20)

 def draw_grid_plot():
     if not has_plt:
         return
     plt.imshow(z_sample, cmap=plt.cm.binary, origin='lower',
                interpolation='nearest',
                extent=(grid_min, grid_max, grid_min, grid_max))
     plt.hold(True)  # XXX: restore plt state
     plt.contour(np.linspace(grid_min, grid_max, grid_n),
                 np.linspace(grid_min, grid_max, grid_n),
                 z_actual)
     plt.savefig(grid_filename)
     plt.close()

    def show(self, rootFilename):
        """
        Plot the solution
        @param rootFilename a wild card expression, eg 'XXX_rk%d.txt'
        """
        import glob
        from matplotlib import pylab
        import re
        import os
        
        dataWindows = {}
        for f in glob.glob(rootFilename):
            pe = int(re.search(r'_rk(\d+)\.txt', f).group(1))
            dataWindows[pe] = numpy.loadtxt(f, unpack=False)
        
        nBigSizes = (self.decomp[0]*self.n,
                     self.decomp[1]*self.n)
        bigData = numpy.zeros(nBigSizes, numpy.float64)
        for j in range(self.decomp[0]):
            begJ = j*self.n
            endJ = begJ + self.n
            for i in range(self.decomp[1]):
                pe = j*self.decomp[1] + i
                begI = i*self.n
                endI = begI + self.n
                minVal = min(dataWindows[pe].flat)
                maxVal = max(dataWindows[pe].flat)
                bigData[begJ:endJ, begI:endI] = dataWindows[pe]

        xs = self.x0s[1] + numpy.array([(i + 0.5)*self.hs[1] for i in \
                          range(1, self.decomp[1]*self.n-1)])
        ys = self.x0s[0] + numpy.array([(j + 0.5)*self.hs[0] for j in \
                          range(1, self.decomp[0]*self.n-1)])
        pylab.contour(xs, ys, bigData[1:-1,1:-1], 21)
        pylab.colorbar()
        # plot the domain decomp
        for m in range(1, self.decomp[0]):
            yVal = self.x0s[0] + m * self.ls[0] - self.hs[0] / 2.0
            pylab.plot([xs[0], xs[-1]], [yVal, yVal], 'k--')
            yVal += self.hs[0]
            pylab.plot([xs[0], xs[-1]], [yVal, yVal], 'k--')
        for n in range(1, self.decomp[1]):
            xVal = self.x0s[1] + n * self.ls[1] - self.hs[1] / 2.0
            pylab.plot([xVal, xVal], [ys[0], ys[-1]], 'k--')
            xVal += self.hs[1]
            pylab.plot([xVal, xVal], [ys[0], ys[-1]], 'k--')
        pylab.show()
        fname = re.sub(r'_rk\*.txt', '', rootFilename) + '.png'
        print 'saving colorplot in file ' + fname
        pylab.savefig(fname)

def boundary(X, y, theta):
    y_ = y.reshape(y.size,)
    classes = (0., 1.)
    colors = ('b', 'r')
    marks = ('o', '+')
    for i, color, mark in zip(classes, colors, marks):
        plt.scatter(np.array(X)[y_ == i, 0], np.array(X)[y_ == i, 1], c=color, marker=mark)

    def fn((ui,vj)):
        return map_feature(np.array(ui), np.array(vj)).dot(theta)
    plt.contour(*calc_contour(X, fn))
    
    plt.legend(['y=%s' % i for i in classes] + ['Decision boundary'])
    plt.show()    

    def plot(self, func, interp=True, plotter='imshow'):
        import matplotlib as mpl
        from matplotlib import pylab as pl
        if interp:
            lpi = self.interpolator(func)
            z = lpi[self.yrange[0]:self.yrange[1]:complex(0, self.nrange),
                    self.xrange[0]:self.xrange[1]:complex(0, self.nrange)]
        else:
            y, x = np.mgrid[
                self.yrange[0]:self.yrange[1]:complex(0, self.nrange),
                self.xrange[0]:self.xrange[1]:complex(0, self.nrange)]
            z = func(x, y)

        z = np.where(np.isinf(z), 0.0, z)

        extent = (self.xrange[0], self.xrange[1],
            self.yrange[0], self.yrange[1])
        pl.ioff()
        pl.clf()
        pl.hot()  # Some like it hot
        if plotter == 'imshow':
            pl.imshow(np.nan_to_num(z), interpolation='nearest', extent=extent,
                      origin='lower')
        elif plotter == 'contour':
            Y, X = np.ogrid[
                self.yrange[0]:self.yrange[1]:complex(0, self.nrange),
                self.xrange[0]:self.xrange[1]:complex(0, self.nrange)]
            pl.contour(np.ravel(X), np.ravel(Y), z, 20)
        x = self.x
        y = self.y
        lc = mpl.collections.LineCollection(
            np.array([((x[i], y[i]), (x[j], y[j]))
                      for i, j in self.tri.edge_db]),
            colors=[(0, 0, 0, 0.2)])
        ax = pl.gca()
        ax.add_collection(lc)

        if interp:
            title = '%s Interpolant' % self.name
        else:
            title = 'Reference'
        if hasattr(func, 'title'):
            pl.title('%s: %s' % (func.title, title))
        else:
            pl.title(title)

        pl.show()
        pl.ion()

def contour(X,Y,Z, 
           extent=None, vrange=None, levels=None, extend='both', 
           inaxis=None,
           cmap=None, addcolorbar=True, clabel=None,
           smooth=True, smoothlen=None):
    '''Build a super fancy contour image'''
    
    # Build up some nice ranges and levels to be plotted 
    XX, YY = _getmesh(X,Y,Z)
    extent = _getextent(extent, XX, YY)
    vrange = _getvrange(vrange, XX,YY,Z, inaxis=inaxis)
    levels = _getlevels(levels, vrange)
    cmap   = _getcmap(cmap)
    
    # Smooth if needed
    if smooth:
        X,Y,Z = _smooth(X,Y,Z, smoothlen)
    
    cs = pylab.contourf(X, Y, Z, levels,
                       vmin=vrange[0], vmax=vrange[1],
                       extent=extent, extend='both',
                       cmap=cmap)
    ccs = pylab.contour(X, Y, Z, levels, vmin=vrange[0], vmax=vrange[1],
                        cmap=cmap)
    
    # setup a colorbar, add in the lines, and then return it all out.
    if addcolorbar:
        cb = colorbar(cs, ccs, levels=levels, clabel=clabel)
        return cs, ccs, cb
    
    return cs, ccs

def plot_contour(z,x,y,title="TITLE",xtitle="",ytitle="",xrange=None,yrange=None,plot_points=0,contour_levels=20,
                 cmap=None,cbar=True,fill=False,cbar_title="",show=1):

    fig = plt.figure()

    if fill:
        fig = plt.contourf(x, y, z.T, contour_levels, cmap=cmap, origin='lower')
    else:
        fig = plt.contour( x, y, z.T, contour_levels, cmap=cmap, origin='lower')

    if cbar:
        cbar = plt.colorbar(fig)
        cbar.ax.set_ylabel(cbar_title)


    plt.title(title)
    plt.xlabel(xtitle)
    plt.ylabel(ytitle)

    # the scatter plot:
    if plot_points:
        axScatter = plt.subplot(111)
        axScatter.scatter( np.outer(x,np.ones_like(y)), np.outer(np.ones_like(x),y))

    # set axes range
    plt.xlim(xrange)
    plt.ylim(yrange)

    if show:
        plt.show()

    return fig

def viz_data(X, mu, sigma_sq):
  X1, X2 = np.meshgrid(np.arange(0,35,0.5), np.arange(0,35,0.5))
  X1_flattened = X1.flatten()
  X2_flattened = X2.flatten()
  new_X = np.empty((X1_flattened.size,2))
  new_X[:, 0] = X1_flattened
  new_X[:, 1] = X2_flattened
  Z = multivariate_gaussian(new_X, mu, sigma_sq)
  Z = Z.reshape(X1.shape)
  plt.hold()
  plt.plot(X[:,0], X[:,1], 'bx')
  plt.contour(X1, X2, Z, levels=10.0**(np.arange(-20.0,0,3)))
  plt.axis([0, 30, 0, 30])
  plt.xlabel('Latency ms')
  plt.ylabel('Tput Mbps') 
  plt.hold()

def plot_model(data_array, mixture, axis=(0, 1), samples=20, contour_lines=20):
    """ Plot the scaterplot of data_array and the contour lines of the
    probability for the mixture.
     
    """
    import matplotlib
    import matplotlib.pylab as plt
    import matplotlib.cm as cm
    
    axis = list(axis)
    
    if isinstance(mixture, GMClusterModel):
        mixture = mixture.norm_model
    
    if isinstance(data_array, Orange.data.Table):
        data_array, _, _ = data_array.to_numpy_MA()
    array = data_array[:, axis]
    
    weights = mixture.weights
    means = mixture.means[:, axis]
    
    covariances = [cov[axis,:][:, axis] for cov in mixture.covariances] # TODO: Need the whole marginal distribution. 
    
    gmm = GMModel(weights, means, covariances)
    
    min = numpy.min(array, 0)
    max = numpy.max(array, 0)
    extent = (min[0], max[0], min[1], max[1])
    
    X = numpy.linspace(min[0], max[0], samples)
    Y = numpy.linspace(min[1], max[1], samples)
    
    Z = numpy.zeros((X.shape[0], Y.shape[0]))
    
    for i, x in enumerate(X):
        for j, y in enumerate(Y):
            Z[i, j] = gmm([x, y])
            
    plt.plot(array[:,0], array[:,1], "ro")
    plt.contour(X, Y, Z.T, contour_lines,
                extent=extent)
    
    im = plt.imshow(Z.T, interpolation='bilinear', origin='lower',
                cmap=cm.gray, extent=extent)
    
    plt.plot(means[:, 0], means[:, 1], "b+")
    plt.show()

def plot(points, col):
    n = 256
    x = np.linspace(-100, 100, n)
    y = np.linspace(-100, 100, n)
    X, Y = np.meshgrid(x, y)
    
    xs = []
    ys = []
    
    pl.contourf(X, Y, fun([X, Y]), 8, alpha=.75, cmap='jet')
    pl.contour(X, Y, fun([X, Y]), 8, colors='black', linewidth=.5) 
    
    for i in range(len(points)):
        xs.append(points[i][0])
        ys.append(points[i][1])
    
    pl.plot(xs, ys, marker='o', linestyle='--', color=str(col), label='Square')

def plot_contour_pair(xi, yi, zi):
    fig = plt.figure(figsize=(2 * w, golden_mean * w * 2))
    ax = fig.add_subplot(1, 2, 1, projection='3d')
    cs = plt.contour(xi, yi, zi, 15, linewidths=0.5, color='k')
    ax = fig.add_subplot(1, 2, 2, projection='3d')

    xig, yig = np.meshgrid(xi, yi)

    surf = ax.plot_surface(xig, yig, zi, linewidth=0)
    return fig

def conto_plots(Data,DataMHD,**kwargs):
    Qu = jana.quantities()
    Cur = Qu.Current(Data)
    CurMHD = Qu.Current(DataMHD)
    Fastsp = Qu.Magspeed(Data)['fast']
    Slowsp = Qu.Magspeed(Data)['slow']
    Alfvsp = Qu.Magspeed(Data)['alfven']
    Vpol = np.sqrt(Data.v1**2 + Data.v2**2)

    FastspMHD = Qu.Magspeed(DataMHD)['fast']
    SlowspMHD = Qu.Magspeed(DataMHD)['slow']
    AlfvspMHD = Qu.Magspeed(DataMHD)['alfven']
    VpolMHD = np.sqrt(DataMHD.v1**2 + DataMHD.v2**2)

    fastrat = Fastsp/Vpol
    slowrat = Slowsp/Vpol
    Alfvrat = Alfvsp/Vpol

    fastratMHD = FastspMHD/VpolMHD
    slowratMHD = SlowspMHD/VpolMHD
    AlfvratMHD = AlfvspMHD/VpolMHD

    

    f1 = plt.figure(num=1)
    ax1 = f1.add_subplot(121)
    currcont=plt.contour(Data.x1,Data.x2,Cur.T,kwargs.get('Currents',[-0.7,-0.6,-0.5,-0.4,-0.3]),colors=kwargs.get('colors','r'),linestyles=kwargs.get('ls','-'),lw=kwargs.get('lw',2.0))
    plt.clabel(currcont,manual=True)
    currmhdcont = plt.contour(DataMHD.x1,DataMHD.x2,CurMHD.T,kwargs.get('Currents',[-0.7,-0.6,-0.5,-0.4,-0.3]),colors=kwargs.get('colors','k'),linestyles=kwargs.get('ls','-'),lw=kwargs.get('lw',2.0))
    plt.clabel(currmhdcont,manual=True)
    plt.xlabel(r'r [AU]')
    plt.ylabel(r'z [AU]')
    plt.title(r'Total Current -- $\int\int J_{z} r dr d\phi$')
    ax2 = f1.add_subplot(122)
    plt.xlabel(r'r [AU]')
    #plt.ylabel(r'$z [AU]$')
    plt.title(r'Critical Surfaces')
    
    fcont=plt.contour(Data.x1,Data.x2,fastrat.T,[1.0],colors='r',linestyles='solid')
    plt.clabel(fcont,inline=1,fmt=r'Fast')
    scont=plt.contour(Data.x1,Data.x2,slowrat.T,[1.0],colors='r',linestyles='dashdot')
    plt.clabel(scont,inline=1,fmt=r'Slow')
    acont=plt.contour(Data.x1,Data.x2,Alfvrat.T,[1.0],colors='r',linestyles='dashed')
    plt.clabel(acont,inline=1,fmt=r'Alfv$\acute{e}$n')

    mfcont=plt.contour(DataMHD.x1,DataMHD.x2,fastratMHD.T,[1.0],colors='k',linestyles='solid')
    plt.clabel(mfcont,manual=1,fmt=r'Fast')
    mscont=plt.contour(DataMHD.x1,DataMHD.x2,slowratMHD.T,[1.0],colors='k',linestyles='dashdot')
    plt.clabel(mscont,manual=1,fmt=r'Slow')
    macont=plt.contour(DataMHD.x1,DataMHD.x2,AlfvratMHD.T,[1.0],colors='k',linestyles='dashed')
    plt.clabel(macont,manual=1,fmt=r'Alfv$\acute{e}$n')