23425: Bug fixes, edited examples

src/sage/dynamics/complex_dynamics/mandel_julia.py

@@ -25,8 +25,8 @@
 #*****************************************************************************
 
 from __future__ import absolute_import, division
-from sage.dynamics.complex_dynamics.mandel_julia_helper import fast_mandelbrot_plot, \
- fast_external_ray, convert_to_pixels, get_line
+from sage.dynamics.complex_dynamics.mandel_julia_helper import \
+ fast_mandelbrot_plot, fast_external_ray, convert_to_pixels, get_line
 from sagenb.notebook.interact import interact, slider, input_box, color_selector
 from sage.plot.colors import Color
 from sage.repl.image import Image
@@ -151,32 +151,11 @@ def external_ray(theta, **kwds):
 
     kwds:
 
-    Adjusting the image of the Mandelbrot set:
+    - ``image`` -- 24-bit RGB image (optional - default: None) user specified image of Mandelbrot set.
 
-    - ``x_center`` -- double (optional - default: ``-1.0``), Real part of center point.
-
-    - ``y_center`` -- double (optional - default: ``0.0``), Imaginary part of center point.
-
-    - ``image_width`` -- double (optional - default: ``4.0``), width of image in the complex plane.
-
-    - ``max_iteration`` -- long (optional - default: ``500``), maximum number of iterations the map ``Q_c(z)``.
-
-    - ``pixel_count`` -- long (optional - default: ``500``), side length of image in number of pixels.
-
-    - ``base_color`` -- RGB color (optional - default: ``[40, 40, 40]``) color used to determine the coloring of set.
-
-    - ``iteration_level`` -- long (optional - default: 1) number of iterations between each color level.
-
-    - ``number_of_colors`` -- long (optional - default: 30) number of colors used to plot image.
-
-    - ``image`` -- 24-bit RGB image (optional - default: ``mandelbrot_plot()`` with parameters given above)  image of Mandelbrot set.
-
-
-    Adjusting the External Ray(s):
+    - ``D`` -- long (optional - default: ``25``) depth of the approximation. As ``D`` increases, the external ray gets closer to the boundary of the Mandelbrot set. If the ray doesn't reach the boundary of the Mandelbrot set, increase ``D``.
 
-    - ``D`` -- long (optional - default: ``25``) depth of the approximation. As ``D`` increases, the external ray gets closer to the boundary of the Mandelbrot set.
-
-    - ``S`` -- long (optional - default: ``10``) sharpness of the approximation. Adjusts the number of points used to approximate the external ray (number of points is equal to ``S*D``).
+    - ``S`` -- long (optional - default: ``10``) sharpness of the approximation. Adjusts the number of points used to approximate the external ray (number of points is equal to ``S*D``). If ray looks jagged, increase ``S``.
 
     - ``R`` -- long (optional - default: ``100``) radial parameter. If ``R`` is sufficiently large, the external ray reaches enough close to infinity.
 
@@ -195,7 +174,12 @@ def external_ray(theta, **kwds):
 
     ::
 
-        sage: external_ray([0, 2/7, 0.1, pi/2]) # long time
+        sage: external_ray(0.6, ray_color=[255, 0, 0])
+        500x500px 24-bit RGB image
+
+    ::
+
+        sage: external_ray([0, 0.2, 0.4, 0.7]) # long time
         500x500px 24-bit RGB image
 
     ::
@@ -205,7 +189,7 @@ def external_ray(theta, **kwds):
 
     WARNING:
 
-    If you are passing in an image, make sure you specify the
+    If you are passing in an image, make sure you specify
     which parameters to use when drawing the external ray.
     For example, the following is incorrect::
 
@@ -224,29 +208,19 @@ def external_ray(theta, **kwds):
     The ``copy()`` function for bitmap images needs to be implemented in Sage.
     """
 
-    # Keywords for adjusting image of Mandelbrot set
-    x_0 = kwds.pop("x_center", -1)
-    y_0 = kwds.pop("y_center", 0)
-    plot_width = kwds.pop("image_width", 4)
-    max_it = kwds.pop("max_iteration", 500)
-    pixel_width = kwds.pop("pixel_count", 500)
-    mandel_color = kwds.pop("base_color", [40, 40, 40])
-    it_lvl = kwds.pop("iteration_level", 2)
-    color_num = kwds.pop("number_of_colors", 30)
-
-    # Key words for adjusting External Ray(s)
-    depth = kwds.pop("D", 25)
-    sharpness = kwds.pop("S", 10)
-    radial_parameter = kwds.pop("R", 100)
-    precision = kwds.pop("prec", 300)
+    x_0 = kwds.get("x_center", -1)
+    y_0 = kwds.get("y_center", 0)
+    plot_width = kwds.get("image_width", 4)
+    pixel_width = kwds.get("pixel_count", 500)
+    depth = kwds.get("D", 25)
+    sharpness = kwds.get("S", 10)
+    radial_parameter = kwds.get("R", 100)
+    precision = kwds.get("prec", 300)
     precision = max(precision, -log(pixel_width * 0.001, 2).round() + 10)
-    ray_color = kwds.pop("ray_color", [255]*3)
-
-    # User can also use their own image of the Mandelbrot set
-    image = kwds.pop("image", mandelbrot_plot(x_center=x_0, y_center= y_0,
-     image_width=plot_width, max_iteration=max_it, pixel_count=pixel_width,
-     base_color=mandel_color, iteration_level=it_lvl, number_of_colors=color_num))
-
+    ray_color = kwds.get("ray_color", [255]*3)
+    image = kwds.get("image", None)
+    if image is None:
+        image = mandelbrot_plot(**kwds)
 
     # Make a copy of the bitmap image.
     # M = copy(image)
@@ -261,26 +235,36 @@ def external_ray(theta, **kwds):
     if type(theta) != list:
         theta = [theta]
 
+    # Check if theta is in the invterval [0,1]
+    for angle in theta:
+        if angle < 0 or angle > 1:
+            raise \
+            ValueError("values for theta must be in the closed interval [0,1].")
+
     # Loop through each value for theta in list and plot the external ray.
     for angle in theta:
         E = fast_external_ray(angle, D=depth, S=sharpness, R=radial_parameter,
-         prec=precision, pixel_count=pixel_width)
+         prec=precision, image_width=plot_width, pixel_count=pixel_width)
 
         # Convert points to pixel coordinates.
         pixel_list = convert_to_pixels(E, x_0, y_0, plot_width, pixel_width)
 
         # Find the pixels between points in pixel_list.
         extra_points = []
         for i in range(len(pixel_list) - 1):
-            for j in get_line(pixel_list[i], pixel_list[i+1]):
-                extra_points.append(j)
+            if min(pixel_list[i+1]) >= 0 and max(pixel_list[i+1]) < pixel_width:
+                for j in get_line(pixel_list[i], pixel_list[i+1]):
+                    extra_points.append(j)
 
         # Add these points to pixel_list to fill in gaps in the ray.
         pixel_list += extra_points
 
+        # Remove duplicates from list.
+        pixel_list = list(set(pixel_list))
+
         # Check if point is in window and if it is, plot it on the image to
         # create an external ray.
         for k in pixel_list:
-            if max(k[0], k[1]) < pixel_width and min(k[0], k[1]) >= 0:
+            if max(k) < pixel_width and min(k) >= 0:
                 pixel[int(k[0]), int(k[1])] = tuple(ray_color)
     return M

src/sage/dynamics/complex_dynamics/mandel_julia_helper.pyx

@@ -25,7 +25,7 @@ from sage.repl.image import Image
 from copy import copy
 from cysignals.signals cimport sig_check
 from sage.rings.complex_field import ComplexField
-from sage.functions.log import exp
+from sage.functions.log import exp, log
 from sage.symbolic.constants import pi
 
 def fast_mandelbrot_plot(double x_center, double y_center, double image_width,
@@ -135,7 +135,8 @@ def fast_mandelbrot_plot(double x_center, double y_center, double image_width,
                     pixel[col,row] = color_list[-1]
     return M
 
-def fast_external_ray(double theta, long D = 30, long S = 10, long R = 100, pixel_count = 500, long prec = 300):
+cpdef fast_external_ray(double theta, long D=30, long S=10, long R=100,
+ long pixel_count=500, double image_width=4, long prec=300):
     r"""
     Returns a list of points that approximate the external ray for a given angle.
 
@@ -151,6 +152,8 @@ def fast_external_ray(double theta, long D = 30, long S = 10, long R = 100, pixe
 
     - ``pixel_count`` -- long (optional - default: ``500``) side length of image in number of pixels.
 
+    - ``image_width`` -- double (optional - default: ``4``) width of the image in the complex plane.
+
     - ``prec`` -- long (optional - default: ``300``) specifies the bits of precision used by the Complex Field when using Newton's method to compute points on the external ray.
 
     OUTPUT:
@@ -163,46 +166,49 @@ def fast_external_ray(double theta, long D = 30, long S = 10, long R = 100, pixe
         sage: fast_external_ray(0,S=1,D=1)
         [(100.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000,
           0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000),
-         (9.51254777713729104959838878130540251731872558593750000000000000000000000000000000000000000,
+         (9.51254777713729174697578576623132297117784691109499464854806785133621315075854778426714908,
           0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000)]
 
+
     ::
 
         sage: from sage.dynamics.complex_dynamics.mandel_julia_helper import fast_external_ray
         sage: fast_external_ray(1/3,S=1,D=1)
         [(-49.9999999999999786837179271969944238662719726562500000000000000000000000000000000000000000,
           86.6025403784438765342201804742217063903808593750000000000000000000000000000000000000000000),
-         (-5.50628047023173028406972662196494638919830322265625000000000000000000000000000000000000000,
-          8.64947510053972479227013536728918552398681640625000000000000000000000000000000000000000000)]
+         (-5.50628047023173006234970878097113901879832542655926629309001652388544528575532346900138516,
+          8.64947510053972513843999918917106032664030380426885745306040284140385975750462108180377187)]
 
     ::
 
         sage: from sage.dynamics.complex_dynamics.mandel_julia_helper import fast_external_ray
         sage: fast_external_ray(0.75234,S=1,D=1)
         [(1.47021239172637052661229972727596759796142578125000000000000000000000000000000000000000000,
           -99.9891917935294287644865107722580432891845703125000000000000000000000000000000000000000000),
-         (-0.352790406744857509835355813265778124332427978515625000000000000000000000000000000000000000,
-          -9.98646630765023601838947797659784555435180664062500000000000000000000000000000000000000000)]
+         (-0.352790406744857508500937144524776555433184352559852962308757189778284058275081335121601384,
+          -9.98646630765023514178761177926164047797465369576787921409326037870837930920646860774032363)]
     """
 
     cdef:
         CF = ComplexField(prec)
         PI = CF.pi()
         I = CF.gen()
-        c_0, r_m, t_m, temp_c, C_k, D_k, old_c, x, y
+        c_0, r_m, t_m, temp_c, C_k, D_k, old_c, x, y, dist
         int k, j, t
         double difference, m
         double error = pixel_count * 0.0001
 
+        double pixel_width = image_width / pixel_count
+
         # initialize list with c_0
         c_list = [CF(R*exp(2*PI*I*theta))]
 
     # Loop through each subinterval and approximate point on external ray.
     for k in range(1,D+1):
         for j in range(1,S+1):
             m = (k-1)*S + j
-            r_m = R**(2**(-m/S))
-            t_m = r_m**(2**k) * exp(2*PI*I*theta * 2**k)
+            r_m = CF(R**(2**(-m/S)))
+            t_m = CF(r_m**(2**k) * exp(2*PI*I*theta * 2**k))
             temp_c = c_list[-1]
             difference = error
 
@@ -211,12 +217,18 @@ def fast_external_ray(double theta, long D = 30, long S = 10, long R = 100, pixe
                 sig_check()
                 old_c = temp_c
                 # Recursive formula for iterates of q(z) = z^2 + c
-                C_k, D_k = old_c, 1.0
+                C_k, D_k = CF(old_c), CF(1)
                 for t in range(k):
-                    C_k, D_k = C_k**2 + old_c, 2*D_k*C_k + 1.0
+                    C_k, D_k = C_k**2 + old_c, CF(2)*D_k*C_k + CF(1)
                 temp_c = old_c - (C_k - t_m) / D_k   # Newton map
                 difference = abs(old_c) - abs(temp_c)
+
+            dist = (2*C_k.abs()*(C_k.abs()).log()) / D_k.abs()
+            if dist < pixel_width:
+                break
             c_list.append(CF(temp_c))
+        if dist < pixel_width:
+            break
 
     # Convert Complex Field elements into tuples.
     for k in range(len(c_list)):
@@ -225,7 +237,8 @@ def fast_external_ray(double theta, long D = 30, long S = 10, long R = 100, pixe
 
     return c_list
 
-def convert_to_pixels(point_list, double x_0, double y_0, double width, long number_of_pixels):
+cpdef convert_to_pixels(point_list, double x_0, double y_0, double width,
+ long number_of_pixels):
     r"""
     Converts cartesian coordinates to pixels within a specified window.
 
@@ -253,7 +266,7 @@ def convert_to_pixels(point_list, double x_0, double y_0, double width, long num
     """
     cdef:
         k, pixel_list, x_corner, y_corner, step_size
-        int x_pixel, y_pixel
+        long x_pixel, y_pixel
     pixel_list = []
 
     # Compute top left corner of window and step size
@@ -269,7 +282,7 @@ def convert_to_pixels(point_list, double x_0, double y_0, double width, long num
         pixel_list.append((x_pixel, y_pixel))
     return pixel_list
 
-def get_line(start, end):
+cpdef get_line(start, end):
     r"""
     Produces a list of pixel coordinates approximating a line from a starting
     point to an ending point using the Bresenham's Line Algorithm.
@@ -327,7 +340,7 @@ def get_line(start, end):
     dx, dy = x2 - x1, y2 - y1
 
     # Calculate error
-    error = int(dx / 2.0) #TODO: Do I want true division here?
+    error = int(dx / 2.0)
     ystep = 1 if y1 < y2 else -1
 
     # Iterate over bounding box generating points between start and end