View Javadoc
1   package org.djutils.math.functions;
2   
3   import java.util.Objects;
4   import java.util.SortedSet;
5   import java.util.TreeSet;
6   
7   import org.djutils.exceptions.Throw;
8   import org.djutils.math.AngleUtil;
9   
10  /**
11   * Sine function.
12   * <p>
13   * Copyright (c) 2023-2025 Delft University of Technology, Jaffalaan 5, 2628 BX Delft, the Netherlands. All rights reserved. See
14   * for project information <a href="https://djutils.org" target="_blank"> https://djutils.org</a>. The DJUTILS project is
15   * distributed under a three-clause BSD-style license, which can be found at
16   * <a href="https://djutils.org/docs/license.html" target="_blank"> https://djutils.org/docs/license.html</a>.
17   * </p>
18   * @author <a href="https://www.tudelft.nl/averbraeck">Alexander Verbraeck</a>
19   */
20  public class Sine implements MathFunction
21  {
22      /** Size multiplier. */
23      private final double amplitude;
24  
25      /** Time scale multiplier. */
26      private final double omega;
27  
28      /** Time scale shift. */
29      private final double shift;
30  
31      /** The function that yields x (may be null). */
32      private final MathFunction chain;
33  
34      /**
35       * Construct a new Sine function <code>factor * sin(omega * x + shift)</code>.
36       * @param chain the MathFunction that yields the <code>x</code> for this power function
37       * @param amplitude multiplication factor for the output
38       * @param omega radial frequency; multiplication factor for the input
39       * @param shift time shift for the input (applied <b>after</b> the <code>omega</code> factor)
40       */
41      public Sine(final MathFunction chain, final double amplitude, final double omega, final double shift)
42      {
43          this.amplitude = amplitude;
44          this.omega = omega;
45          this.shift = shift;
46          this.chain = chain;
47      }
48  
49      /**
50       * Construct a new Sine function <code>factor * sin(omega * x + shift)</code>.
51       * @param amplitude multiplication factor for the output
52       * @param omega radial frequency; multiplication factor for the input
53       * @param shift time shift for the input (applied <b>after</b> the <code>omega</code> factor)
54       */
55      public Sine(final double amplitude, final double omega, final double shift)
56      {
57          this.amplitude = amplitude;
58          this.omega = omega;
59          this.shift = shift;
60          this.chain = null;
61      }
62  
63      /**
64       * Construct a cosine function <code>factor * cos(omega * x + shift)</code>. The result is actually a Sine with the
65       * correctly adjusted <code>shift</code>.
66       * @param amplitude multiplication factor for the output
67       * @param omega radial frequency; multiplication factor for the input
68       * @param shift time shift for the input (applied <b>after</b> the <code>omega</code> factor)
69       * @return Sine with the requested <code>amplitude</code>, <code>omega</code> and adjusted <code>shift</code>
70       */
71      public static Sine cosine(final double amplitude, final double omega, final double shift)
72      {
73          return new Sine(amplitude, omega, AngleUtil.normalizeAroundZero(shift + Math.PI / 2));
74      }
75  
76      @Override
77      public double get(final double x)
78      {
79          if (this.amplitude == 0.0)
80          {
81              return 0.0;
82          }
83          double xValue = this.chain == null ? x : this.chain.get(x);
84          return this.amplitude * Math.sin(this.omega * xValue + this.shift);
85      }
86  
87      @Override
88      public MathFunction getDerivative()
89      {
90          Sine myDerivative = new Sine(this.chain, this.amplitude * this.omega, this.omega,
91                  AngleUtil.normalizeAroundZero(this.shift + Math.PI / 2));
92          if (this.chain == null)
93          {
94              return myDerivative.simplify();
95          }
96          MathFunction myChainDerivative = new Sine(this.chain, myDerivative.amplitude, myDerivative.omega, myDerivative.shift);
97          return new Product(myChainDerivative, this.chain.getDerivative()).simplify();
98      }
99  
100     @Override
101     public MathFunction simplify()
102     {
103         if (this.amplitude == 0.0)
104         {
105             return Constant.ZERO;
106         }
107         if (this.chain != null && this.chain instanceof Constant)
108         {
109             return new Constant(get(0)).simplify();
110         }
111         return this;
112     }
113 
114     @Override
115     public double getScale()
116     {
117         return this.amplitude;
118     }
119 
120     @Override
121     public MathFunction scaleBy(final double scaleFactor)
122     {
123         if (scaleFactor == 0.0)
124         {
125             return Constant.ZERO;
126         }
127         if (scaleFactor == 1.0)
128         {
129             return this;
130         }
131         return new Sine(this.chain, scaleFactor * this.amplitude, this.omega, this.shift);
132     }
133 
134     @Override
135     public int sortPriority()
136     {
137         return 4;
138     }
139 
140     @Override
141     public int compareWithinSubType(final MathFunction other)
142     {
143         Throw.when(!(other instanceof Sine), IllegalArgumentException.class, "other is of wrong type");
144         Sine otherSine = (Sine) other;
145         if (this.omega < otherSine.omega)
146         {
147             return -1;
148         }
149         if (this.omega > otherSine.omega)
150         {
151             return 1;
152         }
153         if (this.amplitude < otherSine.amplitude)
154         {
155             return -1;
156         }
157         if (this.amplitude > otherSine.amplitude)
158         {
159             return 1;
160         }
161         if (this.shift < otherSine.shift)
162         {
163             return -1;
164         }
165         if (this.shift > otherSine.shift)
166         {
167             return 1;
168         }
169         return compareChains(this.chain, otherSine.chain);
170     }
171 
172     /**
173      * Compute the sum of two sines. The sines must have the same <code>omega</code> but have their own amplitude and shift.
174      * @param chain the chained MathFunction
175      * @param amplitude1 amplitude of the first sine
176      * @param amplitude2 amplitude of the second sine
177      * @param omega angular frequency of both sines
178      * @param shift1 phase shift of first sine
179      * @param shift2 phase shift of second sine
180      * @return a new Sine that represents the sum of the supplied sines
181      */
182     private static Sine sumSines(final MathFunction chain, final double amplitude1, final double amplitude2, final double omega,
183             final double shift1, final double shift2)
184     {
185         // There is probably a way to calculate the result that takes fewer CPU cycles.
186         double re = amplitude1 * Math.cos(shift1) + amplitude2 * Math.cos(shift2);
187         double im = amplitude1 * Math.sin(shift1) + amplitude2 * Math.sin(shift2);
188         double resultAmplitude = Math.hypot(re, im);
189         double resultShift = Math.atan2(im, re);
190         return new Sine(chain, resultAmplitude, omega, resultShift);
191     }
192 
193     @Override
194     public MathFunction mergeAdd(final MathFunction other)
195     {
196         if (other instanceof Sine)
197         {
198             Sine otherSine = (Sine) other;
199             if (this.omega == otherSine.omega && (this.chain == null && otherSine.chain == null
200                     || (this.chain != null && this.chain.equals(otherSine.chain))))
201             {
202                 return sumSines(this.chain, this.amplitude, otherSine.amplitude, this.omega, this.shift, otherSine.shift);
203             }
204         }
205         return null;
206     }
207 
208     @Override
209     public MathFunction mergeMultiply(final MathFunction other)
210     {
211         if (other instanceof Sine)
212         {
213             Sine otherSine = (Sine) other;
214             if (this.omega == otherSine.omega && (this.chain == null && otherSine.chain == null
215                     || (this.chain != null && this.chain.equals(otherSine.chain))))
216             {
217                 /*-
218                  * a * sin(x + theta) * b * sin(x + phi)
219                  * == (cos(x + theta - x - phi) - cos(x + theta + x + phi))               * a * b / 2
220                  * == (cos(theta - phi)         - cos(2 * x + theta + phi))               * a * b / 2
221                  * == (cos(theta - phi)         - sin(pi/2 + 2 * x + theta + phi))        * a * b / 2
222                  * == (cos(theta - phi)         + sin(pi + pi / 2 + 2 * x + theta + phi)) * a * b / 2
223                  * == (cos(theta - phi)         + sin(2 * x + theta + phi + pi + pi / 2)) * a * b / 2
224                  */
225                 double scale = this.amplitude * otherSine.amplitude / 2;
226                 double phaseShift = AngleUtil.normalizeAroundZero(this.shift + otherSine.shift + Math.PI + Math.PI / 2);
227                 return new Sum(new Constant(Math.cos(this.shift - otherSine.shift) * scale),
228                         new Sine(this.chain, scale, this.omega * 2, phaseShift));
229             }
230         }
231         return null;
232     }
233 
234     @Override
235     public KnotReport getKnotReport(final Interval<?> interval)
236     {
237         return this.chain == null ? KnotReport.NONE : this.chain.getKnotReport(interval);
238     }
239 
240     @Override
241     public SortedSet<Double> getKnots(final Interval<?> interval)
242     {
243         return this.chain == null ? new TreeSet<Double>() : this.chain.getKnots(interval);
244     }
245 
246     @Override
247     public String toString()
248     {
249         double quadrant = AngleUtil.normalizeAroundPi(this.shift) / (Math.PI / 2);
250         int roundedQuadrant = (int) Math.round(quadrant);
251         double deviation = Math.abs(roundedQuadrant - quadrant);
252         boolean closeTo90 = deviation < 10 * Math.ulp(Math.PI);
253         boolean useCosine = closeTo90 && roundedQuadrant % 2 == 1;
254         boolean useSine = closeTo90 && roundedQuadrant % 2 == 0;
255         // System.out.println("roundedQuadrant=" + roundedQuadrant);
256         boolean negativeSign = (useSine || useCosine) && ((roundedQuadrant >= 2) != (this.amplitude < 0));
257 
258         StringBuilder result = new StringBuilder();
259         if (negativeSign)
260         {
261             result.append("-");
262         }
263         if (this.amplitude != 1.0 && this.amplitude != -1.0)
264         {
265             result.append(printValue(Math.abs(this.amplitude)));
266         }
267         result.append(useCosine ? "cos(" : "sin(");
268         if (this.omega != 1.0)
269         {
270             result.append(printValue(this.omega));
271         }
272         result.append(this.chain == null ? "x" : ("(" + this.chain.toString() + ")"));
273         if (this.shift != 0.0)
274         {
275             if (!closeTo90)
276             {
277                 if (this.shift >= 0.0)
278                 {
279                     result.append("+");
280                 }
281                 result.append(printValue(this.shift));
282             }
283         }
284         result.append(")");
285         return result.toString();
286     }
287 
288     @Override
289     public int hashCode()
290     {
291         return Objects.hash(this.amplitude, this.chain, this.omega, this.shift);
292     }
293 
294     @SuppressWarnings("checkstyle:needbraces")
295     @Override
296     public boolean equals(final Object obj)
297     {
298         if (this == obj)
299             return true;
300         if (obj == null)
301             return false;
302         if (getClass() != obj.getClass())
303             return false;
304         Sine other = (Sine) obj;
305         return Double.doubleToLongBits(this.amplitude) == Double.doubleToLongBits(other.amplitude)
306                 && Objects.equals(this.chain, other.chain)
307                 && Double.doubleToLongBits(this.omega) == Double.doubleToLongBits(other.omega)
308                 && Double.doubleToLongBits(this.shift) == Double.doubleToLongBits(other.shift);
309     }
310 
311 }