Line data    Source code

       1             : /*************************************************/
       2             : /*           DO NOT MODIFY THIS HEADER           */
       3             : /*                                               */
       4             : /*                     MASTODON                  */
       5             : /*                                               */
       6             : /*    (c) 2015 Battelle Energy Alliance, LLC     */
       7             : /*            ALL RIGHTS RESERVED                */
       8             : /*                                               */
       9             : /*   Prepared by Battelle Energy Alliance, LLC   */
      10             : /*     With the U. S. Department of Energy       */
      11             : /*                                               */
      12             : /*     See COPYRIGHT for full restrictions       */
      13             : /*************************************************/
      14             : 
      15             : #include "OrmsbyWavelet.h"
      16             : 
      17             : registerMooseObject("MastodonApp", OrmsbyWavelet);
      18             : 
      19             : InputParameters
      20           6 : OrmsbyWavelet::validParams()
      21             : {
      22           6 :   InputParameters params = Function::validParams();
      23          12 :   params.addRequiredParam<Real>("f1", "First frequency for defining the Ormsby wavelet.");
      24          12 :   params.addRequiredParam<Real>("f2", "Second frequency for defining the Ormsby wavelet.");
      25          12 :   params.addRequiredParam<Real>("f3", "Third frequency for defining the Ormsby wavelet.");
      26          12 :   params.addRequiredParam<Real>("f4", "Fourth frequency for defining the Ormsby wavelet.");
      27          12 :   params.addRequiredParam<Real>("ts", "Time of the peak of the Ormsby wavelet.");
      28          12 :   params.addParam<Real>("scale_factor", 1.0, "Amplitude scale factor to be applied to wavelet.");
      29           6 :   params.addClassDescription(
      30             :       "Calculates an amplitude normalized Ormsby wavelet with the given input parameters.");
      31           6 :   return params;
      32           0 : }
      33             : 
      34           3 : OrmsbyWavelet::OrmsbyWavelet(const InputParameters & parameters)
      35           6 :   : Function(parameters), _scale_factor(getParam<Real>("scale_factor"))
      36             : {
      37           3 : }
      38             : 
      39             : Real
      40        1203 : OrmsbyWavelet::value(Real t, const Point &) const
      41             : {
      42        2406 :   Real f1 = getParam<Real>("f1");
      43        2406 :   Real f2 = getParam<Real>("f2");
      44        2406 :   Real f3 = getParam<Real>("f3");
      45        2406 :   Real f4 = getParam<Real>("f4");
      46        2406 :   Real ts = getParam<Real>("ts");
      47             : 
      48             :   Real c1, c2, c3, c4, c;
      49             : 
      50        2406 :   c1 = libMesh::pi * f1 * f1 / (f2 - f1) * sinc(libMesh::pi * f1 * (t - ts)) *
      51             :        sinc(libMesh::pi * f1 * (t - ts));
      52        2406 :   c2 = libMesh::pi * f2 * f2 / (f2 - f1) * sinc(libMesh::pi * f2 * (t - ts)) *
      53             :        sinc(libMesh::pi * f2 * (t - ts));
      54        2406 :   c3 = libMesh::pi * f3 * f3 / (f3 - f4) * sinc(libMesh::pi * f3 * (t - ts)) *
      55             :        sinc(libMesh::pi * f3 * (t - ts));
      56        2406 :   c4 = libMesh::pi * f4 * f4 / (f3 - f4) * sinc(libMesh::pi * f4 * (t - ts)) *
      57             :        sinc(libMesh::pi * f4 * (t - ts));
      58             : 
      59             :   // c is the max value
      60        1203 :   c = (libMesh::pi * f4 * f4 / (f3 - f4) - libMesh::pi * f3 * f3 / (f3 - f4)) -
      61        1203 :       (libMesh::pi * f2 * f2 / (f2 - f1) - libMesh::pi * f1 * f1 / (f2 - f1));
      62             : 
      63        1203 :   return _scale_factor / c * ((c4 - c3) - (c2 - c1));
      64             : }
      65             : 
      66             : // sinc function (tends to 1 as x -> 0)
      67             : inline Real
      68             : OrmsbyWavelet::sinc(Real x) const
      69             : {
      70        1203 :   return (x == 0) ? 1.0 : sin(x) / x;
      71             : }