Implement SphericalPlot3D

.gitignore

@@ -25,7 +25,7 @@ ChangeLog.orig
 ChangeLog.rej
 Documents/
 Homepage/
-Test/
+#Test/
 _Copies_/
 _Database_/
 build/

mathics/builtin/drawing/plot.py

@@ -22,7 +22,7 @@
 from mathics.core.evaluation import Evaluation
 from mathics.core.expression import Expression
 from mathics.core.list import ListExpression
-from mathics.core.symbols import Symbol
+from mathics.core.symbols import Symbol, SymbolList
 from mathics.core.systemsymbols import (
     SymbolAll,
     SymbolAutomatic,
@@ -413,14 +413,27 @@ class PlotOptions:
     plot_points: list
     maxdepth: int
 
-    def __init__(self, builtin, range_exprs, options, dim, evaluation):
+    def __init__(self, builtin, functions, range_exprs, options, dim, evaluation):
         def error(*args, **kwargs):
             evaluation.message(builtin.get_name(), *args, **kwargs)
             raise ValueError()
 
+        # convert functions to list of lists of exprs
+        def to_list(expr):
+            if isinstance(expr, Expression) and expr.head is SymbolList:
+                return [to_list(e) for e in expr.elements]
+            else:
+                return expr
+
+        functions = to_list(functions)
+        self.functions = functions if isinstance(functions, list) else [functions]
+
         # plot ranges of the form {x,xmin,xmax} etc. (returns Symbol)
         self.ranges = []
-        for range_expr in range_exprs:
+        for i, range_expr in enumerate(range_exprs):
+            if isinstance(range_expr, Symbol) and hasattr(builtin, "default_ranges"):
+                self.ranges.append([range_expr, *builtin.default_ranges[i]])
+                continue
             if not range_expr.has_form("List", 3):
                 error("invrange", range_expr)
             if not isinstance(range_expr.elements[0], Symbol):
@@ -471,18 +484,25 @@ def error(*args, **kwargs):
         self.exclusions = exclusions
 
         # Mesh option (returns Symbol)
-        mesh = builtin.get_option(options, "Mesh", evaluation)
-        if mesh not in (SymbolNone, SymbolFull, SymbolAll):
+        mesh = builtin.get_option(options, "Mesh", evaluation).to_python(
+            preserve_symbols=True
+        )
+        if isinstance(mesh, (list, tuple)) and all(isinstance(m, int) for m in mesh):
+            self.mesh = mesh
+        elif mesh not in (SymbolNone, SymbolFull, SymbolAll):
             evaluation.message("Mesh", "ilevels", mesh)
-            mesh = SymbolFull
-        self.mesh = mesh
+            self.mesh = SymbolFull
+        else:
+            self.mesh = mesh
 
         # PlotPoints option (returns Symbol)
         plot_points_option = builtin.get_option(options, "PlotPoints", evaluation)
         pp = plot_points_option.to_python(preserve_symbols=True)
-        npp = len(self.ranges)
-        if builtin.get_name() in ("System`ComplexPlot3D", "System`ComplexPlot"):
-            npp = 2
+        npp = (
+            builtin.num_plot_points
+            if hasattr(builtin, "num_plot_points")
+            else len(self.ranges)
+        )
         if pp == SymbolNone:
             pp = None
         else:

mathics/builtin/drawing/plot_plot.py

@@ -74,7 +74,9 @@ def eval(self, functions, ranges, evaluation: Evaluation, options: dict):
         # parse options, bailing out if anything is wrong
         try:
             ranges = ranges.elements if ranges.head is SymbolSequence else [ranges]
-            plot_options = plot.PlotOptions(self, ranges, options, 2, evaluation)
+            plot_options = plot.PlotOptions(
+                self, functions, ranges, options, 2, evaluation
+            )
         except ValueError:
             return None
 
@@ -84,10 +86,15 @@ def eval(self, functions, ranges, evaluation: Evaluation, options: dict):
         apply_function = self.apply_function
         if not plot.use_vectorized_plot:
             apply_function = lru_cache(apply_function)
+        plot_options.apply_function = apply_function
 
-        # additional options specific to this class
+        # TODO: PlotOptions has already regularized .functions to be a list
+        # (of lists) of functions, used by the _Plot3d builtins.
+        # But _Plot builtins still need to be reworked to use it,
+        # so we still use the old mechanism here.
         plot_options.functions = self.get_functions_param(functions)
-        plot_options.apply_function = apply_function
+
+        # additional options specific to this class
         plot_options.use_log_scale = self.use_log_scale
         plot_options.expect_list = self.expect_list
         if plot_options.plot_points is None:

mathics/builtin/drawing/plot_plot3d.py

@@ -35,6 +35,7 @@ class _Plot3D(Builtin):
     # Check for correct number of args
     eval_error = Builtin.generic_argument_error
     expected_args = 3
+    is_cartesian = True
 
     messages = {
         "invmaxrec": (
@@ -99,22 +100,27 @@ def eval(
         try:
             dim = 3 if self.graphics_class is Graphics3D else 2
             ranges = ranges.elements if ranges.head is SymbolSequence else [ranges]
-            plot_options = plot.PlotOptions(self, ranges, options, dim, evaluation)
+            plot_options = plot.PlotOptions(
+                self, functions, ranges, options, dim, evaluation
+            )
         except ValueError:
             return None
 
-        # TODO: consult many_functions variable set by subclass and error
-        # if many_functions is False but multiple are supplied
-        if functions.has_form("List", None):
-            plot_options.functions = functions.elements
-        else:
-            plot_options.functions = [functions]
-
-        # supply default value
+        # supply default value for PlotPoints
         if plot_options.plot_points is None:
-            default_plot_points = (200, 200) if plot.use_vectorized_plot else (7, 7)
+            if isinstance(self, ParametricPlot3D) and len(plot_options.ranges) == 1:
+                # ParametricPlot3D with one independent variable generating a curve
+                default_plot_points = (1000,)
+            elif plot.use_vectorized_plot:
+                default_plot_points = (200, 200)
+            else:
+                default_plot_points = (7, 7)
             plot_options.plot_points = default_plot_points
 
+        # supply apply_function which knows how to take the plot parameters
+        # and produce xs, ys, and zs
+        plot_options.apply_function = self.apply_function
+
         # subclass must set eval_function and graphics_class
         eval_function = plot.get_plot_eval_function(self.__class__)
         with np.errstate(all="ignore"):  # suppress numpy warnings
@@ -123,13 +129,16 @@ def eval(
             return
 
         # now we have a list of length dim
-        # handle Automatic ~ {xmin,xmax} etc.
+        # handle Automatic ~ {xmin,xmax} etc., but only if is_cartesion: the independent variables are x and y
         # TODO: dowstream consumers might be happier if we used data range where applicable
-        for i, (pr, r) in enumerate(zip(plot_options.plot_range, plot_options.ranges)):
-            # TODO: this treats Automatic and Full as the same, which isn't quite right
-            if isinstance(pr, (str, Symbol)) and not isinstance(r[1], complex):
-                # extract {xmin,xmax} from {x,xmin,xmax}
-                plot_options.plot_range[i] = r[1:]
+        if self.is_cartesian:
+            for i, (pr, r) in enumerate(
+                zip(plot_options.plot_range, plot_options.ranges)
+            ):
+                # TODO: this treats Automatic and Full as the same, which isn't quite right
+                if isinstance(pr, (str, Symbol)) and not isinstance(r[1], complex):
+                    # extract {xmin,xmax} from {x,xmin,xmax}
+                    plot_options.plot_range[i] = r[1:]
 
         # unpythonize and update PlotRange option
         options[str(SymbolPlotRange)] = to_mathics_list(*plot_options.plot_range)
@@ -140,6 +149,10 @@ def eval(
         )
         return graphics_expr
 
+    def apply_function(self, function, names, us, vs):
+        parms = {str(names[0]): us, str(names[1]): vs}
+        return us, vs, function(**parms)
+
 
 class ComplexPlot3D(_Plot3D):
     """
@@ -161,8 +174,13 @@ class ComplexPlot3D(_Plot3D):
     options = _Plot3D.options3d | {"Mesh": "None"}
 
     many_functions = True
+    num_plot_points = 2  # different from number of ranges
     graphics_class = Graphics3D
 
+    def apply_function(self, function, names, us, vs):
+        parms = {str(names[0]): us + vs * 1j}
+        return us, vs, function(**parms)
+
 
 class ComplexPlot(_Plot3D):
     """
@@ -184,8 +202,13 @@ class ComplexPlot(_Plot3D):
     options = _Plot3D.options2d
 
     many_functions = False
+    num_plot_points = 2  # different from number of ranges
     graphics_class = Graphics
 
+    def apply_function(self, function, names, us, vs):
+        parms = {str(names[0]): us + vs * 1j}
+        return us, vs, function(**parms)
+
 
 class ContourPlot(_Plot3D):
     """
@@ -244,6 +267,47 @@ class DensityPlot(_Plot3D):
     graphics_class = Graphics
 
 
+class ParametricPlot3D(_Plot3D):
+    """
+    <url>:Parametric equation: https://en.wikipedia.org/wiki/Parametric_equation</url>
+    <url>:WMA link: https://reference.wolfram.com/language/ref/ParametricPlot3D.html</url>
+    <dl>
+      <dt>'ParametricPlot3D'[${x(u,v), y(u,v), z(u,v)}$, {$u$, $u_{min}$, $u_{max}$}, {$v$, $v_{min}$, $v_{max}$}]
+      <dd>creates a three-dimensional surface using the functions $x$, $y$, $z$ over the specified ranges for parameters $u$ and $v$.
+
+      <dt>'ParametricPlot3D'[${x(u), y(u), z(u)}$, {$u$, $u_{min}$, $u_{max}$}]
+      <dd>creates a three-dimensional space curve using the functions $x$, $y$, $z$ over the specified range for parameter $u$.
+
+          See <url>:Drawing Option and Option Values:
+    /doc/reference-of-built-in-symbols/graphics-and-drawing/drawing-options-and-option-values
+    </url> for a list of Plot options.
+    </dl>
+
+    >> ParametricPlot3D[{Sin[t] + 2 Sin[2 t], Cos[t] - 2 Cos[2 t], -Sin[3 t]}, {t, 0, 2 Pi}]
+     = ...
+
+    A function of a single parameter $t$ generates a trefoil knot.
+
+    >> ParametricPlot3D[{(2 + Cos[v]) Cos[u], (2 + Cos[v]) Sin[u], Sin[v]}, {u, 0, 2 Pi}, {v, 0, 2 Pi}]
+     = ...
+
+    A function of two parameters $u$ and $v$ generates a torus.
+
+    """
+
+    summary_text = "plot a parametric surface or curve in three dimensions"
+    expected_args = 3
+    options = _Plot3D.options3d
+
+    is_cartesian = False
+    many_functions = True
+    graphics_class = Graphics3D
+
+    def apply_function(self, functions, names, *parms):
+        parms = {str(n): p for n, p in zip(names, parms)}
+        return [f(**parms) for f in functions]
+
+
 class Plot3D(_Plot3D):
     """
     <url>:WMA link: https://reference.wolfram.com/language/ref/Plot3D.html</url>
@@ -279,3 +343,44 @@ class Plot3D(_Plot3D):
 
     many_functions = True
     graphics_class = Graphics3D
+
+
+class SphericalPlot3D(_Plot3D):
+    """
+    <url>:Spherical coordinate system: https://en.wikipedia.org/wiki/Spherical_coordinate_system</url>
+    <url>:WMA link: https://reference.wolfram.com/language/ref/SphericalPlot3D.html</url>
+    <dl>
+      <dt>'SphericalPlot3D'[$r(theta, phi)$, {$theta$, $theta_{min}$, $theta_{max}$}, {$phi$, $phi_{min}$, $phi_{max}$}]
+      <dd>creates a three-dimensional surface at radius $r(theta, phi)$ for spherical angles $theta$ and $phi$ over the specified ranges
+
+      <dt>'SphericalPlot3D'[$r(theta, phi)$, $theta$, $phi$]
+      <dd>creates a three-dimensional surface at radius $r(theta, phi)$ for spherical angles $theta$ and $phi$
+          in the ranges $0 < theta < pi$ and $0 < phi < 2pi$ covering the entire sphere
+
+          See <url>:Drawing Option and Option Values:
+    /doc/reference-of-built-in-symbols/graphics-and-drawing/drawing-options-and-option-values
+    </url> for a list of Plot options.
+    </dl>
+
+    >> SphericalPlot3D[1 + 0.4 Abs[SphericalHarmonicY[10, 4, theta, phi]], theta, phi]
+     = ...
+
+    Spherical harmonics are the canonical use case for spherical plots.
+
+
+    """
+
+    summary_text = "produce a surface plot functions spherical angles theta and phi"
+    expected_args = 3
+    options = _Plot3D.options3d | {"BoxRatios": "{1,1,1}"}
+
+    is_cartesian = False
+    many_functions = True
+    graphics_class = Graphics3D
+    default_ranges = [[0, np.pi], [0, 2 * np.pi]]
+
+    def apply_function(self, function, names, θ, φ):
+        parms = {names[0]: θ, names[1]: φ}
+        r = function(**parms)
+        x, y, z = r * np.sin(θ) * np.cos(φ), r * np.sin(θ) * np.sin(φ), r * np.cos(θ)
+        return x, y, z

mathics/core/systemsymbols.py

@@ -26,6 +26,7 @@
 SymbolAborted = Symbol("System`$Aborted")
 SymbolAbs = Symbol("System`Abs")
 SymbolAbsoluteTime = Symbol("AbsoluteTime")
+SymbolAbsoluteThickness = Symbol("System`AbsoluteThickness")
 SymbolAccuracy = Symbol("System`Accuracy")
 SymbolAlignmentPoint = Symbol("System`AlignmentPoint")
 SymbolAll = Symbol("System`All")

mathics/core/util.py

@@ -134,7 +134,10 @@ def print_expression_tree(
     if file is None:
         file = sys.stdout
 
-    if isinstance(expr, Symbol):
+    if isinstance(expr, (tuple, list)):
+        for e in expr:
+            print_expression_tree(e, indent, marker, file, approximate)
+    elif isinstance(expr, Symbol):
         print(f"{indent}{marker(expr)}{expr}", file=file)
     elif not hasattr(expr, "elements"):
         if isinstance(expr, MachineReal) and approximate:

mathics/eval/drawing/plot3d.py

@@ -513,3 +513,17 @@ def eval_ContourPlot(
     evaluation: Evaluation,
 ):
     return None
+
+
+def eval_ParametricPlot3D(
+    plot_options,
+    evaluation: Evaluation,
+):
+    return None
+
+
+def eval_SphericalPlot3D(
+    plot_options,
+    evaluation: Evaluation,
+):
+    return None