Ticket #11273: trac_11273-rebase-part1.patch

File trac_11273-rebase-part1.patch, 47.1 KB (added by kcrisman, 8 years ago)

Rebased to 5.1.beta0

  • sage/calculus/riemann.pyx

    # HG changeset patch
    # User Ethan Van Andel <evlutte@gmail.com>
    # Date 1304016273 14400
    # Node ID 06289daa915f6f1c0bc4fe054adf8cfa3603864a
    # Parent  ec6a3e8365f9c930df8da62874d84f198ad59ca4
    Major additions: Exterior Map, Ahlfors Spiderweb, Analytic Error Checking, docs
    
    diff --git a/sage/calculus/riemann.pyx b/sage/calculus/riemann.pyx
    a b  
    33
    44AUTHORS:
    55
    6 - Ethan Van Andel (2009): initial version
     6- Ethan Van Andel (2009-2011): initial version and upgrades
    77
    88- Robert Bradshaw (2009): his "complex_plot" was adapted for plot_colored
    99
     
    1111"""
    1212
    1313#*****************************************************************************
    14 #       Copyright (C) 2009 Ethan Van Andel <evlutte@gmail.com>,
     14#       Copyright (C) 2011 Ethan Van Andel <evlutte@gmail.com>,
    1515#
    1616#  Distributed under the terms of the GNU General Public License (GPL)
    1717#
     
    2828include "sage/ext/stdsage.pxi"
    2929include "sage/ext/interrupt.pxi"
    3030
    31 from sage.plot.primitive import GraphicPrimitive
    3231from sage.misc.decorators import options
    3332from sage.plot.all import list_plot, Graphics
    34 from sage.plot.misc import setup_for_eval_on_grid
    3533
    3634from sage.ext.fast_eval import fast_callable
    3735
     
    4341
    4442from sage.plot.complex_plot import ComplexPlot
    4543
     44from sage.gsl.integration import numerical_integral
     45
     46
    4647import numpy as np
    4748cimport numpy as np
    4849
    4950from math import pi
    5051from math import sin
    5152from math import cos
     53from math import sqrt
     54
     55from math import log # used for complex plot lightness
     56from math import atan
    5257
    5358from cmath import exp
    5459from cmath import phase
    5560
    56 cdef double PI = pi
    57 cdef double TWOPI = 2*PI
    58 cdef I = complex(0,1)
     61from random import uniform # used for accuracy tests
    5962
    6063FLOAT = np.float64
    6164ctypedef np.float64_t FLOAT_T
    6265
     66COMPLEX = np.complex128
     67ctypedef np.complex128_t COMPLEX_T
     68
     69cdef FLOAT_T PI = pi
     70cdef FLOAT_T TWOPI = 2*PI
     71cdef COMPLEX_T I = complex(0,1)
    6372
    6473cdef class Riemann_Map:
    6574    """
    66     The ``Riemann_Map`` class computes a Riemann or Ahlfors map from
    67     supplied data. It also has various methods to provide information about
    68     the map. A Riemann map conformally maps a simply connected region in
     75    The ``Riemann_Map`` class computes an interior or exterior Riemann map, or
     76    an Ahlfors map of a region given by the supplied boundary curve(s) and center '
     77    point. The class also provides various methods to evaluate, visualize, or extract
     78    data from the map.
     79   
     80    A Riemann map conformally maps a simply connected region in
    6981    the complex plane to the unit disc. The Ahlfors map does the same thing
    7082    for multiply connected regions.
    71 
    72     Note that all the methods are numeric rather than analytic, for unusual
    73     regions or insufficient collocation points may give very inaccurate
    74     results.
     83   
     84    Note that all the methods are numerical. As a result all answers have some
     85    imprecision. Moreover, maps computed with small number of collocation points, or
     86    for unusually shaped regions may be very inaccurate. Error computations for the
     87    ellipse can be found in the documentation for ``analytic_boundary()`` and
     88    ``analytic_interior()``.
     89   
     90    [BSV] provides an overview of the Riemann map and discusses the research
     91    that lead to the creation of this module.
    7592
    7693    INPUT:
    7794
     
    85102      Must be in the same order as ``fs``.
    86103
    87104    - ``a`` -- Complex, the center of the Riemann map. Will be mapped to
    88       the origin of the unit disc.
     105      the origin of the unit disc. Note that ``a`` MUST be within
     106      the region in order for the results to be mathematically valid.
    89107
    90108    The following inputs must all be passed in as named parameters:
    91109
    92110    - ``N`` -- integer (default: ``500``), the number of collocation points
    93111      used to compute the map. More points will give more accurate results,
    94112      especially near the boundaries, but will take longer to compute.
    95 
     113   
     114    - ``exterior`` -- boolean (default: ``False``), if set to ``True``, the
     115      exterior map will be computed, mapping the exterior of the region to the
     116      exterior of the unit circle.
     117   
     118    The following inputs may be passed as named parameters in unusual circumstances:
     119   
    96120    - ``ncorners`` -- integer (default: ``4``), if mapping a figure with
    97       (equally t-spaced) corners, better results may be obtained by
     121      (equally t-spaced) corners--corners that make a significant change in
     122      the direction of the boundary, better results may be sometimes obtained by
    98123      accurately giving this parameter. Used to add the proper constant to
    99       the theta correspondance function.
     124      the theta correspondence function.
    100125
    101126    - ``opp`` -- boolean (default: ``False``), set to ``True`` in very rare
    102       cases where the theta correspondance function is off by ``pi``, that
     127      cases where the theta correspondence function is off by ``pi``, that
    103128      is, if red is mapped left of the origin in the color plot.
     129     
    104130
    105131    EXAMPLES:
    106132
     
    142168    .. [KT] N. Kerzman and M. R. Trummer. "Numerical Conformal Mapping via
    143169      the Szego kernel". Journal of Computational and Applied Mathematics,
    144170      14(1-2): 111--123, 1986.
     171   
     172    .. [BSV] M. Bolt, S. Snoeyink, E. Van Andel. "Visual representation of
     173      the Riemann map and Ahlfors map via the Kerzman-Stein equation".
     174      Involve 3-4 (2010), 405-420.
    145175    """
    146176    cdef int N, B, ncorners
    147177    cdef f
    148178    cdef opp
    149     cdef double complex a
    150     cdef np.ndarray tk, tk2, cps, dps, szego, p_vector, pre_q_vector
     179    cdef COMPLEX_T a
     180    cdef np.ndarray tk, tk2
     181    cdef np.ndarray cps, dps, szego, p_vector, pre_q_vector
    151182    cdef np.ndarray p_vector_inverse, sinalpha, cosalpha, theta_array
    152183    cdef x_range, y_range
     184    cdef exterior
    153185
    154     def __init__(self, fs, fprimes, a, int N=500, int ncorners=4, opp=False):
     186    def __init__(self, fs, fprimes, COMPLEX_T a, int N=500, int ncorners=4, opp=False, exterior = False):
    155187        """
    156188        Initializes the ``Riemann_Map`` class. See the class ``Riemann_Map``
    157189        for full documentation on the input of this initialization method.
     
    160192
    161193            sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
    162194            sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
    163             sage: m = Riemann_Map([f], [fprime], 0)
     195            sage: m = Riemann_Map([f], [fprime], 0, N = 10)
    164196        """
    165         a = np.complex128(a)
     197        cdef double xmax, xmin, ymax, ymin, space
     198        cdef int i, k
    166199        if N <= 2:
    167200            raise ValueError(
    168201                "The number of collocation points must be > 2.")
     202        if exterior and a!=0:
     203            raise ValueError("The exterior map requires a=0")
    169204        try:
    170205            fs = [fast_callable(f, domain=CDF) for f in fs]
    171206            fprimes = [fast_callable(f, domain=CDF) for f in fprimes]
     
    176211        self.ncorners = ncorners
    177212        self.N = N  # Number of collocation pts
    178213        self.opp = opp
    179         cdef int i, k
     214        self.exterior = exterior
    180215        self.tk = np.array(np.arange(N) * TWOPI / N + 0.001 / N,
    181                            dtype=np.float64)
    182         self.tk2 = np.zeros(N + 1, dtype=np.float64)
     216                           dtype=FLOAT)
     217        self.tk2 = np.zeros(N + 1, dtype=FLOAT)
    183218        for i in xrange(N):
    184219            self.tk2[i] = self.tk[i]
    185220        self.tk2[N] = TWOPI
    186221        self.B = len(fs) # number of boundaries of the figure
    187         self.cps = np.zeros([self.B, N], dtype=np.complex128)
    188         self.dps = np.zeros([self.B, N], dtype=np.complex128)
     222        if exterior and (self.B > 1):
     223            raise ValueError(
     224                "The exterior map is undefined for multiply connected domains")
     225        cdef np.ndarray[COMPLEX_T,ndim=2] cps = np.zeros([self.B, N], dtype=COMPLEX)
     226        cdef np.ndarray[COMPLEX_T,ndim=2] dps = np.zeros([self.B, N], dtype=COMPLEX)
    189227        # Find the points on the boundaries and their derivatives.
    190         for k in xrange(self.B):
    191             for i in xrange(N):
    192                 self.cps[k, i] = np.complex(fs[k](self.tk[i]))
    193                 self.dps[k, i] = np.complex(fprimes[k](self.tk[i]))
    194         cdef double xmax = self.cps.real.max()
    195         cdef double xmin = self.cps.real.min()
    196         cdef double ymax = self.cps.imag.max()
    197         cdef double ymin = self.cps.imag.min()
    198         cdef double space = 0.1 * max(xmax - xmin, ymax - ymin)
    199         #The default plotting window, mainly for color plot.
     228        if self.exterior:
     229            for k in xrange(self.B):
     230                for i in xrange(N):
     231                    fk = fs[k](self.tk[N-i-1])
     232                    cps[k, i] = np.complex(1/fk)
     233                    dps[k, i] = np.complex(1/fk**2*fprimes[k](self.tk[N-i-1]))
     234        else:
     235            for k in xrange(self.B):
     236                for i in xrange(N):
     237                    cps[k, i] = np.complex(fs[k](self.tk[i]))
     238                    dps[k, i] = np.complex(fprimes[k](self.tk[i]))
     239        if exterior:
     240            xmax = (1/cps).real.max()
     241            xmin = (1/cps).real.min()
     242            ymax = (1/cps).imag.max()
     243            ymin = (1/cps).imag.min()
     244        else:
     245            xmax = cps.real.max()
     246            xmin = cps.real.min()
     247            ymax = cps.imag.max()
     248            ymin = cps.imag.min()
     249        space = 0.1 * max(xmax - xmin, ymax - ymin)
     250        #The default plotting window for this map.
     251        self.cps = cps
     252        self.dps = dps
    200253        self.x_range = (xmin - space, xmax + space)
    201254        self.y_range = (ymin - space, ymax + space)
     255        # Now we do some more computation
    202256        self._generate_theta_array()
    203257        self._generate_interior_mapper()
    204258        self._generate_inverse_mapper()
    205259
     260
    206261    def _repr_(self):
    207262        """
    208263        Return a string representation of this ``Riemann_Map`` object.
     
    211266           
    212267            sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
    213268            sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
    214             sage: isinstance(Riemann_Map([f], [fprime], 0)._repr_(), str)  # long time
     269            sage: isinstance(Riemann_Map([f], [fprime], 0, N = 10)._repr_(), str)  # long time
    215270            True
    216271        """
    217         return "A Riemann mapping of a figure to the unit circle."
     272        return "A Riemann or Alfohrs mapping of a figure to the unit circle."
    218273
    219274    cdef _generate_theta_array(self):
    220275        """
    221276        Generates the essential data for the Riemann map, primarily the
    222         Szego kernel and boundary correspondance.
     277        Szego kernel and boundary correspondence.  See [KT]_ for the algorithm.
    223278
    224279        TESTS::
    225280
    226281            sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
    227282            sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
    228             sage: m = Riemann_Map([f], [fprime], 0)
     283            sage: m = Riemann_Map([f], [fprime], 0, N = 10)
    229284        """
    230         cp = self.cps.flatten()
    231         dp = self.dps.flatten()
     285        cdef np.ndarray[COMPLEX_T,ndim =1] cp = self.cps.flatten()
     286        cdef np.ndarray[COMPLEX_T,ndim =1] dp = self.dps.flatten()
    232287        cdef int N = self.N
    233288        cdef int NB = N * self.B
    234289        cdef int B = self.B
    235290        cdef int i, k
    236291        cdef FLOAT_T saa, t0
    237         cdef np.ndarray[FLOAT_T, ndim=1] adp
    238         cdef np.ndarray[FLOAT_T, ndim=1] sadp
    239         cdef np.ndarray h
    240         cdef np.ndarray hconj
    241         cdef np.ndarray g
    242         cdef np.ndarray K
    243         cdef np.ndarray phi
     292        cdef np.ndarray[FLOAT_T, ndim=1] adp, sadp
     293        cdef np.ndarray[COMPLEX_T,ndim =1] h, hconj, g, normalized_dp, C, phi
     294        cdef np.ndarray[COMPLEX_T,ndim =2] K
    244295        cdef np.ndarray[FLOAT_T, ndim=2] theta_array
    245         # Seting things up to use the Nystrom method
     296        # Setting things up to use the Nystrom method
    246297        adp = abs(dp)
    247298        sadp = np.sqrt(adp)
    248299        h = 1 / (TWOPI * I) * ((dp / adp) / (self.a - cp))
    249         hconj = np.array(map(np.complex.conjugate, h), dtype=np.complex128)
     300        hconj = np.array(map(np.complex.conjugate, h), dtype=COMPLEX)
    250301        g = -sadp * hconj
    251302        normalized_dp=dp/adp
    252303        C = I / N * sadp # equivalent to -TWOPI / N * 1 / (TWOPI * I) * sadp
     
    260311        np.seterr(divide=errdivide,invalid=errinvalid) # resets the error handling
    261312        for i in xrange(NB):
    262313            K[i, i] = 1
    263         phi = np.linalg.solve(K, g) / NB * TWOPI  # Nystrom
     314        phi = np.linalg.solve(K, g) / NB * TWOPI  # Nystrom Method for solving 2nd kind integrals
    264315        # the all-important Szego kernel
    265         szego = np.array(phi.flatten() / np.sqrt(dp), dtype=np.complex128)
     316        szego = np.array(phi.flatten() / np.sqrt(dp), dtype=COMPLEX)
    266317        self.szego = szego.reshape([B, N])
    267318        start = 0
    268         # Finding the theta correspondance using phase. Misbehaves for some
     319        # Finding the theta correspondence using phase. Misbehaves for some
    269320        # regions.
    270321        if B != 1:
    271322            theta_array = np.zeros([1, NB])
     
    275326                [theta_array.reshape([B, N]), np.zeros([B, 1])], axis=1)
    276327            for k in xrange(B):
    277328                self.theta_array[k, N] = self.theta_array[k, 0] + TWOPI
    278         # Finding the theta correspondance using ab. Well behaved, but
     329        # Finding the theta correspondence using abs. Well behaved, but
    279330        # doesn't work on multiply connected domains.
    280331        else:
    281332            phi2 = phi.reshape([self.B, N])
     
    342393            sage: s = spline(points)
    343394            sage: s(3*pi / 4)
    344395            0.00121587378429...
     396            sage: plot(s,0,2*pi) # plot the kernel
    345397
    346398        The unit circle with a small hole::
    347399
     
    378430        """
    379431        Returns an array of points of the form
    380432        ``[t value, theta in e^(I*theta)]``, that is, a discretized version
    381         of the theta/boundary correspondence function. For multiply
    382         connected domains, ``get_theta_points`` will list the points for
    383         each boundary in the order that they were supplied.
     433        of the theta/boundary correspondence function. In other words, a point
     434        in this array [t1, t2] represents that the boundary point given by f(t1)
     435        is mapped to a point on the boundary of the unit circle given by e^(I*t2).
     436       
     437        For multiply connected domains, ``get_theta_points`` will list the
     438        points for each boundary in the order that they were supplied.
    384439
    385440        INPUT:
    386441
    387442        The following input must all be passed in as named parameters:
    388443
    389         - ``boundary`` -- integer (default: ``-``1) if < 0,
     444        - ``boundary`` -- integer (default: ``-1``) if < 0,
    390445          ``get_theta_points()`` will return the points for all boundaries.
    391446          If >= 0, ``get_theta_points()`` will return only the points for
    392447          the boundary specified.
     
    409464        Extending the points by a spline::
    410465
    411466            sage: s = spline(points)
    412             sage: s(3*pi / 4)
     467            sage: s(3*pi / 4) 
    413468            1.62766037996...
    414469
    415470        The unit circle with a small hole::
     
    420475            sage: hfprime(t) = 0.5*-I*e^(-I*t)
    421476            sage: m = Riemann_Map([f, hf], [hf, hfprime], 0.5 + 0.5*I)
    422477
    423         Getting the szego for a specifc boundary::
     478        Getting the boundary correspondence for a specifc boundary::
    424479
    425480            sage: tp0 = m.get_theta_points(boundary=0)
    426481            sage: tp1 = m.get_theta_points(boundary=1)
     
    437492
    438493    cdef _generate_interior_mapper(self):
    439494        """
    440         Generates the data necessary to use the ``reimann_map()`` function.
    441         As much setup as possible is done here to minimize what has to be
    442         done in ``riemann_map()``.
     495        Generates the data necessary to use the ``riemann_map`` function.
     496        As much setup as possible is done here to minimize the computation
     497        that must be done in ``riemann_map``.
    443498
    444499        TESTS::
    445500
    446501            sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
    447502            sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
    448             sage: m = Riemann_Map([f], [fprime], 0)
     503            sage: m = Riemann_Map([f], [fprime], 0, N = 10)
    449504        """
    450505        cdef int N = self.N
    451506        cdef double complex coeff = 3*I / (8*N)
     
    488543        cdef np.ndarray[double complex, ndim=1] pq = self.cps[:,list(range(N))+[0]].flatten()
    489544        self.pre_q_vector = pq
    490545
    491     cpdef riemann_map(self, pt):
     546    cpdef riemann_map(self, COMPLEX_T pt):
    492547        """
    493548        Returns the Riemann mapping of a point. That is, given ``pt`` on
    494         the interior of the mapped region, ``riemann_map()`` will return
     549        the interior of the mapped region, ``riemann_map`` will return
    495550        the point on the unit disk that ``pt`` maps to. Note that this
    496         method only works for interior points; it breaks down very close
     551        method only works for interior points; accuracy breaks down very close
    497552        to the boundary. To get boundary corrospondance, use
    498553        ``get_theta_points()``.
    499554
     
    525580            sage: m.riemann_map(np.complex(-3, 0.0001))
    526581            (1.405757...e-05+8.06...e-10j)
    527582        """
    528         pt1 = np.complex(pt)
    529         cdef np.ndarray[double complex, ndim=1] q_vector = 1 / (
    530             self.pre_q_vector - pt1)
    531         return -np.dot(self.p_vector, q_vector)
     583   
     584        cdef COMPLEX_T pt1
     585        cdef np.ndarray[COMPLEX_T, ndim=1] q_vector
     586        if self.exterior:
     587            pt1 = 1/pt
     588            q_vector = 1 / (
     589                self.pre_q_vector - pt1)
     590            return -1/np.dot(self.p_vector, q_vector)
     591        else:
     592            pt1 = pt
     593            q_vector = 1 / (
     594                self.pre_q_vector - pt1)
     595            return -np.dot(self.p_vector, q_vector)
    532596
    533597    cdef _generate_inverse_mapper(self):
    534598        """
    535599        Generates the data necessary to use the
    536         ``inverse_reimann_map()`` function. As much setup as possible is
    537         done here to minimize what has to be done in
     600        ``inverse_riemann_map()`` function. As much setup as possible is
     601        done here to minimize the computation that must be done in
    538602        ``inverse_riemann_map()``.
    539603
    540604        TESTS::
    541605
    542606            sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
    543607            sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
    544             sage: m = Riemann_Map([f], [fprime], 0)
     608            sage: m = Riemann_Map([f], [fprime], 0, N = 10)
    545609        """
    546610        cdef int N = self.N
    547611        cdef int B = self.B
     
    567631            for i in xrange(N):
    568632                self.cosalpha[k, i] = cos(-theta_array[k, i])
    569633
    570     def inverse_riemann_map(self, pt):
     634    cpdef inverse_riemann_map(self, COMPLEX_T pt):
    571635        """
    572636        Returns the inverse Riemann mapping of a point. That is, given ``pt``
    573         on the interior of the unit disc, ``inverse_reimann_map()`` will
     637        on the interior of the unit disc, ``inverse_riemann_map()`` will
    574638        return the point on the original region that would be Riemann
    575         mapped to ``pt``.
    576 
    577         .. NOTE::
    578 
    579             This method does not work for multiply connected domains.
     639        mapped to ``pt``. Note that this method does not work for multiply connected
     640        domains.
    580641
    581642        INPUT:
    582643
     
    604665            sage: m.inverse_riemann_map(np.complex(-0.2, 0.5))
    605666            (-0.156280570579...+0.321819151891...j)
    606667        """
    607         pt = np.complex128(pt)
     668        if self.exterior:
     669            pt = 1/pt
    608670        r = abs(pt)
    609671        if r == 0:
    610672            stheta = 0
     
    613675            stheta = pt.imag / r
    614676            ctheta = pt.real / r
    615677        k = 0
    616         return (1 - r**2) / TWOPI * np.dot(
     678        mapped = (1 - r**2) / TWOPI * np.dot(
    617679                self.p_vector_inverse[k] * self.cps[k],
    618680                1 / (1 + r**2 - 2*r *
    619681                     (ctheta * self.cosalpha[k] - stheta * self.sinalpha[k])))
     682        if self.exterior:
     683            return 1/mapped
     684        else:
     685            return mapped
    620686
    621687    def plot_boundaries(self, plotjoined=True, rgbcolor=[0,0,0], thickness=1):
    622688        """
     
    657723        """
    658724        plots = range(self.B)
    659725        for k in xrange(self.B):
    660             # This should be eliminated when the thickness/pointsize issue
     726            # This conditional should be eliminated when the thickness/pointsize issue
    661727            # is resolved later. Same for the others in plot_spiderweb().
    662728            if plotjoined:
    663729                plots[k] = list_plot(
     
    668734                    comp_pt(self.cps[k], 1), rgbcolor=rgbcolor,
    669735                    pointsize=thickness)
    670736        return sum(plots)
     737       
     738   
     739    cpdef compute_on_grid(self, plot_range, int x_points):
     740        """
     741        Computes the riemann map on a grid of points. Note that these points
     742        are complex of the form z = x + y*i.
     743       
     744        INPUT:
     745       
     746        - ``plot_range`` -- a tuple of the form [xmin, xmax, ymin, ymax].
     747          If the value is ``[]``, the default plotting window of the map will
     748          be used.
     749         
     750        - ``x_points`` -- int, the size of the grid in the x direction
     751          The number of points in the y_direction is scaled accordingly
     752       
     753        OUTPUT:
     754       
     755        - a tuple containing [z_values, xmin, xmax, ymin, ymax] where z_values is
     756          the evaluation of the map on the specified grid.
     757           
     758        EXAMPLES:
    671759
     760        General usage::
     761
     762            sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
     763            sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
     764            sage: m = Riemann_Map([f], [fprime], 0)
     765            sage: data = m.compute_on_grid([],5)
     766            sage: print data[0][8,1]
     767            (-0.0879...+0.9709...j)
     768        """
     769        cdef FLOAT_T xmin, xmax, xstep, ymin, ymax, ystep
     770        cdef int y_points
     771        if plot_range == []:
     772            xmin = self.x_range[0]
     773            xmax = self.x_range[1]
     774            ymin = self.y_range[0]
     775            ymax = self.y_range[1]
     776        else:
     777            xmin = plot_range[0]
     778            xmax = plot_range[1]
     779            ymin = plot_range[2]
     780            ymax = plot_range[3]
     781        xstep = (xmax - xmin) / x_points
     782        ystep = xstep
     783        y_points = int((ymax-ymin)/ ystep)
     784        cdef Py_ssize_t i, j
     785        cdef COMPLEX_T pt
     786        cdef np.ndarray[COMPLEX_T] pre_q_vector = self.pre_q_vector
     787        cdef np.ndarray[COMPLEX_T] p_vector = self.p_vector
     788        cdef np.ndarray[COMPLEX_T, ndim=2] z_values = np.empty(
     789            [y_points, x_points], dtype=np.complex128)
     790        if self.exterior:
     791            for i in xrange(x_points):
     792                for j in xrange(y_points):
     793                    pt = 1/(xmin + 0.5*xstep + i*xstep + I*(ymin + 0.5*ystep + j*ystep))
     794                    z_values[j, i] = 1/(-np.dot(p_vector,1/(pre_q_vector - pt)))
     795        else:
     796            for i in xrange(x_points):
     797                for j in xrange(y_points):
     798                    pt = xmin + 0.5*xstep + i*xstep + I*(ymin + 0.5*ystep + j*ystep)
     799                    z_values[j, i] = -np.dot(p_vector,1/(pre_q_vector - pt))
     800        return z_values, xmin, xmax, ymin, ymax
     801       
     802       
     803    @options(interpolation='catrom')
    672804    def plot_spiderweb(self, spokes=16, circles=4, pts=32, linescale=0.99,
    673                        rgbcolor=[0,0,0], thickness=1, plotjoined=True):
     805                       rgbcolor=[0,0,0], thickness=1, plotjoined=True, withcolor = False,plot_points = 200, **options):
    674806        """
    675         Generates a traditional "spiderweb plot" of the Riemann map. Shows
    676         what concentric circles and radial lines map to. Note that this
    677         method DOES NOT work for multiply connected domains.
     807        Generates a traditional "spiderweb plot" of the Riemann map. Shows
     808        what concentric circles and radial lines map to. The radial lines
     809        may exhibit erratic behavior near the boundary, if this occurs,
     810        decreasing ``linescale`` may mitigate the problem.
     811       
     812        Note that this method requires significantly more computation for
     813        multiply connected domains.
    678814
    679815        INPUT:
    680816
     
    690826          plot. Each radial line is made by ``1*pts`` points, each circle
    691827          has ``2*pts`` points. Note that high values may cause erratic
    692828          behavior of the radial lines near the boundaries.
     829          - only for simply connected domains
    693830
    694831        - ``linescale`` -- float between 0 and 1. Shrinks the radial lines
    695832          away from the boundary to reduce erratic behavior.
     833          - only for simply connected domains
    696834
    697835        - ``rgbcolor`` -- float array (default: ``[0,0,0]``) the
    698836          red-green-blue color of the spiderweb.
     
    703841        - ``plotjoined`` -- boolean (default: ``True``) If ``False``,
    704842          discrete points will be drawn; otherwise they will be connected
    705843          by lines.
     844          - only for simply connected domains
    706845
     846        - ``withcolor`` -- boolean (default: ``False``) If ``True``,
     847          The spiderweb will be overlaid on the basic color plot.
     848         
     849        - ``plot_points`` -- integer (default: ``200``) the size of the grid in the x direction
     850          The number of points in the y_direction is scaled accordingly.
     851          Note that very large values can cause this function to run slowly.
     852          - only for multiply connected domains
    707853        EXAMPLES:
    708854
    709855        General usage::
     
    732878            sage: m = Riemann_Map([f], [fprime], 0, 1000)
    733879            sage: m.plot_spiderweb()
    734880        """
    735         edge = self.plot_boundaries(plotjoined=plotjoined, rgbcolor=rgbcolor,
    736                                     thickness=thickness)
    737         circle_list = range(circles)
    738         theta_array = self.theta_array[0]
    739         s = spline(np.column_stack([self.theta_array[0], self.tk2]).tolist())
    740         tmax = self.theta_array[0, self.N]
    741         tmin = self.theta_array[0, 0]
    742881        cdef int k, i
    743         for k in xrange(circles):
    744             temp = range(pts*2)
    745             for i in xrange(2*pts):
    746                 temp[i] = self.inverse_riemann_map(
    747                     (k + 1) / (circles + 1.0) * exp(I*i * TWOPI / (2*pts)))
    748             if plotjoined:
    749                 circle_list[k] = list_plot(
    750                     comp_pt(temp, 1), rgbcolor=rgbcolor, thickness=thickness,
    751                     plotjoined=True)
     882        if self.B == 1: #The efficient simply connected
     883            edge = self.plot_boundaries(plotjoined=plotjoined, rgbcolor=rgbcolor,
     884                                        thickness=thickness)
     885            circle_list = range(circles)
     886            theta_array = self.theta_array[0]
     887            s = spline(np.column_stack([self.theta_array[0], self.tk2]).tolist())
     888            tmax = self.theta_array[0, self.N]
     889            tmin = self.theta_array[0, 0]
     890            for k in xrange(circles):
     891                temp = range(pts*2)
     892                for i in xrange(2*pts):
     893                    temp[i] = self.inverse_riemann_map(
     894                        (k + 1) / (circles + 1.0) * exp(I*i * TWOPI / (2*pts)))
     895                if plotjoined:
     896                    circle_list[k] = list_plot(
     897                        comp_pt(temp, 1), rgbcolor=rgbcolor, thickness=thickness,
     898                        plotjoined=True)
     899                else:
     900                    circle_list[k] = list_plot(
     901                        comp_pt(temp, 1), rgbcolor=rgbcolor, pointsize=thickness)
     902            line_list = range(spokes)
     903            for k in xrange(spokes):
     904                temp = range(pts)
     905                angle = (k*1.0) / spokes * TWOPI
     906                if angle >= tmax:
     907                    angle -= TWOPI
     908                elif angle <= tmin:
     909                    angle += TWOPI
     910                for i in xrange(pts - 1):
     911                    temp[i] = self.inverse_riemann_map(
     912                        (i * 1.0) / (pts * 1.0) * exp(I * angle) * linescale)
     913                temp[pts - 1] = np.complex(
     914                    self.f(s(angle)) if angle <= tmax else self.f(s(angle-TWOPI)))
     915                if plotjoined:
     916                    line_list[k] = list_plot(
     917                        comp_pt(temp, 0), rgbcolor=rgbcolor, thickness=thickness,
     918                        plotjoined=True)
     919                else:
     920                    line_list[k] = list_plot(
     921                        comp_pt(temp, 0), rgbcolor=rgbcolor, pointsize=thickness)
     922            if withcolor:
     923                return edge + sum(circle_list) + sum(line_list) + self.plot_colored(plot_points=plot_points)
    752924            else:
    753                 circle_list[k] = list_plot(
    754                     comp_pt(temp, 1), rgbcolor=rgbcolor, pointsize=thickness)
    755         line_list = range(spokes)
    756         for k in xrange(spokes):
    757             temp = range(pts)
    758             angle = (k*1.0) / spokes * TWOPI
    759             if angle >= tmax:
    760                 angle -= TWOPI
    761             elif angle <= tmin:
    762                 angle += TWOPI
    763             for i in xrange(pts - 1):
    764                 temp[i] = self.inverse_riemann_map(
    765                     (i * 1.0) / (pts * 1.0) * exp(I * angle) * linescale)
    766             temp[pts - 1] = np.complex(
    767                 self.f(s(angle)) if angle <= tmax else self.f(s(angle-TWOPI)))
    768             if plotjoined:
    769                 line_list[k] = list_plot(
    770                     comp_pt(temp, 0), rgbcolor=rgbcolor, thickness=thickness,
    771                     plotjoined=True)
    772             else:
    773                 line_list[k] = list_plot(
    774                     comp_pt(temp, 0), rgbcolor=rgbcolor, pointsize=thickness)
    775         return edge + sum(circle_list) + sum(line_list)
     925                return edge + sum(circle_list) + sum(line_list)
     926        else: # The more difficult multiply connected
     927            z_values, xmin, xmax, ymin, ymax = self.compute_on_grid([], plot_points)
     928            xstep = (xmax-xmin)/plot_points
     929            ystep = (ymax-ymin)/plot_points
     930            dr, dtheta= get_derivatives(z_values, xstep, ystep) # clean later
     931           
     932            g = Graphics()
     933            g.add_primitive(ComplexPlot(complex_to_spiderweb(z_values,dr,dtheta, spokes, circles, rgbcolor,thickness, withcolor), (xmin, xmax), (ymin, ymax),options))
     934            return g + self.plot_boundaries(thickness = thickness)
     935           
    776936
    777937    @options(interpolation='catrom')
    778938    def plot_colored(self, plot_range=[], int plot_points=100, **options):
    779939        """
    780         Draws a colored plot of the Riemann map. A red point on the
    781         colored plot corresponds to a red point on the unit disc. Note that
    782         this method DOES work for multiply connected domains.
     940        Generates a colored plot of the Riemann map. A red point on the
     941        colored plot corresponds to a red point on the unit disc.
    783942
    784943        INPUT:
    785944
     
    789948          ``(xmin, xmax, ymin, ymax)``. Declare if you do not want the plot
    790949          to use the default range for the figure.
    791950
    792         - ``plot_points`` -- integer (default: ``100``), number of points
    793           to plot in each direction of the grid. Note that very large values
    794           can cause this function to run slowly.
     951        - ``plot_points`` -- integer (default: ``100``), number of points to plot in
     952          the x direction. Points in the y direction are scaled accordingly.
     953          Note that very large values can cause this function to run slowly.
     954         
    795955
    796956        EXAMPLES:
    797957
     
    818978            sage: m = Riemann_Map([f], [fprime], 0, 1000)
    819979            sage: m.plot_colored()
    820980        """
    821         cdef int i, j
    822         cdef double xmin
    823         cdef double xmax
    824         cdef double ymin
    825         cdef double ymax
    826         if plot_range == []:
    827             xmin = self.x_range[0]
    828             xmax = self.x_range[1]
    829             ymin = self.y_range[0]
    830             ymax = self.y_range[1]
    831         else:
    832             xmin = plot_range[0]
    833             xmax = plot_range[1]
    834             ymin = plot_range[2]
    835             ymax = plot_range[3]
    836         xstep = (xmax - xmin) / plot_points
    837         ystep = (ymax - ymin) / plot_points
    838         cdef np.ndarray z_values = np.empty(
    839             [plot_points, plot_points], dtype=np.complex128)
    840         for i in xrange(plot_points):
    841             for j in xrange(plot_points):
    842                 z_values[j, i] = self.riemann_map(
    843                     np.complex(xmin + 0.5*xstep + i*xstep,
    844                                ymin + 0.5*ystep + j*ystep))
     981        z_values, xmin, xmax, ymin, ymax = self.compute_on_grid(plot_range, plot_points)
    845982        g = Graphics()
    846983        g.add_primitive(ComplexPlot(complex_to_rgb(z_values), (xmin, xmax), (ymin, ymax),options))
    847984        return g
    848985
    849986cdef comp_pt(clist, loop=True):
    850987    """
    851     This function converts a list of complex numbers to the plottable
     988    Converts a list of complex numbers xderivs = get_derivatives(z_values, xstep, ystep)[0]to the plottable
    852989    `(x,y)` form. If ``loop=True``, then the first point will be
    853990    added as the last, i.e. to plot a closed circle.
    854991
     
    8741011    if loop:
    8751012        list2[len(clist)] = list2[0]
    8761013    return list2
     1014   
     1015cpdef get_derivatives(np.ndarray[COMPLEX_T, ndim=2]z_values, FLOAT_T xstep, FLOAT_T ystep):
     1016    """
     1017    Computes the r*e^(I*theta) form derivatives from the grid of points.
     1018    The derivatives are computed using quick-and-dirty taylor expansion and assuming analyticy.
     1019    As such ``get_derivatives`` is primarily intended to be used for comparisions in
     1020    ``plot_spiderweb`` and not for applications that require great precision.
     1021   
     1022    INPUT:
     1023    - ``z_values`` -- The values for a complex function evaluated on a grid in the complex
     1024      plane, usually from ``compute_on_grid``.
     1025     
     1026    - ``xstep`` -- float, the spacing of the grid points in the real direction
     1027   
     1028    OUTPUT:
     1029   
     1030    - A tuple of arrays, [``dr``, ``dtheta``], each array is 2 less in both dimensions
     1031      than ``z_values``
     1032      ``dr`` - the absolute value of the derivative of the function in the +r direction
     1033      ``dtheta`` - the rate of accumulation of angle in the +theta direction
     1034     
     1035    EXAMPLES:
     1036       
     1037        Standard usage with compute_on_grid::
     1038            sage: from sage.calculus.riemann import get_derivatives
     1039            sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
     1040            sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
     1041            sage: m = Riemann_Map([f], [fprime], 0)
     1042            sage: data = m.compute_on_grid([],19)
     1043            sage: xstep = (data[2]-data[1])/19
     1044            sage: ystep = (data[4]-data[3])/19
     1045            sage: dr, dtheta = get_derivatives(data[0],xstep,ystep)
     1046            sage: dr[8,8]
     1047            0.241...
     1048            sage: dtheta[5,5]
     1049            5.907...
     1050     """
     1051    cdef np.ndarray[COMPLEX_T, ndim=2] xderiv
     1052    cdef np.ndarray[FLOAT_T, ndim = 2] dr, dtheta, zabs
     1053    imax = len(z_values)-2
     1054    jmax = len(z_values[0])-2
     1055    xderiv = (z_values[1:-1,2:]-z_values[1:-1,:-2])/(2*xstep) #(f(x+delta)-f(x-delta))/2delta
     1056    dr = np.abs(xderiv) #b/c the function is analytic, we know it's abs(derivative) is equal in all directions
     1057    zabs = np.abs(z_values[1:-1,1:-1]) # the abs(derivative) scaled by distance from origin
     1058    dtheta = np.divide(dr,zabs)
     1059    return dr, dtheta
    8771060
    878 cdef inline double mag_to_lightness(double r):
     1061-cdef inline double mag_to_lightness(double r):
     1062-    Tweak this to adjust how the magnitude affects the color.
     1063-    Note this method is customized for riemann plots, the
     1064-    magnitude loops rather than fading to black.
     1065
     1066     INPUT:
     1067
     1068-    - ``r`` -- a non-negative real number.
     1069
     1070     OUTPUT:
     1071 
     1072-    A value between `-1` (black) and `+1` (white), inclusive.
     1073
     1074-    EXAMPLES:
     1075
     1076 
     1077-    This tests it implicitly::
     1078-
     1079-        sage: from sage.calculus.riemann import complex_to_rgb
     1080-        sage: import numpy
     1081-        sage: complex_to_rgb(numpy.array([[0, 1, 1000]],dtype = numpy.complex128))
     1082-        array([[[   1.,    1.,    1.],
     1083-                [   1.,    0.,    0.],
     1084-                [-998.,    0.,    0.]]])
     1085     """
     1086-    return 1 - r
     1087-
     1088-cpdef complex_to_rgb(np.ndarray z_values):
     1089-    r"""
     1090-    Convert from an array of complex numbers to its corresponding matrix of
     1091
     1092cpdef complex_to_spiderweb(np.ndarray[COMPLEX_T, ndim = 2]z_values, np.ndarray[FLOAT_T, ndim = 2] dr, np.ndarray[FLOAT_T, ndim = 2] dtheta, spokes, circles, rgbcolor, thickness, withcolor):
    8791093    """
    880     Tweak this to adjust how the magnitude affects the color.
     1094    Converts a grid of complex numbers into a matrix containing rgb data
     1095    for the Riemann spiderweb plot.
     1096   
     1097     INPUT:
    8811098
    882     INPUT:
     1099    - ``z_values`` -- A grid of complex numbers, as a list of lists.
     1100   
     1101    - ``dr`` -- grid of floats, the r derivative of ``z_values``.
     1102      Used to determine precision.
     1103   
     1104    - ``dtheta`` -- grid of floats, the theta derivative of ``z_values``.
     1105      Used to determine precision.
     1106     
     1107    - ``spokes`` -- integer the number of equallyspaced radial lines to plot.
    8831108
    884     - ``r`` -- a non-negative real number.
     1109    - ``circles`` -- integer the number of equally spaced circles about the center to plot.
    8851110
    886     OUTPUT:
     1111    - ``rgbcolor`` -- float array the red-green-blue color of the lines of the spiderweb.
    8871112
    888     A value between `-1` (black) and `+1` (white), inclusive.
     1113    - ``thickness`` -- positive float the thickness of the lines or points in the spiderweb.
     1114     
     1115    - ``withcolor`` -- boolean If ``True`` the spiderweb will be overlaid on the basic color plot.
    8891116
     1117     OUTPUT:
     1118 
     1119    An `N x M x 3` floating point Numpy array ``X``, where
     1120    ``X[i,j]`` is an (r,g,b) tuple.
     1121   
    8901122    EXAMPLES:
     1123        sage: from sage.calculus.riemann import complex_to_spiderweb
     1124        sage: import numpy
     1125        sage: zval = numpy.array([[0, 1, 1000],[.2+.3j,1,-.3j],[0,0,0]],dtype = numpy.complex128)
     1126        sage: deriv = numpy.array([[.1]],dtype = numpy.float64)
     1127        sage: complex_to_spiderweb(zval, deriv,deriv, 4,4,[0,0,0],1,False)
     1128        array([[[ 1.,  1.,  1.],
     1129                [ 1.,  1.,  1.],
     1130                [ 1.,  1.,  1.]],
     1131        <BLANKLINE>
     1132               [[ 1.,  1.,  1.],
     1133                [ 0.,  0.,  0.],
     1134                [ 1.,  1.,  1.]],
     1135        <BLANKLINE>
     1136               [[ 1.,  1.,  1.],
     1137                [ 1.,  1.,  1.],
     1138                [ 1.,  1.,  1.]]])
    8911139
    892     This tests it implicitly::
     1140        sage: complex_to_spiderweb(zval, deriv,deriv, 4,4,[0,0,0],1,True)
     1141        array([[[ 1.        ,  1.        ,  1.        ],
     1142                [ 1.        ,  0.05558355,  0.05558355],
     1143                [ 0.17301243,  0.        ,  0.        ]],
     1144        <BLANKLINE>
     1145               [[ 1.        ,  0.96804683,  0.48044583],
     1146                [ 0.        ,  0.        ,  0.        ],
     1147                [ 0.77351965,  0.5470393 ,  1.        ]],
     1148        <BLANKLINE>
     1149               [[ 1.        ,  1.        ,  1.        ],
     1150                [ 1.        ,  1.        ,  1.        ],
     1151                [ 1.        ,  1.        ,  1.        ]]])
     1152     """
     1153    cdef Py_ssize_t i, j, imax, jmax
     1154    cdef FLOAT_T x, y, mag, arg, width, target, precision, dmag, darg
     1155    cdef COMPLEX_T z
     1156    cdef FLOAT_T DMAX = 70 # change to adjust rate_of_change cutoff below
     1157    cdef FLOAT_T MMIN = .0001
     1158    precision = thickness/150.0
     1159    imax = len(z_values)
     1160    jmax = len(z_values[0])
     1161    cdef np.ndarray[FLOAT_T, ndim=3, mode="c"] rgb
     1162    if withcolor:
     1163        rgb = complex_to_rgb(z_values)
     1164    else:
     1165        rgb = np.zeros(dtype=FLOAT, shape=(imax, jmax, 3))
     1166        rgb += 1
     1167    if circles != 0:
     1168        circ_radii = srange(0,1.0,1.0/circles)
     1169    else:
     1170        circ_radii = []
     1171    if spokes != 0:
     1172        spoke_angles = srange(-PI,PI+TWOPI/spokes,TWOPI/spokes) # both -pi and pi are included
     1173    else:
     1174        spoke_angles = []
     1175    for i in xrange(imax-2): # the d arrays are 1 smaller on each side
     1176        for j in xrange(jmax-2):
     1177            z = z_values[i+1,j+1]
     1178            mag = abs(z)
     1179            arg = phase(z)
     1180            dmag = dr[i,j]
     1181            darg = dtheta[i,j]
     1182            if darg < DMAX and mag > MMIN:
     1183            #points that change too rapidly are presumed to be borders
     1184            #points that are too small are presumed to be outside
     1185                for target in circ_radii:
     1186                    if abs(mag - target)/dmag < precision:
     1187                        rgb[i+1,j+1] = rgbcolor
     1188                        break
     1189                for target in spoke_angles:
     1190                    if abs(arg - target)/darg < precision:
     1191                        rgb[i+1,j+1] = rgbcolor
     1192                        break
     1193    return rgb
     1194   
    8931195
    894         sage: from sage.calculus.riemann import complex_to_rgb
    895         sage: import numpy
    896         sage: complex_to_rgb(numpy.array([[0, 1, 1000]],dtype = numpy.complex128))
    897         array([[[   1.,    1.,    1.],
    898                 [   1.,    0.,    0.],
    899                 [-998.,    0.,    0.]]])
    900     """
    901     return 1 - r
    902 
    903 cpdef complex_to_rgb(np.ndarray z_values):
     1196cpdef complex_to_rgb(np.ndarray[COMPLEX_T, ndim = 2] z_values):
    9041197    r"""
    9051198    Convert from a (Numpy) array of complex numbers to its corresponding
    9061199    matrix of RGB values.  For internal use of :meth:`~Riemann_Map.plot_colored`
     
    9201213        sage: from sage.calculus.riemann import complex_to_rgb
    9211214        sage: import numpy
    9221215        sage: complex_to_rgb(numpy.array([[0, 1, 1000]],dtype = numpy.complex128))
    923         array([[[   1.,    1.,    1.],
    924                 [   1.,    0.,    0.],
    925                 [-998.,    0.,    0.]]])
     1216        array([[[ 1.        ,  1.        ,  1.        ],
     1217                [ 1.        ,  0.05558355,  0.05558355],
     1218                [ 0.17301243,  0.        ,  0.        ]]])
     1219
    9261220        sage: complex_to_rgb(numpy.array([[0, 1j, 1000j]],dtype = numpy.complex128))
    927         array([[[  1.00000000e+00,   1.00000000e+00,   1.00000000e+00],
    928                 [  5.00000000e-01,   1.00000000e+00,   0.00000000e+00],
    929                 [ -4.99000000e+02,  -9.98000000e+02,   0.00000000e+00]]])
     1221        array([[[ 1.        ,  1.        ,  1.        ],
     1222                [ 0.52779177,  1.        ,  0.05558355],
     1223                [ 0.08650622,  0.17301243,  0.        ]]])
     1224
    9301225
    9311226    TESTS::
    9321227
     
    9351230        ...
    9361231        TypeError: Argument 'z_values' has incorrect type (expected numpy.ndarray, got list)
    9371232    """
    938     cdef unsigned int i, j, imax, jmax
    939     cdef double x, y, mag, arg
    940     cdef double lightness, hue, top, bot
    941     cdef double r, g, b
     1233    cdef Py_ssize_t i, j, imax, jmax
     1234    cdef FLOAT_T x, y, mag, arg
     1235    cdef FLOAT_T lightness, hue, top, bot
     1236    cdef FLOAT_T r, g, b
    9421237    cdef int ihue
    943     cdef z
    944 
    945     from cmath import phase
     1238    cdef COMPLEX_T z
    9461239
    9471240    imax = len(z_values)
    9481241    jmax = len(z_values[0])
    949     cdef np.ndarray[np.float_t, ndim=3, mode="c"] rgb = np.empty(
    950         dtype=np.float, shape=(imax, jmax, 3))
     1242    cdef np.ndarray[FLOAT_T, ndim=3, mode="c"] rgb = np.empty(
     1243        dtype=FLOAT, shape=(imax, jmax, 3))
    9511244
    9521245    sig_on()
    953     for i from 0 <= i < imax:
     1246    for i from 0 <= i < imax: #replace with xrange?
    9541247        row = z_values[i]
    955         for j from 0 <= j < jmax:
     1248        for j from 0 <= j < jmax: #replace with xrange?
    9561249            z = row[j]
    9571250            mag = abs(z)
    9581251            arg = phase(z)
    959             lightness = mag_to_lightness(mag)
     1252            # tweak these levels to adjust how bright/dark the colors appear
     1253            # output can range from -1 (black) to 1 (white)
     1254            lightness = -(atan(log(mag*1.5+1)) * (4/PI) - 1)
    9601255            # in hsv, variable value, full saturation (s=1, v=1+lightness)
    9611256            if lightness < 0:
    9621257                bot = 0
     
    10021297            rgb[i, j, 2] = b
    10031298    sig_off()
    10041299    return rgb
     1300   
     1301cpdef analytic_boundary(FLOAT_T t, int n):
     1302    """
     1303    Provides an exact (for n = infinity) riemann boundary
     1304    correspondence for the ellipse with axes 1 + epsilon and 1 - epsilon. The
     1305    boundary is therefore given by e^(I*t)+epsilon*e^(-I*t). It is primarily useful
     1306    for testing the accuracy of the numerical Riemann Map.
     1307   
     1308    INPUT:
     1309   
     1310    - ``t`` -- The boundary parameter, from 0 to two*pi
     1311   
     1312    - ``n`` -- Integer, The number of terms to include. 10 is fairly accurate,
     1313        20 is very accurate.
     1314       
     1315    OUTPUT:
     1316   
     1317    A theta value from 0 to 2*pi, corresponding to the point on the circle e^(I*theta)
     1318   
     1319    TESTS:
     1320   
     1321    Checking the accuracy for different n values::
     1322        sage: from sage.calculus.riemann import analytic_boundary
     1323        sage: t100 = analytic_boundary(pi/2,100)
     1324        sage: abs(analytic_boundary(pi/2,10) - t100) < 10^-8
     1325        True
     1326        sage: abs(analytic_boundary(pi/2,20) - t100) < 10^-15
     1327        True
     1328       
     1329    Using this to check the accuracy of the Riemann_Map boundary::
     1330        sage: f(t) = e^(I*t)+.3*e^(-I*t)
     1331        sage: fp(t) = I*e^(I*t)-I*.3*e^(-I*t)
     1332        sage: m = Riemann_Map([f], [fp],0,200)
     1333        sage: s = spline(m.get_theta_points())
     1334        sage: test_pt = uniform(0,2*pi)
     1335        sage: s(test_pt) - analytic_boundary(test_pt,20) < 10^-4
     1336        True
     1337    """
     1338    cdef FLOAT_T i
     1339    cdef FLOAT_T result = t
     1340    for i from 1 <= i < n+1:
     1341        result += (2*(-1)**i/i)*(.3**i/(1+.3**(2*i)))*sin(2*i*t)
     1342    return result
    10051343
     1344   
     1345
     1346cpdef cauchy_kernel(t, args):
     1347    """
     1348    Intermediate function for the integration in ``analytic_interior``.
     1349   
     1350    INPUT:
     1351   
     1352    - ``t`` -- The boundary parameter, meant to be integrated over
     1353   
     1354    - ``args`` -- a tuple containing:
     1355        - ``z`` -- Complex, the point to be mapped.
     1356       
     1357        - ``n`` -- Integer, The number of terms to include. 10 is fairly accurate,
     1358            20 is very accurate.
     1359           
     1360        - ``part`` -- will return the real ('r'), imaginary ('i') or complex ('c')
     1361        value of the kernel
     1362   
     1363    TESTS:
     1364        Primarily tested implicitly by analytic_interior,
     1365       
     1366        Simple test::
     1367        sage: from sage.calculus.riemann import cauchy_kernel
     1368        sage: cauchy_kernel(.5,(.1+.2*I, 10,'c'))
     1369        (-0.584136405997...+0.5948650858950...j)
     1370
     1371
     1372    """
     1373    cdef COMPLEX_T result
     1374    cdef COMPLEX_T z = args[0]
     1375    cdef int n = args[1]
     1376    part = args[2]
     1377    result = exp(I*analytic_boundary(t,n))/(exp(I*t)+.3*exp(-I*t)-z)*(I*exp(I*t)-I*.3*exp(-I*t))
     1378    if part == 'c':
     1379        return result
     1380    elif part == 'r':
     1381        return result.real
     1382    elif part == 'i':
     1383        return result.imag
     1384    else: return None
     1385   
     1386cpdef analytic_interior(COMPLEX_T z, int n):
     1387    """
     1388    Provides a nearly exact compuation of the Riemann Map of an interior
     1389    point of the ellipse with axes 1 + epsilon and 1 - epsilon. It is primarily useful
     1390    for testing the accuracy of the numerical Riemann Map.
     1391   
     1392    INPUT:
     1393   
     1394    - ``z`` -- Complex, the point to be mapped.
     1395   
     1396    - ``n`` -- Integer, The number of terms to include. 10 is fairly accurate,
     1397        20 is very accurate.
     1398
     1399    TESTS:
     1400        Testing the accuracy of Riemann_Map::
     1401        sage: from sage.calculus.riemann import analytic_interior
     1402        sage: f(t) = e^(I*t)+.3*e^(-I*t)
     1403        sage: fp(t) = I*e^(I*t)-I*.3*e^(-I*t)
     1404        sage: m = Riemann_Map([f],[fp],0,200)
     1405        sage: abs(m.riemann_map(.5)-analytic_interior(.5,20)) < 10^-4
     1406        True
     1407        sage: m = Riemann_Map([f],[fp],0,2000)
     1408        sage: abs(m.riemann_map(.5)-analytic_interior(.5,20)) < 10^-6
     1409        True
     1410    """
     1411    # evaluates the cauchy integral of the boundary, split into the real
     1412    # and imaginary results because numerical_integral can't handle complex data.
     1413    rp = 1/(TWOPI)*numerical_integral(cauchy_kernel,0,2*pi, params = [z,n,'i'])[0]
     1414    ip = 1/(TWOPI*I)*numerical_integral(cauchy_kernel,0,2*pi, params = [z,n,'r'])[0]
     1415    return rp + ip