Commit c40a97d1 authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

implemented exact transformation for single velocity

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 5126
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 7f6c8a8b
......@@ -101,7 +101,15 @@ siggenx: %x: %.o
-ldatrwxx -lsffxx \
-lgsexx -ltime++ -laff -lgsl -lm -lgslcblas
lisousi deconv foutra sigval: \
lisousi: \
%: %.o
$(CXX) -o $@ $^ -I$(LOCINCLUDEDIR) -lfourierxx -lfftw3 -lm \
-lsffxx -ldatrwxx -llinearxx -lgsl -lgslcblas \
-ltsxx -ltfxx -lsffxx -lgsexx -ltime++ -laff \
-lsffxx $(FOTRANLIB) -lm\
-L$(LOCLIBDIR) $(CXXFLAGS) $(FLAGS) $(LDFLAGS)
deconv foutra sigval: \
%: %.o
$(CXX) -o $@ $^ -I$(LOCINCLUDEDIR) -lfourierxx -lfftw3 -lm \
-lsffxx -ldatrwxx -llinearxx \
......
......@@ -49,11 +49,13 @@
* transformation
* - 19/10/2012 V1.9 disabled tapslo default value
* - 20/10/2012 V1.10 properly implemented reflected wave transformation
* - 12/04/2013 V1.11 implement near-field transformation for constant
* velocity approach
*
* ============================================================================
*/
#define LISOUSI_VERSION \
"LISOUSI V1.10 line source simulation"
"LISOUSI V1.11 line source simulation"
#define LISOUSI_CVSID \
"$Id$"
......@@ -61,6 +63,7 @@
#include <fstream>
#include <string>
#include <sstream>
#include <gsl/gsl_sf_bessel.h>
#include <tfxx/commandline.h>
#include <tfxx/xcmdline.h>
#include <tfxx/error.h>
......@@ -90,7 +93,7 @@ typedef Ttimeseries::Tseries Tseries;
/*----------------------------------------------------------------------*/
struct Options {
bool verbose, debug;
bool verbose, debug, exact;
bool overwrite, taperfirst, sqrttaper, limitlength;
bool fredomain, fdfilter, tdfilter, nointeg, transition, tapsloset;
double velocity, tfac, tshift, tlim, integshift, tapdel;
......@@ -199,7 +202,7 @@ int main(int iargc, char* argv[])
"usage: lisousi [-v] [-o] [-type type] [-Type type]" "\n"
" [-fdfilter] [-tdfilter] [-fredomain]\n"
" [-taperlast] [-sqrttaper] [-limitlength]" "\n"
" [-velocity v]" "\n"
" [-velocity v] [-exact]" "\n"
" [-tapdel t] [-tapslo s]" "\n"
" [-pad n] [-tshift f] [-tfac f] [-tlim f]" "\n"
" [-nointeg] [-integshift f]\n"
......@@ -283,6 +286,8 @@ int main(int iargc, char* argv[])
"Single velocity transformation:\n"
"-fredomain apply transformation in the frequency domain for a\n"
" single direct wave velocity.\n"
"-exact apply exact transformation for single velocity approach\n"
" (default is: far-field approximation)\n"
"\n"
"-velocity v set assumed wave velocity in km/s\n"
"\n"
......@@ -478,6 +483,9 @@ int main(int iargc, char* argv[])
{"tapslo",arg_yes,"0."},
// 22: transition from fredomain to standard approach
{"transition",arg_yes,"0.,1."},
// 23: exact single velocity transformation (instead of far-field
// approximation)
{"exact",arg_no,"-"},
{NULL}
};
......@@ -558,6 +566,7 @@ int main(int iargc, char* argv[])
opt.tapslo=cmdline.double_arg(21);
opt.tapsloset=cmdline.optset(21);
opt.transition=cmdline.optset(22);
opt.exact=cmdline.optset(23);
if (opt.transition)
{
std::istringstream transopt(tfxx::string::patsubst(cmdline.string_arg(22),
......@@ -797,13 +806,38 @@ int main(int iargc, char* argv[])
<< "wave velocity: " << opt.velocity << " km/s" << endl;
}
TFourier::Tspectrum coeff=Fourier(series, dt);
coeff*=(CFTfac*sqrt(offset));
// scale with frequency independent and offset dependent factor
if (opt.exact)
{
coeff*=(-TFourier::Tcoeff(0.,1.)*offset*M_PI);
}
else
{
coeff*=(CFTfac*sqrt(offset));
}
// factor in argument of Hankel function and complex exponential
// in case of exact solution
double argfact=2.e-3*M_PI*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/T;
coeff(i)*=sqrt(1.e3*opt.velocity/f);
if (opt.exact)
{
// std::cout << "exact!" << std::endl;
double argument=f*argfact;
TFourier::Tcoeff
numerator=gsl_sf_bessel_J0(argument)-
TFourier::Tcoeff(0.,1.)*gsl_sf_bessel_Y0(argument);
TFourier::Tcoeff
denominator=exp(-TFourier::Tcoeff(0.,1.)*argument);
coeff(i) *= (numerator/denominator);
}
else
{
coeff(i)*=sqrt(1.e3*opt.velocity/f);
}
}
singlevelocityseries=Fourier(coeff, dt);
} // if (opt.fredomain)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment