sffxx.cc

Go to the documentation of this file.

#define TF_SFFXX_CC_VERSION \
  "TF_SFFXX_CC   V1.19"

#include<sstream>
#include <sffxx.h>
#include <gsexx.h>
#include <cstdio>
#include <cmath>
#include <stdlib.h>

namespace sff {

  namespace helper {
    template<class C>
    bool IDmatch(const std::string& line, const bool& debug=false)
    { 
      if (debug) 
      {
        std::cerr << "DEBUG: IDmatch: "
          << line.substr(0,4) << "=?=" << std::string(C::LINEID)
          << std::endl;
      }
      return(line.substr(0,4)==std::string(C::LINEID)); 
    }
  } // namespace helper


  const double STAT::libversion=2.00;
  const double STAT::decode_libversion=1.09;
  const char* const STAT::LINEID="STAT";
  const char* const FREE::LINEID="FREE";
  const char* const SRCE::LINEID="SRCE";
  const char* const DAST::LINEID="DAST";
  const char* const INFO::LINEID="INFO";

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

  char coosysID(const Ecoosys& csid)
  {
    char retval;
    if (csid == CS_cartesian) retval='C';
    else if (csid == CS_spherical) retval='S';
    else throw
      GSE2::Terror("ERROR (sff::coosysID): library inconsistency!");
    return(retval);
  } // coosysID

  Ecoosys coosysID(const char& csid)
  {
    Ecoosys retval;
    if (csid=='C') { retval=CS_cartesian; }
    else if (csid=='S') { retval=CS_spherical; }
    else throw
      GSE2::Terror("ERROR (sff::coosysID): unknown coordinate system key!");
    return(retval);
  } // coosysID

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

  STAT::STAT(): hasfree(false), hassrce(false) 
  { setstamp(libtime::now()); }

  /*----------------------------------------------------------------------*/
  
  void STAT::setstamp(const libtime::TAbsoluteTime& date) const
  {
    char stamp[14];
    int yeardigits, century;
    century=int(date.year()/100);
    yeardigits=date.year()-100*century;
    sprintf(stamp, "%2.2d%2.2d%2.2d.%2.2d%2.2d%2.2d",
            static_cast<int>(yeardigits), static_cast<int>(date.month()),
            static_cast<int>(date.day()), static_cast<int>(date.hour()),
            static_cast<int>(date.minute()), static_cast<int>(date.second()));
    timestamp=std::string(stamp);
  } // void STAT::setstamp(const libtime::TAbsoluteTime& date) const

  /*----------------------------------------------------------------------*/
  
  std::string STAT::line() const
  {
    this->setstamp(libtime::now());
    char charline[40];
    std::string code("");
    if (this->hasfree) { code.append("F"); }
    if (this->hassrce) { code.append("S"); }
    sprintf(charline, "%-4s  %6.2f %-13s %-10s\n",
            STAT::LINEID, 
            STAT::libversion,
            timestamp.substr(0,13).c_str(),
            code.substr(0,10).c_str());
    std::string retval(charline);
    return(retval);
  } // std::string STAT::line() const

  /*----------------------------------------------------------------------*/
  
  void STAT::read(std::istream& fis, const bool& debug)
  {
    if (debug) { std::cerr << "DEBUG (STAT::read):" << std::endl; }
    // read complete line
    std::string theline;
    std::getline(fis, theline);

    if (debug) { std::cerr << theline << std::endl; }
    std::string lineID=theline.substr(0,4);
    if (debug) { std::cerr << lineID << std::endl; }
    if (!helper::IDmatch<STAT>(lineID)) throw
       GSE2::Terror("ERROR (STAT::read): missing STAT ID!");

    double inlibversion=atof(theline.substr(5,7).c_str());
    if (inlibversion>STAT::libversion) 
    { 
      std::cerr << "ERROR (STAT::read): "
        << "data file has format version " 
        << inlibversion << std::endl;
      std::cerr << "ERROR (STAT::read): "
        << "library can decode formats up to version " 
        << STAT::libversion << std::endl;
      throw 
        GSE2::Terror("ERROR (STAT::read): file library version too large!"); 
    }
    if (inlibversion<STAT::decode_libversion) 
    { 
//      std::cerr << "theline " << theline << std::endl;
//      std::cerr << "theline.lv " << theline.substr(5,6) << std::endl;
//      std::cerr << "inlibversion " << inlibversion << std::endl;
//      std::cerr << "decode_libversion " << decode_libversion << std::endl;
      std::cerr << "ERROR (STAT::read): "
        << "data file has format version " 
        << inlibversion << std::endl;
      std::cerr << "ERROR (STAT::read): "
        << "library can decode formats from version "
        << STAT::decode_libversion 
        << " up to version " 
        << STAT::libversion << std::endl;
      throw 
        GSE2::Terror("ERROR (STAT::read): incompatible SFF version!"); 
    }

    this->timestamp=theline.substr(13,13);

    std::string code=theline.substr(27,10);
//    std::cerr << "code " << code << std::endl;
    this->hasfree=(code.find('F')!=std::string::npos);
    this->hassrce=(code.find('S')!=std::string::npos);
    if (debug)
    {
      if (this->hasfree) 
      { std::cerr << "DEBUG (STAT::read): has FREE block" << std::endl; }
      if (this->hassrce) 
      { std::cerr << "DEBUG (STAT::read): has SRCE line" << std::endl; }
      std::cerr << "DEBUG (STAT::read): finished" << std::endl;
    }
  } // void STAT::read(std::istream& fis, const bool& debug)

/*----------------------------------------------------------------------*/
// SRCE
// ----

  SRCE::SRCE(): 
    type("NSP"), date(libtime::now()),
    cs(CS_cartesian), cx(0.), cy(0.), cz(0.) { }

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

  std::string SRCE::line() const
  {
    char charline[95];
    int yeardigits, century;
    century=int(date.year()/100);
    yeardigits=date.year()-100*century;
    sprintf(charline, "%-4s %-20s %1c %15.6f%15.6f%15.6f "
            "%2.2i%2.2i%2.2i %2.2i%2.2i%2.2i.%3.3i\n",
            SRCE::LINEID, 
            type.substr(0,20).c_str(), 
            coosysID(cs), cx, cy, cz,
            yeardigits, static_cast<int>(date.month()), static_cast<int>(date.day()),
            static_cast<int>(date.hour()), static_cast<int>(date.minute()), 
            static_cast<int>(date.second()), static_cast<int>(date.milsec()));
    std::string retval(charline);
    return(retval);
  } // std::string SRCE::line() const

  /*----------------------------------------------------------------------*/
  
  void SRCE::read(std::istream& fis, const bool& debug)
  {
    std::string theline;
    std::getline(fis, theline);

    std::string lineID=theline.substr(0,4);
    if (!helper::IDmatch<SRCE>(lineID)) throw
       GSE2::Terror("ERROR (SRCE::read): missing SRCE ID!");

    this->type=theline.substr(5,20);
    if (debug)
    { std::cerr << "DEBUG (SRCE::read): type: " << type << std::endl; }

    this->cs=coosysID(*theline.substr(26,1).c_str());
    this->cx=atof(theline.substr(28,15).c_str());
    this->cy=atof(theline.substr(43,15).c_str());
    this->cz=atof(theline.substr(58,15).c_str());
    if (debug)
    { 
      std::cerr << "DEBUG (SRCE::read): cs,cx,cy,cz: " 
        << coosysID(cs) << "," << cx << "," << cy << "," << cz << std::endl; 
    }

    std::string datestring,timestring;
    datestring=theline.substr(74,6);
    timestring=theline.substr(81,10);
    if (debug)
    { 
      std::cerr << "DEBUG (SRCE::read): datestring: " 
        << datestring << std::endl;
      std::cerr << "DEBUG (SRCE::read): timestring: " 
        << timestring << std::endl;
    }
    std::string fulldate("");
    fulldate+=datestring.substr(0,2);
    fulldate+="/";
    fulldate+=datestring.substr(2,2);
    fulldate+="/";
    fulldate+=datestring.substr(4,2);
    fulldate+=" ";
    fulldate+=timestring.substr(0,2);
    fulldate+=":";
    fulldate+=timestring.substr(2,2);
    fulldate+=":";
    fulldate+=timestring.substr(4,6);
    if (debug)
    { 
      std::cerr << "DEBUG (SRCE::read): convert string: " 
        << fulldate << std::endl; 
    }
    this->date=libtime::TAbsoluteTime(fulldate);
    if (debug)
    { 
      std::cerr << "DEBUG (SRCE::read): time: " 
        << date.timestring() << std::endl; 
    }
  } // void SRCE::read(std::istream& fis, const bool& debug)

/*----------------------------------------------------------------------*/
// DAST
// ----

  DAST::DAST(): 
    nchar(-1), ampfac(1.), hasfree(false), 
    hasinfo(false), last(false) { }

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

  std::string DAST::line() const
  {
    char charline[50];
    std::string code("");
    if (this->hasfree) { code.append("F"); }
    if (this->hasinfo) { code.append("I"); }
    if (!this->last) { code.append("D"); }
    // write -1 to nchar field in any case
    sprintf(charline, "%-4s  %10i %16.6E %-10s\n",
            DAST::LINEID, 
            -1, ampfac, code.substr(0,10).c_str());
    std::string retval(charline);
    return(retval);
  } // std::string DAST::line() const

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

  void DAST::read(std::istream& fis, const bool& debug)
  {
    if (debug)
    {
      std::cerr << "DEBUG (DAST::read):" << std::endl;
    }
    std::string theline;
    std::getline(fis, theline);
    if (debug) { std::cerr << theline << std::endl; }

    std::string lineID=theline.substr(0,4);
    if (debug) { std::cerr << lineID << std::endl; }
    if (!helper::IDmatch<DAST>(lineID)) throw
       GSE2::Terror("ERROR (DAST::read): missing DAST ID!");

    this->nchar=atoi(theline.substr(6,10).c_str());
    this->ampfac=atof(theline.substr(17,16).c_str());

    std::string code=theline.substr(34,10);
    if (debug)
    {
      std::cerr << "DEBUG (DAST): read nchar=" << nchar
        << " ampfac=" << ampfac << " code=" << code << std::endl;
    }
    this->hasinfo=(code.find('I')!=std::string::npos);
    this->hasfree=(code.find('F')!=std::string::npos);
    this->last=(code.find('D')==std::string::npos);
  } // DAST::read(std::istream& fis, const bool& debug)

/*----------------------------------------------------------------------*/
// FREE
// ----

  FREE::FREE() { lines.clear(); }

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

  void FREE::write(std::ostream& os) const
  {
    os << FREE::LINEID << " " << std::endl;
    for(Tlines::const_iterator I=lines.begin();
        I != lines.end(); I++)
    { os << *I << std::endl; }
    os << FREE::LINEID << " " << std::endl;
  } // void FREE::write(std::ostream& os) const

  /*----------------------------------------------------------------------*/
    
  void FREE::read(std::istream& is, const bool& debug)
  {
    std::string lineID;
    // getline(is,lineID);
    is >> lineID;
    if (debug) 
    { std::cerr << "DEBUG (FREE): lineID=" << lineID << std::endl; }
    if (!helper::IDmatch<FREE>(lineID.substr(0,4), debug)) throw
       GSE2::Terror("ERROR (FREE::read): missing FREE ID!");
    lines.clear();
    is.ignore(10,'\n');
    getline(is,lineID);
    if (debug) { std::cerr << "FREE " << lineID << std::endl; }
    if (debug) { std::cerr << "ID  \"" << lineID.substr(0,4) << "\"" << std::endl; }
    while (!helper::IDmatch<FREE>(lineID.substr(0,4), debug))
    {
      lines.push_back(lineID);
      getline(is,lineID);
    }
  } // void FREE::read(std::istream& is, const bool& debug)

  /*----------------------------------------------------------------------*/
  
  void FREE::append(const Tlines& newlines)
  {
    Tlines::const_iterator I(newlines.begin());
    while (I != newlines.end())
    {
      this->append(*I);
      ++I;
    }
  } // void FREE::append(const Tlines& lines)

/*----------------------------------------------------------------------*/
// WID2
// ----

  WID2::WID2():
    date(libtime::now()), station("NSP"), channel("NSP"), auxid("NSP"),
    nsamples(-1), dt(-1.), calib(-1.), calper(-1.), instype("NSP"),
    hang(-1.),vang(-1.) { }

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

  std::string WID2::line() const
  {
    GSE2::waveform::TWID2 wid2line;
    // apply proper rounding to nearest millisecond
    libtime::TAbsoluteTime rnddate(date.year(), date.month(),
                                   date.day(), date.hour(),
                                   date.minute(), date.second(),
                                   date.milsec());
    if (date.micsec()>499) { rnddate += libtime::TRelativeTime(0,0,0,0,1); }
    wid2line.Fyear=rnddate.year();
    wid2line.Fmonth=rnddate.month();
    wid2line.Fday=rnddate.day();
    wid2line.Fhour=rnddate.hour();
    wid2line.Fminute=rnddate.minute();
    wid2line.Fmilsec=rnddate.second()*1000+rnddate.milsec();
    wid2line.Fstation=station;
    wid2line.Fchannel=channel;
    wid2line.Fauxid=auxid;
    wid2line.Fsamps=nsamples;
    wid2line.Fsamprate=1./dt;
    wid2line.Fcalib=calib;
    wid2line.Fcalper=calper;
    wid2line.Finstype=instype;
    wid2line.Fhang=hang;
    wid2line.Fvang=vang;
    return(wid2line.line());
  } // std::string WID2::line() const

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

  void WID2::read(std::istream& is)
  {
    std::string line;
    std::getline(is, line);
    if (line.substr(0,4) == WIDXID)
    {
      *this = WIDXline(line);
    }
    else
    {
      std::istringstream iss(line);
      GSE2::waveform::TWID2 wid2line;
      wid2line.read(iss);
      int second=int(wid2line.Fmilsec/1000);
      int milsec=wid2line.Fmilsec-1000*second;
      date=libtime::TAbsoluteTime(wid2line.Fyear,
                                  wid2line.Fmonth,
                                  wid2line.Fday,
                                  wid2line.Fhour,
                                  wid2line.Fminute,
                                  second, milsec);
      this->station=wid2line.Fstation;
      this->channel=wid2line.Fchannel;
      this->auxid=wid2line.Fauxid;
      this->nsamples=wid2line.Fsamps;
      this->dt=1./wid2line.Fsamprate;
      this->calib=wid2line.Fcalib;
      this->calper=wid2line.Fcalper;
      this->instype=wid2line.Finstype;
      this->hang=wid2line.Fhang;
      this->vang=wid2line.Fvang;
    }
  } // void WID2::read(std::istream& is)

/*----------------------------------------------------------------------*/
// INFO
// ----

  INFO::INFO():
    cs(CS_cartesian), cx(0.), cy(0.), cz(0.), nstacks(0) { }

  /*----------------------------------------------------------------------*/
  
  std::string INFO::line() const
  {
    char charline[60];
    sprintf(charline, "%-4s %1c %15.6f%15.6f%15.6f %4i\n",
            INFO::LINEID, 
            coosysID(cs), cx, cy, cz, nstacks);
    std::string retval(charline);
    return(retval);
  } // std::string INFO::line() const

  /*----------------------------------------------------------------------*/
  
  void INFO::read(std::istream& fis)
  {
    std::string theline;
    std::getline(fis, theline);

    std::string lineID=theline.substr(0,4);
    if (!helper::IDmatch<INFO>(lineID)) throw
       GSE2::Terror("ERROR (INFO::read): missing INFO ID!");

    this->cs=coosysID(*theline.substr(5,1).c_str());
    this->cx=atof(theline.substr(7,15).c_str());
    this->cy=atof(theline.substr(22,15).c_str());
    this->cz=atof(theline.substr(37,15).c_str());
    this->nstacks=atoi(theline.substr(53,4).c_str());
  } // void INFO::read(std::istream& fis)

  /*----------------------------------------------------------------------*/
  
  bool INFO::operator==(const INFO& info) const
  {
    return ((info.cs==this->cs) &&
            (info.cx==this->cx) &&
            (info.cy==this->cy) &&
            (info.cz==this->cz) &&
            (info.nstacks==this->nstacks));
  } // bool INFO::operator==(const INFO& info) const

/*----------------------------------------------------------------------*/
// FileHeader
// ----------

  void FileHeader::write(std::ostream& os) const
  {
    os << Mstat.line();
    if (Mstat.hasfree) { Mfree.write(os); }
    if (Mstat.hassrce) { os << Msrce.line(); }
  }

  /*----------------------------------------------------------------------*/
  
  void FileHeader::read(std::istream& is, const bool& debug)
  {
    // reading the FREE block first matches the order of header elements as
    // written by sff_WOpenFS in libstuff.f
    Mstat.read(is, debug);
    if (Mstat.hasfree) 
    { 
      Mfree.read(is, debug); 
      if (debug)
      { std::cerr << "DEBUG (FileHeader::read): file FREE read" << std::endl; }
    }
    if (Mstat.hassrce) 
    { 
      Msrce.read(is, debug); 
      if (debug)
      { std::cerr << "DEBUG (FileHeader::read): SRCE line read" << std::endl; }
    }
    if (debug)
    { std::cerr << "DEBUG (FileHeader::read): finished" << std::endl; }
  }

/*----------------------------------------------------------------------*/
// TraceHeader
// -----------

  void TraceHeader::writeheader(std::ostream& os) const
  {
    if (Mdebug) { std::cerr << "DEBUG: write DAST line" << std::endl; } 
    os << Mdast.line();
    if (Mdebug) { std::cerr << "DEBUG: write WID2 line" << std::endl; } 
    os << Mwid2.line();
  }

  /*----------------------------------------------------------------------*/
  
  void TraceHeader::writetrailer(std::ostream& os) const
  {
    if (Mdast.hasfree) 
    {
      if (Mdebug) { std::cerr << "DEBUG: write FREE block" << std::endl; } 
      Mfree.write(os);
    }
    if (Mdast.hasinfo) 
    {
      if (Mdebug) { std::cerr << "DEBUG: write INFO line" << std::endl; } 
      os << Minfo.line();
    }
  }

  /*----------------------------------------------------------------------*/
  
  void TraceHeader::readheader(std::istream& is) 
  {
    if (Mdebug) { std::cerr << "DEBUG: read DAST line" << std::endl; } 
    Mdast.read(is,Mdebug);
    if (Mdebug) { std::cerr << "DEBUG: read WID2 line" << std::endl; } 
    Mwid2.read(is);
  }

  /*----------------------------------------------------------------------*/
  
  void TraceHeader::readtrailer(std::istream& is) 
  {
    if (Mdebug) { std::cerr << "DEBUG: read trace trailer" << std::endl; } 
    if (Mdast.hasfree) 
    {
      if (Mdebug) { std::cerr << "DEBUG: read FREE block" << std::endl; } 
      Mfree.read(is, Mdebug); 
      if (Mdebug) { Mfree.write(std::cerr); }
    }
    if (Mdast.hasinfo) 
    {
      if (Mdebug) { std::cerr << "DEBUG: read INFO line" << std::endl; } 
      Minfo.read(is);
    }
  }

/*----------------------------------------------------------------------*/
// WaveformNormalizer
// ------------------

  const int WaveformNormalizer::limit=0x1ffffff;

  WaveformNormalizer::WaveformNormalizer(const Enormmode& nm, 
                                         const double& maxval):
    Mmaxval(maxval), Mnorm(nm)
  {
    if (Mnorm == NM_one) 
    { 
      Mampfac=1.; 
      Mscale=false;
      if (Mmaxval > double(WaveformNormalizer::limit)) throw
        GSE2::Terror("ERROR (sff::WaveformNormalizer::scan): "
                     "dynamic range to large for non-normalizing mode");
    }
    else if (Mnorm == NM_ifneeded)
    { 
      Mampfac=1.;
      Mscale=false;
      if (Mmaxval > double(WaveformNormalizer::limit))
      {
        Mampfac=Mmaxval/double(WaveformNormalizer::limit); 
        Mscale=true;
      }
    }
    else if (Mnorm == NM_maxdyn)
    { 
      Mampfac=Mmaxval/double(WaveformNormalizer::limit); 
      Mscale=true;
    }
    else throw
      GSE2::Terror("ERROR (sff::WaveformNormalizer::scan): "
                   "library inconsistency!");
  }

/*----------------------------------------------------------------------*/
// SkipWaveform
// ------------

    void SkipWaveform::read(std::istream& is) 
    {
      Mheader.readheader(is);;
      int nsamples=Mheader.wid2().nsamples;
      GSE2::waveform::TDAT2readCM6 freader(nsamples);
      int idata;
      for (int i=0; i<nsamples; i++)
      { idata=freader(is); }
      Mheader.readtrailer(is);
      Mvalid=true;
    } // SkipWaveform::read

} // namespace sff

/* ----- END OF sffxx.cc ----- */