1 package org.djutils.complex;
2
3 import org.djutils.base.AngleUtil;
4
5 /**
6 * ComplexMath.java. Math with complex operands and results. <br>
7 * Copyright (c) 2021-2023 Delft University of Technology, Jaffalaan 5, 2628 BX Delft, the Netherlands. All rights reserved. See
8 * for project information <a href="https://djutils.org" target="_blank"> https://djutils.org</a>. The DJUTILS project is
9 * distributed under a three-clause BSD-style license, which can be found at
10 * <a href="https://djutils.org/docs/license.html" target="_blank"> https://djutils.org/docs/license.html</a>. <br>
11 * @author <a href="https://www.tudelft.nl/averbraeck">Alexander Verbraeck</a>
12 * @author <a href="https://www.tudelft.nl/pknoppers">Peter Knoppers</a>
13 */
14 public final class ComplexMath
15 {
16
17 /**
18 * Do not instantiate.
19 */
20 private ComplexMath()
21 {
22 // Do not instantiate.
23 }
24
25 /**
26 * Principal square root of a Complex operand. The principal square root of a complex number has a non-negative real
27 * component.
28 * @param z Complex; the operand
29 * @return Complex; the principal square root of the operand
30 */
31 public static Complex sqrt(final Complex z)
32 {
33 double norm = z.norm();
34 return new Complex(Math.sqrt((z.re + norm) / 2), (z.im >= 0 ? 1 : -1) * Math.sqrt((-z.re + norm) / 2));
35 }
36
37 /**
38 * Principal cube root of a Complex operand. The principal cube root of a complex number has <i>the most positive</i>. real
39 * component.
40 * @param z Complex; the operand
41 * @return Complex; the principal cube root of the operand
42 */
43 public static Complex cbrt(final Complex z)
44 {
45 double cbrtOfNorm = Math.cbrt(z.norm());
46 double phi = z.phi() / 3;
47 return new Complex(cbrtOfNorm * Math.cos(phi), cbrtOfNorm * Math.sin(phi));
48 }
49
50 /**
51 * Exponential function of a Complex operand.
52 * @param z Complex; the operand
53 * @return Complex; the result of the exponential function applied to the operand
54 */
55 public static Complex exp(final Complex z)
56 {
57 double factor = Math.exp(z.re);
58 return new Complex(factor * Math.cos(z.im), factor * Math.sin(z.im));
59 }
60
61 /**
62 * Principal value of the natural logarithm of a Complex operand. See
63 * <a href="https://en.wikipedia.org/wiki/Complex_logarithm">Wikipedia Complex logarithm</a>.
64 * @param z Complex; the operand
65 * @return Complex; the principal value of the natural logarithm of the Complex operand
66 */
67 public static Complex ln(final Complex z)
68 {
69 return new Complex(Math.log(z.norm()), z.phi());
70 }
71
72 /**
73 * Sine function of a Complex operand. See <a href="https://proofwiki.org/wiki/Sine_of_Complex_Number">ProofWiki Sine of
74 * Complex Number</a>.
75 * @param z Complex; the operand
76 * @return Complex; the result of the sine function applied to the operand
77 */
78 public static Complex sin(final Complex z)
79 {
80 return new Complex(Math.sin(z.re) * Math.cosh(z.im), Math.cos(z.re) * Math.sinh(z.im));
81 }
82
83 /**
84 * Cosine function of Complex operand. See <a href="https://proofwiki.org/wiki/Cosine_of_Complex_Number">ProofWiki Cosine of
85 * Complex Number</a>.
86 * @param z Complex; the operand
87 * @return Complex; the result of the cosine function applied to the operand
88 */
89 public static Complex cos(final Complex z)
90 {
91 return new Complex(Math.cos(z.re) * Math.cosh(z.im), -Math.sin(z.re) * Math.sinh(z.im));
92 }
93
94 /**
95 * Tangent function of a Complex operand. See <a href="https://proofwiki.org/wiki/Tangent_of_Complex_Number">ProofWiki
96 * Tangent of Complex Number</a>.
97 * @param z Complex; the operand
98 * @return Complex; the result of the tangent function applied to the operand
99 */
100 public static Complex tan(final Complex z)
101 {
102 // Using Formulation 4 of the reference as it appears to need the fewest trigonometric operations
103 double divisor = Math.cos(2 * z.re) + Math.cosh(2 * z.im);
104 return new Complex(Math.sin(2 * z.re) / divisor, Math.sinh(2 * z.im) / divisor);
105 }
106
107 /**
108 * Hyperbolic sine function of a Complex operand.
109 * @param z Complex; the operand
110 * @return Complex; the result of the sinh function applied to the operand
111 */
112 public static Complex sinh(final Complex z)
113 {
114 return new Complex(Math.sinh(z.re) * Math.cos(z.im), Math.cosh(z.re) * Math.sin(z.im));
115 }
116
117 /**
118 * Hyperbolic cosine function of Complex operand.
119 * @param z Complex; the operand
120 * @return Complex; the result of the cosh function applied to the operand
121 */
122 public static Complex cosh(final Complex z)
123 {
124 return new Complex(Math.cosh(z.re) * Math.cos(z.im), Math.sinh(z.re) * Math.sin(z.im));
125 }
126
127 /**
128 * Hyperbolic tangent function of a Complex operand. Based on
129 * <a href="https://proofwiki.org/wiki/Hyperbolic_Tangent_of_Complex_Number">ProofWiki: Hyperbolic Tangent of Complex
130 * Number, Formulation 4</a>.
131 * @param z Complex; the operand
132 * @return Complex; the result of the tanh function applied to the operand
133 */
134 public static Complex tanh(final Complex z)
135 {
136 double re2 = z.re * 2;
137 double im2 = z.im * 2;
138 double divisor = Math.cosh(re2) + Math.cos(im2);
139 return new Complex(Math.sinh(re2) / divisor, Math.sin(im2) / divisor);
140 }
141
142 /**
143 * Inverse sine function of a Complex operand. Derived from <a href=
144 * "http://cvsweb.netbsd.org/bsdweb.cgi/~checkout~/src/lib/libm/complex/casin.c?rev=1.1&content-type=text/plain">NetBSD
145 * Complex casin.c</a>.
146 * @param z Complex; the operand
147 * @return Complex; the result of the asin function applied to the operand
148 */
149 public static Complex asin(final Complex z)
150 {
151 Complex ct = z.times(Complex.I);
152 Complex zz = new Complex((z.re - z.im) * (z.re + z.im), (2.0 * z.re * z.im));
153
154 zz = new Complex(1.0 - zz.re, -zz.im);
155 zz = ct.plus(sqrt(zz));
156 zz = ln(zz);
157 /* multiply by 1/i = -i */
158 return zz.times(Complex.MINUS_I);
159 }
160
161 /**
162 * Inverse cosine function of a Complex operand. Derived from <a href=
163 * "http://cvsweb.netbsd.org/bsdweb.cgi/~checkout~/src/lib/libm/complex/cacos.c?rev=1.1&content-type=text/plain">NetBSD
164 * Complex cacos.c</a>.
165 * @param z Complex; the operand
166 * @return Complex; the result of the acos function applied to the operand
167 */
168 public static Complex acos(final Complex z)
169 {
170 Complex asin = asin(z);
171 return new Complex(Math.PI / 2 - asin.re, -asin.im);
172 }
173
174 /** Maximum value of a Complex. May be returned by the atan function. */
175 private static final Complex MAXIMUM = new Complex(Double.MAX_VALUE, Double.MAX_VALUE);
176
177 /**
178 * Inverse tangent function of a Complex operand. Derived from <a href=
179 * "http://cvsweb.netbsd.org/bsdweb.cgi/~checkout~/src/lib/libm/complex/catan.c?rev=1.2&content-type=text/plain">NetBSD
180 * Complex catan.c</a>.
181 * @param z Complex; the operand
182 * @return Complex; the result of the atan function applied to the operand
183 */
184 public static Complex atan(final Complex z)
185 {
186 if ((z.re == 0.0) && (z.im > 1.0))
187 {
188 return MAXIMUM;
189 }
190
191 double x2 = z.re * z.re;
192 double a = 1.0 - x2 - (z.im * z.im);
193 if (a == 0.0)
194 {
195 return MAXIMUM;
196 }
197
198 double re = AngleUtil.normalizeAroundZero(Math.atan2(2.0 * z.re, a)) / 2;
199 double t = z.im - 1.0;
200 a = x2 + (t * t);
201 if (a == 0.0)
202 {
203 return MAXIMUM;
204 }
205
206 t = z.im + 1.0;
207 a = (x2 + (t * t)) / a;
208 return new Complex(re, 0.25 * Math.log(a));
209 }
210
211 /**
212 * Inverse hyperbolic sine of a Complex operand.
213 * @param z Complex; the operand
214 * @return Complex; the result of the asinh function applied to the operand
215 */
216 public static Complex asinh(final Complex z)
217 {
218 return Complex.MINUS_I.times(asin(z.times(Complex.I)));
219 }
220
221 /**
222 * Inverse hyperbolic cosine of a Complex operand.
223 * @param z Complex; the operand
224 * @return Complex; the result of the acosh function applied to the operand
225 */
226 public static Complex acosh(final Complex z)
227 {
228 return ln(z.plus(sqrt(z.plus(Complex.ONE)).times(sqrt(z.minus(Complex.ONE)))));
229 }
230
231 /**
232 * Inverse hyperbolic tangent of a Complex operand. Derived from
233 * <a href="http://cvsweb.netbsd.org/bsdweb.cgi/~checkout~/src/lib/libm/complex/catanh.c?rev=1.1&content-type=text/plain">
234 * NetBSD Complex catanh.c</a>.
235 * @param z Complex; the operand
236 * @return Complex; the result of the atanh function applied to the operand
237 */
238 public static Complex atanh(final Complex z)
239 {
240 return Complex.MINUS_I.times(atan(z.times(Complex.I)));
241 }
242
243 }