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    
018    package org.apache.commons.math3.analysis.function;
019    
020    import org.apache.commons.math3.analysis.DifferentiableUnivariateFunction;
021    import org.apache.commons.math3.analysis.FunctionUtils;
022    import org.apache.commons.math3.analysis.ParametricUnivariateFunction;
023    import org.apache.commons.math3.analysis.UnivariateFunction;
024    import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
025    import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction;
026    import org.apache.commons.math3.exception.DimensionMismatchException;
027    import org.apache.commons.math3.exception.NullArgumentException;
028    import 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     */
037    public 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    }