sffxx.h

Go to the documentation of this file.

// include guard
#ifndef TF_SFFXX_H_VERSION

#define TF_SFFXX_H_VERSION \
  "TF_SFFXX_H   V1.10"

#include<string>
#include<cmath>
#include<list>
#include<iostream>
#include<sstream>
#include<libtime++.h>
#include<aff/iterator.h>
//#include<tfxx/commandline.h>
#include<gsexx.h>

namespace sff {

/*======================================================================*/
// basic error handling
// --------------------

  class Terror: public GSE2::Terror 
  {
    public:
      Terror(const std::string& message): GSE2::Terror(message) { }
  }; // class Terror

#define SFF_assert( C , M ) \
  if (!(C)) { \
    std::ostringstream os; \
    os << "Condition " << #C << " at line #" << __LINE__  \
      << " in " << __FILE__ << " is false: " << M; \
    throw(Terror(os.str())); \
  }

#define SFF_abort( M ) \
  { \
    std::ostringstream os; \
    os << "ABORT at line #" << __LINE__  \
      << " in " << __FILE__ << " because of " << M; \
    throw(Terror(os.str())); \
  }

/*======================================================================*/
// enum
// ----
  
  enum Ecoosys {
    CS_cartesian,
    CS_spherical
  }; // enum Ecoosys

  char coosysID(const Ecoosys& csid);
  Ecoosys coosysID(const char& csid);

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

  enum Enormmode {
    NM_one,       
    NM_maxdyn,    
    NM_ifneeded   
  }; // enum Enormmode

/*======================================================================*/
// constants
// ---------

  const char *const WIDXID = "WIDX";

/*======================================================================*/
// SFF structs
// -----------

  struct STAT {
    static const double libversion;
    static const double decode_libversion;
    static const char* const LINEID;
    STAT();
    STAT(std::istream& is, const bool& debug=false) { read(is, debug); }
    std::string line() const;
    void read(std::istream& is, const bool& debug=false);
    void setstamp(const libtime::TAbsoluteTime& date) const;
    public:
      mutable std::string timestamp;
      bool hasfree;
      bool hassrce;
  }; // struct STAT

  struct FREE {
    typedef std::list<std::string> Tlines;
    static const char* const LINEID;
    FREE();
    FREE(std::istream& is) { read(is); }
    void write(std::ostream& os) const;
    void read(std::istream& is, const bool& debug=false);
    void append(const std::string& line) { lines.push_back(line); }
    void append(const Tlines& newlines);
    void append(const FREE& free) { this->append(free.lines); }
    Tlines lines;
  }; // struct FREE

  struct SRCE {
    static const char* const LINEID;
    SRCE();
    SRCE(std::istream& is) { read(is); }
    std::string line() const;
    void read(std::istream& is, const bool& debug=false);
    public:
      std::string type;
      libtime::TAbsoluteTime date;        
      Ecoosys cs;
      double cx, cy, cz;
  }; // struct SRCE

  struct DAST {
    static const char* const LINEID;
    DAST();
    DAST(std::istream& is) { read(is); }
    std::string line() const;
    void read(std::istream& is, const bool& debug=false);
    public:
      int nchar;
      double ampfac;
      bool hasfree;
      bool hasinfo;
      bool last;
  }; // struct DAST

  struct INFO {
    static const char* const LINEID;
    INFO();
    INFO(std::istream& is) { read(is); }
    std::string line() const;
    void read(std::istream& is);
    bool operator==(const INFO& info) const;
    public:
      Ecoosys cs;
      double cx, cy, cz;
      int nstacks;
  }; // struct INFO

  struct WID2 {
    WID2();
    WID2(std::istream& is) { read(is); }
    std::string line() const;
    void read(std::istream& is);
    public:
      libtime::TAbsoluteTime date;       
      std::string           station;     
      std::string           channel;     
      std::string           auxid;       
      int                   nsamples;    
      double                dt;          
      double                calib;       
      double                calper;      
      std::string           instype;     
      double                hang;        
      double                vang;        
  }; // struct WID2

/*======================================================================*/
// file I/O classes
// ----------------

  class FileHeader {
    public:
      FileHeader()
      { Mstat.hasfree=false; Mstat.hassrce=false; }
      FileHeader(const FREE& free): 
        Mfree(free) 
      { Mstat.hasfree=true; Mstat.hassrce=false; }
      FileHeader(const SRCE& srce): 
        Msrce(srce)
      { Mstat.hasfree=false; Mstat.hassrce=true; }
      FileHeader(const SRCE& srce, const FREE& free):
        Mfree(free), Msrce(srce)
      { Mstat.hasfree=true; Mstat.hassrce=true; }
      FileHeader(std::istream& is, const bool& debug=false) 
      { read(is, debug); }
      void write(std::ostream&) const;
      void read(std::istream&, const bool& debug=false);
      const STAT& stat() const { return(Mstat); }
      const FREE& free() const { return(Mfree); }
      const SRCE& srce() const { return(Msrce); }
      const bool& hassrce() const { return(Mstat.hassrce); }
      const bool& hasfree() const { return(Mstat.hasfree); }
      void appendfree(const std::string& line)
      { Mstat.hasfree=true; Mfree.lines.push_back(line); }
      void setfree(const FREE& free)
      { Mstat.hasfree=true; Mfree=free; }
      void setsrce(const SRCE& srce)
      { Mstat.hassrce=true; Msrce=srce; }
    private:
      STAT Mstat;
      FREE Mfree;
      SRCE Msrce;
  }; // class FileHeader

  class WaveformNormalizer {
    public:
      static const int limit;
      WaveformNormalizer(const Enormmode& nm, const double& maxval); 
      const double& maxval() const { return(Mmaxval); }
      const double& ampfac() const { return(Mampfac); }
      const bool& scale() const { return(Mscale); }
    private:
      double Mampfac;
      double Mmaxval;
      bool Mscale;
      Enormmode Mnorm;
  }; // class WaveformNormalizer

  class TraceHeader {
    public:
      TraceHeader(): Mdebug(false) { }
      TraceHeader(const WID2& wid2, const bool& last=false): 
        Mwid2(wid2), Mdebug(false)
      { Mdast.hasfree=false; Mdast.hasinfo=false; Mdast.last=last; }
      TraceHeader(const WID2& wid2, const FREE& free, const bool& last=false): 
        Mwid2(wid2), Mfree(free), Mdebug(false)
      { Mdast.hasfree=true; Mdast.hasinfo=false; Mdast.last=last; }
      TraceHeader(const WID2& wid2, const INFO& info, const bool& last=false): 
        Mwid2(wid2), Minfo(info), Mdebug(false)
      { Mdast.hasfree=false; Mdast.hasinfo=true; Mdast.last=last; }
      TraceHeader(const WID2& wid2, const INFO& info, 
                  const FREE& free, const bool& last=false):
        Mwid2(wid2), Mfree(free), Minfo(info), Mdebug(false)
      { Mdast.hasfree=true; Mdast.hasinfo=true; Mdast.last=last; }
      void writeheader(std::ostream&) const;
      void writetrailer(std::ostream&) const;
      void readheader(std::istream&);
      void readtrailer(std::istream&);
      bool last() const { return(Mdast.last); }
      template<class C> void scanseries(const C&, 
                                        const Enormmode& nm=NM_maxdyn);
      const WID2& wid2() const { return(Mwid2); }
      const DAST& dast() const { return(Mdast); }
      const FREE& free() const { return(Mfree); }
      const INFO& info() const { return(Minfo); }
      const bool& scale() const { return(Mscale); }
      const bool& hasinfo() const { return(Mdast.hasinfo); }
      const bool& hasfree() const { return(Mdast.hasfree); }
      void setlast(const bool& flag) { Mdast.last=flag; }
      void setdebug(const bool& flag) { Mdebug=flag; }
      void setwid2(const WID2& wid2line) { Mwid2=wid2line; }
      void setnsamples(const long int& n) { Mwid2.nsamples=n; }
      void setinfo(const INFO& infoline) 
      { Mdast.hasinfo=true; Minfo=infoline; }
      void setfree(const FREE& free)
      { Mdast.hasfree=true; Mfree=free; }
      void appendfree(const std::string& line)
      { Mdast.hasfree=true; Mfree.lines.push_back(line); }
    private:
      WID2 Mwid2;
      DAST Mdast;
      FREE Mfree;
      INFO Minfo;
      bool Mscale;
      bool Mdebug;
  }; // class TraceHeader

  template<class C>
    class OutputWaveform {
      public:
        typedef typename C::Tcoc Tcoc;
        OutputWaveform(const Tcoc& c, const TraceHeader& th,
                       const Enormmode& nm=NM_maxdyn):
          Mseries(c), Mheader(th) 
          { Mheader.scanseries(Mseries,nm); }
        void write(std::ostream& os) const;
      private:
        Tcoc Mseries;
        mutable TraceHeader Mheader;
    }; // class Waveform

  template<class C>
    class InputWaveform {
      public:
        typedef C Tcontainer;
        InputWaveform(const bool& debug=false): 
          Mvalid(false), Mdebug(debug) { }
        InputWaveform(std::istream& is, const bool& debug=false):
          Mdebug(debug) { this->read(is); }
        void read(std::istream& is);
        const bool& valid() const { return(Mvalid); }
        Tcontainer series() const { return(Mseries); }
        TraceHeader header() const { return(Mheader); }
        const bool& last() const { return(Mheader.dast().last); }
      private:
        bool Mvalid;
        Tcontainer Mseries;
        TraceHeader Mheader;
        bool Mdebug;
    }; // class InputWaveform

  class SkipWaveform {
    public:
      SkipWaveform(): Mvalid(false) { }
      SkipWaveform(std::istream& is) { this->read(is); }
      void read(std::istream& is);
      const bool& valid() const { return(Mvalid); }
      TraceHeader header() const { return(Mheader); }
      const bool& last() const { return(Mheader.dast().last); }
    private:
      bool Mvalid;
      TraceHeader Mheader;
  }; // class SkipWaveform

/*======================================================================*/
// I/O operators
// -------------

  inline std::istream& operator >> (std::istream& is, FileHeader& fh)
  { fh.read(is); return(is); }
  template<class C>
  inline std::istream& operator >> (std::istream& is, InputWaveform<C>& wf)
  { wf.read(is); return(is); }
  inline std::istream& operator >> (std::istream& is, SkipWaveform& swf)
  { swf.read(is); return(is); }
  inline std::ostream& operator << (std::ostream& os, const FileHeader& fh)
  { fh.write(os); return(os); }
  inline std::ostream& operator << (std::ostream& os, const TraceHeader& th)
  { 
    th.writeheader(os);
    os << "DATA!" << std::endl;
    th.writetrailer(os);
    return(os);
  }
  template<class C> 
    inline std::ostream& operator << (std::ostream& os, const
                                      OutputWaveform<C>& wf)
  { wf.write(os); return(os); }

/*======================================================================*/
// WaveformNormalizer template functions
// -------------------------------------

  template<class C>
    inline
    void TraceHeader::scanseries(const C& c,
                                 const Enormmode& nm)
    {
      Mwid2.nsamples=0;
      double maxval=0.;
      double absval, value;
      double null(0);
      for(aff::Browser<C> i(c); i.valid(); ++i)
      { 
        Mwid2.nsamples++;
        value= *i;
        absval= (value < null) ? -value : value;
        maxval= (maxval < absval) ? absval : maxval;
      }
      double dmaxval=double((maxval==0) ? WaveformNormalizer::limit : maxval);
      WaveformNormalizer normalizer(nm, dmaxval);
      Mdast.ampfac=normalizer.ampfac();
      Mscale=normalizer.scale();
    } // TraceHeader::scanseries

/*----------------------------------------------------------------------*/
// template OutputWaveform functions
// ---------------------------------

  template<class C>
    inline
    void OutputWaveform<C>::write(std::ostream& os) const
    {
      Mheader.setnsamples(Mseries.size());
      Mheader.writeheader(os);;
      GSE2::waveform::TDAT2writeCM6 fwriter(Mheader.wid2().nsamples);
      int idata;
      typename C::Tvalue data;
      for(aff::Browser<Tcoc> i(Mseries); i.valid(); ++i)
      { 
        data= *i;
        if (Mheader.scale())
        { idata=int(round(data/Mheader.dast().ampfac)); }
        else
        { idata=int(round(data)); }
        os << fwriter(idata);
      }
      Mheader.writetrailer(os);
    } // OutputWaveform::write

/*----------------------------------------------------------------------*/
// template InputWaveform functions
// --------------------------------

  template<class C>
    inline
    void InputWaveform<C>::read(std::istream& is) 
    {
      typedef typename C::Tvalue Tvalue;
      if (Mdebug)
      {
        std::cerr << "DEBUG (InputWaveform<C>::read): " << std::endl
          << "  calling function readheader() of member Mheader." 
          << std::endl;
      }
      Mheader.setdebug(Mdebug);
      Mheader.readheader(is);
      int nsamples=Mheader.wid2().nsamples;
      GSE2::waveform::TDAT2readCM6 freader(nsamples);
      try {
        Mseries=C(nsamples);
      }
      catch(...) {
        std::cerr << "ERROR (InputWaveform::read): "
          << "allocating series for " << nsamples << " samples!" << std::endl;
        throw;
      }
      for(aff::Iterator<C> i(Mseries); i.valid(); ++i)
      { (*i) = Tvalue(freader(is)*Mheader.dast().ampfac); }
      Mheader.readtrailer(is);
      Mvalid=true;
    } // InputWaveform::read

/*======================================================================*/
// some utilities
// --------------


  /*----------------------------------------------------------------------*/
  // compare WID2 headers
  // --------------------
  enum Ewid2field {
    Fdate    =1<<0,
    Fstation =1<<1,
    Fchannel =1<<2,
    Fauxid   =1<<3,
    Fnsamples=1<<4,
    Fdt      =1<<5,
    Fcalib   =1<<6,
    Fcalper  =1<<7,
    Finstype =1<<8,
    Fhang    =1<<9,
    Fvang    =1<<10
  }; // enum Ewid2field

  class WID2compare {
    public:
      WID2compare(const int& flags=(Fstation | Fchannel | Fdt)):
        Mflags(flags), Mdttolerance(0.), Mdatetolerance(0.) { }
      void set(const int& flags) { Mflags=Mflags | flags; }
      void clear(const int& flags) { Mflags=Mflags & (0xffffffff ^ flags); }
      void setdttolerance(const double& tol) { Mdttolerance=tol; }
      void setdatetolerance(const double& tol) { Mdatetolerance=tol; }
      bool operator()(const WID2& hd1, const WID2& hd2) const;
      int flags() const { return(Mflags); }
      double dttolerance() const { return(Mdttolerance); }
      double datetolerance() const { return(Mdatetolerance); }
    private:
      int Mflags;
      double Mdttolerance;
      double Mdatetolerance; 
  }; // class WID2compare

  /*======================================================================*/
  // functions
  // ---------
  
  libtime::TAbsoluteTime wid2lastsample(const WID2& wid2);
  libtime::TAbsoluteTime wid2nextdate(const WID2& wid2);
  long int wid2isample(const WID2& wid2, 
                       const libtime::TAbsoluteTime& idate);
  libtime::TAbsoluteTime wid2isample(const WID2& wid2, 
                                     const long int& i);
  libtime::TRelativeTime wid2isamplerest(const WID2& wid2, 
                                         const libtime::TAbsoluteTime& idate);

  std::string srce_reference_ID();
  sff::SRCE srce_reference();

  double offset(const SRCE& srce, const INFO& info,
                const double& radius=6371.);
  double sourcedistance(const SRCE& srce, const INFO& info);
  double offsetdeg(const SRCE& srce, const INFO& info, 
                   const double& radius=6371.);

  std::string WIDXline(const sff::WID2& wid2, const bool& debug=false);

  sff::WID2 WIDXline(const std::string& line);

  /*======================================================================*/
  // verbose output
  // --------------

  void verbose(std::ostream& os, const WID2& wid2);
  void verbose(std::ostream& os, const SRCE& srce);
  void verbose(std::ostream& os, const DAST& dast);
  void verbose(std::ostream& os, const INFO& info);
  void verbose(std::ostream& os, const FREE& free);
  void verbose(std::ostream& os, const STAT& stat);
  void verbose(std::ostream& os, const FileHeader& fh);
  void verbose(std::ostream& os, const TraceHeader& th);

} // namespace sff

#endif // TF_SFFXX_H_VERSION (includeguard)

/* ----- END OF sffxx.h ----- */