ovtaper.cc

Go to the documentation of this file.

#define TSXX_OVTAPER_CC_VERSION \
  "TSXX_OVTAPER_CC   V1.2"

#include <fstream>
#include <sstream>
#include <tsxx/tsxx.h>
#include <tsxx/debug.h>
#include <tsxx/ovtaper.h>
#include <tsxx/error.h>

namespace ts {

  namespace tapers {
    
    namespace ovtaper {
        
      void Picks::read(std::istream& is)
      {
        Mpicks.clear();
        std::string line;
        do {
          TSXX_assert(getline(is, line), 
                      "Picks::read: error while reading taper definition");
        } while (line.substr(0,9)!="taper set");

        TSXX_debug(Mdebug, "Picks::read",
                   "found data line:\n" << line);

        std::istringstream iss(line.substr(17,4));
        int n;
        iss >> n;
        TSXX_debug(Mdebug, "Picks::read",
                   "read " << n << " picks");

        // skip one line
        TSXX_assert(getline(is, line), 
                    "Picks::read: error while reading taper definition");
        for (int i=0; i<n; ++i)
        {
          Pick p;
          TSXX_assert(is >> p.x >> p.t,
                      "Picks::read: error while reading taper definition");
          Mpicks.push_back(p);
          TSXX_debug(Mdebug, "Picks::read",
                     "#" << i << " x=" << p.x << ", t=" << p.t);
        }

        Mpicks.sort();
      } // 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
      {
        Pick p=this->pick(offset);
        return(p.t);
      } // double Picks::time(const double& offset) const

    } // namespace ovtaper

    /*======================================================================*/

    ts::tapers::FourPoint OffsetVariableTaper::taper(const double& offset,
                                                     const double& T0,
                                                     const double& T) const
    {
      TSXX_assert(Mvalid,
                  "OffsetVariableTaper::taper: "
                  "taper is undefined");
      // calculate times in units of the duration of the time series
      double t1=(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;
      // keep values in bounds of time series
      t1 = t1 >=  0. ? t1 : 0.; 
      t2 = t2 >=  1.e-6 ? t2 : 1.e-6; 
      t3 = t3 <=  .99999 ? t3 : .99999; 
      t4 = t4 <=  1. ? t4 : 1.; 
      // initialize taper
      ts::tapers::FourPoint fpt(t1, t2, t3, t4);
      return(fpt);
    } // ts::tapers::FourPoint OffsetVariableTaper::taper(...)

    /*----------------------------------------------------------------------*/

    void OffsetVariableTaper::read(std::istream& is)
    {
      Mt1.read(is);
      Mt2.read(is);
      Mt3.read(is);
      Mt4.read(is);
      Mvalid=true;
    } // void OffsetVariableTaper::read(std::istream& is)

    /*----------------------------------------------------------------------*/

    void OffsetVariableTaper::read(const std::string& filename)
    {
      std::ifstream is(filename.c_str());
      TSXX_assert(is.good(), 
                  "OffsetVariableTaper::read(): "
                  "error when opening taper file");
      this->read(is);
    } // void OffsetVariableTaper::read(const std::string& filename)

  } // namespace tapers

} // namespace ts

/* ----- END OF ovtaper.cc ----- */