1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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
41
42 public class GnuplotComplexPlotter {
43
44
45 private Complex n;
46
47
48 private Complex m;
49
50
51 private FieldJacobiElliptic<Complex> jacobi;
52
53
54 private List<Predefined> functions;
55
56
57 private File output;
58
59
60 private DomainColoring coloring;
61
62
63 private int width;
64
65
66 private int height;
67
68
69 private double xMin;
70
71
72 private double xMax;
73
74
75 private double yMin;
76
77
78 private double yMax;
79
80
81 private double zMax;
82
83
84 private double viewXRot;
85
86
87 private double viewZRot;
88
89
90 private boolean use3D;
91
92
93 private int maxEval;
94
95
96 final ComplexUnivariateIntegrator integrator;
97
98
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
124
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
232
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
249
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
313 private interface Evaluator {
314
315
316
317
318
319
320 Complex value(GnuplotComplexPlotter plotter, Complex z);
321
322 }
323
324
325 private enum Predefined {
326
327
328 id((plotter, z) -> z),
329
330
331 sn((plotter, z) -> plotter.jacobi.valuesN(z).sn()),
332
333
334 cn((plotter, z) -> plotter.jacobi.valuesN(z).cn()),
335
336
337 dn((plotter, z) -> plotter.jacobi.valuesN(z).dn()),
338
339
340 cs((plotter, z) -> plotter.jacobi.valuesS(z).cs()),
341
342
343 ds((plotter, z) -> plotter.jacobi.valuesS(z).ds()),
344
345
346 ns((plotter, z) -> plotter.jacobi.valuesS(z).ns()),
347
348
349 dc((plotter, z) -> plotter.jacobi.valuesC(z).dc()),
350
351
352 nc((plotter, z) -> plotter.jacobi.valuesC(z).nc()),
353
354
355 sc((plotter, z) -> plotter.jacobi.valuesC(z).sc()),
356
357
358 nd((plotter, z) -> plotter.jacobi.valuesD(z).nd()),
359
360
361 sd((plotter, z) -> plotter.jacobi.valuesD(z).sd()),
362
363
364 cd((plotter, z) -> plotter.jacobi.valuesD(z).cd()),
365
366
367 arcsn((plotter, z) -> plotter.jacobi.arcsn(z)),
368
369
370 arccn((plotter, z) -> plotter.jacobi.arccn(z)),
371
372
373 arcdn((plotter, z) -> plotter.jacobi.arcdn(z)),
374
375
376 arccs((plotter, z) -> plotter.jacobi.arccs(z)),
377
378
379 arcds((plotter, z) -> plotter.jacobi.arcds(z)),
380
381
382 arcns((plotter, z) -> plotter.jacobi.arcns(z)),
383
384
385 arcdc((plotter, z) -> plotter.jacobi.arcdc(z)),
386
387
388 arcnc((plotter, z) -> plotter.jacobi.arcnc(z)),
389
390
391 arcsc((plotter, z) -> plotter.jacobi.arcsc(z)),
392
393
394 arcnd((plotter, z) -> plotter.jacobi.arcnd(z)),
395
396
397 arcsd((plotter, z) -> plotter.jacobi.arcsd(z)),
398
399
400 arccd((plotter, z) -> plotter.jacobi.arccd(z)),
401
402
403 K((plotter, z) -> LegendreEllipticIntegral.bigK(z)),
404
405
406 KPrime((plotter, z) -> LegendreEllipticIntegral.bigKPrime(z)),
407
408
409 Fzm((plotter, z) -> LegendreEllipticIntegral.bigF(z, plotter.m)),
410
411
412 integratedFzm((plotter, z) -> LegendreEllipticIntegral.bigF(z, plotter.m, plotter.integrator, plotter.maxEval)),
413
414
415 E((plotter, z) -> LegendreEllipticIntegral.bigE(z)),
416
417
418 Ezm((plotter, z) -> LegendreEllipticIntegral.bigE(z, plotter.m)),
419
420
421 integratedEzm((plotter, z) -> LegendreEllipticIntegral.bigE(z, plotter.m, plotter.integrator, plotter.maxEval)),
422
423
424 Pi((plotter, z) -> LegendreEllipticIntegral.bigPi(plotter.n, z)),
425
426
427 Pizm((plotter, z) -> LegendreEllipticIntegral.bigPi(plotter.n, z, plotter.m)),
428
429
430 integratedPizm((plotter, z) -> LegendreEllipticIntegral.bigPi(plotter.n, z, plotter.m, plotter.integrator, plotter.maxEval)),
431
432
433 sin((plotter, z) -> FastMath.sin(z)),
434
435
436 cos((plotter, z) -> FastMath.cos(z)),
437
438
439 tan((plotter, z) -> FastMath.tan(z)),
440
441
442 asin((plotter, z) -> FastMath.asin(z)),
443
444
445 acos((plotter, z) -> FastMath.acos(z)),
446
447
448 atan((plotter, z) -> FastMath.atan(z)),
449
450
451 sinh((plotter, z) -> FastMath.sinh(z)),
452
453
454 cosh((plotter, z) -> FastMath.cosh(z)),
455
456
457 tanh((plotter, z) -> FastMath.tanh(z)),
458
459
460 asinh((plotter, z) -> FastMath.asinh(z)),
461
462
463 acosh((plotter, z) -> FastMath.acosh(z)),
464
465
466 atanh((plotter, z) -> FastMath.atanh(z));
467
468
469 private final Evaluator evaluator;
470
471
472
473
474 Predefined(final Evaluator evaluator) {
475 this.evaluator = evaluator;
476 }
477
478
479
480
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 }