Ticket #11273: trac_11273rebasepart1.patch
File trac_11273rebasepart1.patch, 47.1 KB (added by , 8 years ago) 


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 3 3 4 4 AUTHORS: 5 5 6  Ethan Van Andel (2009 ): initial version6  Ethan Van Andel (20092011): initial version and upgrades 7 7 8 8  Robert Bradshaw (2009): his "complex_plot" was adapted for plot_colored 9 9 … … 11 11 """ 12 12 13 13 #***************************************************************************** 14 # Copyright (C) 20 09Ethan Van Andel <evlutte@gmail.com>,14 # Copyright (C) 2011 Ethan Van Andel <evlutte@gmail.com>, 15 15 # 16 16 # Distributed under the terms of the GNU General Public License (GPL) 17 17 # … … 28 28 include "sage/ext/stdsage.pxi" 29 29 include "sage/ext/interrupt.pxi" 30 30 31 from sage.plot.primitive import GraphicPrimitive32 31 from sage.misc.decorators import options 33 32 from sage.plot.all import list_plot, Graphics 34 from sage.plot.misc import setup_for_eval_on_grid35 33 36 34 from sage.ext.fast_eval import fast_callable 37 35 … … 43 41 44 42 from sage.plot.complex_plot import ComplexPlot 45 43 44 from sage.gsl.integration import numerical_integral 45 46 46 47 import numpy as np 47 48 cimport numpy as np 48 49 49 50 from math import pi 50 51 from math import sin 51 52 from math import cos 53 from math import sqrt 54 55 from math import log # used for complex plot lightness 56 from math import atan 52 57 53 58 from cmath import exp 54 59 from cmath import phase 55 60 56 cdef double PI = pi 57 cdef double TWOPI = 2*PI 58 cdef I = complex(0,1) 61 from random import uniform # used for accuracy tests 59 62 60 63 FLOAT = np.float64 61 64 ctypedef np.float64_t FLOAT_T 62 65 66 COMPLEX = np.complex128 67 ctypedef np.complex128_t COMPLEX_T 68 69 cdef FLOAT_T PI = pi 70 cdef FLOAT_T TWOPI = 2*PI 71 cdef COMPLEX_T I = complex(0,1) 63 72 64 73 cdef class Riemann_Map: 65 74 """ 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 69 81 the complex plane to the unit disc. The Ahlfors map does the same thing 70 82 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. 75 92 76 93 INPUT: 77 94 … … 85 102 Must be in the same order as ``fs``. 86 103 87 104  ``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. 89 107 90 108 The following inputs must all be passed in as named parameters: 91 109 92 110  ``N``  integer (default: ``500``), the number of collocation points 93 111 used to compute the map. More points will give more accurate results, 94 112 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 96 120  ``ncorners``  integer (default: ``4``), if mapping a figure with 97 (equally tspaced) corners, better results may be obtained by 121 (equally tspaced) cornerscorners that make a significant change in 122 the direction of the boundary, better results may be sometimes obtained by 98 123 accurately giving this parameter. Used to add the proper constant to 99 the theta correspond ance function.124 the theta correspondence function. 100 125 101 126  ``opp``  boolean (default: ``False``), set to ``True`` in very rare 102 cases where the theta correspond ance function is off by ``pi``, that127 cases where the theta correspondence function is off by ``pi``, that 103 128 is, if red is mapped left of the origin in the color plot. 129 104 130 105 131 EXAMPLES: 106 132 … … 142 168 .. [KT] N. Kerzman and M. R. Trummer. "Numerical Conformal Mapping via 143 169 the Szego kernel". Journal of Computational and Applied Mathematics, 144 170 14(12): 111123, 1986. 171 172 .. [BSV] M. Bolt, S. Snoeyink, E. Van Andel. "Visual representation of 173 the Riemann map and Ahlfors map via the KerzmanStein equation". 174 Involve 34 (2010), 405420. 145 175 """ 146 176 cdef int N, B, ncorners 147 177 cdef f 148 178 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 151 182 cdef np.ndarray p_vector_inverse, sinalpha, cosalpha, theta_array 152 183 cdef x_range, y_range 184 cdef exterior 153 185 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): 155 187 """ 156 188 Initializes the ``Riemann_Map`` class. See the class ``Riemann_Map`` 157 189 for full documentation on the input of this initialization method. … … 160 192 161 193 sage: f(t) = e^(I*t)  0.5*e^(I*t) 162 194 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) 164 196 """ 165 a = np.complex128(a) 197 cdef double xmax, xmin, ymax, ymin, space 198 cdef int i, k 166 199 if N <= 2: 167 200 raise ValueError( 168 201 "The number of collocation points must be > 2.") 202 if exterior and a!=0: 203 raise ValueError("The exterior map requires a=0") 169 204 try: 170 205 fs = [fast_callable(f, domain=CDF) for f in fs] 171 206 fprimes = [fast_callable(f, domain=CDF) for f in fprimes] … … 176 211 self.ncorners = ncorners 177 212 self.N = N # Number of collocation pts 178 213 self.opp = opp 179 cdef int i, k214 self.exterior = exterior 180 215 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) 183 218 for i in xrange(N): 184 219 self.tk2[i] = self.tk[i] 185 220 self.tk2[N] = TWOPI 186 221 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) 189 227 # 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[Ni1]) 232 cps[k, i] = np.complex(1/fk) 233 dps[k, i] = np.complex(1/fk**2*fprimes[k](self.tk[Ni1])) 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 200 253 self.x_range = (xmin  space, xmax + space) 201 254 self.y_range = (ymin  space, ymax + space) 255 # Now we do some more computation 202 256 self._generate_theta_array() 203 257 self._generate_interior_mapper() 204 258 self._generate_inverse_mapper() 205 259 260 206 261 def _repr_(self): 207 262 """ 208 263 Return a string representation of this ``Riemann_Map`` object. … … 211 266 212 267 sage: f(t) = e^(I*t)  0.5*e^(I*t) 213 268 sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(I*t) 214 sage: isinstance(Riemann_Map([f], [fprime], 0 )._repr_(), str) # long time269 sage: isinstance(Riemann_Map([f], [fprime], 0, N = 10)._repr_(), str) # long time 215 270 True 216 271 """ 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." 218 273 219 274 cdef _generate_theta_array(self): 220 275 """ 221 276 Generates the essential data for the Riemann map, primarily the 222 Szego kernel and boundary correspond ance.277 Szego kernel and boundary correspondence. See [KT]_ for the algorithm. 223 278 224 279 TESTS:: 225 280 226 281 sage: f(t) = e^(I*t)  0.5*e^(I*t) 227 282 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) 229 284 """ 230 c p = 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() 232 287 cdef int N = self.N 233 288 cdef int NB = N * self.B 234 289 cdef int B = self.B 235 290 cdef int i, k 236 291 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 244 295 cdef np.ndarray[FLOAT_T, ndim=2] theta_array 245 # Set ing things up to use the Nystrom method296 # Setting things up to use the Nystrom method 246 297 adp = abs(dp) 247 298 sadp = np.sqrt(adp) 248 299 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) 250 301 g = sadp * hconj 251 302 normalized_dp=dp/adp 252 303 C = I / N * sadp # equivalent to TWOPI / N * 1 / (TWOPI * I) * sadp … … 260 311 np.seterr(divide=errdivide,invalid=errinvalid) # resets the error handling 261 312 for i in xrange(NB): 262 313 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 264 315 # the allimportant 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) 266 317 self.szego = szego.reshape([B, N]) 267 318 start = 0 268 # Finding the theta correspond ance using phase. Misbehaves for some319 # Finding the theta correspondence using phase. Misbehaves for some 269 320 # regions. 270 321 if B != 1: 271 322 theta_array = np.zeros([1, NB]) … … 275 326 [theta_array.reshape([B, N]), np.zeros([B, 1])], axis=1) 276 327 for k in xrange(B): 277 328 self.theta_array[k, N] = self.theta_array[k, 0] + TWOPI 278 # Finding the theta correspond ance using ab. Well behaved, but329 # Finding the theta correspondence using abs. Well behaved, but 279 330 # doesn't work on multiply connected domains. 280 331 else: 281 332 phi2 = phi.reshape([self.B, N]) … … 342 393 sage: s = spline(points) 343 394 sage: s(3*pi / 4) 344 395 0.00121587378429... 396 sage: plot(s,0,2*pi) # plot the kernel 345 397 346 398 The unit circle with a small hole:: 347 399 … … 378 430 """ 379 431 Returns an array of points of the form 380 432 ``[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. 384 439 385 440 INPUT: 386 441 387 442 The following input must all be passed in as named parameters: 388 443 389  ``boundary``  integer (default: `` ``1) if < 0,444  ``boundary``  integer (default: ``1``) if < 0, 390 445 ``get_theta_points()`` will return the points for all boundaries. 391 446 If >= 0, ``get_theta_points()`` will return only the points for 392 447 the boundary specified. … … 409 464 Extending the points by a spline:: 410 465 411 466 sage: s = spline(points) 412 sage: s(3*pi / 4) 467 sage: s(3*pi / 4) 413 468 1.62766037996... 414 469 415 470 The unit circle with a small hole:: … … 420 475 sage: hfprime(t) = 0.5*I*e^(I*t) 421 476 sage: m = Riemann_Map([f, hf], [hf, hfprime], 0.5 + 0.5*I) 422 477 423 Getting the szegofor a specifc boundary::478 Getting the boundary correspondence for a specifc boundary:: 424 479 425 480 sage: tp0 = m.get_theta_points(boundary=0) 426 481 sage: tp1 = m.get_theta_points(boundary=1) … … 437 492 438 493 cdef _generate_interior_mapper(self): 439 494 """ 440 Generates the data necessary to use the ``r eimann_map()`` function.441 As much setup as possible is done here to minimize what has to be442 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``. 443 498 444 499 TESTS:: 445 500 446 501 sage: f(t) = e^(I*t)  0.5*e^(I*t) 447 502 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) 449 504 """ 450 505 cdef int N = self.N 451 506 cdef double complex coeff = 3*I / (8*N) … … 488 543 cdef np.ndarray[double complex, ndim=1] pq = self.cps[:,list(range(N))+[0]].flatten() 489 544 self.pre_q_vector = pq 490 545 491 cpdef riemann_map(self, pt):546 cpdef riemann_map(self, COMPLEX_T pt): 492 547 """ 493 548 Returns the Riemann mapping of a point. That is, given ``pt`` on 494 the interior of the mapped region, ``riemann_map ()`` will return549 the interior of the mapped region, ``riemann_map`` will return 495 550 the point on the unit disk that ``pt`` maps to. Note that this 496 method only works for interior points; itbreaks down very close551 method only works for interior points; accuracy breaks down very close 497 552 to the boundary. To get boundary corrospondance, use 498 553 ``get_theta_points()``. 499 554 … … 525 580 sage: m.riemann_map(np.complex(3, 0.0001)) 526 581 (1.405757...e05+8.06...e10j) 527 582 """ 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) 532 596 533 597 cdef _generate_inverse_mapper(self): 534 598 """ 535 599 Generates the data necessary to use the 536 ``inverse_r eimann_map()`` function. As much setup as possible is537 done here to minimize what has to be done in600 ``inverse_riemann_map()`` function. As much setup as possible is 601 done here to minimize the computation that must be done in 538 602 ``inverse_riemann_map()``. 539 603 540 604 TESTS:: 541 605 542 606 sage: f(t) = e^(I*t)  0.5*e^(I*t) 543 607 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) 545 609 """ 546 610 cdef int N = self.N 547 611 cdef int B = self.B … … 567 631 for i in xrange(N): 568 632 self.cosalpha[k, i] = cos(theta_array[k, i]) 569 633 570 def inverse_riemann_map(self,pt):634 cpdef inverse_riemann_map(self, COMPLEX_T pt): 571 635 """ 572 636 Returns the inverse Riemann mapping of a point. That is, given ``pt`` 573 on the interior of the unit disc, ``inverse_r eimann_map()`` will637 on the interior of the unit disc, ``inverse_riemann_map()`` will 574 638 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. 580 641 581 642 INPUT: 582 643 … … 604 665 sage: m.inverse_riemann_map(np.complex(0.2, 0.5)) 605 666 (0.156280570579...+0.321819151891...j) 606 667 """ 607 pt = np.complex128(pt) 668 if self.exterior: 669 pt = 1/pt 608 670 r = abs(pt) 609 671 if r == 0: 610 672 stheta = 0 … … 613 675 stheta = pt.imag / r 614 676 ctheta = pt.real / r 615 677 k = 0 616 return(1  r**2) / TWOPI * np.dot(678 mapped = (1  r**2) / TWOPI * np.dot( 617 679 self.p_vector_inverse[k] * self.cps[k], 618 680 1 / (1 + r**2  2*r * 619 681 (ctheta * self.cosalpha[k]  stheta * self.sinalpha[k]))) 682 if self.exterior: 683 return 1/mapped 684 else: 685 return mapped 620 686 621 687 def plot_boundaries(self, plotjoined=True, rgbcolor=[0,0,0], thickness=1): 622 688 """ … … 657 723 """ 658 724 plots = range(self.B) 659 725 for k in xrange(self.B): 660 # This should be eliminated when the thickness/pointsize issue726 # This conditional should be eliminated when the thickness/pointsize issue 661 727 # is resolved later. Same for the others in plot_spiderweb(). 662 728 if plotjoined: 663 729 plots[k] = list_plot( … … 668 734 comp_pt(self.cps[k], 1), rgbcolor=rgbcolor, 669 735 pointsize=thickness) 670 736 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: 671 759 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((ymaxymin)/ 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') 672 804 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): 674 806 """ 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. 678 814 679 815 INPUT: 680 816 … … 690 826 plot. Each radial line is made by ``1*pts`` points, each circle 691 827 has ``2*pts`` points. Note that high values may cause erratic 692 828 behavior of the radial lines near the boundaries. 829  only for simply connected domains 693 830 694 831  ``linescale``  float between 0 and 1. Shrinks the radial lines 695 832 away from the boundary to reduce erratic behavior. 833  only for simply connected domains 696 834 697 835  ``rgbcolor``  float array (default: ``[0,0,0]``) the 698 836 redgreenblue color of the spiderweb. … … 703 841  ``plotjoined``  boolean (default: ``True``) If ``False``, 704 842 discrete points will be drawn; otherwise they will be connected 705 843 by lines. 844  only for simply connected domains 706 845 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 707 853 EXAMPLES: 708 854 709 855 General usage:: … … 732 878 sage: m = Riemann_Map([f], [fprime], 0, 1000) 733 879 sage: m.plot_spiderweb() 734 880 """ 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]742 881 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(angleTWOPI))) 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) 752 924 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(angleTWOPI))) 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 = (xmaxxmin)/plot_points 929 ystep = (ymaxymin)/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 776 936 777 937 @options(interpolation='catrom') 778 938 def plot_colored(self, plot_range=[], int plot_points=100, **options): 779 939 """ 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. 783 942 784 943 INPUT: 785 944 … … 789 948 ``(xmin, xmax, ymin, ymax)``. Declare if you do not want the plot 790 949 to use the default range for the figure. 791 950 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 795 955 796 956 EXAMPLES: 797 957 … … 818 978 sage: m = Riemann_Map([f], [fprime], 0, 1000) 819 979 sage: m.plot_colored() 820 980 """ 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) 845 982 g = Graphics() 846 983 g.add_primitive(ComplexPlot(complex_to_rgb(z_values), (xmin, xmax), (ymin, ymax),options)) 847 984 return g 848 985 849 986 cdef comp_pt(clist, loop=True): 850 987 """ 851 This function converts a list of complex numbersto the plottable988 Converts a list of complex numbers xderivs = get_derivatives(z_values, xstep, ystep)[0]to the plottable 852 989 `(x,y)` form. If ``loop=True``, then the first point will be 853 990 added as the last, i.e. to plot a closed circle. 854 991 … … 874 1011 if loop: 875 1012 list2[len(clist)] = list2[0] 876 1013 return list2 1014 1015 cpdef 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 quickanddirty 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(xdelta))/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 877 1060 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 nonnegative 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 1092 cpdef 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): 879 1093 """ 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: 881 1098 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. 883 1108 884  `` r``  a nonnegative real number.1109  ``circles``  integer the number of equally spaced circles about the center to plot. 885 1110 886 OUTPUT:1111  ``rgbcolor``  float array the redgreenblue color of the lines of the spiderweb. 887 1112 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. 889 1116 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 890 1122 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.]]]) 891 1139 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(imax2): # the d arrays are 1 smaller on each side 1176 for j in xrange(jmax2): 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 893 1195 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): 1196 cpdef complex_to_rgb(np.ndarray[COMPLEX_T, ndim = 2] z_values): 904 1197 r""" 905 1198 Convert from a (Numpy) array of complex numbers to its corresponding 906 1199 matrix of RGB values. For internal use of :meth:`~Riemann_Map.plot_colored` … … 920 1213 sage: from sage.calculus.riemann import complex_to_rgb 921 1214 sage: import numpy 922 1215 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 926 1220 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.00000000e01, 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 930 1225 931 1226 TESTS:: 932 1227 … … 935 1230 ... 936 1231 TypeError: Argument 'z_values' has incorrect type (expected numpy.ndarray, got list) 937 1232 """ 938 cdef unsigned int i, j, imax, jmax939 cdef doublex, y, mag, arg940 cdef doublelightness, hue, top, bot941 cdef doubler, g, b1233 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 942 1237 cdef int ihue 943 cdef z 944 945 from cmath import phase 1238 cdef COMPLEX_T z 946 1239 947 1240 imax = len(z_values) 948 1241 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)) 951 1244 952 1245 sig_on() 953 for i from 0 <= i < imax: 1246 for i from 0 <= i < imax: #replace with xrange? 954 1247 row = z_values[i] 955 for j from 0 <= j < jmax: 1248 for j from 0 <= j < jmax: #replace with xrange? 956 1249 z = row[j] 957 1250 mag = abs(z) 958 1251 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) 960 1255 # in hsv, variable value, full saturation (s=1, v=1+lightness) 961 1256 if lightness < 0: 962 1257 bot = 0 … … 1002 1297 rgb[i, j, 2] = b 1003 1298 sig_off() 1004 1299 return rgb 1300 1301 cpdef 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 1005 1343 1344 1345 1346 cpdef 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 1386 cpdef 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