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


sage/calculus/riemann.pyx
# HG changeset patch # User KarlDieter Crisman <kcrisman@gmail.com> # Date 1337993071 14400 # Node ID 4253cee9001cbd6d8476d6a6293a12a98d6ccc3b # Parent 1f9aeedfb653ac789d45ce55c14704d8ae97e3e9 Trac 11273 reviewer patch Lots of whitespace removal and fixing documentation issues diff git a/sage/calculus/riemann.pyx b/sage/calculus/riemann.pyx
a b 72 72 73 73 cdef class Riemann_Map: 74 74 """ 75 76 The ``Riemann_Map`` class computes an interior or exterior Riemann map, 77 or an Ahlfors map of a region given by the supplied boundary curve(s) 78 and center ' point. The class also provides various methods to75 76 The ``Riemann_Map`` class computes an interior or exterior Riemann map, 77 or an Ahlfors map of a region given by the supplied boundary curve(s) 78 and center point. The class also provides various methods to 79 79 evaluate, visualize, or extract data from the map. 80 80 81 81 A Riemann map conformally maps a simply connected region in 82 82 the complex plane to the unit disc. The Ahlfors map does the same thing 83 83 for multiply connected regions. 84 85 Note that all the methods are numerical. As a result all answers have 86 some imprecision. Moreover, maps computed with small number of 87 collocation points, or for unusually shaped regions may be very88 inaccurate. Error computations for the ellipse can be found in the 89 documentation for ``analytic_boundary()`` and ``analytic_interior()``.90 91 [BSV] provides an overview of the Riemann map and discusses the research84 85 Note that all the methods are numerical. As a result all answers have 86 some imprecision. Moreover, maps computed with small number of 87 collocation points, or for unusually shaped regions, may be very 88 inaccurate. Error computations for the ellipse can be found in the 89 documentation for :meth:`analytic_boundary` and :meth:`analytic_interior`. 90 91 [BSV]_ provides an overview of the Riemann map and discusses the research 92 92 that lead to the creation of this module. 93 93 94 94 INPUT: … … 111 111  ``N``  integer (default: ``500``), the number of collocation points 112 112 used to compute the map. More points will give more accurate results, 113 113 especially near the boundaries, but will take longer to compute. 114 114 115 115  ``exterior``  boolean (default: ``False``), if set to ``True``, the 116 116 exterior map will be computed, mapping the exterior of the region to the 117 117 exterior of the unit circle. 118 119 The following inputs may be passed as named parameters in unusual 118 119 The following inputs may be passed as named parameters in unusual 120 120 circumstances: 121 121 122 122  ``ncorners``  integer (default: ``4``), if mapping a figure with 123 (equally tspaced) corners corners that make a significant change in124 the direction of the boundary ,better results may be sometimes obtained by123 (equally tspaced) corners  corners that make a significant change in 124 the direction of the boundary  better results may be sometimes obtained by 125 125 accurately giving this parameter. Used to add the proper constant to 126 126 the theta correspondence function. 127 127 128 128  ``opp``  boolean (default: ``False``), set to ``True`` in very rare 129 129 cases where the theta correspondence function is off by ``pi``, that 130 130 is, if red is mapped left of the origin in the color plot. 131 131 132 132 133 133 EXAMPLES: 134 134 … … 145 145 sage: hf(t) = 0.5*e^(I*t) 146 146 sage: hfprime(t) = 0.5*I*e^(I*t) 147 147 sage: m = Riemann_Map([f, hf], [fprime, hfprime], 0.5 + 0.5*I) 148 sage: m.plot_colored() + m.plot_spiderweb() # long time 148 149 149 150 A square:: 150 151 151 sage: ps = polygon_spline([(1, 1), (1, 1), (1, 1), (1, 1)]) # long time152 sage: f = lambda t: ps.value(real(t)) # long time153 sage: fprime = lambda t: ps.derivative(real(t)) # long time154 sage: m = Riemann_Map([f], [fprime], 0.25, ncorners=4) # long time152 sage: ps = polygon_spline([(1, 1), (1, 1), (1, 1), (1, 1)]) 153 sage: f = lambda t: ps.value(real(t)) 154 sage: fprime = lambda t: ps.derivative(real(t)) 155 sage: m = Riemann_Map([f], [fprime], 0.25, ncorners=4) 155 156 sage: m.plot_colored() + m.plot_spiderweb() # long time 156 157 157 158 Compute rough error for this map:: … … 170 171 .. [KT] N. Kerzman and M. R. Trummer. "Numerical Conformal Mapping via 171 172 the Szego kernel". Journal of Computational and Applied Mathematics, 172 173 14(12): 111123, 1986. 173 174 .. [BSV] M. Bolt, S. Snoeyink, E. Van Andel. "Visual representation of 175 the Riemann map and Ahlfors map via the KerzmanStein equation". 174 175 .. [BSV] M. Bolt, S. Snoeyink, E. Van Andel. "Visual representation of 176 the Riemann map and Ahlfors map via the KerzmanStein equation". 176 177 Involve 34 (2010), 405420. 178 177 179 """ 178 180 cdef int N, B, ncorners 179 181 cdef f … … 185 187 cdef x_range, y_range 186 188 cdef exterior 187 189 188 def __init__(self, fs, fprimes, COMPLEX_T a, int N=500, int ncorners=4, 190 def __init__(self, fs, fprimes, COMPLEX_T a, int N=500, int ncorners=4, 189 191 opp=False, exterior = False): 190 192 191 193 """ 192 Initializes the ``Riemann_Map`` class. See the class ``Riemann_Map``194 Initializes the ``Riemann_Map`` class. See the class :class:`Riemann_Map` 193 195 for full documentation on the input of this initialization method. 194 196 195 197 TESTS:: … … 226 228 if exterior and (self.B > 1): 227 229 raise ValueError( 228 230 "The exterior map is undefined for multiply connected domains") 229 cdef np.ndarray[COMPLEX_T,ndim=2] cps = np.zeros([self.B, N], 231 cdef np.ndarray[COMPLEX_T,ndim=2] cps = np.zeros([self.B, N], 230 232 dtype=COMPLEX) 231 cdef np.ndarray[COMPLEX_T,ndim=2] dps = np.zeros([self.B, N], 233 cdef np.ndarray[COMPLEX_T,ndim=2] dps = np.zeros([self.B, N], 232 234 dtype=COMPLEX) 233 235 # Find the points on the boundaries and their derivatives. 234 236 if self.exterior: … … 266 268 267 269 def _repr_(self): 268 270 """ 269 Return a string representation of this ``Riemann_Map`` object.271 Return a string representation of this :class:`Riemann_Map` object. 270 272 271 273 TESTS:: 272 274 273 275 sage: f(t) = e^(I*t)  0.5*e^(I*t) 274 276 sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(I*t) 275 277 sage: isinstance(Riemann_Map([f], [fprime], 0, N = 10)._repr_(), str) # long time 276 278 True 277 279 """ 278 return "A Riemann or A lfohrs mapping of a figure to the unit circle."280 return "A Riemann or Ahlfors mapping of a figure to the unit circle." 279 281 280 282 cdef _generate_theta_array(self): 281 283 """ … … 317 319 for i in xrange(NB): 318 320 K[i, i] = 1 319 321 # Nystrom Method for solving 2nd kind integrals 320 phi = np.linalg.solve(K, g) / NB * TWOPI # Nystrom Method for solving 2nd kind integrals322 phi = np.linalg.solve(K, g) / NB * TWOPI 321 323 # the allimportant Szego kernel 322 324 szego = np.array(phi.flatten() / np.sqrt(dp), dtype=COMPLEX) 323 325 self.szego = szego.reshape([B, N]) … … 370 372 The following inputs must all be passed in as named parameters: 371 373 372 374  ``boundary``  integer (default: ``1``) if < 0, 373 ``get_theta_points()`` will return the points for all boundaries.374 If >= 0, ``get_theta_points()`` will return only the points for375 :meth:`get_theta_points` will return the points for all boundaries. 376 If >= 0, :meth:`get_theta_points` will return only the points for 375 377 the boundary specified. 376 378 377 379  ``absolute_value``  boolean (default: ``False``) if ``True``, will … … 437 439 Returns an array of points of the form 438 440 ``[t value, theta in e^(I*theta)]``, that is, a discretized version 439 441 of the theta/boundary correspondence function. In other words, a point 440 in this array [t1, t2] represents that the boundary point given by f(t1) 442 in this array [t1, t2] represents that the boundary point given by f(t1) 441 443 is mapped to a point on the boundary of the unit circle given by e^(I*t2). 442 443 For multiply connected domains, ``get_theta_points`` will list the 444 445 For multiply connected domains, ``get_theta_points`` will list the 444 446 points for each boundary in the order that they were supplied. 445 447 446 448 INPUT: … … 470 472 Extending the points by a spline:: 471 473 472 474 sage: s = spline(points) 473 sage: s(3*pi / 4) 475 sage: s(3*pi / 4) 474 476 1.627660... 475 477 476 478 The unit circle with a small hole:: … … 498 500 499 501 cdef _generate_interior_mapper(self): 500 502 """ 501 Generates the data necessary to use the ``riemann_map`` function.502 As much setup as possible is done here to minimize the computation 503 Generates the data necessary to use the :meth:`riemann_map` function. 504 As much setup as possible is done here to minimize the computation 503 505 that must be done in ``riemann_map``. 504 506 505 507 TESTS:: 506 508 507 509 sage: f(t) = e^(I*t)  0.5*e^(I*t) 508 510 sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(I*t) 509 sage: m = Riemann_Map([f], [fprime], 0, N = 10) 511 sage: m = Riemann_Map([f], [fprime], 0, N = 10) # indirect doctest 510 512 """ 511 513 cdef int N = self.N 512 514 cdef double complex coeff = 3*I / (8*N) … … 556 558 the point on the unit disk that ``pt`` maps to. Note that this 557 559 method only works for interior points; accuracy breaks down very close 558 560 to the boundary. To get boundary corrospondance, use 559 ``get_theta_points()``.561 :meth:`get_theta_points`. 560 562 561 563 INPUT: 562 564 … … 577 579 sage: m = Riemann_Map([f], [fprime], 0) 578 580 sage: m.riemann_map(0.25 + sqrt(0.5)) 579 581 (0.137514...+0.876696...j) 582 sage: I = CDF.gen() 580 583 sage: m.riemann_map(1.3*I) 581 584 (1.56...e05+0.989694...j) 582 sage: I = CDF.gen()583 585 sage: m.riemann_map(0.4) 584 586 (0.73324...+3.2...e06j) 585 587 sage: import numpy as np 586 588 sage: m.riemann_map(np.complex(3, 0.0001)) 587 589 (1.405757...e05+8.06...e10j) 588 590 """ 589 591 590 592 cdef COMPLEX_T pt1 591 593 cdef np.ndarray[COMPLEX_T, ndim=1] q_vector 592 594 if self.exterior: … … 603 605 cdef _generate_inverse_mapper(self): 604 606 """ 605 607 Generates the data necessary to use the 606 ``inverse_riemann_map()`` function. As much setup as possible is607 done here to minimize the computation that must be done in 608 ``inverse_riemann_map ()``.608 :meth:`inverse_riemann_map` function. As much setup as possible is 609 done here to minimize the computation that must be done in 610 ``inverse_riemann_map``. 609 611 610 612 TESTS:: 611 613 612 614 sage: f(t) = e^(I*t)  0.5*e^(I*t) 613 615 sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(I*t) 614 sage: m = Riemann_Map([f], [fprime], 0, N = 10) 616 sage: m = Riemann_Map([f], [fprime], 0, N = 10) # indirect doctest 615 617 """ 616 618 cdef int N = self.N 617 619 cdef int B = self.B … … 640 642 cpdef inverse_riemann_map(self, COMPLEX_T pt): 641 643 """ 642 644 Returns the inverse Riemann mapping of a point. That is, given ``pt`` 643 on the interior of the unit disc, ``inverse_riemann_map()`` will 644 return the point on the original region that would be Riemann 645 mapped to ``pt``. Note that this method does not work for multiply 645 on the interior of the unit disc, ``inverse_riemann_map()`` will 646 return the point on the original region that would be Riemann 647 mapped to ``pt``. Note that this method does not work for multiply 646 648 connected domains. 647 649 648 650 INPUT: … … 729 731 """ 730 732 plots = range(self.B) 731 733 for k in xrange(self.B): 732 # This conditional should be eliminated when the thickness/pointsize 734 # This conditional should be eliminated when the thickness/pointsize 733 735 # issue is resolved later. Same for the others in plot_spiderweb(). 734 736 if plotjoined: 735 737 plots[k] = list_plot( … … 740 742 comp_pt(self.cps[k], 1), rgbcolor=rgbcolor, 741 743 pointsize=thickness) 742 744 return sum(plots) 743 744 745 746 745 747 cpdef compute_on_grid(self, plot_range, int x_points): 746 748 """ 747 749 Computes the riemann map on a grid of points. Note that these points 748 750 are complex of the form z = x + y*i. 749 751 750 752 INPUT: 751 752  ``plot_range``  a tuple of the form [xmin, xmax, ymin, ymax].753 If the value is ``[]``, the default plotting window of the map will 753 754  ``plot_range``  a tuple of the form ``[xmin, xmax, ymin, ymax]``. 755 If the value is ``[]``, the default plotting window of the map will 754 756 be used. 755 757 756 758  ``x_points``  int, the size of the grid in the x direction 757 759 The number of points in the y_direction is scaled accordingly 758 760 759 761 OUTPUT: 760 761  a tuple containing [z_values, xmin, xmax, ymin, ymax] where z_values762 is the evaluation of the map on the specified grid.763 762 763  a tuple containing ``[z_values, xmin, xmax, ymin, ymax]`` 764 where ``z_values`` is the evaluation of the map on the specified grid. 765 764 766 EXAMPLES: 765 767 766 768 General usage:: … … 804 806 pt = xmin + 0.5*xstep + i*xstep + I*(ymin + 0.5*ystep + j*ystep) 805 807 z_values[j, i] = np.dot(p_vector,1/(pre_q_vector  pt)) 806 808 return z_values, xmin, xmax, ymin, ymax 807 808 809 810 809 811 @options(interpolation='catrom') 810 812 def plot_spiderweb(self, spokes=16, circles=4, pts=32, linescale=0.99, 811 813 rgbcolor=[0,0,0], thickness=1, plotjoined=True, withcolor = False, 812 814 plot_points = 200, **options): 813 815 """ 814 Generates a traditional "spiderweb plot" of the Riemann map. Shows 815 what concentric circles and radial lines map to. The radial lines 816 may exhibit erratic behavior near the boundary , if this occurs,816 Generates a traditional "spiderweb plot" of the Riemann map. Shows 817 what concentric circles and radial lines map to. The radial lines 818 may exhibit erratic behavior near the boundary; if this occurs, 817 819 decreasing ``linescale`` may mitigate the problem. 818 819 Note that this method requires significantly more computation for 820 821 Note that this method requires significantly more computation for 820 822 multiply connected domains. 821 823 822 824 INPUT: … … 852 854 853 855  ``withcolor``  boolean (default: ``False``) If ``True``, 854 856 The spiderweb will be overlaid on the basic color plot. 855 857 856 858  ``plot_points``  integer (default: ``200``) the size of the grid in the x direction 857 The number of points in the y_direction is scaled accordingly. 859 The number of points in the y_direction is scaled accordingly. 858 860 Note that very large values can cause this function to run slowly. 859 861  only for multiply connected domains 862 860 863 EXAMPLES: 861 864 862 865 General usage:: … … 886 889 sage: m.plot_spiderweb() 887 890 """ 888 891 cdef int k, i 889 if self.B == 1: #The efficient simply connected 890 edge = self.plot_boundaries(plotjoined=plotjoined, 892 if self.B == 1: #The efficient simply connected 893 edge = self.plot_boundaries(plotjoined=plotjoined, 891 894 rgbcolor=rgbcolor, thickness=thickness) 892 895 circle_list = range(circles) 893 896 theta_array = self.theta_array[0] … … 900 903 temp[i] = self.inverse_riemann_map( 901 904 (k + 1) / (circles + 1.0) * exp(I*i * TWOPI / (2*pts))) 902 905 if plotjoined: 903 circle_list[k] = list_plot(comp_pt(temp, 1), 906 circle_list[k] = list_plot(comp_pt(temp, 1), 904 907 rgbcolor=rgbcolor, thickness=thickness, plotjoined=True) 905 908 else: 906 circle_list[k] = list_plot(comp_pt(temp, 1), 909 circle_list[k] = list_plot(comp_pt(temp, 1), 907 910 rgbcolor=rgbcolor, pointsize=thickness) 908 911 line_list = range(spokes) 909 912 for k in xrange(spokes): … … 931 934 else: 932 935 return edge + sum(circle_list) + sum(line_list) 933 936 else: # The more difficult multiply connected 934 z_values, xmin, xmax, ymin, ymax = self.compute_on_grid([], 937 z_values, xmin, xmax, ymin, ymax = self.compute_on_grid([], 935 938 plot_points) 936 939 xstep = (xmaxxmin)/plot_points 937 940 ystep = (ymaxymin)/plot_points 938 941 dr, dtheta= get_derivatives(z_values, xstep, ystep) # clean later 939 942 940 943 g = Graphics() 941 g.add_primitive(ComplexPlot(complex_to_spiderweb(z_values,dr,dtheta, 942 spokes, circles, rgbcolor,thickness, withcolor), 944 g.add_primitive(ComplexPlot(complex_to_spiderweb(z_values,dr,dtheta, 945 spokes, circles, rgbcolor,thickness, withcolor), 943 946 (xmin, xmax), (ymin, ymax),options)) 944 947 return g + self.plot_boundaries(thickness = thickness) 945 948 946 949 947 950 @options(interpolation='catrom') 948 951 def plot_colored(self, plot_range=[], int plot_points=100, **options): … … 958 961 ``(xmin, xmax, ymin, ymax)``. Declare if you do not want the plot 959 962 to use the default range for the figure. 960 963 961  ``plot_points``  integer (default: ``100``), number of points to 962 plot in the x direction. Points in the y direction are scaled 963 accordingly. Note that very large values can cause this function to 964  ``plot_points``  integer (default: ``100``), number of points to 965 plot in the x direction. Points in the y direction are scaled 966 accordingly. Note that very large values can cause this function to 964 967 run slowly. 965 968 966 969 967 970 EXAMPLES: 968 971 … … 983 986 984 987 To generate the unit circle map, it's helpful to see what the 985 988 colors correspond to:: 986 989 987 990 sage: f(t) = e^(I*t) 988 991 sage: fprime(t) = I*e^(I*t) 989 992 sage: m = Riemann_Map([f], [fprime], 0, 1000) 990 993 sage: m.plot_colored() 991 994 """ 992 z_values, xmin, xmax, ymin, ymax = self.compute_on_grid(plot_range, 995 z_values, xmin, xmax, ymin, ymax = self.compute_on_grid(plot_range, 993 996 plot_points) 994 997 g = Graphics() 995 g.add_primitive(ComplexPlot(complex_to_rgb(z_values), (xmin, xmax), 998 g.add_primitive(ComplexPlot(complex_to_rgb(z_values), (xmin, xmax), 996 999 (ymin, ymax),options)) 997 1000 return g 998 1001 999 1002 cdef comp_pt(clist, loop=True): 1000 1003 """ 1001 Converts a list of complex numbers xderivs = get_derivatives(z_values, xstep, ystep)[0]to the plottable 1004 Utility function to convert the list of complex numbers 1005 ``xderivs = get_derivatives(z_values, xstep, ystep)[0]`` to the plottable 1002 1006 `(x,y)` form. If ``loop=True``, then the first point will be 1003 1007 added as the last, i.e. to plot a closed circle. 1004 1008 … … 1024 1028 if loop: 1025 1029 list2[len(clist)] = list2[0] 1026 1030 return list2 1027 1028 cpdef get_derivatives(np.ndarray[COMPLEX_T, ndim=2] z_values, FLOAT_T xstep,1031 1032 cpdef get_derivatives(np.ndarray[COMPLEX_T, ndim=2] z_values, FLOAT_T xstep, 1029 1033 FLOAT_T ystep): 1030 1034 """ 1031 Computes the r*e^(I*theta) form derivatives from the grid of points. The1032 derivatives are computed using quickanddirty taylor expansion and 1033 assuming analytic y. As such ``get_derivatives`` is primarily intended1034 to be used for comparisions in ``plot_spiderweb`` and not for 1035 Computes the r*e^(I*theta) form of derivatives from the grid of points. The 1036 derivatives are computed using quickanddirty taylor expansion and 1037 assuming analyticity. As such ``get_derivatives`` is primarily intended 1038 to be used for comparisions in ``plot_spiderweb`` and not for 1035 1039 applications that require great precision. 1036 1040 1037 1041 INPUT: 1038  ``z_values``  The values for a complex function evaluated on a grid 1039 in the complexplane, usually from ``compute_on_grid``. 1040 1042 1043  ``z_values``  The values for a complex function evaluated on a grid 1044 in the complex plane, usually from ``compute_on_grid``. 1045 1041 1046  ``xstep``  float, the spacing of the grid points in the real direction 1042 1047 1043 1048 OUTPUT: 1044 1045  A tuple of arrays, [``dr``, ``dtheta``], each array is 2 less in both1049 1050  A tuple of arrays, [``dr``, ``dtheta``], with each array 2 less in both 1046 1051 dimensions than ``z_values`` 1047 ``dr``  the abs of the derivative of the function in the +r direction 1048 ``dtheta``  the rate of accumulation of angle in the +theta direction 1049 1052 1053  ``dr``  the abs of the derivative of the function in the +r direction 1054  ``dtheta``  the rate of accumulation of angle in the +theta direction 1055 1050 1056 EXAMPLES: 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1057 1058 Standard usage with compute_on_grid:: 1059 1060 sage: from sage.calculus.riemann import get_derivatives 1061 sage: f(t) = e^(I*t)  0.5*e^(I*t) 1062 sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(I*t) 1063 sage: m = Riemann_Map([f], [fprime], 0) 1064 sage: data = m.compute_on_grid([],19) 1065 sage: xstep = (data[2]data[1])/19 1066 sage: ystep = (data[4]data[3])/19 1067 sage: dr, dtheta = get_derivatives(data[0],xstep,ystep) 1068 sage: dr[8,8] 1069 0.241... 1070 sage: dtheta[5,5] 1071 5.907... 1072 """ 1067 1073 cdef np.ndarray[COMPLEX_T, ndim=2] xderiv 1068 1074 cdef np.ndarray[FLOAT_T, ndim = 2] dr, dtheta, zabs 1069 1075 imax = len(z_values)2 1070 1076 jmax = len(z_values[0])2 1071 1077 #(f(x+delta)f(xdelta))/2delta 1072 xderiv = (z_values[1:1,2:]z_values[1:1,:2])/(2*xstep) 1073 #b/c the function is analytic, we know it 's abs(derivative) is equal1078 xderiv = (z_values[1:1,2:]z_values[1:1,:2])/(2*xstep) 1079 #b/c the function is analytic, we know its abs(derivative) is equal 1074 1080 #in all directions 1075 1081 dr = np.abs(xderiv) 1076 1082 # the abs(derivative) scaled by distance from origin … … 1078 1084 dtheta = np.divide(dr,zabs) 1079 1085 return dr, dtheta 1080 1086 1081 cpdef complex_to_spiderweb(np.ndarray[COMPLEX_T, ndim = 2] z_values,1082 np.ndarray[FLOAT_T, ndim = 2] dr, np.ndarray[FLOAT_T, ndim = 2] dtheta, 1087 cpdef complex_to_spiderweb(np.ndarray[COMPLEX_T, ndim = 2] z_values, 1088 np.ndarray[FLOAT_T, ndim = 2] dr, np.ndarray[FLOAT_T, ndim = 2] dtheta, 1083 1089 spokes, circles, rgbcolor, thickness, withcolor): 1084 1090 """ 1085 Converts a grid of complex numbers into a matrix containing rgb data 1091 Converts a grid of complex numbers into a matrix containing rgb data 1086 1092 for the Riemann spiderweb plot. 1087 1088 1093 1094 INPUT: 1089 1095 1090 1096  ``z_values``  A grid of complex numbers, as a list of lists. 1091 1092  ``dr``  grid of floats, the r derivative of ``z_values``. 1097 1098  ``dr``  grid of floats, the r derivative of ``z_values``. 1093 1099 Used to determine precision. 1094 1095  ``dtheta``  grid of floats, the theta derivative of ``z_values``. 1100 1101  ``dtheta``  grid of floats, the theta derivative of ``z_values``. 1096 1102 Used to determine precision. 1097 1103 1098 1104  ``spokes``  integer  the number of equally spaced radial lines to plot. 1099 1105 1100 1106  ``circles``  integer  the number of equally spaced circles about the … … 1105 1111 1106 1112  ``thickness``  positive float  the thickness of the lines or points 1107 1113 in the spiderweb. 1108 1114 1109 1115  ``withcolor``  boolean  If ``True`` the spiderweb will be overlaid 1110 1116 on the basic color plot. 1111 1117 1112 1113 1118 OUTPUT: 1119 1114 1120 An `N x M x 3` floating point Numpy array ``X``, where 1115 1121 ``X[i,j]`` is an (r,g,b) tuple. 1116 1122 1117 1123 EXAMPLES:: 1118 1124 1119 1125 sage: from sage.calculus.riemann import complex_to_spiderweb … … 1154 1160 precision = thickness/150.0 1155 1161 imax = len(z_values) 1156 1162 jmax = len(z_values[0]) 1157 cdef np.ndarray[FLOAT_T, ndim=3, mode="c"] rgb 1158 if withcolor: 1163 cdef np.ndarray[FLOAT_T, ndim=3, mode="c"] rgb 1164 if withcolor: 1159 1165 rgb = complex_to_rgb(z_values) 1160 1166 else: 1161 1167 rgb = np.zeros(dtype=FLOAT, shape=(imax, jmax, 3)) … … 1166 1172 circ_radii = [] 1167 1173 if spokes != 0: 1168 1174 # both pi and pi are included 1169 spoke_angles = srange(PI,PI+TWOPI/spokes,TWOPI/spokes) 1175 spoke_angles = srange(PI,PI+TWOPI/spokes,TWOPI/spokes) 1170 1176 else: 1171 1177 spoke_angles = [] 1172 1178 for i in xrange(imax2): # the d arrays are 1 smaller on each side … … 1178 1184 darg = dtheta[i,j] 1179 1185 #points that change too rapidly are presumed to be borders 1180 1186 #points that are too small are presumed to be outside 1181 if darg < DMAX and mag > MMIN: 1187 if darg < DMAX and mag > MMIN: 1182 1188 for target in circ_radii: 1183 1189 if abs(mag  target)/dmag < precision: 1184 1190 rgb[i+1,j+1] = rgbcolor 1185 1191 break 1186 1192 for target in spoke_angles: 1187 if abs(arg  target)/darg < precision: 1193 if abs(arg  target)/darg < precision: 1188 1194 rgb[i+1,j+1] = rgbcolor 1189 1195 break 1190 1196 return rgb 1191 1197 1192 1198 1193 1199 cpdef complex_to_rgb(np.ndarray[COMPLEX_T, ndim = 2] z_values): 1194 1200 r""" … … 1294 1300 rgb[i, j, 2] = b 1295 1301 sig_off() 1296 1302 return rgb 1297 1303 1298 1304 cpdef analytic_boundary(FLOAT_T t, int n): 1299 1305 """ 1300 Provides an exact (for n = infinity) riemann boundary1301 correspondence for the ellipse with axes 1 + epsilon and 1  epsilon. The 1302 boundary is therefore given by e^(I*t)+epsilon*e^(I*t). It is primarily 1303 useful for testing the accuracy of the numerical Riemann Map.1304 1306 Provides an exact (for n = infinity) Riemann boundary 1307 correspondence for the ellipse with axes 1 + epsilon and 1  epsilon. The 1308 boundary is therefore given by e^(I*t)+epsilon*e^(I*t). It is primarily 1309 useful for testing the accuracy of the numerical :class:`Riemann_Map`. 1310 1305 1311 INPUT: 1306 1307  ``t``  The boundary parameter, from 0 to two*pi1308 1309  ``n``  Integer, The number of terms to include. 10 is fairly accurate,1310 1311 1312 1313  ``t``  The boundary parameter, from 0 to 2*pi 1314 1315  ``n``  integer  the number of terms to include. 1316 10 is fairly accurate, 20 is very accurate. 1317 1312 1318 OUTPUT: 1313 1314 A theta value from 0 to 2*pi, corresponding to the point on the 1319 1320 A theta value from 0 to 2*pi, corresponding to the point on the 1315 1321 circle e^(I*theta) 1316 1322 1317 1323 TESTS: 1318 1319 Checking the accuracy for different n values::1320 1324 1325 Checking the accuracy of this function for different n values:: 1326 1321 1327 sage: from sage.calculus.riemann import analytic_boundary 1322 1328 sage: t100 = analytic_boundary(pi/2,100) 1323 1329 sage: abs(analytic_boundary(pi/2,10)  t100) < 10^8 1324 1330 True 1325 1331 sage: abs(analytic_boundary(pi/2,20)  t100) < 10^15 1326 1332 True 1327 1333 1328 1334 Using this to check the accuracy of the Riemann_Map boundary:: 1329 1335 1330 1336 sage: f(t) = e^(I*t)+.3*e^(I*t) 1331 1337 sage: fp(t) = I*e^(I*t)I*.3*e^(I*t) 1332 1338 sage: m = Riemann_Map([f], [fp],0,200) … … 1341 1347 result += (2*(1)**i/i)*(.3**i/(1+.3**(2*i)))*sin(2*i*t) 1342 1348 return result 1343 1349 1344 1350 1345 1351 1346 1352 cpdef cauchy_kernel(t, args): 1347 1353 """ 1348 Intermediate function for the integration in ``analytic_interior``.1349 1354 Intermediate function for the integration in :meth:`~Riemann_Map.analytic_interior`. 1355 1350 1356 INPUT: 1351 1357 1352 1358  ``t``  The boundary parameter, meant to be integrated over 1353 1359 1354 1360  ``args``  a tuple containing: 1355 1356  ``z``  Complex,the point to be mapped.1357 1358  ``n``  Integer, The number of terms to include. 10 is fairly1359 1360 1361  ``part``  will return the real ('r'), imaginary ('i') or1362 1363 1361 1362  ``z``  complex  the point to be mapped. 1363 1364  ``n``  integer  the number of terms to include. 1365 10 is fairly accurate, 20 is very accurate. 1366 1367  ``part``  will return the real ('r'), imaginary ('i') or 1368 complex ('c') value of the kernel 1369 1364 1370 TESTS: 1365 Primarily tested implicitly by analytic_interior 1366 1367 Simple test::1368 1371 1372 This is primarily tested implicitly by :meth:`~Riemann_Map.analytic_interior`. 1373 Here is a simple test:: 1374 1369 1375 sage: from sage.calculus.riemann import cauchy_kernel 1370 1376 sage: cauchy_kernel(.5,(.1+.2*I, 10,'c')) 1371 1377 (0.584136405997...+0.5948650858950...j) 1372 1373 1374 1378 """ 1375 cdef COMPLEX_T result 1379 cdef COMPLEX_T result 1376 1380 cdef COMPLEX_T z = args[0] 1377 1381 cdef int n = args[1] 1378 1382 part = args[2] … … 1385 1389 elif part == 'i': 1386 1390 return result.imag 1387 1391 else: return None 1388 1392 1389 1393 cpdef analytic_interior(COMPLEX_T z, int n): 1390 1394 """ 1391 Provides a nearly exact compuation of the Riemann Map of an interior 1392 point of the ellipse with axes 1 + epsilon and 1  epsilon. It is 1395 Provides a nearly exact compuation of the Riemann Map of an interior 1396 point of the ellipse with axes 1 + epsilon and 1  epsilon. It is 1393 1397 primarily useful for testing the accuracy of the numerical Riemann Map. 1394 1398 1395 1399 INPUT: 1396 1397  ``z``  Complex,the point to be mapped.1398 1399  ``n``  Integer, The number of terms to include. 10 is fairly accurate,1400 1400 1401  ``z``  complex  the point to be mapped. 1402 1403  ``n``  integer  the number of terms to include. 1404 10 is fairly accurate, 20 is very accurate. 1401 1405 1402 1406 TESTS: 1403 1404 Testing the accuracy of Riemann_Map::1405 1407 1408 Testing the accuracy of :class:`Riemann_Map`:: 1409 1406 1410 sage: from sage.calculus.riemann import analytic_interior 1407 1411 sage: f(t) = e^(I*t)+.3*e^(I*t) 1408 1412 sage: fp(t) = I*e^(I*t)I*.3*e^(I*t) … … 1415 1419 """ 1416 1420 # evaluates the cauchy integral of the boundary, split into the real 1417 1421 # and imaginary results because numerical_integral can't handle complex data. 1418 rp = 1/(TWOPI)*numerical_integral(cauchy_kernel,0,2*pi, 1422 rp = 1/(TWOPI)*numerical_integral(cauchy_kernel,0,2*pi, 1419 1423 params = [z,n,'i'])[0] 1420 ip = 1/(TWOPI*I)*numerical_integral(cauchy_kernel,0,2*pi, 1424 ip = 1/(TWOPI*I)*numerical_integral(cauchy_kernel,0,2*pi, 1421 1425 params = [z,n,'r'])[0] 1422 1426 return rp + ip