setplot.py.html
Source file:   setplot.py
Directory:   /var/www/html/clawpack/links/shockvacuum10/shockvacuum
Converted:   Mon Aug 16 2010 at 14:13:32   using clawcode2html
This documentation file will not reflect any later changes in the source file.

 

"""
Set up the plot figures, axes, and items to be done for each frame.

This module is imported by the plotting routines and then the
function setplot is called to set the plot parameters.

This version requires Clawpack 4.4.
   $ svn co  http://kingkong.amath.washington.edu/svn/claw4/trunk claw4
Or download the tarfile (see documentation).

For documentation, see ...
   http://kingkong.amath.washington.edu/claw4/www/users/

"""

import os
import numpy as np
import pylab
from pyclaw.data import Data
import similarity_solution as ss
reload(ss)  # in case it has been changed by user

# Global parameters used in functions below:
#-------------------------------------------

# Specify where the data comes from for plotting.  
#outdir = '_output_poly_n1_1200'
#outdir = '_output_poly_n1_2400'
#outdir = '_output_poly_n1_4800'
#outdir = '_output_poly_n1_9600'

#outdir = '_output_euler_n1_alpha1_1200'
#outdir = '_output_euler_n1_alpha1_4800'

#outdir = '_output_euler_n1_alpha1.5_1200'
#outdir = '_output_euler_n1_alpha1.5_4800'

#outdir = '_output_euler_n1.5_alpha1.5_1200'
#outdir = '_output_euler_n1.5_alpha1.5_4800'

#outdir = '_output'

# Read the values from setprob.data, which should reflect the values set
# at the top of setrun.py.
setprob_file = os.path.join(outdir, 'setprob.data')
params = Data(setprob_file)

ieos = params.ieos       # = 1 for Euler, 2 for polytropic
n = params.npoly
nstar = params.nstar     # = alpha
K = params.sK
grav = params.grav
gamma = 1. + 1./n

gamma_star = 1. + 1./nstar

tfontsize = 20  # fontsize for some titles

plot_simsol = True   # Set to False to suppress similarity solution plots.


if plot_simsol & (ieos==1):
    if (n==1) & (nstar==1.):
        fname_simsol = 'simsol_euler_n1.dat'
    elif (n==1.1) & (nstar==1.1):
        fname_simsol = 'simsol_euler_n1.1.dat'
    elif (n==1.5) & (nstar==1.5):
        fname_simsol = 'simsol_euler_n3o2.dat'
        #fname_simsol = 'simsol_euler_n1.5_alpha1.5.dat'
    elif (n==1.) & (nstar==1.5):
        fname_simsol = 'simsol_euler_n1_alpha1.5.dat'
    elif (n==1.5) & (nstar==3.):
        fname_simsol = 'simsol_euler_n3o2_alpha3.dat'
    elif n==3:
        fname_simsol = 'simsol_euler_n3.dat'
    else:
        print '*** Cannot find a data file for similarity solution'
        plot_simsol = False

print "Assuming computation was done with ieos = %i, n = %g" % (ieos,n)
if ieos==1:
    print " and nstar = ",nstar


#------------------------------------------


# Read in t and xs from fort.xs
try:
    fname = outdir+'/fort.xs'
    txs_array = np.loadtxt(fname, dtype=float)
except:
    print "*** Error reading ",fname
    txs_array = None
    print "exists? ",os.path.isfile(fname)
    raise


if plot_simsol:

    if ieos==2:
        # compute the similarity solution for polytropic:
        #vstar = 4.8905  # for n=2, from numerics with order=2, mx=2400
        vstar = 3.7418 # for n=1, mx=1200
        vstar = 4.89  # for n=1. mx=4800
        vstar = 4.5  # for n=3. mx=4800
        vstar = 4.89  # for n=1. mx=2400
        tfinal = 2.59125913623e-07  # for polytropic npoly=1, mx=2400
        tfinal = 2.4477500128e-07  # for polytropic npoly=1, mx=4800
        tfinal = 2.641111388e-07  # for polytropic npoly=3, mx=4800
        tfinal = 2.4967677e-07  # for polytropic npoly=1, mx=2400
        tfinal = 2.447750014889652e-07  # for polytropic npoly=1, mx=4800
        # xreset = -1e-4:
        tfinal = 4.0066713534085611e-06  # for polytropic npoly=1, mx=4800,

        tfinal, vstar = ss.estimate_params_polytropic(txs_array)
        print 'tfinal = ',tfinal
        #vstar = 4.54  # for n=1.5 mx=2400
        #vstar = 4.16  # for n=3.0 mx=1200
        #vstar = 4.89  # for n=1.0 mx=1200
        #vstar = 4.92  # for n=1.0 mx=1200

        # Reset vstar based on Riemann invariant to better match...
        vstar = 4.893
        print 'Reset vstar = ',vstar

        pylab.figure(23)
        pylab.clf()
        xs = txs_array[1:,1]
        t_xs = txs_array[1:,0] - tfinal
        pylab.loglog(-t_xs, -xs/vstar, 'b', linewidth=3)
        pylab.title('Computed shock location')
        pylab.xlabel(r'$t_f-t_c$', fontsize=tfontsize)
        pylab.ylabel(r'$-x_s/v_*$', fontsize=tfontsize)
        pylab.axis('scaled')
        pylab.xticks(np.logspace(-15,-5,6), fontsize=tfontsize)
        pylab.yticks(np.logspace(-15,-5,6), fontsize=tfontsize)

        nypts = 1000  # number of points to evaluate sim. soln.
        ymax = 1000.
        print "Computing similarity solution using tfinal = ",tfinal
        simsol = ss.compute_simsol_polytropic(ymax, nypts, grav, n, vstar, tfinal)

    else:
        # Read in the tabulated similarity solution for Euler:
        print "Reading similarity solution from file ",fname_simsol
        #print "  using tfinal = ",tfinal
        simsol = ss.read_tabular_simsol(fname_simsol)
        simsol.grav = grav

        tfinal, C_y = ss.estimate_params(simsol.beta, simsol.ys, txs_array)
        simsol.C_y = C_y
        simsol.tfinal = tfinal
        print "Estimated tfinal = %22.15e" % tfinal
        print "Estimated    C_y = %22.15e" % C_y

    simsol.K = K
    simsol.gamma = gamma
    C_rho = (grav / ((nstar+1)*K))**nstar
    simsol.C_rho = C_rho
    print "C_rho = ",C_rho


ftreset = os.path.join(outdir, 'fort.treset')
if os.path.isfile(ftreset):
    treset = float(open(ftreset,'r').readline().split()[0])
else:
    treset = 0.



#---------------------
def setplot(plotdata):
#---------------------

    """
    Specify what is to be plotted at each frame.
    Input:  plotdata, an instance of pyclaw.plotters.data.ClawPlotData.
    Output: a modified version of plotdata.
    """

    from pyclaw.plotters.data import ClawPlotData
    from pyclaw.plotters.frametools import var_limits

    if not isinstance(plotdata,ClawPlotData):
        print '*** Error in setplot: plotdata must be a ClawPlotData object'
        return

    plotdata.clearfigures()  # clear any old figures,axes,items data
    plotdata.outdir = outdir # where to find solution

    plotdata.beforeframe = beforeframe_simsol
    #plotdata.afterframe = None
    plotdata.afterframe = afterframe_savefigs

    if plotdata.mode() == 'printframes':
        # call function defined below to set html and latex parameters:
        plotdata = setplot_print(plotdata)


    # limits of x-axis, set to 'auto' or set to an interval [a,b]
    xlim = 'auto'   

    # determine min and max of each component of q over all frames
    # and use this to set limits on y-axis for plots so they
    # don't change from one frame to the next:

    # set of variables that will be plotted:
    vars = [2, velocity, soundspeed, vp2nc, vm2nc, h25, h30]

    use_global_limits = False

    if use_global_limits:
        # scan all frames to determine limits:
        varmin,varmax,varlim = var_limits(plotdata,vars,padding=0.1)

        # v+2nc and v-2nc appear on same plot, and we always have v-2nc < v+2nc:
        vpm2nc_lim = [varmin[vm2nc], varmax[vp2nc]]

        # limits for powers of h:
        hpowers_lim = 'auto'
    else:
        # choose all axes limits automatically (may change each frame):
        varlim = {}
        for var in vars:
            varlim[var] = 'auto'
        vpm2nc_lim = 'auto'
        hpowers_lim = 'auto'


    plotstyle = 'o'
    color = 'b'

    def add_xshock(current_data):
        import pylab
        xshock = current_data.user.xshock
        if xshock is not None:
            pylab.xlabel("shock location = %8.2e" % xshock, fontsize=15)


    #=== Figure 1 ====
    title = 'rho'

    plotfigure = plotdata.new_plotfigure(name=title, figno=1)
    #plotfigure.kwargs = {'figsize':(12,6)}
    plotfigure.kwargs = {'figsize':(8,4)}
    plotfigure.show = True

    plotaxes = plotfigure.new_plotaxes(name=title)
    plotaxes.title = title
    plotaxes.xlimits = xlim
    plotaxes.ylimits = 'auto'
    if plot_simsol:
        plotaxes.afteraxes = plot_rho_simsol   # add plot of similarity soln

    plotitem = plotaxes.new_plotitem('rho','1d_plot')       
    plotitem.plot_var = 0
    plotitem.plotstyle = 'b-'
    plotitem.kwargs = {'linewidth': 3}

    #=== Figure 2 ====
    title = 'velocity'

    plotfigure = plotdata.new_plotfigure(name=title, figno=2)
    plotfigure.kwargs = {'figsize':(8,4)}
    plotfigure.show = True

    plotaxes = plotfigure.new_plotaxes(name=title)
    plotaxes.title = title
    plotaxes.xlimits = xlim
    plotaxes.ylimits = 'auto'

    plotitem = plotaxes.new_plotitem('velocity','1d_plot')       
    plotitem.plot_var = velocity
    plotitem.plotstyle = 'b-'
    plotitem.kwargs = {'linewidth': 3}
    if plot_simsol:
        plotaxes.afteraxes = plot_v_simsol   # add plot of similarity soln

    #plotitem = plotaxes.new_plotitem('polytropic','1d_plot')       
    #plotitem.outdir = '_output_npoly1'
    #plotitem.plot_var = velocity
    #plotitem.plotstyle = 'g-'

    #=== Figure 22 ====
    title = 'v and c'

    plotfigure = plotdata.new_plotfigure(name=title, figno=22)
    plotfigure.kwargs = {'figsize':(8,4)}
    plotfigure.show = True

    plotaxes = plotfigure.new_plotaxes(name=title)
    plotaxes.title = title
    plotaxes.xlimits = xlim
    plotaxes.ylimits = [-1., 5.]
    def add_legend(current_data):
        import pylab
        tsim = current_data.user.tsim
        pylab.legend(('sound speed','velocity'))
        pylab.title("v and c at time t = %14.8e" % tsim, fontsize=tfontsize)
    plotaxes.afteraxes = add_legend

    plotitem = plotaxes.new_plotitem('soundspeed','1d_plot')       
    plotitem.plot_var = soundspeed
    plotitem.plotstyle = 'b-'
    plotitem.kwargs = {'linewidth': 3}
    plotitem = plotaxes.new_plotitem('velocity','1d_plot')       
    plotitem.plot_var = velocity
    plotitem.plotstyle = 'k-'
    plotitem.kwargs = {'linewidth': 2}

    #=== Figure 3 ====
    title = 'internal energy per unit mass'

    plotfigure = plotdata.new_plotfigure(name=title, figno=3)
    plotfigure.kwargs = {'figsize':(8,4)}
    plotfigure.show = True

    plotaxes = plotfigure.new_plotaxes(name=title)
    plotaxes.title = title
    plotaxes.xlimits = xlim
    plotaxes.ylimits = 'auto'

    plotitem = plotaxes.new_plotitem('energy','1d_plot')       
    plotitem.plot_var = energy
    plotitem.plotstyle = 'bo-'



    #=== Figure 4 ====
    title = 'Riemann invariants'

    plotfigure = plotdata.new_plotfigure(name=title, figno=4)
    plotfigure.kwargs = {'figsize':(8,4)}
    plotfigure.show = True

    plotaxes = plotfigure.new_plotaxes(name=title)
    plotaxes.title = title
    plotaxes.xlimits = xlim
    plotaxes.ylimits = vpm2nc_lim
    if plot_simsol:
        plotaxes.afteraxes = plot_vpm2nc

    plotitem = plotaxes.new_plotitem('vp2nc','1d_plot')       
    plotitem.plot_var = vp2nc                     # function defined below
    plotitem.plotstyle = plotstyle
    plotitem.color = 'g'                          # green for all grids


    plotitem = plotaxes.new_plotitem('vm2nc','1d_plot')       
                                                    # needs distinct name
    plotitem.plot_var = vm2nc                     # function defined below
    plotitem.plotstyle = plotstyle
    plotitem.color = 'r'                          # red for all grids


    #=== Figure 5 ====
    title = 'Sound speed'

    plotfigure = plotdata.new_plotfigure(name=title, figno=5)
    plotfigure.kwargs = {'figsize':(8,4)}
    plotfigure.show = True

    plotaxes = plotfigure.new_plotaxes(name=title)
    plotaxes.title = title
    plotaxes.xlimits = xlim
    plotaxes.ylimits = 'auto'

    plotitem = plotaxes.new_plotitem('soundspeed','1d_plot')       
    plotitem.plot_var = soundspeed
    plotitem.plotstyle = 'b-'
    plotitem.kwargs = {'linewidth': 3}
    if plot_simsol:
        plotaxes.afteraxes = plot_c_simsol   # add plot of similarity soln

    #=== Figure 8 ====
    title = "v vs. y"

    plotfigure = plotdata.new_plotfigure(name='v vs. y', figno=8)
    plotfigure.kwargs = {'figsize':(8,4)}
    plotfigure.show = False

    plotaxes = plotfigure.new_plotaxes(name=title, type='empty')
    plotaxes.title = title
    # plot is done in afteraxes, plotting against y instead of x:
    if plot_simsol:
        plotaxes.afteraxes = plot_v_vs_y

    #=== Figure 9 ====
    title = "c vs. y"

    plotfigure = plotdata.new_plotfigure(name='c vs. y', figno=9)
    plotfigure.kwargs = {'figsize':(8,4)}
    plotfigure.show = True

    plotaxes = plotfigure.new_plotaxes(name=title, type='empty')
    plotaxes.title = title
    # plot is done in afteraxes, plotting against y instead of x:
    if plot_simsol:
        plotaxes.afteraxes = plot_c_vs_y

    return plotdata


#------------------------------------------------------

def afterframe_savefigs(current_data):
    """
    Save the figures needed for the paper...  may need to change the
    filenames below to be appropriate for this run.
    """

    from pylab import figure, savefig
    n = current_data.frameno
    if ieos==1:
        nlist = [30, 40, 50]  # for Euler
        fname = 'euler_n%3.1f_nstar%3.1f_' % (simsol.npoly,simsol.nstar)
    else:
        nlist = [0, 1, 2, 25, 30, 35, 40, 45, 50]  # for polytropic
        fname = 'poly_n%3.1f_' % simsol.npoly

    if n in nlist:
        figure(1)  # rho
        savefig(fname + 'rho%s.eps' % n)
        figure(2)  # velocity
        savefig(fname + 'v%s.eps' % n)
        figure(5)  # sound speed
        savefig(fname + 'c%s.eps' % n)
        if n < 10:
            figure(22)  # v and c
            savefig(fname + 'vc%s.eps' % n)
        print 'Saved figures for frame ',n
    return

#------------------------------------------------------

#  Define the functions to be plotted.

#   q is the full q array (all columns, i.e. what was called data in Matlab) 
#   on a single grid.  This is called from the plotting routine for 
#   each grid in an AMR calculation.

#  Note that indexing starts at 0 in Python so e.g.  q[:,0]= density


#------------------------------------------------------
def surface(cd):
    grid = cd.grid
    x = cd.x
    B = x  # need to fix for general polytropic!
    surface = cd.q[:,0] + B
    return surface


#------------------------------------------------------
def energy(cd):
    from numpy import where
    h = cd.q[:,0]
    hu = cd.q[:,1]
    E = cd.q[:,2]
    energy = where(h>0, (E-0.5*hu**2/h)/h, 0) 
    return energy

#------------------------------------------------------
def velocity(cd):
    from numpy import where
    h = cd.q[:,0]
    hu = cd.q[:,1]
    velocity = where(h>0, hu/h, 0)  # h = 0 ==> v = 0
    return velocity


#------------------------------------------------------

def soundspeed(cd):
    from numpy import sqrt, where
    if ieos==2:
        rho = cd.q[:,0]
        #c = sqrt(grav*rho)
        c = sqrt(K*gamma*rho**(gamma-1))
    else:
        e = energy(cd)
        c = sqrt(gamma*(gamma-1)*e)
    return c

#------------------------------------------------------

def vp2nc(cd):
    c = soundspeed(cd)
    v = velocity(cd)
    vp2nc = v + 2.*n*c
    return vp2nc

#------------------------------------------------------

def vm2nc(cd):
    c = soundspeed(cd)
    v = velocity(cd)
    vm2nc = v - 2.*n*c
    return vm2nc

#------------------------------------------------------

def vmc(cd):
    c = soundspeed(cd)
    v = velocity(cd)
    vmc = v - c
    return vmc

#------------------------------------------------------

def h25(cd):
    h25 = cd.q[:,0] ** 2.5
    return h25

#------------------------------------------------------

def h30(cd):
    h30 = cd.q[:,0] ** 3.0
    return h30


#------------------------------------------------------

def beforeframe_simsol(current_data):
    # Compute the x,c,u from the similarity solution at the current time.
    #  for adding to plots of computed soln vs. x.
    #  e.g. in plot_v_simsol.

    from numpy import hstack

    t = current_data.t  # time from this frame
    if (t==0.) or (t>1e-5):
        t = t - treset

    if ieos==2:
        (x,c,v,rho) = ss.compute_xcv_polytropic(t,simsol)
    else:
        (x,c,v,rho) = ss.compute_xcv_euler(t,simsol)

    current_data.user.xshocksim = x[0]
    print "from simsol, xshock, cshock = ",x[0],c[0]
    if hasattr(simsol,'C_rho'):
        rhos = simsol.C_rho * (-x[0])**nstar
    else:
        print "*** simsol lacks attribute C_rho"
        rhos = 0.  # need to fix
    print "+++ C_rho, x[0] = ",simsol.C_rho, x[0]
    print "+++ rhos = ",rhos
    #drho = (simsol.gamma - 1.) / (simsol.gamma + 1.)
    #drho2 = y[0]**(-nstar) * rhobar[0]
    #print "Delta rho = ", drho, drho2

    x = hstack((0., x[0], x))
    c = hstack((0., 0., c))   # not correct since c !=0 to right of shock,
                              # but ok for plotting.
    v = hstack((0., 0., v))
    rho = hstack((0., rhos, rho)) 
    current_data.user.xsim = x
    current_data.user.csim = c
    current_data.user.vsim = v
    current_data.user.rhosim = rho

    tsim = t - tfinal  # shift by treset has already been done
    current_data.user.tsim = tsim
    print "+++ t, treset, tfinal, tsim: ",t,treset,tfinal,tsim

    # compute xshock, current shock location:
    if (t>0.) and (t< 1e-5):
        xshock = np.interp(t, txs_array[:,0], txs_array[:,1])
        current_data.user.xshock = xshock
    else:
        current_data.user.xshock = None

    print "xshock = ",current_data.user.xshock

    return current_data


#------------------------------------------------------

def plot_v_vs_y(cd):
    # plot velocity vs. y instead of x
    # Similarity solution assumed to be precomputed in top of this module.

    from pylab import figure, plot,axes, hold, array, clf, isnan, title,xlim

    t = cd.t
    if ieos==2:
        ycell   = ss.mapx2y_polytropic(cd.x, t, simsol)
    else:
        ycell   = ss.mapx2y_euler(cd.x, t, simsol)
    v = velocity(cd)   # from computed results

    clf()
    plot(ycell,v,'bo')
    title("velocity vs. y at t = %14.8f " % cd.t)

    # add plot of similarity solution:
    hold(True)
    ysim = simsol.y
    vsim = cd.user.vsim
    plot(ysim, vsim, 'r-', linewidth=2)
    xlim([ycell.min(),ycell.max()])

# --------------------------------------------------

def plot_c_vs_y(cd):
    # plot soundspeed vs. y instead of x
    # Similarity solution assumed to be precomputed in top of this module.

    from pylab import figure, plot,axes, hold, array, clf, isnan, title,xlim

    t = cd.t
    if ieos==2:
        ycell   = ss.mapx2y_polytropic(cd.x, t, simsol)
    else:
        ycell   = ss.mapx2y_euler(cd.x, t, simsol)
    c = soundspeed(cd)   # from computed results

    clf()
    plot(ycell,c,'bo')
    title("soundspeed vs. y at t = %14.8f " % cd.t)

    # add plot of similarity solution:
    hold(True)
    ysim = simsol.y
    csim = cd.user.csim[2:]
    plot(ysim, csim, 'r-', linewidth=2)
    xlim([ycell.min(),ycell.max()])


# --------------------------------------------------

def plot_rho_simsol(current_data):
    # Plot rho(x,t) on top of density plot.
    # Similarity solution assumed to be precomputed in top of this module.
    # Also assumes x,u have been computed from simsol in beforeframe.

    from pylab import plot,xlim, legend,title

    # Only plot for times after t was reset:
    t = current_data.t
    tsim = current_data.user.tsim
    x = current_data.user.xsim
    rho = current_data.user.rhosim
    plot(x,rho,'r--', linewidth=2)
    xlim([current_data.xlower, 0.])

    #rhoss = simsol.C_rho * (-x)**nstar
    #plot(x,0.5*rhoss,'g--', linewidth=2)   
                
    if current_data.user.xshock is not None:
        xs = current_data.user.xshock
        xssim = current_data.user.xshocksim
        legend(('computed soln, xshock = %9.3e' % xs, \
                'similarity soln,  xshock = %9.3e' % xssim),loc='upper right')
    title("density at time t = %14.8e" % tsim, fontsize=tfontsize)



# --------------------------------------------------

def plot_v_simsol(current_data):
    # Plot v(x,t) on top of velocity plot.
    # Similarity solution assumed to be precomputed in top of this module.
    # Also assumes x,u have been computed from simsol in beforeframe.

    from pylab import plot,xlim, legend,title

    # Only plot for times after t was reset:
    t = current_data.t
    tsim = current_data.user.tsim
    x = current_data.user.xsim
    v = current_data.user.vsim
    plot(x,v,'r--', linewidth=2)
    xlim([current_data.xlower, 0.])
    if current_data.user.xshock is not None:
        xs = current_data.user.xshock
        xssim = current_data.user.xshocksim
        legend(('computed soln, xshock = %9.3e' % xs, \
                'similarity soln,  xshock = %9.3e' % xssim),loc='lower left')
    title("velocity at time t = %14.8e" % tsim, fontsize=tfontsize)



#------------------------------------------------------

def plot_c_simsol(current_data):
    # Plot c(x,t) on top of velocity plot.
    # Similarity solution assumed to be precomputed in top of this module.
    # Also assumes x,c have been computed from simsol in beforeframe.

    from pylab import plot,xlim, legend,title, axes,ylim,xticks

    # Only plot for times after t was reset:
    t = current_data.t
    tsim = current_data.user.tsim
    x = current_data.user.xsim
    c = current_data.user.csim
    plot(x,c,'r--', linewidth=2)
    xlim([current_data.xlower, 0.])
    if current_data.user.xshock is not None:
    
        xs = current_data.user.xshock
        xssim = current_data.user.xshocksim
        legend(('computed soln, xshock = %9.3e' % xs, \
                'similarity soln,  xshock = %9.3e' % xssim),loc='lower left')
    title("sound speed at time t = %14.8e" % tsim, fontsize=tfontsize)

    if 0:
        ylims = ylim()        
        # sub-plot with zoomed view:
        # set up for polytropic with n=1
        a = axes([.2,.2,.3,.3])
        c1 = current_data.var
        x1 = current_data.x
        plot(x1,c1,'b-', linewidth=2)
        plot(x,c,'r--', linewidth=2)
        xlim([1.05*xs, 0.95*xs])
        ylim([0., 0.5*ylims[1]])
        title("Zoom near shock")
        xticks([1.05*xs, xs, 0.95*xs])



#------------------------------------------------------

def plot_vpm2nc(current_data):
    # Plot v+2*n*c, v-2*n*c from similarity solution
    # Similarity solution assumed to be precomputed in top of this module.
    # Also assumes x,c,u have been computed from simsol in beforeframe.

    from pylab import plot,xlim

    # Only plot for times after t was reset:
    if current_data.t < 1e-5:
        n = simsol.npoly
        x = current_data.user.xsim
        c = current_data.user.csim
        v = current_data.user.vsim
        plot(x,v+2*n*c,'g', linewidth=2)
        plot(x,v-2*n*c,'r', linewidth=2)
        xlim([current_data.xlower, 0.])

    return


#-------------------
def setplot_print(plotdata):
#-------------------
    """
    Parameters used only when creating html and/or latex hardcopy
    e.g., via pyclaw.plotters.frametools.printframes 
    """

    plotdata.printfigs = True                # print figures
    plotdata.print_format = 'png'            # file format

    #plotdata.print_framenos = [0,10,20,50,70,80] + range(90,101)
    #plotdata.print_framenos = range(0,121,3) 
    plotdata.print_framenos = 'all'

    plotdata.print_fignos = 'all'
    plotdata.html = True                     # create html files of plots?
    plotdata.latex = True                    # create latex file of plots?
    plotdata.latex_figsperline = 2           # layout of plots
    plotdata.latex_framesperline = 1         # layout of plots
    plotdata.latex_makepdf = False           # also run pdflatex?
    return plotdata