001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018package org.apache.commons.math4.analysis.function;
019
020import org.apache.commons.math4.analysis.ParametricUnivariateFunction;
021import org.apache.commons.math4.analysis.differentiation.DerivativeStructure;
022import org.apache.commons.math4.analysis.differentiation.UnivariateDifferentiableFunction;
023import org.apache.commons.math4.exception.DimensionMismatchException;
024import org.apache.commons.math4.exception.NullArgumentException;
025import org.apache.commons.math4.util.FastMath;
026
027/**
028 * <a href="http://en.wikipedia.org/wiki/Harmonic_oscillator">
029 *  simple harmonic oscillator</a> function.
030 *
031 * @since 3.0
032 */
033public class HarmonicOscillator implements UnivariateDifferentiableFunction {
034    /** Amplitude. */
035    private final double amplitude;
036    /** Angular frequency. */
037    private final double omega;
038    /** Phase. */
039    private final double phase;
040
041    /**
042     * Harmonic oscillator function.
043     *
044     * @param amplitude Amplitude.
045     * @param omega Angular frequency.
046     * @param phase Phase.
047     */
048    public HarmonicOscillator(double amplitude,
049                              double omega,
050                              double phase) {
051        this.amplitude = amplitude;
052        this.omega = omega;
053        this.phase = phase;
054    }
055
056    /** {@inheritDoc} */
057    @Override
058    public double value(double x) {
059        return value(omega * x + phase, amplitude);
060    }
061
062    /**
063     * Parametric function where the input array contains the parameters of
064     * the harmonic oscillator function, ordered as follows:
065     * <ul>
066     *  <li>Amplitude</li>
067     *  <li>Angular frequency</li>
068     *  <li>Phase</li>
069     * </ul>
070     */
071    public static class Parametric implements ParametricUnivariateFunction {
072        /**
073         * Computes the value of the harmonic oscillator at {@code x}.
074         *
075         * @param x Value for which the function must be computed.
076         * @param param Values of norm, mean and standard deviation.
077         * @return the value of the function.
078         * @throws NullArgumentException if {@code param} is {@code null}.
079         * @throws DimensionMismatchException if the size of {@code param} is
080         * not 3.
081         */
082        @Override
083        public double value(double x, double ... param)
084            throws NullArgumentException,
085                   DimensionMismatchException {
086            validateParameters(param);
087            return HarmonicOscillator.value(x * param[1] + param[2], param[0]);
088        }
089
090        /**
091         * Computes the value of the gradient at {@code x}.
092         * The components of the gradient vector are the partial
093         * derivatives of the function with respect to each of the
094         * <em>parameters</em> (amplitude, angular frequency and phase).
095         *
096         * @param x Value at which the gradient must be computed.
097         * @param param Values of amplitude, angular frequency and phase.
098         * @return the gradient vector at {@code x}.
099         * @throws NullArgumentException if {@code param} is {@code null}.
100         * @throws DimensionMismatchException if the size of {@code param} is
101         * not 3.
102         */
103        @Override
104        public double[] gradient(double x, double ... param)
105            throws NullArgumentException,
106                   DimensionMismatchException {
107            validateParameters(param);
108
109            final double amplitude = param[0];
110            final double omega = param[1];
111            final double phase = param[2];
112
113            final double xTimesOmegaPlusPhase = omega * x + phase;
114            final double a = HarmonicOscillator.value(xTimesOmegaPlusPhase, 1);
115            final double p = -amplitude * FastMath.sin(xTimesOmegaPlusPhase);
116            final double w = p * x;
117
118            return new double[] { a, w, p };
119        }
120
121        /**
122         * Validates parameters to ensure they are appropriate for the evaluation of
123         * the {@link #value(double,double[])} and {@link #gradient(double,double[])}
124         * methods.
125         *
126         * @param param Values of norm, mean and standard deviation.
127         * @throws NullArgumentException if {@code param} is {@code null}.
128         * @throws DimensionMismatchException if the size of {@code param} is
129         * not 3.
130         */
131        private void validateParameters(double[] param)
132            throws NullArgumentException,
133                   DimensionMismatchException {
134            if (param == null) {
135                throw new NullArgumentException();
136            }
137            if (param.length != 3) {
138                throw new DimensionMismatchException(param.length, 3);
139            }
140        }
141    }
142
143    /**
144     * @param xTimesOmegaPlusPhase {@code omega * x + phase}.
145     * @param amplitude Amplitude.
146     * @return the value of the harmonic oscillator function at {@code x}.
147     */
148    private static double value(double xTimesOmegaPlusPhase,
149                                double amplitude) {
150        return amplitude * FastMath.cos(xTimesOmegaPlusPhase);
151    }
152
153    /** {@inheritDoc}
154     * @since 3.1
155     */
156    @Override
157    public DerivativeStructure value(final DerivativeStructure t)
158        throws DimensionMismatchException {
159        final double x = t.getValue();
160        double[] f = new double[t.getOrder() + 1];
161
162        final double alpha = omega * x + phase;
163        f[0] = amplitude * FastMath.cos(alpha);
164        if (f.length > 1) {
165            f[1] = -amplitude * omega * FastMath.sin(alpha);
166            final double mo2 = - omega * omega;
167            for (int i = 2; i < f.length; ++i) {
168                f[i] = mo2 * f[i - 2];
169            }
170        }
171
172        return t.compose(f);
173
174    }
175
176}