def crystal_embed(e,eta):
    """
    Evaluates crystal embedding potential by solving the one-dimensional
    Schroedinger equation through one unit cell in both directions.
    """
    y0=[1.0,0.0,0.0,0.0]
    z1=np.linspace(-z_left,-z_left+alat,101)
    sol1=odeint(schroed,y0,z1,args=(e,eta))
    z2=np.linspace(-z_left+alat,-z_left,101)
    sol2=odeint(schroed,y0,z2,args=(e,eta))
    psi1=complex(sol1[50,0],sol1[50,1])
    psi1_prime=complex(sol1[50,2],sol1[50,3])
    psi2=complex(sol2[50,0],sol2[50,1])
    psi2_prime=complex(sol2[50,2],sol2[50,3])
    wronskian=psi1*psi2_prime-psi2*psi1_prime
    psi1=complex(sol1[100,0],sol1[100,1])
    psi2=complex(sol2[100,0],sol2[100,1])
    cos_ka=0.5*(psi1+psi2)
    ka=arccos(cos_ka)
    if ka.imag<0: ka=np.conj(ka)
    exp_ka=cos_ka+1.0j*sqrt(1.0-cos_ka**2)
    emb_cryst=0.5*wronskian/(exp_ka-psi2)
    if emb_cryst.imag>0.0:
        exp_ka=cos_ka-1.0j*sqrt(1.0-cos_ka**2)
        emb_cryst=0.5*wronskian/(exp_ka-psi2)
    return emb_cryst

def getEta(data,t,xi):
	global samp, ETA, time, agent, XI
	if t < len(data):
		return data[t]
	else:
		XI   = xi # update input function from paramteter
		if len(data) == 0:
			ETA0 = getInitialEta(agent.beta,agent.gamma,XI)
			data.append(ETA0[:])

		for T in range(len(data),t+1):
			# TODO: should this be samp*t so that accuracy is not lost far from 0???
			logging.info('solving ode @ t='+str(T)+', using '+str(samp)+' sub-samples')
			time = linspace(0,T,samp) #(start,end,nSamples)
			etadot_0 = [0,0,0,0,0]	#assumption of 1st order model
			#get arrays of data len=samp*t
			ETA[0] = integrate.odeint(eta1Func,[data[0][0],etadot_0[0]],time)
			ETA[1] = integrate.odeint(eta2Func,[data[0][1],etadot_0[1]],time)
			ETA[2] = integrate.odeint(eta3Func,[data[0][2],etadot_0[2]],time)
			ETA[3] = integrate.odeint(eta4Func,[data[0][3],etadot_0[3]],time)
			ETA[4] = integrate.odeint(eta5Func,[data[0][4],etadot_0[4]],time)
			logging.debug('len(result)='+str(len(ETA[0][:,0])))
			# restructure ETA using [eta#][time , eta_or_dEta] )
			E = [ETA[0][-1,0],\
			     ETA[1][-1,0],\
			     ETA[2][-1,0],\
			     ETA[3][-1,0],\
			     ETA[4][-1,0]]
			data.append(E)
		return data[t]

def ode_wrapper(fun,x0,tsp):
    # array of times at which dynamic uncertainty happens
    tdisc = np.arange(tsp[0],tsp[-1]+DT,DT)

    # determine (ad hoc) whether the called function takes 2 or 3 args
    flag_stoch = False
    try:
        dy = fun(x0,tsp[0],0.0)
        flag_stoch = True
    except TypeError:
        dy = fun(x0,tsp[0])
        flag_stoch = False

    yt = np.zeros((len(tdisc),len(x0)))
    y0 = x0.copy()
    for k in range(len(tdisc)-1):
        if flag_stoch:
            # compute the noise
            v = np.random.normal(scale=q_w)
            yp = sp.odeint(fun,y0,tdisc[k:k+2],args=(v,))
        else:
            yp = sp.odeint(fun,y0,tdisc[k:k+2])
        yt[k,:] = y0.copy()
        y0 = yp[-1,:].copy()
    yt[-1,:] = y0.copy()
    # interpolate in the new history to match the passed-in history tsp
    ysp = np.zeros((len(tsp),len(x0)))
    for k in range(len(x0)):
        ysp[:,k] = np.interp(tsp,tdisc,yt[:,k])
    return ysp

def compute_jacobian(star, differences, surface_guesses, core_guesses, core_masses, surface_masses, mass_step):
    jacobian =(np.zeros((4,4)))

    for i in range(0,4):
        guess_star = copy(star)
        if i == 0:
            step_size = core_guesses[0]*0.01
            guess_star.core_pressure = core_guesses[0] + step_size
        elif i == 1:
            step_size = core_guesses[1]*0.01
            guess_star.core_temp     = core_guesses[1] + step_size
        elif i == 2:
            step_size = surface_guesses[2]*0.01
            guess_star.total_radius  = surface_guesses[2] + step_size
        elif i == 3:
            step_size = surface_guesses[3]*0.01
            guess_star.total_lum     = surface_guesses[3] + step_size

        new_surface = inward_start(guess_star)
        new_core = outward_start(guess_star, mass_step)

        new_differences = difference_is(odeint(derivatives, new_core, core_masses, args=(guess_star,)),
                            odeint(derivatives, new_surface, surface_masses, args=(guess_star,)))
        jacobian[:,i] = np.asarray((new_differences - differences)/step_size)

    return np.linalg.inv(jacobian)

    def __init__(self, f, u0, s, t, dfdu=None):
        self.f = f
        self.t = np.array(t, float).copy()
        self.s = np.array(s, float).copy()

        if self.s.ndim == 0:
            self.s = self.s[np.newaxis]

        if dfdu is None:
            dfdu = ddu(f)
        self.dfdu = dfdu

        u0 = np.array(u0, float)
        if u0.ndim == 1:
            # run up to t[0]
            f = lambda u, t : self.f(u, s)
            assert t[0] >= 0 and t.size > 1
            N0 = int(t[0] / (t[-1] - t[0]) * t.size)
            u0 = odeint(f, u0, np.linspace(0, t[0], N0+1))[-1]

            # compute a trajectory
            self.u = odeint(f, u0, t - t[0])
        else:
            assert (u0.shape[0],) == t.shape
            self.u = u0.copy()

        self.dt = t[1:] - t[:-1]
        self.uMid = 0.5 * (self.u[1:] + self.u[:-1])
        self.dudt = (self.u[1:] - self.u[:-1]) / self.dt[:,np.newaxis]

    def Advect(self, flow_type = None, t0 = 0, dt = 0, extra_args = None):

        if flow_type is not None:
    
            # Number of particles
            num_particles = len(self.Particles);
        
            # Get the positions of all of the particles
            x0, y0, z0 = self.GetCoordinates();
        
            # Time vector
            tf = t0 + dt;
        
            # Time vector
            t = np.array([t0, tf]);
        
            # Initial positions as a numpy vector
            # in the form compatible with the 
            # velocity functions
            xy0 = np.array(x0 + y0);
            # xy0 = np.array([x0, y0]);
    
            # Choose between fields
            if "hama" in flow_type.lower():
                xy = (odeint(HamaVelocity, y0 = xy0, t = t, args = (extra_args,)))[-1, :];
                
            elif "cpipe" in flow_type.lower():
                xy = (odeint(cpipe, y0 = xy0, t = t, args = (extra_args,)))[-1, :];
                
            # Extract the new positions
            x_new, y_new = Parse_Vector_2d(xy);
            
            # Set the new coordinates
            self.SetCoordinates(x = x_new, y = y_new);

def exercise3part1(T=200, Deltat=0.1, epsilon=0.5):
    
    "In this exercise we want to perturb the parameters of the Lorenz system to see how the shape of the Lorenz attractor change"
    
    y0=np.random.randn(3)
    t = np.arange(0.,T, Deltat)
    data = integrate.odeint(vectorfield, y0, t=t,args=(10,28,8/3))
    data_sigma=integrate.odeint(vectorfield, y0, t=t,args=(10+epsilon,28,8/3))
    data_rho=integrate.odeint(vectorfield, y0, t=t,args=(10,28+epsilon,8/3))
    data_beta=integrate.odeint(vectorfield, y0, t=t,args=(10,28,8/3+epsilon))
    fig = plt.figure()
    ax0 = fig.add_subplot(2, 2, 1, projection='3d')
    ax1 = fig.add_subplot(2, 2, 2, projection='3d')
    ax2 = fig.add_subplot(2, 2, 3, projection='3d')
    ax3 = fig.add_subplot(2, 2, 4, projection='3d')
    
    ax0.plot(data[:,0],data[:,1],data[:,2])
    ax1.plot(data_sigma[:,0],data_sigma[:,1],data_sigma[:,2])
    ax2.plot(data_rho[:,0],data_rho[:,1],data_rho[:,2])
    ax3.plot(data_beta[:,0],data_beta[:,1],data_beta[:,2])

    ax0.set_title('Standard parameters')
    ax1.set_title('Sigma perturbed')
    ax2.set_title('Rho perturbed')
    ax3.set_title('Beta perturbed')

    pylab.show()
    print 'From the plot, the system seems to be more sensibile in beta perturbations'

def test_repeated_t_values():
    """Regression test for gh-8217."""

    def func(x, t):
        return -0.25*x

    t = np.zeros(10)
    sol = odeint(func, [1.], t)
    assert_array_equal(sol, np.ones((len(t), 1)))

    tau = 4*np.log(2)
    t = [0]*9 + [tau, 2*tau, 2*tau, 3*tau]
    sol = odeint(func, [1, 2], t, rtol=1e-12, atol=1e-12)
    expected_sol = np.array([[1.0, 2.0]]*9 +
                            [[0.5, 1.0],
                             [0.25, 0.5],
                             [0.25, 0.5],
                             [0.125, 0.25]])
    assert_allclose(sol, expected_sol)

    # Edge case: empty t sequence.
    sol = odeint(func, [1.], [])
    assert_array_equal(sol, np.array([], dtype=np.float64).reshape((0, 1)))

    # t values are not monotonic.
    assert_raises(ValueError, odeint, func, [1.], [0, 1, 0.5, 0])
    assert_raises(ValueError, odeint, func, [1, 2, 3], [0, -1, -2, 3])

def computeTrajectories(func, E0=np.zeros(3), **keywords):
  """Movement of electron and ion under a constant magnetic field.
  
  Positional arguments:
  func -- the name of the function computing dy/dt at time t0
  Keyword arguments:
  E0 -- Constant component of the electric field
  All other keyword arguments are collected in a 'keywords' dictionary 
  and are specific to each func."""
    
  from scipy.integrate import odeint
    
  # Processes the func specific arguments
  re0, rp0 = "ri" in keywords.keys() and keywords["ri"] or [np.zeros(3),np.zeros(3)]
  if "vi" in keywords.keys():   # Initial velocity
    v0 = keywords["vi"]
  else:
    v0 = np.array([0,1,0])
  wce, wcp = "wc" in keywords.keys() and keywords["wc"] or [0,0]

  tf = 350; NPts = 10*tf; t = np.linspace(0,tf,NPts)  # Time points

  # Integration of the equations of movement
  Q0 = np.concatenate((re0,v0))                  # Initial values
  if "wc" in keywords.keys():
    keywords["wc"] = wce
  Qe = odeint(func, Q0, t, args=(-q/me,E0,B0,keywords))  # Trajectory for the "electron"
  Q0 = np.concatenate((rp0,v0))                  # Initial values
  if "wc" in keywords.keys():
    keywords["wc"] = wcp
  Qp = odeint(func, Q0, t, args=(q/Mp,E0,B0,keywords))   # Trajectory for the "ion"
  
  return Qe, Qp

def Rhalf(Rp,Rn,a,b):
    Rn /= Rp
    c = (a*a + b*b)**0.5 / Rp
    r = linspace(0.001,1,100)

    def numer(f,r):
        lo,hi = theta_lohi(Rn,c,0,r)
        return r*r * (sin(lo)-sin(hi))

    def denom(f,r):
        lo,hi = theta_lohi(Rn,c,0,r)
        return r * (hi - lo)

    def arc(f,r):
        lo,hi = theta_lohi(Rn,c,d,r)
        return r * (hi - lo)

    up = odeint(numer,[0],r)
    dn = odeint(denom,[0],r)

    d = up[-1,0]/dn[-1,0]
    area = odeint(arc,[0],r)

    for i in range(1,len(r)):
        if area[i-1] < 0.5*area[-1] and area[i] > 0.5*area[-1]:
            return Rp * r[i]

def integrate(state, N_t=20, T=1., pts=None):
    """
    flow forward integration

    Points pts are carried along the flow without affecting it
    """
    assert(N)
    assert(DIM)
    assert(SIGMA)
    gaussian.N = N
    gaussian.DIM = DIM
    gaussian.SIGMA = SIGMA

    t_span = np.linspace(0. ,T , N_t )
    #print 'forward integration: SIGMA = ' + str(SIGMA) + ', N = ' + str(N) + ', DIM = ' + str(DIM) + ', N_t = ' + str(N_t) + ', T = ' + str(T)

    if parallel:
        odef = ode_function
    else:
        odef = ode_function_single
    
    if pts == None:
        y_span = odeint( odef , state , t_span)
        return (t_span, y_span)
    else:
        y_span = odeint( odef , np.hstack((state,pts.flatten())) , t_span)
        return (t_span, y_span[:,0:get_dim_state()], y_span[:,get_dim_state():])

def poly(n,zi,dz,filename):
    '''
    Integrates the L-E equation and returns and writes to a file the important values.
    n is the index, zi is the interval to use the polynomial solution, dz is the step size, filename is the prefix for the output file (filename+.fits)
    '''
    #set up equation and variables
    LEeq = LEeqn(n)
    zs = np.arange(0.,zi,dz)
    nz = np.size(zs)

    #integrate over initial interval using polynomial
    y = intg.odeint(thPN,np.array([1.,0.]),zs)

    #integrate until th<0
    while y[-1,0]>0.:
        zs2 = np.arange(zs[-1],zs[-1]+200.*dz,dz)
        zs = np.append(zs,zs2)
        y = np.append(y,intg.odeint(LEeq,y[-1],zs2),axis=0)
        #write progress
        sys.stdout.write("\rz={0},th={1}".format(zs[-1],y[-1,0]))
        sys.stdout.flush()
    sys.stdout.write("\n")

    #only keep values where th>0
    ngood= np.size(np.where(y[:,0]>0.))
    zs = zs[:ngood]
    y = y[:ngood]

    #prepare output array
    data = np.array([zs,y[:,0],-(zs**2)*y[:,1],(-3./zs)*y[:,1]])

    #write and return data
    fits.writeto(filename+'.fits',data)
    return data

	def delta_r(self):
		"""
		returns the coeficient delta indicating the RADIAL stability of the movement
		"""
		time = np.linspace(0, 5, 1000)
		def sm1(x, t):
			if x[0] < 0:
				return np.array([x[1], 0]) 
			else:
				return np.array([x[1], self.eta*self.Ez0(x[0])])
		traj = odeint(sm1,[0,1.],time)
		if traj[-1,0] < 0:
			ind = np.where(traj[:,0] > 0)[0][-1]
			tm = [time[ind+1],time[ind]]
			zm = [traj[ind+1,0],traj[ind,0]]
			period = np.interp(0,zm,tm)
			def sm2(x, t):
				t = t - np.floor(t/period)*period
				pos = np.interp(t, time, traj[:,0])
				psi = self.alpha + self.b**2/4. + (self.eta/2.)*self.d2V(pos)
				return np.array([x[1], psi*x[0], x[3], psi*x[2]])
			x0 = [1, 0, 0, 1]
			sol = odeint(sm2, x0, np.linspace(0, period, 1000))
			delta = np.abs((sol[-1,0] + sol[-1,3])/2)
			beta  = np.arccos((sol[-1,0] + sol[-1,3])/2)/np.pi
			if self.b == 0.:
				return delta
			else:
				out = np.log(delta + np.sqrt(np.abs(delta**2-1)))/period - self.b
				return out
		else:
			return 1.0

def runClock(dut, t_start, init_cond):
        global ts, dt, sig_jit
        global p0, z0
        global t_ab, t_ba

        alpha = random.normal(0, sig_jit)
        beta  = random.normal(0, sig_jit)
        print (alpha, beta) 
        t1 = arange(0, t_ab + alpha, dt)
        t2 = arange(0, t_ba + beta,  dt)

        concatenate((t1, [t_ab + alpha])) # 0 ~ t_ab + jitter
        concatenate((t2, [t_ba + beta ])) # 0 ~ t_ba + jitter
        t2 = map(lambda x: x + t1[-1], t2) # 0 ~ t_ab+jitter, t_ab+jitter ~ (t_ab+t_ba+2*jitter)
        
        state_half = odeint(dut, init_cond,      t1, args=(t_start, p0, z0))
        eventAtA(state_half[-1,:]) # comparator / DAC
        
        if(t_ba > dt): # when t_ba = 0 sec (NRZ DAC)
                state_full = odeint(dut, state_half[-1], t2, args=(t_start, p0, z0))
                eventAtB(state_full[-1,:]) # DAC (RZ)
                t     = concatenate((t1,         t2        ))
                state = concatenate((state_half, state_full))
        else:
                t     = t1 
                state = state_half 

        return (t, state)

def Gate(x,t):
  global pars
  xinf = ActivationCurve(V[:,1], pars[0], pars[1])
  tau = Kinetic(V[:,1], pars[2], pars[3], pars[4], pars[5])
  return (xinf-x)/tau

  # the following function is the optimisation part of the algorithm. The two chanels are modeled with initial conditions 
  # and the result is compared to the data, taking account the parameter physiological ranges for identificationdef FullTrace(p,time,y):
  pars = p[0:6]
  rinit1 = ActivationCurve(V[:,0], pars[0], pars[1]) # initial values
  r1 = odeint(Gate,rinit1,time)
  pars = p[7:13]
  rinit2 = ActivationCurve(V[:,0], pars[0], pars[1]) # initial values
  r2 = odeint(Gate,rinit2,time)
  I = p[6]*r1*(V[:,1]-E)+p[13]*r2*(V[:,1]-E)
  
  
  # this plot is to see in real time the reconstructed curves trying to fit the data (very fun)
  plot(time,y,'b'); hold(True)
  plot(time,I,'--r');hold(False)
  draw()  
  
  # this part is added to constrain the algorithm with physiological search ranges
  A = 0
  for a in arange(len(p)):
    if p[a] > HB[a] :
      A = A + (p[a] - HB[a])*1e8
    if p[a] < LB[a] :
      A = A + (LB[a]-p[a])*1e8
  return sum(array(y - I)*array(y - I))+A

    def solve_leg_transient(self):

        """Solves leg based on array of transient BC's."""

        self.delta_x = self.x[1] - self.x[0]

        self.y0 = self.T_x

        try: 
            self.T_xt

        except AttributeError:
            self.odeint_output = odeint(
                self.get_dTx_dt, y0=self.y0, t=self.t_array,
                full_output=1 
                )
            self.T_xt = self.odeint_output[0]

        else:
            self.y0 = self.T_xt[-1,:]
            self.odeint_output = odeint(
                self.get_dTx_dt, y0=self.y0, t=self.t_array,
                full_output=1 
                )
            self.T_xt = np.concatenate((self.T_xt, self.odeint_output[0]))

 def modelrun(self):
     """
     Method to integrate a specific size-based model variant.
     The integration it is done with the module "integrate.odeint",
     from the library scipy.
     """
     try:
         if self.Model == "Imm":
             outarray = odeint(self.sizemodel_imm, self.initcond, self.timedays)
             return outarray      
         elif self.Model == "TraitDif":
             outarray = odeint(self.sizemodel_traitdif, self.initcond, self.timedays)
             return outarray      
         elif self.Model == "FixVar":
             outarray = odeint(self.sizemodel_fixvar, self.initcond, self.timedays)
             return outarray  
         elif self.Model == "UnsustVar":
             outarray = odeint(self.sizemodel_unvar, self.initcond, self.timedays)
             return outarray      
         elif self.Model == "FullModel":
             outarray = odeint(self.fullmodel, self.initcond, self.timedays, rtol=1e-12, atol=1e-12)
             return outarray
         raise Exception("InputError:", "it is not a valid size model variant. Please specify one model variant "
                                        "from: Imm, TraitDif, Fixvar, UnsustVar or FullModel.\n")
     except Exception as expt:
         print expt.args[0], self.Model, expt.args[1]
         raise 

def newton_ang_func(L_val,c2,m):

   L = L_val

   slope = (m*(m+1)+c2-L)/(m+1)/2.
   z0 = [1+step*slope,slope]

   z = odeint(f,z0,x_ang,args=(c2,L,m))

   temp=1-pow(x_ang,2.0)
   temp=pow(temp,m/2.)

   zz=temp*z[:,0]

   first_zz = np.array([1])
   zz=np.append(first_zz, zz)

   sloper = -(m*(m+1)+c2-L)/(m+1)/2.
   z0r = [1-step*sloper,sloper]

   zr = odeint(f,z0r,x_angr,args=(c2,L,m))

   zzr=temp*zr[:,0]
   zzr=np.append(first_zz, zzr)
 
   return  z[:,1][-1]

def get_pressure(Planet,layers):
    radii   = Planet.get('radius')
    density = Planet.get('density')
    gravity = Planet.get('gravity')

    num_mantle_layers, num_core_layers, number_h2o_layers = layers


    # convert radii to depths
    depths       = radii[-1] - radii
    depths_core  = depths[:num_core_layers]
    gravity_core = gravity[:num_core_layers]
    density_core = density[:num_core_layers]

    depths_mant  = depths[num_core_layers:(num_core_layers+num_mantle_layers)]
    gravity_mant = gravity[num_core_layers:(num_core_layers+num_mantle_layers)]
    density_mant = density[num_core_layers:(num_core_layers+num_mantle_layers)]



    # Make a spline fit of density as a function of depth
    rhofunc_mant = interpolate.UnivariateSpline(depths_mant[::-1], density_mant[::-1])
    # Make a spline fit of gravity as a function of depth
    gfunc_mant   = interpolate.UnivariateSpline(depths_mant[::-1], gravity_mant[::-1])

    rhofunc_core = interpolate.UnivariateSpline(depths_core[::-1], density_core[::-1])
    gfunc_core   = interpolate.UnivariateSpline(depths_core[::-1], gravity_core[::-1])

    if number_h2o_layers > 0:
        depths_water  = depths[(num_core_layers + num_mantle_layers):]
        gravity_water = gravity[(num_core_layers + num_mantle_layers):]
        density_water = density[(num_core_layers + num_mantle_layers):]

        rhofunc_water = interpolate.UnivariateSpline(depths_water[::-1], density_water[::-1])
        gfunc_water   = interpolate.UnivariateSpline(depths_water[::-1], gravity_water[::-1])

        #integrate from 1 bar
        pressure_water = np.ravel(odeint((lambda p, x: gfunc_water(x) * rhofunc_water(x)),(1./10000.)*1.e9, depths_water[::-1]))
        WMB_pres       = pressure_water[-1] #does not account for very small water layers and ends up breaking code
        WMB_pres = 5.e8


    else:
        WMB_pres = 5.e8



    # integrate the hydrostatic equation
    pressure_mant = np.ravel(odeint((lambda p, x: gfunc_mant(x) * rhofunc_mant(x)),WMB_pres, depths_mant[::-1]))
    CMB_pres      = pressure_mant[-1]
    pressure_core = np.ravel(odeint((lambda p, x: gfunc_core(x) * rhofunc_core(x)),CMB_pres, depths_core[::-1]))

    if number_h2o_layers > 0:
        pressure = np.concatenate((pressure_water,pressure_mant,pressure_core),axis=0)
        pressure = [((i/1.e9)*10000.) for i in pressure]
        return np.asarray(pressure[::-1])
    else:
        pressure = np.concatenate((pressure_mant,pressure_core),axis=0)
        pressure = [((i/1.e9)*10000.) for i in pressure]
        return np.asarray(pressure[::-1])

    def draw(self,file=None):
        """
        Trace l'ensemble des lignes de champ passant par les points de départ 
        stockés dans Startpoints. Si un nom de fichier est donné, enregistre 
        la figure dans le fichier mais n'affiche rien à l'écran
        """

        def fun(P,t):
            B = self.B(P)
            Bx = B[0]
            By = B[1]
            B = np.sqrt(Bx*Bx+By*By)
            return [Bx/pow(B,4./3.),By/pow(B,4./3.)]

        t = np.linspace(0,self.k*self.maxint,self.numpoints/2)
        t2 = - t
        for P0 in self.startpoints:
            sol = odeint(fun,P0,t)
            x = sol[:,0]
            y = sol[:,1]
            pl.plot(x,y,'-',color='k')
            sol = odeint(fun,P0,t2)
            x = sol[1:,0]
            y = sol[1:,1]
            pl.plot(x,y,'-',color='k')
            pl.arrow(x[1],y[1],x[0]-x[1],y[0]-y[1],color='k')
        pl.title(self.title)
        pl.xlim([-self.size,self.size])
        pl.ylim([-self.size,self.size])
        if file:
            pl.savefig(file)
            pl.close()
        else:
            pl.show()