Ticket #11273: trac-11273-riemann-upgrade-part1.patch

File trac-11273-riemann-upgrade-part1.patch, 46.8 KB (added by evanandel, 10 years ago)
  • sage/calculus/riemann.pyx

    # HG changeset patch
    # User Ethan Van Andel <evlutte@gmail.com>
    # Date 1304016273 14400
    # Node ID 1b058b7dfc19837eb2889a3c98df3864fb9be23d
    # Parent  d305a40fdd8e8f574a6a5314b5d362265e56c207
    Major additions: Exterior Map, Ahlfors Spiderweb, Analytic Error Checking, docs
    
    diff -r d305a40fdd8e -r 1b058b7dfc19 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 "../ext/stdsage.pxi"
    2929include "../ext/interrupt.pxi"
    3030
    31 from sage.plot.primitive import GraphicPrimitive
    3231from sage.misc.decorators import options
    33 from sage.plot.plot import list_plot, Graphics, setup_for_eval_on_grid
     32from sage.plot.plot import list_plot, Graphics
    3433
    3534from sage.ext.fast_eval import fast_callable
    3635
     
    4241
    4342from sage.plot.complex_plot import ComplexPlot
    4443
     44from sage.gsl.integration import numerical_integral
     45
     46
    4547import numpy as np
    4648cimport numpy as np
    4749
    4850from math import pi
    4951from math import sin
    5052from math import cos
     53from math import sqrt
     54
     55from math import log # used for complex plot lightness
     56from math import atan
    5157
    5258from cmath import exp
    5359from cmath import phase
    5460
    55 cdef double PI = pi
    56 cdef double TWOPI = 2*PI
    57 cdef I = complex(0,1)
     61from random import uniform # used for accuracy tests
    5862
    5963FLOAT = np.float64
    6064ctypedef np.float64_t FLOAT_T
    6165
     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)
    6272
    6373cdef class Riemann_Map:
    6474    """
    65     The ``Riemann_Map`` class computes a Riemann or Ahlfors map from
    66     supplied data. It also has various methods to provide information about
    67     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
    6881    the complex plane to the unit disc. The Ahlfors map does the same thing
    6982    for multiply connected regions.
    70 
    71     Note that all the methods are numeric rather than analytic, for unusual
    72     regions or insufficient collocation points may give very inaccurate
    73     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.
    7492
    7593    INPUT:
    7694
     
    84102      Must be in the same order as ``fs``.
    85103
    86104    - ``a`` -- Complex, the center of the Riemann map. Will be mapped to
    87       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.
    88107
    89108    The following inputs must all be passed in as named parameters:
    90109
    91110    - ``N`` -- integer (default: ``500``), the number of collocation points
    92111      used to compute the map. More points will give more accurate results,
    93112      especially near the boundaries, but will take longer to compute.
    94 
     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   
    95120    - ``ncorners`` -- integer (default: ``4``), if mapping a figure with
    96       (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
    97123      accurately giving this parameter. Used to add the proper constant to
    98       the theta correspondance function.
     124      the theta correspondence function.
    99125
    100126    - ``opp`` -- boolean (default: ``False``), set to ``True`` in very rare
    101       cases where the theta correspondance function is off by ``pi``, that
     127      cases where the theta correspondence function is off by ``pi``, that
    102128      is, if red is mapped left of the origin in the color plot.
     129     
    103130
    104131    EXAMPLES:
    105132
     
    140167    .. [KT] N. Kerzman and M. R. Trummer. "Numerical Conformal Mapping via
    141168      the Szego kernel". Journal of Computational and Applied Mathematics,
    142169      14(1-2): 111--123, 1986.
     170   
     171    .. [BSV] M. Bolt, S. Snoeyink, E. Van Andel. "Visual representation of
     172      the Riemann map and Ahlfors map via the Kerzman-Stein equation".
     173      Involve 3-4 (2010), 405-420.
    143174    """
    144175    cdef int N, B, ncorners
    145176    cdef f
    146177    cdef opp
    147     cdef double complex a
    148     cdef np.ndarray tk, tk2, cps, dps, szego, p_vector, pre_q_vector
     178    cdef COMPLEX_T a
     179    cdef np.ndarray tk, tk2
     180    cdef np.ndarray cps, dps, szego, p_vector, pre_q_vector
    149181    cdef np.ndarray p_vector_inverse, sinalpha, cosalpha, theta_array
    150182    cdef x_range, y_range
     183    cdef exterior
    151184
    152     def __init__(self, fs, fprimes, a, int N=500, int ncorners=4, opp=False):
     185    def __init__(self, fs, fprimes, COMPLEX_T a, int N=500, int ncorners=4, opp=False, exterior = False):
    153186        """
    154187        Initializes the ``Riemann_Map`` class. See the class ``Riemann_Map``
    155188        for full documentation on the input of this initialization method.
     
    158191
    159192            sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
    160193            sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
    161             sage: m = Riemann_Map([f], [fprime], 0)
     194            sage: m = Riemann_Map([f], [fprime], 0, N = 10)
    162195        """
    163         a = np.complex128(a)
     196        cdef double xmax, xmin, ymax, ymin, space
     197        cdef int i, k
    164198        if N <= 2:
    165199            raise ValueError(
    166200                "The number of collocation points must be > 2.")
     201        if exterior and a!=0:
     202            raise ValueError("The exterior map requires a=0")
    167203        try:
    168204            fs = [fast_callable(f, domain=CDF) for f in fs]
    169205            fprimes = [fast_callable(f, domain=CDF) for f in fprimes]
     
    174210        self.ncorners = ncorners
    175211        self.N = N  # Number of collocation pts
    176212        self.opp = opp
    177         cdef int i, k
     213        self.exterior = exterior
    178214        self.tk = np.array(np.arange(N) * TWOPI / N + 0.001 / N,
    179                            dtype=np.float64)
    180         self.tk2 = np.zeros(N + 1, dtype=np.float64)
     215                           dtype=FLOAT)
     216        self.tk2 = np.zeros(N + 1, dtype=FLOAT)
    181217        for i in xrange(N):
    182218            self.tk2[i] = self.tk[i]
    183219        self.tk2[N] = TWOPI
    184220        self.B = len(fs) # number of boundaries of the figure
    185         self.cps = np.zeros([self.B, N], dtype=np.complex128)
    186         self.dps = np.zeros([self.B, N], dtype=np.complex128)
     221        if exterior and (self.B > 1):
     222            raise ValueError(
     223                "The exterior map is undefined for multiply connected domains")
     224        cdef np.ndarray[COMPLEX_T,ndim=2] cps = np.zeros([self.B, N], dtype=COMPLEX)
     225        cdef np.ndarray[COMPLEX_T,ndim=2] dps = np.zeros([self.B, N], dtype=COMPLEX)
    187226        # Find the points on the boundaries and their derivatives.
    188         for k in xrange(self.B):
    189             for i in xrange(N):
    190                 self.cps[k, i] = np.complex(fs[k](self.tk[i]))
    191                 self.dps[k, i] = np.complex(fprimes[k](self.tk[i]))
    192         cdef double xmax = self.cps.real.max()
    193         cdef double xmin = self.cps.real.min()
    194         cdef double ymax = self.cps.imag.max()
    195         cdef double ymin = self.cps.imag.min()
    196         cdef double space = 0.1 * max(xmax - xmin, ymax - ymin)
    197         #The default plotting window, mainly for color plot.
     227        if self.exterior:
     228            for k in xrange(self.B):
     229                for i in xrange(N):
     230                    fk = fs[k](self.tk[N-i-1])
     231                    cps[k, i] = np.complex(1/fk)
     232                    dps[k, i] = np.complex(1/fk**2*fprimes[k](self.tk[N-i-1]))
     233        else:
     234            for k in xrange(self.B):
     235                for i in xrange(N):
     236                    cps[k, i] = np.complex(fs[k](self.tk[i]))
     237                    dps[k, i] = np.complex(fprimes[k](self.tk[i]))
     238        if exterior:
     239            xmax = (1/cps).real.max()
     240            xmin = (1/cps).real.min()
     241            ymax = (1/cps).imag.max()
     242            ymin = (1/cps).imag.min()
     243        else:
     244            xmax = cps.real.max()
     245            xmin = cps.real.min()
     246            ymax = cps.imag.max()
     247            ymin = cps.imag.min()
     248        space = 0.1 * max(xmax - xmin, ymax - ymin)
     249        #The default plotting window for this map.
     250        self.cps = cps
     251        self.dps = dps
    198252        self.x_range = (xmin - space, xmax + space)
    199253        self.y_range = (ymin - space, ymax + space)
     254        # Now we do some more computation
    200255        self._generate_theta_array()
    201256        self._generate_interior_mapper()
    202257        self._generate_inverse_mapper()
    203258
     259
    204260    def _repr_(self):
    205261        """
    206262        Return a string representation of this ``Riemann_Map`` object.
     
    209265           
    210266            sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
    211267            sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
    212             sage: isinstance(Riemann_Map([f], [fprime], 0)._repr_(), str)  # long time
     268            sage: isinstance(Riemann_Map([f], [fprime], 0, N = 10)._repr_(), str)  # long time
    213269            True
    214270        """
    215         return "A Riemann mapping of a figure to the unit circle."
     271        return "A Riemann or Ahlfors mapping of a figure to the unit circle."
    216272
    217273    cdef _generate_theta_array(self):
    218274        """
    219275        Generates the essential data for the Riemann map, primarily the
    220         Szego kernel and boundary correspondance.
     276        Szego kernel and boundary correspondence. See [KT] for the algorithm.
    221277
    222278        TESTS::
    223279
    224280            sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
    225281            sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
    226             sage: m = Riemann_Map([f], [fprime], 0)
     282            sage: m = Riemann_Map([f], [fprime], 0, N = 10)
    227283        """
    228         cp = self.cps.flatten()
    229         dp = self.dps.flatten()
     284        cdef np.ndarray[COMPLEX_T,ndim =1] cp = self.cps.flatten()
     285        cdef np.ndarray[COMPLEX_T,ndim =1] dp = self.dps.flatten()
    230286        cdef int N = self.N
    231287        cdef int NB = N * self.B
    232288        cdef int B = self.B
    233289        cdef int i, k
    234290        cdef FLOAT_T saa, t0
    235         cdef np.ndarray[FLOAT_T, ndim=1] adp
    236         cdef np.ndarray[FLOAT_T, ndim=1] sadp
    237         cdef np.ndarray h
    238         cdef np.ndarray hconj
    239         cdef np.ndarray g
    240         cdef np.ndarray K
    241         cdef np.ndarray phi
     291        cdef np.ndarray[FLOAT_T, ndim=1] adp, sadp
     292        cdef np.ndarray[COMPLEX_T,ndim =1] h, hconj, g, normalized_dp,C, phi
     293        cdef np.ndarray[COMPLEX_T,ndim =2] K
    242294        cdef np.ndarray[FLOAT_T, ndim=2] theta_array
    243295        # Seting things up to use the Nystrom method
    244296        adp = abs(dp)
    245297        sadp = np.sqrt(adp)
    246298        h = 1 / (TWOPI * I) * ((dp / adp) / (self.a - cp))
    247         hconj = np.array(map(np.complex.conjugate, h), dtype=np.complex128)
     299        hconj = np.array(map(np.complex.conjugate, h), dtype=COMPLEX)
    248300        g = -sadp * hconj
    249301        normalized_dp=dp/adp
    250302        C = I / N * sadp # equivalent to -TWOPI / N * 1 / (TWOPI * I) * sadp
    251303        olderr = np.geterr()['invalid'] # checks the current error handling
    252         np.seterr(invalid='ignore')
     304        np.seterr(invalid='ignore') #there are some 0 divides that get over written later
    253305        K = np.array(
    254306            [C * sadp[t] *
    255307             (normalized_dp/(cp-cp[t]) - (normalized_dp[t]/(cp-cp[t])).conjugate())
    256               for t in np.arange(NB)], dtype=np.complex128)
     308              for t in np.arange(NB)], dtype=COMPLEX)
    257309        np.seterr(invalid=olderr) # resets the error handling
    258310        for i in xrange(NB):
    259311            K[i, i] = 1
    260         phi = np.linalg.solve(K, g) / NB * TWOPI  # Nystrom
     312        phi = np.linalg.solve(K, g) / NB * TWOPI  # Nystrom Method for solvin 2nd kind integrals
    261313        # the all-important Szego kernel
    262         szego = np.array(phi.flatten() / np.sqrt(dp), dtype=np.complex128)
     314        szego = np.array(phi.flatten() / np.sqrt(dp), dtype=COMPLEX)
    263315        self.szego = szego.reshape([B, N])
    264316        start = 0
    265         # Finding the theta correspondance using phase. Misbehaves for some
     317        # Finding the theta correspondence using phase. Misbehaves for some
    266318        # regions.
    267319        if B != 1:
    268320            theta_array = np.zeros([1, NB])
     
    272324                [theta_array.reshape([B, N]), np.zeros([B, 1])], axis=1)
    273325            for k in xrange(B):
    274326                self.theta_array[k, N] = self.theta_array[k, 0] + TWOPI
    275         # Finding the theta correspondance using ab. Well behaved, but
     327        # Finding the theta correspondence using abs. Well behaved, but
    276328        # doesn't work on multiply connected domains.
    277329        else:
    278330            phi2 = phi.reshape([self.B, N])
     
    339391            sage: s = spline(points)
    340392            sage: s(3*pi / 4)
    341393            0.00121587378429...
     394            sage: plot(s,0,2*pi) # plot the kernel
    342395
    343396        The unit circle with a small hole::
    344397
     
    375428        """
    376429        Returns an array of points of the form
    377430        ``[t value, theta in e^(I*theta)]``, that is, a discretized version
    378         of the theta/boundary correspondence function. For multiply
    379         connected domains, ``get_theta_points`` will list the points for
    380         each boundary in the order that they were supplied.
     431        of the theta/boundary correspondence function. In other words, a point
     432        in this array [t1, t2] represents that the boundary point given by f(t1)
     433        is mapped to a point on the boundary of the unit circle given by e^(I*t2).
     434       
     435        For multiply connected domains, ``get_theta_points`` will list the
     436        points for each boundary in the order that they were supplied.
    381437
    382438        INPUT:
    383439
    384440        The following input must all be passed in as named parameters:
    385441
    386         - ``boundary`` -- integer (default: ``-``1) if < 0,
     442        - ``boundary`` -- integer (default: ``-1``) if < 0,
    387443          ``get_theta_points()`` will return the points for all boundaries.
    388444          If >= 0, ``get_theta_points()`` will return only the points for
    389445          the boundary specified.
     
    406462        Extending the points by a spline::
    407463
    408464            sage: s = spline(points)
    409             sage: s(3*pi / 4)
     465            sage: s(3*pi / 4) 
    410466            1.62766037996...
    411467
    412468        The unit circle with a small hole::
     
    417473            sage: hfprime(t) = 0.5*-I*e^(-I*t)
    418474            sage: m = Riemann_Map([f, hf], [hf, hfprime], 0.5 + 0.5*I)
    419475
    420         Getting the szego for a specifc boundary::
     476        Getting the boundary correspondence for a specifc boundary::
    421477
    422478            sage: tp0 = m.get_theta_points(boundary=0)
    423479            sage: tp1 = m.get_theta_points(boundary=1)
     
    434490
    435491    cdef _generate_interior_mapper(self):
    436492        """
    437         Generates the data necessary to use the ``reimann_map()`` function.
    438         As much setup as possible is done here to minimize what has to be
    439         done in ``riemann_map()``.
     493        Generates the data necessary to use the ``riemann_map`` function.
     494        As much setup as possible is done here to minimize the computation
     495        that must be done in ``riemann_map``.
    440496
    441497        TESTS::
    442498
    443499            sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
    444500            sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
    445             sage: m = Riemann_Map([f], [fprime], 0)
     501            sage: m = Riemann_Map([f], [fprime], 0, N = 10)
    446502        """
    447503        cdef int N = self.N
    448504        cdef double complex coeff = 3*I / (8*N)
     
    485541        cdef np.ndarray[double complex, ndim=1] pq = self.cps[:,list(range(N))+[0]].flatten()
    486542        self.pre_q_vector = pq
    487543
    488     cpdef riemann_map(self, pt):
     544    cpdef riemann_map(self, COMPLEX_T pt):
    489545        """
    490546        Returns the Riemann mapping of a point. That is, given ``pt`` on
    491         the interior of the mapped region, ``riemann_map()`` will return
     547        the interior of the mapped region, ``riemann_map`` will return
    492548        the point on the unit disk that ``pt`` maps to. Note that this
    493         method only works for interior points; it breaks down very close
     549        method only works for interior points; accuracy breaks down very close
    494550        to the boundary. To get boundary corrospondance, use
    495551        ``get_theta_points()``.
    496552
     
    522578            sage: m.riemann_map(np.complex(-3, 0.0001))
    523579            (1.405757...e-05+8.06...e-10j)
    524580        """
    525         pt1 = np.complex(pt)
    526         cdef np.ndarray[double complex, ndim=1] q_vector = 1 / (
    527             self.pre_q_vector - pt1)
    528         return -np.dot(self.p_vector, q_vector)
     581   
     582        cdef COMPLEX_T pt1
     583        cdef np.ndarray[COMPLEX_T, ndim=1] q_vector
     584        if self.exterior:
     585            pt1 = 1/pt
     586            q_vector = 1 / (
     587                self.pre_q_vector - pt1)
     588            return -1/np.dot(self.p_vector, q_vector)
     589        else:
     590            pt1 = pt
     591            q_vector = 1 / (
     592                self.pre_q_vector - pt1)
     593            return -np.dot(self.p_vector, q_vector)
    529594
    530595    cdef _generate_inverse_mapper(self):
    531596        """
    532597        Generates the data necessary to use the
    533         ``inverse_reimann_map()`` function. As much setup as possible is
    534         done here to minimize what has to be done in
     598        ``inverse_riemann_map()`` function. As much setup as possible is
     599        done here to minimize the computation that must be done in
    535600        ``inverse_riemann_map()``.
    536601
    537602        TESTS::
    538603
    539604            sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
    540605            sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
    541             sage: m = Riemann_Map([f], [fprime], 0)
     606            sage: m = Riemann_Map([f], [fprime], 0, N = 10)
    542607        """
    543608        cdef int N = self.N
    544609        cdef int B = self.B
     
    564629            for i in xrange(N):
    565630                self.cosalpha[k, i] = cos(-theta_array[k, i])
    566631
    567     def inverse_riemann_map(self, pt):
     632    cpdef inverse_riemann_map(self, COMPLEX_T pt):
    568633        """
    569634        Returns the inverse Riemann mapping of a point. That is, given ``pt``
    570         on the interior of the unit disc, ``inverse_reimann_map()`` will
     635        on the interior of the unit disc, ``inverse_riemann_map()`` will
    571636        return the point on the original region that would be Riemann
    572         mapped to ``pt``.
    573 
    574         .. NOTE::
    575 
    576             This method does not work for multiply connected domains.
     637        mapped to ``pt``. Note that this method does not work for multiply connected
     638        domains.
    577639
    578640        INPUT:
    579641
     
    601663            sage: m.inverse_riemann_map(np.complex(-0.2, 0.5))
    602664            (-0.156280570579...+0.321819151891...j)
    603665        """
    604         pt = np.complex128(pt)
     666        if self.exterior:
     667            pt = 1/pt
    605668        r = abs(pt)
    606669        if r == 0:
    607670            stheta = 0
     
    610673            stheta = pt.imag / r
    611674            ctheta = pt.real / r
    612675        k = 0
    613         return (1 - r**2) / TWOPI * np.dot(
     676        mapped = (1 - r**2) / TWOPI * np.dot(
    614677                self.p_vector_inverse[k] * self.cps[k],
    615678                1 / (1 + r**2 - 2*r *
    616679                     (ctheta * self.cosalpha[k] - stheta * self.sinalpha[k])))
     680        if self.exterior:
     681            return 1/mapped
     682        else:
     683            return mapped
    617684
    618685    def plot_boundaries(self, plotjoined=True, rgbcolor=[0,0,0], thickness=1):
    619686        """
     
    654721        """
    655722        plots = range(self.B)
    656723        for k in xrange(self.B):
    657             # This should be eliminated when the thickness/pointsize issue
     724            # This conditional should be eliminated when the thickness/pointsize issue
    658725            # is resolved later. Same for the others in plot_spiderweb().
    659726            if plotjoined:
    660727                plots[k] = list_plot(
     
    665732                    comp_pt(self.cps[k], 1), rgbcolor=rgbcolor,
    666733                    pointsize=thickness)
    667734        return sum(plots)
     735       
     736   
     737    cpdef compute_on_grid(self, plot_range, int x_points):
     738        """
     739        Computes the riemann map on a grid of points. Note that these points
     740        are complex of the form z = x + y*i.
     741       
     742        INPUT:
     743       
     744        - ``plot_range`` -- a tuple of the form [xmin, xmax, ymin, ymax].
     745          If the value is ``[]``, the default plotting window of the map will
     746          be used.
     747         
     748        - ``x_points`` -- int, the size of the grid in the x direction
     749          The number of points in the y_direction is scaled accordingly
     750       
     751        OUTPUT:
     752       
     753        - a tuple containing [z_values, xmin, xmax, ymin, ymax] where z_values is
     754          the evaluation of the map on the specified grid.
     755           
     756        EXAMPLES:
    668757
     758        General usage::
     759
     760            sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
     761            sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
     762            sage: m = Riemann_Map([f], [fprime], 0)
     763            sage: data = m.compute_on_grid([],5)
     764            sage: print data[0][8,1]
     765            (-0.0879...+0.9709...j)
     766        """
     767        cdef FLOAT_T xmin, xmax, xstep, ymin, ymax, ystep
     768        cdef int y_points
     769        if plot_range == []:
     770            xmin = self.x_range[0]
     771            xmax = self.x_range[1]
     772            ymin = self.y_range[0]
     773            ymax = self.y_range[1]
     774        else:
     775            xmin = plot_range[0]
     776            xmax = plot_range[1]
     777            ymin = plot_range[2]
     778            ymax = plot_range[3]
     779        xstep = (xmax - xmin) / x_points
     780        ystep = xstep
     781        y_points = int((ymax-ymin)/ ystep)
     782        cdef Py_ssize_t i, j
     783        cdef COMPLEX_T pt
     784        cdef np.ndarray[COMPLEX_T] pre_q_vector = self.pre_q_vector
     785        cdef np.ndarray[COMPLEX_T] p_vector = self.p_vector
     786        cdef np.ndarray[COMPLEX_T, ndim=2] z_values = np.empty(
     787            [y_points, x_points], dtype=np.complex128)
     788        if self.exterior:
     789            for i in xrange(x_points):
     790                for j in xrange(y_points):
     791                    pt = 1/(xmin + 0.5*xstep + i*xstep + I*(ymin + 0.5*ystep + j*ystep))
     792                    z_values[j, i] = 1/(-np.dot(p_vector,1/(pre_q_vector - pt)))
     793        else:
     794            for i in xrange(x_points):
     795                for j in xrange(y_points):
     796                    pt = xmin + 0.5*xstep + i*xstep + I*(ymin + 0.5*ystep + j*ystep)
     797                    z_values[j, i] = -np.dot(p_vector,1/(pre_q_vector - pt))
     798        return z_values, xmin, xmax, ymin, ymax
     799       
     800       
     801    @options(interpolation='catrom')
    669802    def plot_spiderweb(self, spokes=16, circles=4, pts=32, linescale=0.99,
    670                        rgbcolor=[0,0,0], thickness=1, plotjoined=True):
     803                       rgbcolor=[0,0,0], thickness=1, plotjoined=True, withcolor = False,plot_points = 200, **options):
    671804        """
    672         Generates a traditional "spiderweb plot" of the Riemann map. Shows
    673         what concentric circles and radial lines map to. Note that this
    674         method DOES NOT work for multiply connected domains.
     805        Generates a traditional "spiderweb plot" of the Riemann map. Shows
     806        what concentric circles and radial lines map to. The radial lines
     807        may exhibit erratic behavior near the boundary, if this occurs,
     808        decreasing ``linescale`` may mitigate the problem.
     809       
     810        Note that this method requires significantly more computation for
     811        multiply connected domains.
    675812
    676813        INPUT:
    677814
     
    687824          plot. Each radial line is made by ``1*pts`` points, each circle
    688825          has ``2*pts`` points. Note that high values may cause erratic
    689826          behavior of the radial lines near the boundaries.
     827          - only for simply connected domains
    690828
    691829        - ``linescale`` -- float between 0 and 1. Shrinks the radial lines
    692830          away from the boundary to reduce erratic behavior.
     831          - only for simply connected domains
    693832
    694833        - ``rgbcolor`` -- float array (default: ``[0,0,0]``) the
    695834          red-green-blue color of the spiderweb.
     
    700839        - ``plotjoined`` -- boolean (default: ``True``) If ``False``,
    701840          discrete points will be drawn; otherwise they will be connected
    702841          by lines.
     842          - only for simply connected domains
    703843
     844        - ``withcolor`` -- boolean (default: ``False``) If ``True``,
     845          The spiderweb will be overlaid on the basic color plot.
     846         
     847        - ``plot_points`` -- integer (default: ``200``) the size of the grid in the x direction
     848          The number of points in the y_direction is scaled accordingly.
     849          Note that very large values can cause this function to run slowly.
     850          - only for multiply connected domains
    704851        EXAMPLES:
    705852
    706853        General usage::
     
    729876            sage: m = Riemann_Map([f], [fprime], 0, 1000)
    730877            sage: m.plot_spiderweb()
    731878        """
    732         edge = self.plot_boundaries(plotjoined=plotjoined, rgbcolor=rgbcolor,
    733                                     thickness=thickness)
    734         circle_list = range(circles)
    735         theta_array = self.theta_array[0]
    736         s = spline(np.column_stack([self.theta_array[0], self.tk2]).tolist())
    737         tmax = self.theta_array[0, self.N]
    738         tmin = self.theta_array[0, 0]
    739879        cdef int k, i
    740         for k in xrange(circles):
    741             temp = range(pts*2)
    742             for i in xrange(2*pts):
    743                 temp[i] = self.inverse_riemann_map(
    744                     (k + 1) / (circles + 1.0) * exp(I*i * TWOPI / (2*pts)))
    745             if plotjoined:
    746                 circle_list[k] = list_plot(
    747                     comp_pt(temp, 1), rgbcolor=rgbcolor, thickness=thickness,
    748                     plotjoined=True)
     880        if self.B == 1: #The efficient simply connected
     881            edge = self.plot_boundaries(plotjoined=plotjoined, rgbcolor=rgbcolor,
     882                                        thickness=thickness)
     883            circle_list = range(circles)
     884            theta_array = self.theta_array[0]
     885            s = spline(np.column_stack([self.theta_array[0], self.tk2]).tolist())
     886            tmax = self.theta_array[0, self.N]
     887            tmin = self.theta_array[0, 0]
     888            for k in xrange(circles):
     889                temp = range(pts*2)
     890                for i in xrange(2*pts):
     891                    temp[i] = self.inverse_riemann_map(
     892                        (k + 1) / (circles + 1.0) * exp(I*i * TWOPI / (2*pts)))
     893                if plotjoined:
     894                    circle_list[k] = list_plot(
     895                        comp_pt(temp, 1), rgbcolor=rgbcolor, thickness=thickness,
     896                        plotjoined=True)
     897                else:
     898                    circle_list[k] = list_plot(
     899                        comp_pt(temp, 1), rgbcolor=rgbcolor, pointsize=thickness)
     900            line_list = range(spokes)
     901            for k in xrange(spokes):
     902                temp = range(pts)
     903                angle = (k*1.0) / spokes * TWOPI
     904                if angle >= tmax:
     905                    angle -= TWOPI
     906                elif angle <= tmin:
     907                    angle += TWOPI
     908                for i in xrange(pts - 1):
     909                    temp[i] = self.inverse_riemann_map(
     910                        (i * 1.0) / (pts * 1.0) * exp(I * angle) * linescale)
     911                temp[pts - 1] = np.complex(
     912                    self.f(s(angle)) if angle <= tmax else self.f(s(angle-TWOPI)))
     913                if plotjoined:
     914                    line_list[k] = list_plot(
     915                        comp_pt(temp, 0), rgbcolor=rgbcolor, thickness=thickness,
     916                        plotjoined=True)
     917                else:
     918                    line_list[k] = list_plot(
     919                        comp_pt(temp, 0), rgbcolor=rgbcolor, pointsize=thickness)
     920            if withcolor:
     921                return edge + sum(circle_list) + sum(line_list) + self.plot_colored(plot_points=plot_points)
    749922            else:
    750                 circle_list[k] = list_plot(
    751                     comp_pt(temp, 1), rgbcolor=rgbcolor, pointsize=thickness)
    752         line_list = range(spokes)
    753         for k in xrange(spokes):
    754             temp = range(pts)
    755             angle = (k*1.0) / spokes * TWOPI
    756             if angle >= tmax:
    757                 angle -= TWOPI
    758             elif angle <= tmin:
    759                 angle += TWOPI
    760             for i in xrange(pts - 1):
    761                 temp[i] = self.inverse_riemann_map(
    762                     (i * 1.0) / (pts * 1.0) * exp(I * angle) * linescale)
    763             temp[pts - 1] = np.complex(
    764                 self.f(s(angle)) if angle <= tmax else self.f(s(angle-TWOPI)))
    765             if plotjoined:
    766                 line_list[k] = list_plot(
    767                     comp_pt(temp, 0), rgbcolor=rgbcolor, thickness=thickness,
    768                     plotjoined=True)
    769             else:
    770                 line_list[k] = list_plot(
    771                     comp_pt(temp, 0), rgbcolor=rgbcolor, pointsize=thickness)
    772         return edge + sum(circle_list) + sum(line_list)
     923                return edge + sum(circle_list) + sum(line_list)
     924        else: # The more difficult multiply connected
     925            z_values, xmin, xmax, ymin, ymax = self.compute_on_grid([], plot_points)
     926            xstep = (xmax-xmin)/plot_points
     927            ystep = (ymax-ymin)/plot_points
     928            dr, dtheta= get_derivatives(z_values, xstep, ystep) # clean later
     929           
     930            g = Graphics()
     931            g.add_primitive(ComplexPlot(complex_to_spiderweb(z_values,dr,dtheta, spokes, circles, rgbcolor,thickness, withcolor), (xmin, xmax), (ymin, ymax),options))
     932            return g + self.plot_boundaries(thickness = thickness)
     933           
    773934
    774935    @options(interpolation='catrom')
    775936    def plot_colored(self, plot_range=[], int plot_points=100, **options):
    776937        """
    777         Draws a colored plot of the Riemann map. A red point on the
    778         colored plot corresponds to a red point on the unit disc. Note that
    779         this method DOES work for multiply connected domains.
     938        Generates a colored plot of the Riemann map. A red point on the
     939        colored plot corresponds to a red point on the unit disc.
    780940
    781941        INPUT:
    782942
     
    786946          ``(xmin, xmax, ymin, ymax)``. Declare if you do not want the plot
    787947          to use the default range for the figure.
    788948
    789         - ``plot_points`` -- integer (default: ``100``), number of points
    790           to plot in each direction of the grid. Note that very large values
    791           can cause this function to run slowly.
     949        - ``plot_points`` -- integer (default: ``100``), number of points to plot in
     950          the x direction. Points in the y direction are scaled accordingly.
     951          Note that very large values can cause this function to run slowly.
     952         
    792953
    793954        EXAMPLES:
    794955
     
    815976            sage: m = Riemann_Map([f], [fprime], 0, 1000)
    816977            sage: m.plot_colored()
    817978        """
    818         cdef int i, j
    819         cdef double xmin
    820         cdef double xmax
    821         cdef double ymin
    822         cdef double ymax
    823         if plot_range == []:
    824             xmin = self.x_range[0]
    825             xmax = self.x_range[1]
    826             ymin = self.y_range[0]
    827             ymax = self.y_range[1]
    828         else:
    829             xmin = plot_range[0]
    830             xmax = plot_range[1]
    831             ymin = plot_range[2]
    832             ymax = plot_range[3]
    833         xstep = (xmax - xmin) / plot_points
    834         ystep = (ymax - ymin) / plot_points
    835         cdef np.ndarray z_values = np.empty(
    836             [plot_points, plot_points], dtype=np.complex128)
    837         for i in xrange(plot_points):
    838             for j in xrange(plot_points):
    839                 z_values[j, i] = self.riemann_map(
    840                     np.complex(xmin + 0.5*xstep + i*xstep,
    841                                ymin + 0.5*ystep + j*ystep))
     979        z_values, xmin, xmax, ymin, ymax = self.compute_on_grid(plot_range, plot_points)
    842980        g = Graphics()
    843981        g.add_primitive(ComplexPlot(complex_to_rgb(z_values), (xmin, xmax), (ymin, ymax),options))
    844982        return g
    845983
    846984cdef comp_pt(clist, loop=True):
    847985    """
    848     This function converts a list of complex numbers to the plottable
     986    Converts a list of complex numbers xderivs = get_derivatives(z_values, xstep, ystep)[0]to the plottable
    849987    `(x,y)` form. If ``loop=True``, then the first point will be
    850988    added as the last, i.e. to plot a closed circle.
    851989
     
    8711009    if loop:
    8721010        list2[len(clist)] = list2[0]
    8731011    return list2
     1012   
     1013cpdef get_derivatives(np.ndarray[COMPLEX_T, ndim=2]z_values, FLOAT_T xstep, FLOAT_T ystep):
     1014    """
     1015    Computes the r*e^(I*theta) form derivatives from the grid of points.
     1016    The derivatives are computed using quick-and-dirty taylor expansion and assuming analyticy.
     1017    As such ``get_derivatives`` is primarily intended to be used for comparisions in
     1018    ``plot_spiderweb`` and not for applications that require great precision.
     1019   
     1020    INPUT:
     1021    - ``z_values`` -- The values for a complex function evaluated on a grid in the complex
     1022      plane, usually from ``compute_on_grid``.
     1023     
     1024    - ``xstep`` -- float, the spacing of the grid points in the real direction
     1025   
     1026    OUTPUT:
     1027   
     1028    - A tuple of arrays, [``dr``, ``dtheta``], each array is 2 less in both dimensions
     1029      than ``z_values``
     1030      ``dr`` - the absolute value of the derivative of the function in the +r direction
     1031      ``dtheta`` - the rate of accumulation of angle in the +theta direction
     1032     
     1033    EXAMPLES:
     1034       
     1035        Standard usage with compute_on_grid::
     1036            sage: from sage.calculus.riemann import get_derivatives
     1037            sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
     1038            sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
     1039            sage: m = Riemann_Map([f], [fprime], 0)
     1040            sage: data = m.compute_on_grid([],19)
     1041            sage: xstep = (data[2]-data[1])/19
     1042            sage: ystep = (data[4]-data[3])/19
     1043            sage: dr, dtheta = get_derivatives(data[0],xstep,ystep)
     1044            sage: dr[8,8]
     1045            0.241...
     1046            sage: dtheta[5,5]
     1047            5.907...
    8741048
    875 cdef inline double mag_to_lightness(double r):
    8761049    """
    877     Tweak this to adjust how the magnitude affects the color.
    878     Note this method is customized for riemann plots, the
    879     magnitude loops rather than fading to black.
     1050    cdef np.ndarray[COMPLEX_T, ndim=2] xderiv
     1051    cdef np.ndarray[FLOAT_T, ndim = 2] dr, dtheta, zabs
     1052    imax = len(z_values)-2
     1053    jmax = len(z_values[0])-2
     1054    xderiv = (z_values[1:-1,2:]-z_values[1:-1,:-2])/(2*xstep) #(f(x+delta)-f(x-delta))/2delta
     1055    dr = np.abs(xderiv) #b/c the function is analytic, we know it's abs(derivative) is equal in all directions
     1056    zabs = np.abs(z_values[1:-1,1:-1]) # the abs(derivative) scaled by distance from origin
     1057    dtheta = np.divide(dr,zabs)
     1058    return dr, dtheta
    8801059
     1060cpdef 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):
     1061    """
     1062    Converts a grid of complex numbers into a matrix containing rgb data
     1063    for the Riemann spiderweb plot.
     1064   
    8811065    INPUT:
    8821066
    883     - ``r`` -- a non-negative real number.
     1067    - ``z_values`` -- A grid of complex numbers, as a list of lists.
     1068   
     1069    - ``dr`` -- grid of floats, the r derivative of ``z_values``.
     1070      Used to determine precision.
     1071   
     1072    - ``dtheta`` -- grid of floats, the theta derivative of ``z_values``.
     1073      Used to determine precision.
     1074     
     1075    - ``spokes`` -- integer the number of equallyspaced radial lines to plot.
    8841076
     1077    - ``circles`` -- integer the number of equally spaced circles about the center to plot.
     1078
     1079    - ``rgbcolor`` -- float array the red-green-blue color of the lines of the spiderweb.
     1080
     1081    - ``thickness`` -- positive float the thickness of the lines or points in the spiderweb.
     1082     
     1083    - ``withcolor`` -- boolean If ``True`` the spiderweb will be overlaid on the basic color plot.
    8851084    OUTPUT:
    8861085
    887     A value between `-1` (black) and `+1` (white), inclusive.
     1086    An `N x M x 3` floating point Numpy array ``X``, where
     1087    ``X[i,j]`` is an (r,g,b) tuple.
     1088   
     1089    EXAMPLES:
     1090        sage: from sage.calculus.riemann import complex_to_spiderweb
     1091        sage: import numpy
     1092        sage: zval = numpy.array([[0, 1, 1000],[.2+.3j,1,-.3j],[0,0,0]],dtype = numpy.complex128)
     1093        sage: deriv = numpy.array([[.1]],dtype = numpy.float64)
     1094        sage: complex_to_spiderweb(zval, deriv,deriv, 4,4,[0,0,0],1,False)
     1095        array([[[ 1.,  1.,  1.],
     1096                [ 1.,  1.,  1.],
     1097                [ 1.,  1.,  1.]],
     1098        <BLANKLINE>
     1099               [[ 1.,  1.,  1.],
     1100                [ 0.,  0.,  0.],
     1101                [ 1.,  1.,  1.]],
     1102        <BLANKLINE>
     1103               [[ 1.,  1.,  1.],
     1104                [ 1.,  1.,  1.],
     1105                [ 1.,  1.,  1.]]])
    8881106
    889     EXAMPLES:
     1107        sage: complex_to_spiderweb(zval, deriv,deriv, 4,4,[0,0,0],1,True)
     1108        array([[[ 1.        ,  1.        ,  1.        ],
     1109                [ 1.        ,  0.05558355,  0.05558355],
     1110                [ 0.17301243,  0.        ,  0.        ]],
     1111        <BLANKLINE>
     1112               [[ 1.        ,  0.96804683,  0.48044583],
     1113                [ 0.        ,  0.        ,  0.        ],
     1114                [ 0.77351965,  0.5470393 ,  1.        ]],
     1115        <BLANKLINE>
     1116               [[ 1.        ,  1.        ,  1.        ],
     1117                [ 1.        ,  1.        ,  1.        ],
     1118                [ 1.        ,  1.        ,  1.        ]]])
    8901119
    891     This tests it implicitly::
    892 
    893         sage: from sage.calculus.riemann import complex_to_rgb
    894         sage: import numpy
    895         sage: complex_to_rgb(numpy.array([[0, 1, 1000]],dtype = numpy.complex128))
    896         array([[[   1.,    1.,    1.],
    897                 [   1.,    0.,    0.],
    898                 [-998.,    0.,    0.]]])
    8991120    """
    900     return 1 - r
    901 
    902 cpdef complex_to_rgb(np.ndarray z_values):
    903     r"""
    904     Convert from an array of complex numbers to its corresponding matrix of
     1121    cdef Py_ssize_t i, j, imax, jmax
     1122    cdef FLOAT_T x, y, mag, arg, width, target, precision, dmag, darg
     1123    cdef COMPLEX_T z
     1124    cdef FLOAT_T DMAX = 70 # change to adjust rate_of_change cutoff below
     1125    cdef FLOAT_T MMIN = .0001
     1126    precision = thickness/150.0
     1127    imax = len(z_values)
     1128    jmax = len(z_values[0])
     1129    cdef np.ndarray[FLOAT_T, ndim=3, mode="c"] rgb
     1130    if withcolor:
     1131        rgb = complex_to_rgb(z_values)
     1132    else:
     1133        rgb = np.zeros(dtype=FLOAT, shape=(imax, jmax, 3))
     1134        rgb += 1
     1135    if circles != 0:
     1136        circ_radii = srange(0,1.0,1.0/circles)
     1137    else:
     1138        circ_radii = []
     1139    if spokes != 0:
     1140        spoke_angles = srange(-PI,PI+TWOPI/spokes,TWOPI/spokes) # both -pi and pi are included
     1141    else:
     1142        spoke_angles = []
     1143    for i in xrange(imax-2): # the d arrays are 1 smaller on each side
     1144        for j in xrange(jmax-2):
     1145            z = z_values[i+1,j+1]
     1146            mag = abs(z)
     1147            arg = phase(z)
     1148            dmag = dr[i,j]
     1149            darg = dtheta[i,j]
     1150            if darg < DMAX and mag > MMIN: #points that change too rapidly are presumed to be borders, points that are too small are presumed to be outside
     1151                for target in circ_radii:
     1152                    if abs(mag - target)/dmag < precision:
     1153                        rgb[i+1,j+1] = rgbcolor
     1154                        break
     1155                for target in spoke_angles:
     1156                    if abs(arg - target)/darg < precision:
     1157                        rgb[i+1,j+1] = rgbcolor
     1158                        break
     1159    return rgb
     1160   
     1161cpdef complex_to_rgb(np.ndarray[COMPLEX_T, ndim = 2] z_values):
     1162    """
     1163    Converts an array of complex numbers to its corresponding matrix of
    9051164    RGB values.
    9061165
    9071166    INPUT:
     
    9101169
    9111170    OUTPUT:
    9121171
    913     An `N \times M \times 3` floating point Numpy array ``X``, where
     1172    An `N x M x 3` floating point Numpy array ``X``, where
    9141173    ``X[i,j]`` is an (r,g,b) tuple.
    9151174
    9161175    EXAMPLES::
    9171176        sage: from sage.calculus.riemann import complex_to_rgb
    9181177        sage: import numpy
    9191178        sage: complex_to_rgb(numpy.array([[0, 1, 1000]],dtype = numpy.complex128))
    920         array([[[   1.,    1.,    1.],
    921                 [   1.,    0.,    0.],
    922                 [-998.,    0.,    0.]]])
     1179        array([[[ 1.        ,  1.        ,  1.        ],
     1180                [ 1.        ,  0.05558355,  0.05558355],
     1181                [ 0.17301243,  0.        ,  0.        ]]])
     1182
    9231183
    9241184        sage: complex_to_rgb(numpy.array([[0, 1j, 1000j]],dtype = numpy.complex128))
    925         array([[[  1.00000000e+00,   1.00000000e+00,   1.00000000e+00],
    926                 [  5.00000000e-01,   1.00000000e+00,   0.00000000e+00],
    927                 [ -4.99000000e+02,  -9.98000000e+02,   0.00000000e+00]]])
     1185        array([[[ 1.        ,  1.        ,  1.        ],
     1186                [ 0.52779177,  1.        ,  0.05558355],
     1187                [ 0.08650622,  0.17301243,  0.        ]]])
     1188
    9281189
    9291190    """
    930     cdef unsigned int i, j, imax, jmax
    931     cdef double x, y, mag, arg
    932     cdef double lightness, hue, top, bot
    933     cdef double r, g, b
     1191    cdef Py_ssize_t i, j, imax, jmax
     1192    cdef FLOAT_T x, y, mag, arg
     1193    cdef FLOAT_T lightness, hue, top, bot
     1194    cdef FLOAT_T r, g, b
    9341195    cdef int ihue
    935     cdef z
    936 
    937     from cmath import phase
    938 
     1196    cdef COMPLEX_T z
    9391197    imax = len(z_values)
    9401198    jmax = len(z_values[0])
    941     cdef np.ndarray[np.float_t, ndim=3, mode="c"] rgb = np.empty(
    942         dtype=np.float, shape=(imax, jmax, 3))
    943 
     1199    cdef np.ndarray[FLOAT_T, ndim=3, mode="c"] rgb = np.empty(
     1200        dtype=FLOAT, shape=(imax, jmax, 3))
    9441201    sig_on()
    945     for i from 0 <= i < imax:
     1202    for i from 0 <= i < imax: #replace with xrange?
    9461203        row = z_values[i]
    947         for j from 0 <= j < jmax:
     1204        for j from 0 <= j < jmax: # replace with xrange?
    9481205            z = row[j]
    9491206            mag = abs(z)
    9501207            arg = phase(z)
    951             lightness = mag_to_lightness(mag)
     1208            # tweak these levels to adjust how bright/dark the colors appear
     1209            # output can range from -1 (black) to 1 (white)
     1210            lightness = -(atan(log(mag*1.5+1)) * (4/PI) - 1)
    9521211            # in hsv, variable value, full saturation (s=1, v=1+lightness)
    9531212            if lightness < 0:
    9541213                bot = 0
     
    9941253            rgb[i, j, 2] = b
    9951254    sig_off()
    9961255    return rgb
     1256   
     1257cpdef analytic_boundary(FLOAT_T t, int n):
     1258    """
     1259    Provides an exact (for n = infinity) riemann boundary
     1260    correspondence for the ellipse with axes 1 + epsilon and 1 - epsilon. The
     1261    boundary is therefore given by e^(I*t)+epsilon*e^(-I*t). It is primarily useful
     1262    for testing the accuracy of the numerical Riemann Map.
     1263   
     1264    INPUT:
     1265   
     1266    - ``t`` -- The boundary parameter, from 0 to two*pi
     1267   
     1268    - ``n`` -- Integer, The number of terms to include. 10 is fairly accurate,
     1269        20 is very accurate.
     1270       
     1271    OUTPUT:
     1272   
     1273    A theta value from 0 to 2*pi, corresponding to the point on the circle e^(I*theta)
     1274   
     1275    TESTS:
     1276   
     1277    Checking the accuracy for different n values::
     1278        sage: from sage.calculus.riemann import analytic_boundary
     1279        sage: t100 = analytic_boundary(pi/2,100)
     1280        sage: abs(analytic_boundary(pi/2,10) - t100) < 10^-8
     1281        True
     1282        sage: abs(analytic_boundary(pi/2,20) - t100) < 10^-15
     1283        True
     1284       
     1285    Using this to check the accuracy of the Riemann_Map boundary::
     1286        sage: f(t) = e^(I*t)+.3*e^(-I*t)
     1287        sage: fp(t) = I*e^(I*t)-I*.3*e^(-I*t)
     1288        sage: m = Riemann_Map([f], [fp],0,200)
     1289        sage: s = spline(m.get_theta_points())
     1290        sage: test_pt = uniform(0,2*pi)
     1291        sage: s(test_pt) - analytic_boundary(test_pt,20) < 10^-4
     1292        True
     1293    """
     1294    cdef FLOAT_T i
     1295    cdef FLOAT_T result = t
     1296    for i from 1 <= i < n+1:
     1297        result += (2*(-1)**i/i)*(.3**i/(1+.3**(2*i)))*sin(2*i*t)
     1298    return result
    9971299
     1300   
     1301
     1302cpdef cauchy_kernel(t, args):
     1303    """
     1304    Intermediate function for the integration in ``analytic_interior``.
     1305   
     1306    INPUT:
     1307   
     1308    - ``t`` -- The boundary parameter, meant to be integrated over
     1309   
     1310    - ``args`` -- a tuple containing:
     1311        - ``z`` -- Complex, the point to be mapped.
     1312       
     1313        - ``n`` -- Integer, The number of terms to include. 10 is fairly accurate,
     1314            20 is very accurate.
     1315           
     1316        - ``part`` -- will return the real ('r'), imaginary ('i') or complex ('c')
     1317        value of the kernel
     1318   
     1319    TESTS:
     1320        Primarily tested implicitly by analytic_interior,
     1321       
     1322        Simple test::
     1323        sage: from sage.calculus.riemann import cauchy_kernel
     1324        sage: cauchy_kernel(.5,(.1+.2*I, 10,'c'))
     1325        (-0.584136405997...+0.5948650858950...j)
     1326
     1327
     1328    """
     1329    cdef COMPLEX_T result
     1330    cdef COMPLEX_T z = args[0]
     1331    cdef int n = args[1]
     1332    part = args[2]
     1333    result = exp(I*analytic_boundary(t,n))/(exp(I*t)+.3*exp(-I*t)-z)*(I*exp(I*t)-I*.3*exp(-I*t))
     1334    if part == 'c':
     1335        return result
     1336    elif part == 'r':
     1337        return result.real
     1338    elif part == 'i':
     1339        return result.imag
     1340    else: return None
     1341   
     1342cpdef analytic_interior(COMPLEX_T z, int n):
     1343    """
     1344    Provides a nearly exact compuation of the Riemann Map of an interior
     1345    point of the ellipse with axes 1 + epsilon and 1 - epsilon. It is primarily useful
     1346    for testing the accuracy of the numerical Riemann Map.
     1347   
     1348    INPUT:
     1349   
     1350    - ``z`` -- Complex, the point to be mapped.
     1351   
     1352    - ``n`` -- Integer, The number of terms to include. 10 is fairly accurate,
     1353        20 is very accurate.
     1354
     1355    TESTS:
     1356        Testing the accuracy of Riemann_Map::
     1357        sage: from sage.calculus.riemann import analytic_interior
     1358        sage: f(t) = e^(I*t)+.3*e^(-I*t)
     1359        sage: fp(t) = I*e^(I*t)-I*.3*e^(-I*t)
     1360        sage: m = Riemann_Map([f],[fp],0,200)
     1361        sage: abs(m.riemann_map(.5)-analytic_interior(.5,20)) < 10^-4
     1362        True
     1363        sage: m = Riemann_Map([f],[fp],0,2000)
     1364        sage: abs(m.riemann_map(.5)-analytic_interior(.5,20)) < 10^-6
     1365        True
     1366    """
     1367    # evaluates the cauchy integral of the boundary, split into the real
     1368    # and imaginary results because numerical_integral can't handle complex data.
     1369    rp = 1/(TWOPI)*numerical_integral(cauchy_kernel,0,2*pi, params = [z,n,'i'])[0]
     1370    ip = 1/(TWOPI*I)*numerical_integral(cauchy_kernel,0,2*pi, params = [z,n,'r'])[0]
     1371    return rp + ip