def main():

    print(__doc__)

    x = symbols('x')

    # a numpy array we can apply the ufuncs to
    grid = np.linspace(-1, 1, 1000)

    # set mpmath precision to 20 significant numbers for verification
    mpmath.mp.dps = 20

    print("Compiling legendre ufuncs and checking results:")

    # Let's also plot the ufunc's we generate
    plot1 = Plot(visible=False)
    for n in range(6):

        # Setup the SymPy expression to ufuncify
        expr = legendre(n, x)
        print("The polynomial of degree %i is" % n)
        pprint(expr)

        # This is where the magic happens:
        binary_poly = ufuncify(x, expr)

        # It's now ready for use with numpy arrays
        polyvector = binary_poly(grid)

        # let's check the values against mpmath's legendre function
        maxdiff = 0
        for j in range(len(grid)):
            precise_val = mpmath.legendre(n, grid[j])
            diff = abs(polyvector[j] - precise_val)
            if diff > maxdiff:
                maxdiff = diff
        print("The largest error in applied ufunc was %e" % maxdiff)
        assert maxdiff < 1e-14

        # We can also attach the autowrapped legendre polynomial to a sympy
        # function and plot values as they are calculated by the binary function
        g = implemented_function('g', binary_poly)
        plot1[n] = g(x), [200]

    print("Here's a plot with values calculated by the wrapped binary functions")
    plot1.show()

def plot_beam(beam, **kwargs):
    """Plot the beam radius as it propagates in space.
    Uses pyglet and ctype libraries

    Parameters
    ==========
    
    beam : BeamParameter for gaussian beam
    z_range : plot range for the beam propagation coordinate

    See Also
    ========    
    
    BeamParameter
    plot_beam2
    
    Examples
    ========
    
    >>> from sympy.physics.gaussopt import BeamParameter, plot_beam
    >>> b = BeamParameter(530e-9, 1, w=1e-5)
    >>> plot_beam(b,z_range=2*b.z_r)    
    """
    
    if len(kwargs) != 1:
        raise ValueError("The function expects only one named argument")    
    elif 'z_range' in kwargs :
        z_range = sympify(kwargs['z_range'])
    else :
        raise ValueError(filldedent('''
            The functions expects the z_range as a named argument'''))

    x = symbols('x')
    # TODO beam.w_0 * 
    # pyglet needs to have better zoom adjustment for geting a normal view
    weist_d = sqrt(1+(x/beam.z_r)**2)
    
    p = Plot(visible=False)
    
    p[1] = weist_d, 'color=black', [x, -z_range, z_range, int(beam.w_0**-1)]
    p[2] = -weist_d, 'color=black', [x, -z_range, z_range, int(beam.w_0**-1)]
    
    p.adjust_all_bounds()
    p.show()

 def _test_plot_log(self):
     from sympy import Plot
     p=Plot(log(x), [x,0,6.282,4], 'mode=polar', visible=False)
     p.wait_for_calculations()

 def test_plot_3d_parametric(self):
     from sympy import Plot
     p=Plot(sin(x), cos(x), x/5.0, [x, 0, 6.282, 4], visible=False)
     p.wait_for_calculations()

 def test_plot_3d_spherical(self):
     from sympy import Plot
     p=Plot(1, [x,0,6.282,4], [y,0,3.141,4], 'mode=spherical;style=wireframe', visible=False)
     p.wait_for_calculations()

 def test_plot_3d_cylinder(self):
     from sympy import Plot
     p=Plot(1/y, [x,0,6.282,4], [y,-1,1,4], 'mode=polar;style=solid', visible=False)
     p.wait_for_calculations()

 def test_plot_2d_polar(self):
     from sympy import Plot
     p=Plot(1/x, [x,-1,1,4], 'mode=polar', visible=False)
     p.wait_for_calculations()

 def test_plot_3d_discontinuous(self):
     from sympy import Plot
     p=Plot(1/x, [x, -3, 3, 6], [y, -1, 1, 1], visible=False)
     p.wait_for_calculations()

 def test_plot_3d(self):
     from sympy import Plot
     p=Plot(x*y, [x, -5, 5, 5], [y, -5, 5, 5], visible=False)
     p.wait_for_calculations()

def test_plot_2d_discontinuous():
    from sympy import Plot
    p=Plot(1/x, [x, -1, 1, 2], visible=False)
    p.wait_for_calculations()

def test_plot_2d():
    from sympy import Plot
    p=Plot(x, [x, -5, 5, 4], visible=False)
    p.wait_for_calculations()

def test_plot_integral():
    # Make sure it doesn't treat x as an independent variable
    from sympy import Plot, Integral
    p = Plot(Integral(z*x, (x, 1, z), (z, 1, y)), visible=False)
    p.wait_for_calculations()