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

soutifu compiles but is not yet tested

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: 3996
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 431a15af
......@@ -25,7 +25,7 @@
PROGRAMS=tsfilt stufi rotate coro xyz2uvw susei evelo tesiff teswf \
phasedsignals hamres siggen smoos dise foutra teseco resaseda gatherdiff \
autocorr cross tidofi fredofi sigfit noisymize sigval fidasexx
autocorr cross tidofi fredofi sigfit noisymize sigval fidasexx soutifu
.PHONY: all
all: install doc
......
......@@ -27,6 +27,7 @@
#include <tfxx/commandline.h>
#include <tsxx/sffheaders.h>
#include <datrwxx/readany.h>
#include <datrwxx/writeany.h>
#include <stfinv/stfinvany.h>
/*----------------------------------------------------------------------*/
......@@ -35,6 +36,15 @@ using std::cout;
using std::cerr;
using std::endl;
/*----------------------------------------------------------------------*/
// function to compare double
// is true if relative residual between a and b is smaller than eps
bool sameineps(const double &a, const double& b, const double& eps=1.e-8)
{
double reldif=std::abs(a-b);
return(reldif<(std::abs(b*eps)));
}
/*----------------------------------------------------------------------*/
struct Options {
......@@ -79,6 +89,10 @@ int main(int iargc, char* argv[])
"--type f use file format type \"f\" for file I/O\n"
"-wc f write convolved synthetics to file \"f\"\n"
"-ws f write source time function to file \"f\"\n"
"\n"
"Upon input the consistency of recorded and synthetic data are\n"
"checked. For coordinates only receiver positions relative to\n"
"source coordinates are checked, not absolute locations.\n"
};
// define commandline options
......@@ -150,7 +164,7 @@ int main(int iargc, char* argv[])
// read input data
// ---------------
typedef aff::Series<double> Tseries;
typedef aff::Series<float> Tseries;
typedef ts::sff::File<Tseries> Tfile;
Tfile recordeddata;
Tfile syntheticdata;
......@@ -170,6 +184,185 @@ int main(int iargc, char* argv[])
datrw::ianystream is(ifs, opt.datafileformat);
syntheticdata.read(is.idatstream(), opt.verbose);
}
}
// check input data consistency and prepare output files
// -----------------------------------------------------
if (opt.verbose)
{
cout << "check input data consistency,\n"
<< "prepare output data containers, and\n"
<< "prepare waveform data workspace" << endl;
}
// number of traces per file
TFXX_assert(recordeddata.size() == syntheticdata.size(),
"ERROR: inconsitent number of traces");
const unsigned int ntraces=recordeddata.size();
// output data
Tfile convolvedsynthetics;
Tfile stffile;
// waveform data workspace
stfinv::Tvectoroftriples vectoroftriples;
stfinv::Waveform stfwaveform;
convolvedsynthetics.fileheader=syntheticdata.fileheader;
// cycle through all traces
for (unsigned int itrace=0; itrace<ntraces; ++itrace)
{
if (opt.verbose)
{
cout << " check trace #" << itrace+1 << endl;
}
// create a reference to the time series with trace header for thsi
// specific trace as read from file
const Tfile::Ttracevector::Ttimeseries& rdtseries=recordeddata[itrace];
const Tfile::Ttracevector::Ttimeseries& sdtseries=syntheticdata[itrace];
// create a time series to be used in the output of convolved synthetics
// based on the input synthetics; copy trace header
Tfile::Ttracevector::Ttimeseries cstseries(rdtseries.shape());
cstseries.Mtraceindex=sdtseries.Mtraceindex;
cstseries.header=sdtseries.header;
// place reference in a waveform triple
stfinv::WaveformTriple tracetriple;
tracetriple.data=rdtseries;
tracetriple.synthetics=sdtseries;
tracetriple.convolvedsynthetics=cstseries;
tracetriple.header.sampling.dt=rdtseries.header.wid2().dt;
tracetriple.header.sampling.n=rdtseries.header.wid2().nsamples;
TFXX_assert((rdtseries.header.wid2().nsamples
== sdtseries.header.wid2().nsamples),
"ERROR: inconsistent trace size");
TFXX_assert(sameineps(rdtseries.header.wid2().dt,
sdtseries.header.wid2().dt),
"ERROR: inconsistent sampling interval");
// read out source coordinates
tracetriple.header.sx=recordeddata.fileheader.srce().cx;
tracetriple.header.sy=recordeddata.fileheader.srce().cy;
tracetriple.header.sz=recordeddata.fileheader.srce().cz;
// read out receiver coordinates
tracetriple.header.rx=rdtseries.header.info().cx;
tracetriple.header.ry=rdtseries.header.info().cy;
tracetriple.header.rz=rdtseries.header.info().cz;
// check coordinate consistency
double ddx=rdtseries.header.info().cx-recordeddata.fileheader.srce().cx;
double ddy=rdtseries.header.info().cy-recordeddata.fileheader.srce().cy;
double ddz=rdtseries.header.info().cz-recordeddata.fileheader.srce().cz;
double sdx=sdtseries.header.info().cx-syntheticdata.fileheader.srce().cx;
double sdy=sdtseries.header.info().cy-syntheticdata.fileheader.srce().cy;
double sdz=sdtseries.header.info().cz-syntheticdata.fileheader.srce().cz;
TFXX_assert((sameineps(ddx,sdx) && sameineps(ddy,sdy) &&
sameineps(ddz,sdz)),
"ERROR: inconsistent receiver positions");
// check time consistency
TFXX_assert((rdtseries.header.wid2().date
== recordeddata.fileheader.srce().date),
"ERROR: trigger delay not supported");
TFXX_assert((sdtseries.header.wid2().date
== syntheticdata.fileheader.srce().date),
"ERROR: trigger delay not supported");
// add this triple to collection
vectoroftriples.push_back(tracetriple);
// add this trace to convolved synthetics
convolvedsynthetics.push_back(cstseries);
// catch values for stf waveform
if (itrace==0)
{
Tfile::Ttracevector::Ttimeseries stftseries(tracetriple.data.shape());
stfwaveform.sampling.dt=tracetriple.header.sampling.dt;
stfwaveform.sampling.n=tracetriple.header.sampling.n;
stfwaveform.series=stftseries;
stftseries.header.wid2().nsamples=stfwaveform.sampling.n;
stftseries.header.wid2().dt=stfwaveform.sampling.dt;
stftseries.header.wid2().date=libtime::now();
stffile.fileheader.srce().date=stftseries.header.wid2().date;
stffile.push_back(stftseries);
} // if (itrace==0)
} // for (unsigned int itrace=0; itrace<ntraces; ++itrace)
// create STF engine
// -----------------
if (opt.verbose)
{
parameters += ":verbose";
cout << "create STF engine" << endl;
}
stfinv::STFEngine engine(vectoroftriples, stfwaveform, parameters);
// run STF engine
// --------------
if (opt.verbose) { cout << "run engine" << endl; }
engine.run();
// write results
// -------------
if (opt.writestf)
{
if (opt.verbose)
{
cout << "write STF to file " << opt.stffilename << endl;
}
if (!opt.overwrite)
{
std::ifstream file(opt.stffilename.c_str(),std::ios_base::in);
TFXX_assert((!file.good()),"ERROR: output file exists!");
}
std::ios_base::openmode oopenmode
=datrw::oanystream::openmode(opt.datafileformat);
std::ofstream ofs(opt.stffilename.c_str(), oopenmode);
datrw::oanystream os(ofs, opt.datafileformat, opt.debug);
os << stffile.fileheader.srce();
for (unsigned int itrace=0; itrace<stffile.size(); ++itrace)
{
os << stffile[itrace].header.wid2();
os << stffile[itrace].header.info();
os << stffile[itrace].series();
}
} // if (opt.writestf)
if (opt.writeconvolved)
{
if (opt.verbose)
{
cout << "write convolved synthetics to file "
<< opt.convolvedfilename << endl;
}
if (!opt.overwrite)
{
std::ifstream file(opt.convolvedfilename.c_str(),std::ios_base::in);
TFXX_assert((!file.good()),"ERROR: output file exists!");
}
std::ios_base::openmode oopenmode
=datrw::oanystream::openmode(opt.datafileformat);
std::ofstream ofs(opt.convolvedfilename.c_str(), oopenmode);
datrw::oanystream os(ofs, opt.datafileformat, opt.debug);
os << convolvedsynthetics.fileheader.srce();
for (unsigned int itrace=0; itrace<convolvedsynthetics.size(); ++itrace)
{
os << convolvedsynthetics[itrace].header.wid2();
os << convolvedsynthetics[itrace].header.info();
os << convolvedsynthetics[itrace].series();
}
} // if (opt.writeconvolved)
} // main
/* ----- END OF soutifu.cc ----- */
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