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