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.math3.analysis.function;
019
020import org.apache.commons.math3.analysis.DifferentiableUnivariateFunction;
021import org.apache.commons.math3.analysis.FunctionUtils;
022import org.apache.commons.math3.analysis.ParametricUnivariateFunction;
023import org.apache.commons.math3.analysis.UnivariateFunction;
024import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
025import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction;
026import org.apache.commons.math3.exception.DimensionMismatchException;
027import org.apache.commons.math3.exception.NullArgumentException;
028import org.apache.commons.math3.util.FastMath;
029
030/**
031 * <a href="http://en.wikipedia.org/wiki/Harmonic_oscillator">
032 *  simple harmonic oscillator</a> function.
033 *
034 * @since 3.0
035 */
036public class HarmonicOscillator implements UnivariateDifferentiableFunction, DifferentiableUnivariateFunction {
037    /** Amplitude. */
038    private final double amplitude;
039    /** Angular frequency. */
040    private final double omega;
041    /** Phase. */
042    private final double phase;
043
044    /**
045     * Harmonic oscillator function.
046     *
047     * @param amplitude Amplitude.
048     * @param omega Angular frequency.
049     * @param phase Phase.
050     */
051    public HarmonicOscillator(double amplitude,
052                              double omega,
053                              double phase) {
054        this.amplitude = amplitude;
055        this.omega = omega;
056        this.phase = phase;
057    }
058
059    /** {@inheritDoc} */
060    public double value(double x) {
061        return value(omega * x + phase, amplitude);
062    }
063
064    /** {@inheritDoc}
065     * @deprecated as of 3.1, replaced by {@link #value(DerivativeStructure)}
066     */
067    @Deprecated
068    public UnivariateFunction derivative() {
069        return FunctionUtils.toDifferentiableUnivariateFunction(this).derivative();
070    }
071
072    /**
073     * Parametric function where the input array contains the parameters of
074     * the harmonic oscillator function, ordered as follows:
075     * <ul>
076     *  <li>Amplitude</li>
077     *  <li>Angular frequency</li>
078     *  <li>Phase</li>
079     * </ul>
080     */
081    public static class Parametric implements ParametricUnivariateFunction {
082        /**
083         * Computes the value of the harmonic oscillator at {@code x}.
084         *
085         * @param x Value for which the function must be computed.
086         * @param param Values of norm, mean and standard deviation.
087         * @return the value of the function.
088         * @throws NullArgumentException if {@code param} is {@code null}.
089         * @throws DimensionMismatchException if the size of {@code param} is
090         * not 3.
091         */
092        public double value(double x, double ... param)
093            throws NullArgumentException,
094                   DimensionMismatchException {
095            validateParameters(param);
096            return HarmonicOscillator.value(x * param[1] + param[2], param[0]);
097        }
098
099        /**
100         * Computes the value of the gradient at {@code x}.
101         * The components of the gradient vector are the partial
102         * derivatives of the function with respect to each of the
103         * <em>parameters</em> (amplitude, angular frequency and phase).
104         *
105         * @param x Value at which the gradient must be computed.
106         * @param param Values of amplitude, angular frequency and phase.
107         * @return the gradient vector at {@code x}.
108         * @throws NullArgumentException if {@code param} is {@code null}.
109         * @throws DimensionMismatchException if the size of {@code param} is
110         * not 3.
111         */
112        public double[] gradient(double x, double ... param)
113            throws NullArgumentException,
114                   DimensionMismatchException {
115            validateParameters(param);
116
117            final double amplitude = param[0];
118            final double omega = param[1];
119            final double phase = param[2];
120
121            final double xTimesOmegaPlusPhase = omega * x + phase;
122            final double a = HarmonicOscillator.value(xTimesOmegaPlusPhase, 1);
123            final double p = -amplitude * FastMath.sin(xTimesOmegaPlusPhase);
124            final double w = p * x;
125
126            return new double[] { a, w, p };
127        }
128
129        /**
130         * Validates parameters to ensure they are appropriate for the evaluation of
131         * the {@link #value(double,double[])} and {@link #gradient(double,double[])}
132         * methods.
133         *
134         * @param param Values of norm, mean and standard deviation.
135         * @throws NullArgumentException if {@code param} is {@code null}.
136         * @throws DimensionMismatchException if the size of {@code param} is
137         * not 3.
138         */
139        private void validateParameters(double[] param)
140            throws NullArgumentException,
141                   DimensionMismatchException {
142            if (param == null) {
143                throw new NullArgumentException();
144            }
145            if (param.length != 3) {
146                throw new DimensionMismatchException(param.length, 3);
147            }
148        }
149    }
150
151    /**
152     * @param xTimesOmegaPlusPhase {@code omega * x + phase}.
153     * @param amplitude Amplitude.
154     * @return the value of the harmonic oscillator function at {@code x}.
155     */
156    private static double value(double xTimesOmegaPlusPhase,
157                                double amplitude) {
158        return amplitude * FastMath.cos(xTimesOmegaPlusPhase);
159    }
160
161    /** {@inheritDoc}
162     * @since 3.1
163     */
164    public DerivativeStructure value(final DerivativeStructure t)
165        throws DimensionMismatchException {
166        final double x = t.getValue();
167        double[] f = new double[t.getOrder() + 1];
168
169        final double alpha = omega * x + phase;
170        f[0] = amplitude * FastMath.cos(alpha);
171        if (f.length > 1) {
172            f[1] = -amplitude * omega * FastMath.sin(alpha);
173            final double mo2 = - omega * omega;
174            for (int i = 2; i < f.length; ++i) {
175                f[i] = mo2 * f[i - 2];
176            }
177        }
178
179        return t.compose(f);
180
181    }
182
183}