singlevelocitytransformation()

TFourier::Tseries singlevelocitytransformation	(	const TFourier::Tseries &	series,
		const Parameters &	par,
		const Exco &	ec,
		const Options &	opt
	)

Definition at line 41 of file singlevelocity.cc.

References CFTfac, Options::debug, Parameters::dt, Fbessel0, Fbessel1, Fcos, Options::fdtype, Ffdexplosion, Ffdfarfield, Ffdlamb, Ffdwizforce, Ffdzforce, Fourier, Fsin, hankel(), IME, Parameters::offset, Options::radial, Parameters::T, Options::velocity, Options::verbose, Options::vpvsratio, wnintegration(), zflinefc(), and zfpointfc().

Referenced by main().

{
  TFXX_debug(opt.debug, "singlevelocitytransformation",
             "frequency domain single velocity transformation");
  if (opt.verbose)
  {
    cout << "  apply 2D/3D Greens function ratio for single "
      << "wave velocity: " << opt.velocity << " km/s" << endl;
  }

  TFXX_debug(opt.debug, "singlevelocity", "  series size: " << series.size());
  TFourier::Tspectrum coeff=Fourier(series, par.dt);
  TFXX_debug(opt.debug, "singlevelocity", "  coefficients size: "
             << coeff.size());
  // scale with frequency independent and offset dependent factor
  if (opt.fdtype==Ffdfarfield)
  {
    coeff*=(CFTfac*sqrt(par.offset));
  }

  // prepare Fourier coefficients of single force line and point source
  // Green's function in full space, if required
  TFourier::Tspectrum zfpoint, zfline;
  if (opt.fdtype==Ffdzforce)
  {
    zfpoint=zfpointfc(series.size(), par.dt, par.offset, 
                      1.e3*opt.velocity, 1.e3*opt.velocity*opt.vpvsratio,
                      opt.debug);
    zfline=zflinefc(series.size(), par.dt, par.offset, 
                    1.e3*opt.velocity, 1.e3*opt.velocity*opt.vpvsratio,
                    opt.debug);
    TFXX_assert(zfpoint.f()==coeff.f(), "shape mismatch (programming error)");
    TFXX_assert(zfline.f()==coeff.f(), "shape mismatch (programming error)");
    TFXX_assert(zfpoint.l()==coeff.l(), "shape mismatch (programming error)");
    TFXX_assert(zfline.l()==coeff.l(), "shape mismatch (programming error)");
  }

  // factor in argument of Hankel function and complex exponential
  // in case of exact solution
  double argfact=2.e-3*M_PI*par.offset/opt.velocity;

  for (unsigned int i=coeff.f(); i<=coeff.l(); ++i)
  {
    int ifre=i-coeff.f();
    if (ifre==0) { ifre = 1; }
    double f=ifre/par.T;
    if (opt.fdtype==Ffdexplosion)
    {
      // std::cout << "exact!" << std::endl;
      double argument=f*argfact;
      TFourier::Tcoeff numerator= -IME*M_PI*hankel(argument)*par.offset;
      TFourier::Tcoeff denominator=exp(-IME*argument);
      coeff(i) *= (numerator/denominator);
    }
    else if (opt.fdtype==Ffdzforce)
    {
      coeff(i) *= (zfline(i)/zfpoint(i));
    }
    else if ((opt.fdtype==Ffdwizforce) || (opt.fdtype==Ffdlamb))
    {
      /*
       * the angular frequency appears as a factor in the denominator;
       * this essentiall is an integration and produces a singularity at
       * f=0Hz; alternatively to using this factor here, we could discard it
       * and do the integration in the time domain...
       */
      TFourier::Tcoeff numerator;
      TFourier::Tcoeff denominator;
      if ((opt.fdtype==Ffdlamb) && (opt.radial))
      {
        numerator = 2.*wnintegration(ec, f, par.offset, Fsin);
        denominator = 2.*M_PI*f*wnintegration(ec, f, par.offset, Fbessel1);
      }
      else
      {
        numerator = 2.*wnintegration(ec, f, par.offset, Fcos);
        denominator = 2.*M_PI*f*wnintegration(ec, f, par.offset, Fbessel0);
      }
      coeff(i) *= (numerator/denominator);
    }
    else
    {
      coeff(i)*=sqrt(1.e3*opt.velocity/f);
    }
  }
  return(Fourier(coeff, par.dt));
}

Here is the call graph for this function:

Here is the caller graph for this function: