View Javadoc
1   package org.djutils.complex;
2   
3   import static org.junit.Assert.assertEquals;
4   
5   import org.djutils.base.AngleUtil;
6   import org.junit.Test;
7   
8   /**
9    * TestComplexMath.java. <br>
10   * <br>
11   * Copyright (c) 2020-2020 Delft University of Technology, Jaffalaan 5, 2628 BX Delft, the Netherlands. All rights reserved. See
12   * for project information <a href="https://djutils.org" target="_blank"> https://djutils.org</a>. The DJUTILS project is
13   * distributed under a three-clause BSD-style license, which can be found at
14   * <a href="https://djutils.org/docs/license.html" target="_blank"> https://djutils.org/docs/license.html</a>. <br>
15   * @author <a href="https://www.tudelft.nl/averbraeck">Alexander Verbraeck</a>
16   * @author <a href="https://www.tudelft.nl/pknoppers">Peter Knoppers</a>
17   */
18  public class TestComplexMath
19  {
20  
21      /**
22       * Test the square root function.
23       */
24      @Test
25      public void testSqrt()
26      {
27          // Start by testing the bloody obvious
28          assertEquals("square root of -1 is I", Complex.I, ComplexMath.sqrt(new Complex(-1)));
29          Complex in = new Complex(0, -4);
30          Complex c = ComplexMath.sqrt(in);
31          assertEquals("square root of " + in + " norm", 2, c.norm(), 0.000001);
32          assertEquals("square root of " + in + " phi", -Math.PI / 4, c.phi(), 0.000001);
33          double[] values = new double[] {0, 1, 0.01, 100, Math.PI, -Math.E};
34          for (double re : values)
35          {
36              for (double im : values)
37              {
38                  in = new Complex(re, im);
39                  c = ComplexMath.sqrt(in);
40                  assertEquals("square root of " + in + " norm", Math.sqrt(in.norm()), c.norm(), 0.0001);
41                  assertEquals("square root of " + in + " phi", in.phi() / 2, c.phi(), 0.0000001);
42                  Complex c2 = c.times(c);
43                  assertEquals("square of square root re", in.re, c2.re, 0.0001);
44                  assertEquals("square of square root im", in.im, c2.im, 0.0001);
45              }
46          }
47      }
48  
49      /**
50       * Test the cube root function.
51       */
52      @Test
53      public void testCbrt()
54      {
55          double[] values = new double[] {0, 1, 0.01, 100, Math.PI, -Math.E};
56          for (double re : values)
57          {
58              for (double im : values)
59              {
60                  Complex in = new Complex(re, im);
61                  Complex c = ComplexMath.cbrt(in);
62                  assertEquals("cube root of " + in + " norm", Math.cbrt(in.norm()), c.norm(), 0.0001);
63                  assertEquals("cube root of " + in + " phi", in.phi() / 3, c.phi(), 0.0000001);
64                  Complex c3 = c.times(c).times(c);
65                  assertEquals("cube of cube root re", in.re, c3.re, 0.0001);
66                  assertEquals("cube of cube root im", in.im, c3.im, 0.0001);
67              }
68          }
69      }
70  
71      /**
72       * Test the exponential function.
73       */
74      @Test
75      public void testExp()
76      {
77          assertEquals("exp of 1 is e; re", Math.E, ComplexMath.exp(Complex.ONE).re, 0.000001);
78          assertEquals("exp of 1 is e; im", 0, ComplexMath.exp(Complex.ONE).im, 0.000001);
79          for (double re : new double[] {0, 1, Math.PI, -Math.E, 10, -10})
80          {
81              for (double im : new double[] {0, 0.1, Math.PI / 2, Math.PI, 5, -1, -Math.E, -50})
82              {
83                  Complex in = new Complex(re, im);
84                  Complex out = ComplexMath.exp(in);
85                  assertEquals("exp(" + in + ") re", Math.exp(re) * Math.cos(im), out.re, 0.01);
86                  assertEquals("exp(" + in + ") im", Math.exp(re) * Math.sin(im), out.im, 0.01);
87              }
88          }
89      }
90  
91      /**
92       * Test the natural logarithm function.
93       */
94      @Test
95      public void testLog()
96      {
97          assertEquals("ln(ONE) is ZERO", Complex.ZERO, ComplexMath.ln(Complex.ONE));
98          Complex in = new Complex(Math.E);
99          Complex out = ComplexMath.ln(in);
100         assertEquals("ln(e) is ONE re", 1, out.re, 0.00000001);
101         assertEquals("ln(e) is ONE im", 0, out.im, 0.00000001);
102         for (double re : new double[] {0, 1, Math.PI, 10, -Math.E, -10})
103         {
104             for (double im : new double[] {0, 0.1, Math.PI / 2, Math.PI, 5, -1, -Math.E, -50})
105             {
106                 in = new Complex(re, im);
107                 out = ComplexMath.ln(in);
108                 assertEquals("ln(" + in + ") re", Math.log(in.norm()), out.re, 0.01);
109                 assertEquals("ln(" + in + ") im", Math.atan2(im, re), out.im, 0.00001);
110             }
111         }
112     }
113 
114     /**
115      * Test the sine, cosine and tangent functions.
116      */
117     @Test
118     public void testSinCosTan()
119     {
120         assertEquals("sin(ZERO) is ZERO", Complex.ZERO, ComplexMath.sin(Complex.ZERO));
121         Complex c = ComplexMath.cos(Complex.ZERO);
122         assertEquals("cos(ZERO) is ONE: re", 1, c.re, 0.00001);
123         assertEquals("cos(ZERO) is ONE: im", 0, c.im, 0.00001);
124         assertEquals("tan(ZERO) is ZERO", Complex.ZERO, ComplexMath.tan(Complex.ZERO));
125         double[] values = new double[] {0, 1, Math.PI, 10, -Math.E, -10};
126         for (double re : values)
127         {
128             for (double im : values)
129             {
130                 Complex in = new Complex(re, im);
131                 Complex sin = ComplexMath.sin(in);
132                 assertEquals("sin(" + in + ") re", Math.sin(re) * Math.cosh(im), sin.re, 0.0001);
133                 assertEquals("sin(" + in + ") im", Math.cos(re) * Math.sinh(im), sin.im, 0.0001);
134                 Complex cos = ComplexMath.cos(in);
135                 assertEquals("cos(" + in + ") re", Math.cos(re) * Math.cosh(im), cos.re, 0.0001);
136                 assertEquals("cos(" + in + ") im", -Math.sin(re) * Math.sinh(im), cos.im, 0.0001);
137                 Complex tan = ComplexMath.tan(in);
138                 Complex div = sin.divideBy(cos);
139                 assertEquals("div norm", sin.norm() / cos.norm(), div.norm(), 0.00001);
140                 assertEquals("div phi", AngleUtil.normalizeAroundZero(sin.phi() - cos.phi()), div.phi(), 0.0001);
141                 assertEquals("tan(" + in + ") re", div.re, tan.re, 0.0001);
142                 assertEquals("tan(" + in + ") im", div.im, tan.im, 0.0001);
143                 Complex sin2plusCos2 = sin.times(sin).plus(cos.times(cos));
144                 assertEquals("sin^2 + cos^2 re", 1, sin2plusCos2.re, 0.00001);
145                 assertEquals("sin^2 + cos^2 im", 0, sin2plusCos2.im, 0.00001);
146             }
147         }
148     }
149 
150     /**
151      * Test the sinh, cosh, tanh functions.
152      */
153     @Test
154     public void testSinhCoshTanH()
155     {
156         double[] values = new double[] {0, 1, Math.PI, 10, -Math.E, -10};
157         for (double re : values)
158         {
159             for (double im : values)
160             {
161                 Complex in = new Complex(re, im);
162                 Complex sinh = ComplexMath.sinh(in);
163                 Complex cosh = ComplexMath.cosh(in);
164                 Complex tanh = ComplexMath.tanh(in);
165 
166                 // System.out.println(" in=" + printComplex(in) + "\ntanh=" + printComplex(tanh));
167                 assertEquals("sinh re", Math.sinh(re) * Math.cos(im), sinh.re, 0.0001);
168                 assertEquals("sinh im", Math.cosh(re) * Math.sin(im), sinh.im, 0.0001);
169                 assertEquals("cosh re", Math.cosh(re) * Math.cos(im), cosh.re, 0.0001);
170                 assertEquals("cosh im", Math.sinh(re) * Math.sin(im), cosh.im, 0.0001);
171                 assertEquals("tanh re", Math.sinh(2 * re) / (Math.cosh(2 * re) + Math.cos(2 * im)), tanh.re, 0.0001);
172                 assertEquals("tanh im", Math.sin(2 * im) / (Math.cosh(2 * re) + Math.cos(2 * im)), tanh.im, 0.0001);
173                 // Alternate way to compute tanh
174                 Complex alternateTanh = sinh.divideBy(cosh);
175                 assertEquals("alternate tanh re", tanh.re, alternateTanh.re, 0.0001);
176                 assertEquals("alternate tanh im", tanh.im, alternateTanh.im, 0.0001);
177                 if (im == 0)
178                 {
179                     // Extra checks
180                     assertEquals("sinh of real re", Math.sinh(re), sinh.re, 0.0001);
181                     assertEquals("sinh of real im", 0, sinh.im, 0.0001);
182                     assertEquals("cosh of real re", Math.cosh(re), cosh.re, 0.0001);
183                     assertEquals("cosh of real im", 0, cosh.im, 0.0001);
184                     assertEquals("tanh of real re", Math.tanh(re), tanh.re, 0.0001);
185                     assertEquals("tahh of real im", 0, tanh.im, 0.0001);
186                 }
187             }
188         }
189     }
190 
191     /**
192      * Test the asin, acos and atan functions.
193      */
194     @Test
195     public void testAsinAcosAtan()
196     {
197         double[] values = new double[] {0, 0.2, 0, 8, 1, -1, -0.2, -0.8, Math.PI, 10, -Math.E, -10};
198         for (double re : values)
199         {
200             for (double im : values)
201             {
202                 Complex in = new Complex(re, im);
203                 Complex asin = ComplexMath.asin(in);
204                 Complex acos = ComplexMath.acos(in);
205                 Complex atan = ComplexMath.atan(in);
206                 // This is a lousy test; we only verify that asin(sin(asin(z)) roughly equals asin(z)
207                 Complex asinOfSinOfAsin = ComplexMath.asin(ComplexMath.sin(asin));
208                 assertEquals("asin re", asinOfSinOfAsin.re, asin.re, 0.0001);
209                 assertEquals("asin im", asinOfSinOfAsin.im, asin.im, 0.0001);
210                 Complex acosOfCosOfAcos = ComplexMath.acos(ComplexMath.cos(acos));
211                 assertEquals("acos re", acosOfCosOfAcos.re, acos.re, 0.0001);
212                 assertEquals("acos im", acosOfCosOfAcos.im, acos.im, 0.0001);
213                 Complex atanOfTanOfAtan = ComplexMath.atan(ComplexMath.tan(atan));
214                 if (Math.abs(atan.re) < 100)
215                 {
216                     assertEquals("atan re", atanOfTanOfAtan.re, atan.re, 0.0001);
217                     assertEquals("atan im", atanOfTanOfAtan.im, atan.im, 0.0001);
218                 }
219                 if (im == 0 && re >= -1 && re <= 1)
220                 {
221                     // Extra checks
222                     assertEquals("asin of real in range -1, 1 re", Math.asin(re), asin.re, 0.00001);
223                     assertEquals("asin of real in range -1, 1 im", 0, asin.im, 0.00001);
224                     assertEquals("acos of real in range -1, 1 re", Math.acos(re), acos.re, 0.00001);
225                     assertEquals("acos of real in range -1, 1 im", 0, acos.im, 0.00001);
226                 }
227                 else if (im == 0)
228                 {
229                     assertEquals("atan of real re", Math.atan(re), atan.re, 0.00001);
230                     assertEquals("atan of real, 1 im", 0, atan.im, 0.00001);
231 
232                 }
233             }
234         }
235     }
236 
237     /**
238      * Test the asinh function.
239      */
240     @Test
241     public void testAsinh()
242     {
243         double[] values = new double[] {0, 0.2, 0, 8, 1, -1, -0.2, -0.8, Math.PI, 10, -Math.E, -10};
244         for (double re : values)
245         {
246             for (double im : values)
247             {
248                 Complex in = new Complex(re, im);
249                 Complex asinh = ComplexMath.asinh(in);
250                 Complex acosh = ComplexMath.acosh(in);
251                 Complex atanh = ComplexMath.atanh(in);
252                 // This is a lousy test; we only verify that asinh(sinh(asinh(z)) roughly equals asinh(z)
253                 Complex asinhOfSinhOfAsinh = ComplexMath.asinh(ComplexMath.sinh(asinh));
254                 assertEquals("asinh re", asinhOfSinhOfAsinh.re, asinh.re, 0.0001);
255                 assertEquals("asinh im", asinhOfSinhOfAsinh.im, asinh.im, 0.0001);
256                 Complex acoshOfCoshOfAcosh = ComplexMath.acosh(ComplexMath.cosh(acosh));
257                 if (im != 0 || re > 1.0)
258                 {
259                     // acosh is unstable around im == 0 && re <= 1.0; see <a
260                     // href="https://mathworld.wolfram.com/InverseHyperbolicCosine.html">Wolfram mathWorld: Inverse Hyperbolic
261                     // Cosine<//a> so we can't use this test there.
262                     assertEquals("acosh re", acoshOfCoshOfAcosh.re, acosh.re, 0.0001);
263                     assertEquals("acosh im", acoshOfCoshOfAcosh.im, acosh.im, 0.0001);
264                 }
265                 Complex atanhOfTanhOfAtanh = ComplexMath.atanh(ComplexMath.tanh(atanh));
266                 if (im != 0 || re > -1.0 && re < 1.0)
267                 {
268                     // atanh is unstable around im == 0 && re <= -1 && re >= 1; see <a
269                     // "https://mathworld.wolfram.com/InverseHyperbolicTangent.html">Wolfram mathWorld: Inverse Hyperbolic
270                     // Tangent</a>, so we can't use this test there.
271                     // System.out.println(" in=" + printComplex(in) + "\natanh=" + printComplex(atanh));
272                     if (im != 1 && im != -1 || re != 0)
273                     {
274                         // Also unstable around i and minus i as the atan function is unstable around -1
275                         assertEquals("atanh re", atanhOfTanhOfAtanh.re, atanh.re, 0.0001);
276                         assertEquals("atanh im", atanhOfTanhOfAtanh.im, atanh.im, 0.0001);
277                     }
278                 }
279 
280                 if (im == 0)
281                 {
282                     // Extra checks
283                     assertEquals("asinh of real re", doubleAsinh(re), asinh.re, 0.00001);
284                     assertEquals("asinh of real im", 0, asinh.im, 0.00001);
285                     if (re >= 1.0)
286                     {
287                         assertEquals("acosh of real re", doubleAcosh(re), acosh.re, 0.00001);
288                         assertEquals("acosh of real im", 0, acosh.im, 0.00001);
289                     }
290                     if (re > -1.0 && re < 1.0)
291                     {
292                         assertEquals("atanh of real re", (Math.log(1 + re) - Math.log(1 - re)) / 2, atanh.re, 0.00001);
293                         assertEquals("acosh of real im", 0, atanh.im, 0.00001);
294                     }
295                 }
296             }
297         }
298     }
299 
300     /**
301      * Copied from <a href="https://forgetcode.com/java/1746-asinh-return-the-hyperbolic-sine-of-value-as-a-argument">Forget
302      * Code asinh</a>.
303      * @param x double; the argument
304      * @return double; the inverse hyperbolic cosine of x
305      */
306     public static double doubleAsinh(final double x)
307     {
308         return Math.log(x + Math.sqrt(x * x + 1.0));
309     }
310 
311     /**
312      * Copied from <a href="https://forgetcode.com/Java/1747-acosh-Return-the-hyperbolic-Cosine-of-value-as-a-Argument">Forget
313      * Code acosh</a>.
314      * @param x double; the argument
315      * @return double; the inverse hyperbolic cosine of x
316      */
317     public static double doubleAcosh(final double x)
318     {
319         return Math.log(x + Math.sqrt(x * x - 1.0));
320     }
321 
322     /**
323      * @param c Complex;
324      * @return String
325      */
326     public static String printComplex(final Complex c)
327     {
328         return String.format("re=%10.6f, im=%10.6f, norm=%10.6f, phi=%10.6f(=%10.6f\u00b0)", c.re, c.im, c.norm(), c.phi(),
329                 Math.toDegrees(c.phi()));
330     }
331 
332 }