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