View Javadoc
1   /*
2    * Licensed to the Hipparchus project under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /** Plotter for complex functions using gnuplot.
19   */
20  package org.hipparchus.samples.complex;
21  
22  import java.io.File;
23  import java.io.IOException;
24  import java.io.PrintStream;
25  import java.nio.charset.StandardCharsets;
26  import java.util.ArrayList;
27  import java.util.List;
28  import java.util.Locale;
29  
30  import org.hipparchus.analysis.integration.IterativeLegendreGaussIntegrator;
31  import org.hipparchus.complex.Complex;
32  import org.hipparchus.complex.ComplexUnivariateIntegrator;
33  import org.hipparchus.exception.MathIllegalStateException;
34  import org.hipparchus.special.elliptic.jacobi.FieldJacobiElliptic;
35  import org.hipparchus.special.elliptic.jacobi.JacobiEllipticBuilder;
36  import org.hipparchus.special.elliptic.legendre.LegendreEllipticIntegral;
37  import org.hipparchus.util.FastMath;
38  import org.hipparchus.util.RyuDouble;
39  
40  /** Program plotting complex functions with domain coloring.
41   */
42  public class GnuplotComplexPlotter {
43  
44      /** Elliptic integrals characteristic. */
45      private Complex n;
46  
47      /** Elliptic integrals parameter. */
48      private Complex m;
49  
50      /** Jacobi functions computer. */
51      private FieldJacobiElliptic<Complex> jacobi;
52  
53      /** Functions to plot. */
54      private List<Predefined> functions;
55  
56      /** Output directory (display to terminal if null). */
57      private File output;
58  
59      /** Domain coloring. */
60      private DomainColoring coloring;
61  
62      /** Output width. */
63      private int width;
64  
65      /** Output height. */
66      private int height;
67  
68      /** Min x. */
69      private double xMin;
70  
71      /** Max x. */
72      private double xMax;
73  
74      /** Min y. */
75      private double yMin;
76  
77      /** Max y. */
78      private double yMax;
79  
80      /** Max z. */
81      private double zMax;
82  
83      /** View X rotation. */
84      private double viewXRot;
85  
86      /** View Z rotation. */
87      private double viewZRot;
88  
89      /** Indicator for 3D surfaces. */
90      private boolean use3D;
91  
92      /** Maximum number of integrands evaluations for each integral evaluation. */
93      private int maxEval;
94  
95      /** Integrator for numerically integrated elliptical integrals. */
96      final ComplexUnivariateIntegrator integrator;
97  
98      /** Default constructor.
99       */
100     private GnuplotComplexPlotter() {
101         n          = new Complex(3.4, 1.3);
102         m          = new Complex(0.64);
103         jacobi     = JacobiEllipticBuilder.build(m);
104         functions  = new ArrayList<>();
105         output     = null;
106         width      = 800;
107         height     = 800;
108         coloring   = new SawToothPhaseModuleValue(1.0, 0.7, 1.0, 15);
109         xMin       = -7;
110         xMax       = +7;
111         yMin       = -7;
112         yMax       = +7;
113         zMax       = +7;
114         viewXRot   = 60;
115         viewZRot   = 30;
116         use3D      = false;
117         maxEval    = 100000;
118         integrator = new ComplexUnivariateIntegrator(new IterativeLegendreGaussIntegrator(24,
119                                                                                           1.0e-6,
120                                                                                           1.0e-6));
121     }
122 
123     /** Main program.
124      * @param args program arguments
125      */
126     public static void main(String[] args) {
127         final GnuplotComplexPlotter plotter = new GnuplotComplexPlotter();
128         try {
129             for (int i = 0; i < args.length; ++i) {
130                 switch (args[i]) {
131                     case "--help" :
132                         usage(0);
133                         break;
134                     case "--output-dir" :
135                         plotter.output = new File(args[++i]);
136                         if (!(plotter.output.exists() && plotter.output.isDirectory() && plotter.output.canWrite())) {
137                             System.err.format(Locale.US, "cannot generate output file in %s%n",
138                                               plotter.output.getAbsolutePath());
139                             System.exit(1);
140                         }
141                         break;
142                     case "--width" :
143                         plotter.width = Integer.parseInt(args[++i]);
144                         break;
145                     case "--height" :
146                         plotter.height = Integer.parseInt(args[++i]);
147                         break;
148                     case "--color" :
149                         switch (args[++i]) {
150                             case "classical" :
151                                 plotter.coloring = new ContinuousModuleValue(1.0);
152                                 break;
153                             case "enhanced-module" :
154                                 plotter.coloring = new SawToothModuleValue(1.0);
155                                 break;
156                             case "enhanced-phase-module" :
157                                 plotter.coloring = new SawToothPhaseModuleValue(1.0, 0.7, 1.0, 15);
158                                 break;
159                             default :
160                                 usage(1);
161                         }
162                         break;
163                     case "--3d" :
164                         plotter.use3D = true;
165                         break;
166                     case "--view" :
167                         plotter.viewXRot = Double.parseDouble(args[++i]);
168                         plotter.viewZRot = Double.parseDouble(args[++i]);
169                         break;
170                     case "--xmin" :
171                         plotter.xMin = Double.parseDouble(args[++i]);
172                         break;
173                     case "--xmax" :
174                         plotter.xMax = Double.parseDouble(args[++i]);
175                         break;
176                     case "--ymin" :
177                         plotter.yMin = Double.parseDouble(args[++i]);
178                         break;
179                     case "--ymax" :
180                         plotter.yMax = Double.parseDouble(args[++i]);
181                         break;
182                     case "--zmax" :
183                         plotter.zMax = Double.parseDouble(args[++i]);
184                         break;
185                     case "--m" : {
186                         plotter.m      = new Complex(Double.parseDouble(args[++i]), Double.parseDouble(args[++i]));
187                         plotter.jacobi = JacobiEllipticBuilder.build(plotter.m);
188                         break;
189                     }
190                     case "--n" : {
191                         plotter.n      = new Complex(Double.parseDouble(args[++i]), Double.parseDouble(args[++i]));
192                         break;
193                     }
194                     case "--maxeval" : {
195                         plotter.maxEval = Integer.parseInt(args[++i]);
196                         break;
197                     }
198                     case "--function" :
199                         try {
200                             plotter.functions.add(Predefined.valueOf(args[++i]));
201                         } catch (IllegalArgumentException iae) {
202                             System.err.format(Locale.US, "unknown function %s, known functions:%n", args[i]);
203                             for (final Predefined predefined : Predefined.values()) {
204                                 System.err.format(Locale.US, " %s", predefined.name());
205                             }
206                             System.err.format(Locale.US, "%n");
207                             System.exit(1);
208                         }
209                         break;
210                     default :
211                         usage(1);
212                 }
213             }
214             if (plotter.functions.isEmpty()) {
215                 usage(1);
216             }
217 
218             plotter.plot();
219 
220         } catch (IndexOutOfBoundsException iobe) {
221             usage(1);
222         } catch (IOException ioe) {
223             System.err.println(ioe.getLocalizedMessage());
224             System.exit(1);
225         }
226 
227         System.exit(0);
228 
229     }
230 
231     /** Display usage.
232      * @param status exit code
233      */
234     private static void usage(final int status) {
235         System.err.println("usage: java org.hipparchus.samples.complex.GnuplotComplexPlotter" +
236                            " [--help]" +
237                            " [--output-dir directory]" +
238                            " [--3d]" +
239                            " [--view xRot zRot]" +
240                            " [--color {classical|enhanced-module|enhanced-phase-module}]" +
241                            " [--xmin xMin] [--xmax xMax] [--ymin yMin] [--ymax yMax] [--zmax zMax]" +
242                            " [--m mRe mIm] [--n nRe nIm] [--maxeval maxEval]" +
243                            " --function {id|sn|cn|dn|cs|...|sin|cos|...} [--function ...]");
244         System.exit(status);
245 
246     }
247 
248     /** Plot the function.
249      * @throws IOException if gnuplot process cannot be run
250      */
251     public void plot() throws IOException {
252         for (final Predefined predefined : functions) {
253             final ProcessBuilder pb = new ProcessBuilder("gnuplot").
254                             redirectOutput(ProcessBuilder.Redirect.INHERIT).
255                             redirectError(ProcessBuilder.Redirect.INHERIT);
256             pb.environment().remove("XDG_SESSION_TYPE");
257             final Process gnuplot = pb.start();
258             try (PrintStream out = new PrintStream(gnuplot.getOutputStream(), false, StandardCharsets.UTF_8.name())) {
259                 if (output == null) {
260                     out.format(Locale.US, "set terminal qt size %d, %d title 'complex plotter'%n", width, height);
261                 } else {
262                     out.format(Locale.US, "set terminal pngcairo size %d, %d%n", width, height);
263                     out.format(Locale.US, "set output '%s'%n", new File(output, predefined.name() + ".png").getAbsolutePath());
264                 }
265                 out.format(Locale.US, "set xrange [%f : %f]%n", xMin, xMax);
266                 out.format(Locale.US, "set yrange [%f : %f]%n", yMin, yMax);
267                 if (use3D) {
268                     out.format(Locale.US, "set zrange [%f : %f]%n", 0.0, zMax);
269                 }
270                 out.format(Locale.US, "set xlabel 'Re(z)'%n");
271                 out.format(Locale.US, "set ylabel 'Im(z)'%n");
272                 out.format(Locale.US, "set key off%n");
273                 out.format(Locale.US, "unset colorbox%n");
274                 out.format(Locale.US, "set title '%s'%n", predefined.title(m));
275                 out.format(Locale.US, "$data <<EOD%n");
276 
277                 for (int i = 0; i < width; ++i) {
278                     final double x = xMin + i * (xMax - xMin) / (width - 1);
279                     for (int j = 0; j < height; ++j) {
280                         final double y = yMin + j * (yMax - yMin) / (height - 1);
281                         Complex z  = Complex.valueOf(x, y);
282                         Complex fz;
283                         try {
284                             fz = predefined.evaluator.value(this, z);
285                         } catch (MathIllegalStateException e) {
286                             fz = Complex.NaN;
287                         }
288                         out.format(Locale.US, "%12.9f %12.9f %12.9f %12.9f %12.9f %12.9f%n",
289                                    z.getRealPart(), z.getImaginaryPart(), fz.norm(),
290                                    coloring.hue(fz), coloring.saturation(fz), coloring.value(fz));
291                     }
292                     out.format(Locale.US, "%n");
293                 }
294                 out.format(Locale.US, "EOD%n");
295                 if (use3D) {
296                     out.format(Locale.US, "set view %f, %f%n", viewXRot, viewZRot);
297                     out.format(Locale.US, "splot $data using 1:2:3:(hsv2rgb($4,$5,$6)) with pm3d lc rgb variable%n");
298                 } else {
299                     out.format(Locale.US, "set view map scale 1%n");
300                     out.format(Locale.US, "splot $data using 1:2:(hsv2rgb($4,$5,$6)) with pm3d lc rgb variable%n");
301                 }
302                 if (output == null) {
303                     out.format(Locale.US, "pause mouse close%n");
304                 } else {
305                     System.out.format(Locale.US, "output written to %s%n",
306                                       new File(output, predefined.name() + ".png").getAbsolutePath());
307                 }
308             }
309         }
310     }
311 
312     /** Interface for evaluating complex functions. */
313     private interface Evaluator {
314 
315         /** Evaluate complex function.
316          * @param plotter associated plotter
317          * @param z free variable
318          * @return value of the complex function
319          */
320         Complex value(GnuplotComplexPlotter plotter, Complex z);
321 
322     }
323 
324     /** Predefined complex functions for plotting. */
325     private enum Predefined {
326 
327         /** id. */
328         id((plotter, z)             -> z),
329 
330         /** sn. */
331         sn((plotter, z)             -> plotter.jacobi.valuesN(z).sn()),
332 
333         /** cn. */
334         cn((plotter, z)             -> plotter.jacobi.valuesN(z).cn()),
335 
336         /** dn. */
337         dn((plotter, z)             -> plotter.jacobi.valuesN(z).dn()),
338 
339         /** cs. */
340         cs((plotter, z)             -> plotter.jacobi.valuesS(z).cs()),
341 
342         /** ds. */
343         ds((plotter, z)             -> plotter.jacobi.valuesS(z).ds()),
344 
345         /** ns. */
346         ns((plotter, z)             -> plotter.jacobi.valuesS(z).ns()),
347 
348         /** dc. */
349         dc((plotter, z)             -> plotter.jacobi.valuesC(z).dc()),
350 
351         /** nc. */
352         nc((plotter, z)             -> plotter.jacobi.valuesC(z).nc()),
353 
354         /** sc. */
355         sc((plotter, z)             -> plotter.jacobi.valuesC(z).sc()),
356 
357         /** nd. */
358         nd((plotter, z)             -> plotter.jacobi.valuesD(z).nd()),
359 
360         /** sd. */
361         sd((plotter, z)             -> plotter.jacobi.valuesD(z).sd()),
362 
363         /** cd. */
364         cd((plotter, z)             -> plotter.jacobi.valuesD(z).cd()),
365 
366         /** arcsn. */
367         arcsn((plotter, z)          -> plotter.jacobi.arcsn(z)),
368 
369         /** arccn. */
370         arccn((plotter, z)          -> plotter.jacobi.arccn(z)),
371 
372         /** arcdn. */
373         arcdn((plotter, z)          -> plotter.jacobi.arcdn(z)),
374 
375         /** arccs. */
376         arccs((plotter, z)          -> plotter.jacobi.arccs(z)),
377 
378         /** arcds. */
379         arcds((plotter, z)          -> plotter.jacobi.arcds(z)),
380 
381         /** arcns. */
382         arcns((plotter, z)          -> plotter.jacobi.arcns(z)),
383 
384         /** arcdc. */
385         arcdc((plotter, z)          -> plotter.jacobi.arcdc(z)),
386 
387         /** arcnc. */
388         arcnc((plotter, z)          -> plotter.jacobi.arcnc(z)),
389 
390         /** arcsc. */
391         arcsc((plotter, z)          -> plotter.jacobi.arcsc(z)),
392 
393         /** arcnd. */
394         arcnd((plotter, z)          -> plotter.jacobi.arcnd(z)),
395 
396         /** arcsd. */
397         arcsd((plotter, z)          -> plotter.jacobi.arcsd(z)),
398 
399         /** arccd. */
400         arccd((plotter, z)          -> plotter.jacobi.arccd(z)),
401 
402         /** K. */
403         K((plotter, z)              -> LegendreEllipticIntegral.bigK(z)),
404 
405         /** KPrime. */
406         KPrime((plotter, z)         -> LegendreEllipticIntegral.bigKPrime(z)),
407 
408         /** Fzm. */
409         Fzm((plotter, z)            -> LegendreEllipticIntegral.bigF(z, plotter.m)),
410 
411         /** integratedFzm. */
412         integratedFzm((plotter, z)  -> LegendreEllipticIntegral.bigF(z, plotter.m, plotter.integrator, plotter.maxEval)),
413 
414         /** E. */
415         E((plotter, z)              -> LegendreEllipticIntegral.bigE(z)),
416 
417         /** Ezm. */
418         Ezm((plotter, z)            -> LegendreEllipticIntegral.bigE(z, plotter.m)),
419 
420         /** integratedEzm. */
421         integratedEzm((plotter, z)  -> LegendreEllipticIntegral.bigE(z, plotter.m, plotter.integrator, plotter.maxEval)),
422 
423         /** Pi. */
424         Pi((plotter, z)             -> LegendreEllipticIntegral.bigPi(plotter.n, z)),
425 
426         /** Pizm. */
427         Pizm((plotter, z)           -> LegendreEllipticIntegral.bigPi(plotter.n, z, plotter.m)),
428 
429         /** integratedPizm. */
430         integratedPizm((plotter, z) -> LegendreEllipticIntegral.bigPi(plotter.n, z, plotter.m, plotter.integrator, plotter.maxEval)),
431 
432         /** sin. */
433         sin((plotter, z)            -> FastMath.sin(z)),
434 
435         /** cos. */
436         cos((plotter, z)            -> FastMath.cos(z)),
437 
438         /** tan. */
439         tan((plotter, z)            -> FastMath.tan(z)),
440 
441         /** asin. */
442         asin((plotter, z)           -> FastMath.asin(z)),
443 
444         /** acos. */
445         acos((plotter, z)           -> FastMath.acos(z)),
446 
447         /** atan. */
448         atan((plotter, z)           -> FastMath.atan(z)),
449 
450         /** sinh. */
451         sinh((plotter, z)           -> FastMath.sinh(z)),
452 
453         /** cosh. */
454         cosh((plotter, z)           -> FastMath.cosh(z)),
455 
456         /** tanh. */
457         tanh((plotter, z)           -> FastMath.tanh(z)),
458 
459         /** asinh. */
460         asinh((plotter, z)          -> FastMath.asinh(z)),
461 
462         /** acosh. */
463         acosh((plotter, z)          -> FastMath.acosh(z)),
464 
465         /** atanh. */
466         atanh((plotter, z)          -> FastMath.atanh(z));
467 
468         /** Function evaluator. */
469         private final Evaluator evaluator;
470 
471         /** Simple constructor.
472          * @param evaluator function evaluator
473          */
474         Predefined(final Evaluator evaluator) {
475             this.evaluator = evaluator;
476         }
477 
478         /** Get plot title.
479          * @param ep elliptic parameter
480          * @return plot title
481          */
482         public String title(final Complex ep) {
483             if (name().endsWith("zm")) {
484                 return name().substring(0, name().length() - 2) +
485                        "(z, m = " +
486                        RyuDouble.doubleToString(ep.getRealPart()) +
487                        (ep.getImaginary() >= 0 ? " + " : " - ") +
488                        RyuDouble.doubleToString(FastMath.abs(ep.getImaginaryPart())) +
489                        "i)";
490             } else {
491                 return name() + "(z)";
492             }
493         }
494 
495     }
496 
497 }