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

FourPoint, OffsetVariableTaper, Picks and Pick are implemented and 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: 4489
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 99df12ec
......@@ -53,6 +53,7 @@
- anyfilter.h
- convolve.h
- correlate.h
- debug.h
- dropcontainer.h
- filterbase.h
- filter.h
......@@ -60,6 +61,7 @@
- innerproduct.h
- ipo.h
- ipolin.h
- ovtaper.h
- random.h
- seifeclass.h
- seifexx.h
......
......@@ -53,6 +53,7 @@ namespace ts {
//! read from file in refract taper file format
void Picks::read(std::istream& is)
{
Mpicks.clear();
std::string line;
do {
TSXX_assert(getline(is, line),
......@@ -82,16 +83,47 @@ namespace ts {
}
Mpicks.sort();
this->Mp1=this->Mpicks.begin();
} // void Picks::read(std::istream& is)
/*----------------------------------------------------------------------*/
Pick Picks::pick(const double& offset) const
{
Pick p(offset);
if (!Mpicks.empty())
{
Tlistofpick::const_iterator P=Mpicks.begin();
if ((*P) > p)
{
p.t=P->t;
}
else
{
Tlistofpick::const_iterator Pp=P;
while ((P!=Mpicks.end()) && (*P<p)) { Pp=P; P++; }
if (P==Pp)
{
p.t=P->t;
}
else if (P==Mpicks.end())
{
p.t=Pp->t;
}
else
{
p.t=Pp->t+(P->t-Pp->t)*(p.x-Pp->x)/(P->x-Pp->x);
}
}
}
return p;
} // Pick Picks::pick(const double& offset) const
/*----------------------------------------------------------------------*/
double Picks::time(const double& offset) const
{
double retval=0.;
TSXX_abort("not yet implemented");
return retval;
Pick p=this->pick(offset);
return(p.t);
} // double Picks::time(const double& offset) const
} // namespace ovtaper
......@@ -106,9 +138,9 @@ namespace ts {
"OffsetVariableTaper::taper: "
"taper is undefined");
double t1=(Mt1.time(offset)-T0)/T;
double t2=(Mt1.time(offset)-T0)/T;
double t3=(Mt1.time(offset)-T0)/T;
double t4=(Mt1.time(offset)-T0)/T;
double t2=(Mt2.time(offset)-T0)/T;
double t3=(Mt3.time(offset)-T0)/T;
double t4=(Mt4.time(offset)-T0)/T;
ts::tapers::FourPoint fpt(t1, t2, t3, t4);
return(fpt);
} // ts::tapers::FourPoint OffsetVariableTaper::taper(...)
......
......@@ -56,13 +56,18 @@ namespace ts {
//! a single pick
struct Pick {
//! time and offset
double t,x;
Pick(): t(0), x(0) { }
Pick(const double& offset): t(0), x(offset) { }
double t; //!< time
double x; //!< offset
}; // struct Pick
/*----------------------------------------------------------------------*/
inline
bool operator<(const Pick& p1, const Pick& p2) { return (p1.x<p2.x); }
inline
bool operator>(const Pick& p1, const Pick& p2) { return (p1.x>p2.x); }
/*----------------------------------------------------------------------*/
......@@ -73,15 +78,15 @@ namespace ts {
Picks(const bool& debug): Mdebug(debug) { }
//! read from file in refract taper file format
void read(std::istream& is);
//! return time for given offset (with interpolation)
//! return interpolated pick for given offset
Pick pick(const double& offset) const;
//! return time for interpolated pick at given offset
double time(const double& offset) const;
private:
//! produce debug output if true
bool Mdebug;
//! picks
Tlistofpick Mpicks;
//! iterators used while seeking entries
mutable Tlistofpick::const_iterator Mp1, Mp2;
}; // class Picks
} // namespace ovtaper
......@@ -124,6 +129,11 @@ namespace ts {
/*! read taper definition from file
*/
void read(const std::string& filename);
ovtaper::Picks t1() const { return Mt1; }
ovtaper::Picks t2() const { return Mt2; }
ovtaper::Picks t3() const { return Mt3; }
ovtaper::Picks t4() const { return Mt4; }
private:
//! produce debug output if true
bool Mdebug;
......
......@@ -40,10 +40,15 @@
#define TF_TAPERS_CC_CVSID \
"$Id$"
//#include <iostream>
#include <tsxx/tapers.h>
#include <tsxx/tsxx.h>
//#include <tsxx/debug.h>
#include <cmath>
//using std::cout;
//using std::endl;
namespace ts {
//! \brief Provides signal tapers.
......@@ -195,43 +200,68 @@ namespace ts {
void FourPoint::init(const int& f, const int& l) const
{
/*
TSXX_debug(true, "FourPoint::init",
TSXX_value(f) << " " <<
TSXX_value(l) << " " <<
TSXX_value(Mt1) << " " <<
TSXX_value(Mt2) << " " <<
TSXX_value(Mt3) << " " <<
TSXX_value(Mt4));
*/
double F=static_cast<double>(f);
double L=static_cast<double>(l);
double T=F-L;
double T=L-F;
Mti1=F+T*Mt1;
Mti2=F+T*Mt2;
Mti3=F+T*Mt3;
Mti4=F+T*Mt4;
Mfac1=M_PI/(Mti2-Mti1);
Mfac2=M_PI/(Mti4-Mti3);
/*
TSXX_debug(true, "FourPoint::init",
TSXX_value(Mti1) << " " <<
TSXX_value(Mti2) << " " <<
TSXX_value(Mti3) << " " <<
TSXX_value(Mti4) << " " <<
TSXX_value(Mfac1) << " " <<
TSXX_value(Mfac2));
*/
} // void FourPoint::init(const int& f, const int& l) const
/*----------------------------------------------------------------------*/
double FourPoint::value(const int& i) const
{
//cout << i << " ";
double retval=0.;
double t=static_cast<double>(i);
if (t<Mti1)
{
//cout << "t<Mti1" << " ";
retval=0.;
}
else if (t<=Mti2)
{
//cout << "t<=Mti2" << " ";
retval=0.5*(1.-cos(Mfac1*(t-Mti1)));
}
else if (t<Mti3)
{
//cout << "t<Mti3" << " ";
retval=1.;
}
else if (t<=Mti4)
{
//cout << "t<=Mti4" << " ";
retval=0.5*(1.+cos(Mfac2*(t-Mti3)));
}
else
{
//cout << "else" << " ";
retval=0.;
}
//cout << "retval=" << retval << endl;
return (retval);
} // double FourPoint::value(const int& i) const
......
......@@ -247,6 +247,24 @@ int main(int iargc, char* argv[])
{
ts::tapers::OffsetVariableTaper ovt(true);
ovt.read(opt.ovtaperfile);
const int n=30;
const double xmax=90.;
const double dx=xmax/n;
ts::tapers::ovtaper::Picks picks=ovt.t3();
for (int i=0; i<n; ++i)
{
double x=dx*i-9.;
cout << "x=" << x << " t=" << picks.time(x) << "\n";
}
ts::tapers::FourPoint taper=ovt.taper(30., 0., 1.);
aff::Series<double> b(-50,50);
b=1;
taper.apply(b);
for (int i=b.f(); i<=b.l(); ++i)
{
cout << i << " " << b(i) << endl;
}
}
}
......
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