mspass::seismic::CoreSeismogram Class Reference

Vector (three-component) seismogram data object. More...

#include <CoreSeismogram.h>

Inheritance diagram for mspass::seismic::CoreSeismogram:

Collaboration diagram for mspass::seismic::CoreSeismogram:

Public Member Functions
	CoreSeismogram ()
	CoreSeismogram (const size_t nsamples)
	CoreSeismogram (const std::vector< mspass::seismic::CoreTimeSeries > &ts, const unsigned int component_to_clone=0)
	CoreSeismogram (const mspass::utility::Metadata &md, const bool load_data=true)
	Construct from Metadata definition that includes data path.
	CoreSeismogram (const CoreSeismogram &)
void	set_dt (const double sample_interval)
	Set the sample interval.
void	set_npts (const size_t npts)
	Set the number of samples attribute for data.
void	sync_npts ()
	Sync the number of samples attribute with actual data size.
void	set_t0 (const double t0in)
	Set the data start time.
CoreSeismogram &	operator= (const CoreSeismogram &)
CoreSeismogram &	operator+= (const CoreSeismogram &d)
	Summation operator.
const CoreSeismogram	operator+ (const CoreSeismogram &other) const
CoreSeismogram &	operator*= (const double)
CoreSeismogram &	operator-= (const CoreSeismogram &d)
	Subtraction operator.
const CoreSeismogram	operator- (const CoreSeismogram &other) const
std::vector< double >	operator[] (const int sample) const
std::vector< double >	operator[] (const double time) const
	Overloaded version of operator[] for time.
virtual	~CoreSeismogram ()
void	rotate_to_standard ()
void	rotate (mspass::utility::SphericalCoordinate &sc)
void	rotate (const double nu[3])
void	rotate (const double phi)
	Rotate horizontals by a simple angle in degrees.
void	transform (const double a[3][3])
void	free_surface_transformation (const mspass::seismic::SlownessVector u, const double vp0, const double vs0)
mspass::utility::dmatrix	get_transformation_matrix () const
bool	set_transformation_matrix (const mspass::utility::dmatrix &A)
	Define the transformaton matrix.
bool	set_transformation_matrix (const double a[3][3])
	Define the transformaton matrix with a C style 3x3 matrix.
bool	set_transformation_matrix (pybind11::object a)
	Define the transformaton matrix with a python object.
bool	cardinal () const
bool	orthogonal () const
double	endtime () const noexcept
Public Member Functions inherited from mspass::seismic::BasicTimeSeries
	BasicTimeSeries ()
	BasicTimeSeries (const BasicTimeSeries &)
virtual	~BasicTimeSeries ()
	Virtual destructor.
double	time (const int i) const
int	sample_number (double t) const
double	endtime () const noexcept
bool	shifted () const
double	get_t0shift () const
double	time_reference () const
void	force_t0_shift (const double t)
	Force a t0 shift value on data.
virtual void	ator (const double tshift)
virtual void	rtoa ()
virtual void	shift (const double dt)
bool	live () const
bool	dead () const
void	kill ()
void	set_live ()
double	dt () const
bool	time_is_UTC () const
bool	time_is_relative () const
TimeReferenceType	timetype () const
double	samprate () const
size_t	npts () const
double	t0 () const
std::vector< double >	time_axis () const
void	set_tref (const TimeReferenceType newtref)
	Force the time standard.
BasicTimeSeries &	operator= (const BasicTimeSeries &parent)
Public Member Functions inherited from mspass::utility::Metadata
	Metadata ()
	Metadata (std::ifstream &ifs, const std::string form=std::string("pf"))
	Metadata (const Metadata &mdold)
virtual	~Metadata ()
Metadata &	operator= (const Metadata &mdold)
Metadata &	operator+= (const Metadata &rhs) noexcept
const Metadata	operator+ (const Metadata &other) const
double	get_double (const std::string key) const override
int	get_int (const std::string key) const override
long	get_long (const std::string key) const
std::string	get_string (const std::string key) const override
bool	get_bool (const std::string key) const override
template<typename T >
T	get (const std::string key) const
template<typename T >
T	get (const char *key) const
	Generic get interface for C char array.
boost::any	get_any (const std::string key) const
std::string	type (const std::string key) const
template<typename T >
void	put (const std::string key, T val) noexcept
template<typename T >
void	put (const char *key, T val) noexcept
void	put (const std::string key, const double val) override
void	put (const std::string key, const int val) override
void	put (const std::string key, const bool val) override
void	put (const std::string key, const std::string val) override
void	put (const char *key, const char *val)
void	put (std::string key, const char *val)
void	put_object (const std::string key, const pybind11::object val)
void	put_int (const std::string key, const int val)
void	put_string (const std::string key, const std::string val)
void	put_bool (const std::string key, const bool val)
void	put_double (const std::string key, const double val)
void	put_long (const std::string key, const long val)
void	append_chain (const std::string key, const std::string val, const std::string separator=std::string(":"))
std::set< std::string >	modified () const
void	clear_modified ()
	Mark all data as unmodified.
std::set< std::string >	keys () const noexcept
bool	is_defined (const std::string key) const noexcept
void	erase (const std::string key)
std::size_t	size () const noexcept
std::map< std::string, boost::any >::const_iterator	begin () const noexcept
std::map< std::string, boost::any >::const_iterator	end () const noexcept
void	change_key (const std::string oldkey, const std::string newkey)
	Change the keyword to access an attribute.
template<>
double	get (const std::string key) const
template<>
int	get (const std::string key) const
template<>
long	get (const std::string key) const
template<>
bool	get (const std::string key) const
template<>
float	get (const std::string key) const
template<>
double	get (const string key) const
template<>
int	get (const string key) const
template<>
long	get (const string key) const
template<>
bool	get (const string key) const
template<>
float	get (const string key) const
Public Member Functions inherited from mspass::utility::BasicMetadata
virtual	~BasicMetadata ()

Public Attributes
mspass::utility::dmatrix	u

Protected Attributes
bool	components_are_orthogonal
bool	components_are_cardinal
double	tmatrix [3][3]
Protected Attributes inherited from mspass::seismic::BasicTimeSeries
bool	mlive
double	mdt
double	mt0
size_t	nsamp
TimeReferenceType	tref
bool	t0shift_is_valid
double	t0shift
Protected Attributes inherited from mspass::utility::Metadata
std::map< std::string, boost::any >	md
std::set< std::string >	changed_or_set

Detailed Description

Vector (three-component) seismogram data object.

A three-component seismogram is a common concept in seismology. The concept used here is that a three-component seismogram is a time series with a 3-vector as the data at each time step. As a result the data are stored internally as a matrix with row defining the component number (C indexing 0,1,2) and the column defining the time variable. The object inherits common concepts of a time series through the BasicTimeSeries object. Auxiliary parameters are defined for the object through inheritance of a Metadata object.

Author

Gary L. Pavlis

Constructor & Destructor Documentation

CoreSeismogram() [1/5]

mspass::seismic::CoreSeismogram::CoreSeismogram

(

)

Default constructor.

Sets ns to zero and builds an empty data matrix. The live variable in BasicTimeSeries is also set false.

                               : BasicTimeSeries(), Metadata() {
  /* mlive and tref are set in BasicTimeSeries so we don't use putters for
  them here.   These three initialize Metadata properly for these attributes*/
  this->set_dt(1.0);
  this->set_t0(0.0);
  this->set_npts(0);
  components_are_orthogonal = true;
  components_are_cardinal = true;
  for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j)
      if (i == j)
        tmatrix[i][i] = 1.0;
      else
        tmatrix[i][j] = 0.0;
}

References components_are_cardinal, components_are_orthogonal, mspass::utility::Metadata::get(), set_dt(), set_npts(), set_t0(), and tmatrix.

CoreSeismogram() [2/5]

mspass::seismic::CoreSeismogram::CoreSeismogram

(

const size_t

nsamples

)

Simplest parameterized constructor.

Initializes data and sets aside memory for matrix of size 3xnsamples. The data matrix is not initialized and the object is marked as not live.

Parameters

nsamples

number of samples expected for holding data.

    : BasicTimeSeries(), Metadata() {
  /* IMPORTANT:  this constructor assumes BasicTimeSeries initializes the
  equivalent of:
  set_dt(1.0)
  set_t0(0.0)
  set_tref(TimeReferenceType::Relative)
  this->kill() - i.e. marked dead
  */
  this->set_npts(
      nsamples); // Assume this is an allocator of the 3xnsamples matrix
  components_are_orthogonal = true;
  components_are_cardinal = true;
  for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j)
      if (i == j)
        tmatrix[i][i] = 1.0;
      else
        tmatrix[i][j] = 0.0;
}

References components_are_cardinal, components_are_orthogonal, mspass::utility::Metadata::get(), set_npts(), and tmatrix.

CoreSeismogram() [3/5]

mspass::seismic::CoreSeismogram::CoreSeismogram	(	const std::vector< mspass::seismic::CoreTimeSeries > &	ts,
		const unsigned int	component_to_clone = 0
	)

Construct a three component seismogram from three TimeSeries objects.

A three component seismogram is commonly assembled from individual single channel components. This constructor does the process taking reasonable care to deal with (potentially) irregular start and end times of the individual components. If the start and end times are all the same it uses a simple copy operation. Otherwise it runs a more complicated (read much slower) algorithm that handles the ragged start and stop times by adding a marked gap.

If start or end times are not constant the algorithm shortens the output to the latest start time and earliest end time respectively.

Note this constructor requires variables hang and vang, which are orientation angles defined in the CSS3.0 schema (NOT spherical coordinates by the way), by set for each component. This is used to construct the transformation matrix for the object that allows, for example, removing raw data orientation errors using rotate_to_standard. The constructor will throw an exception if any component does not have these attributes set in their Metadata area.

Exceptions

SeisppError

exception can be throw for a variety of serious problems.

Parameters

ts	vector of 3 TimeSeries objects to be used to assemble this Seismogram. Input vector order could be arbitrary because a transformation matrix is computed, but for efficiency standard order (E,N,Z) is advised.
component_to_clone	the auxiliary parameters (Metadata and BasicTimeSeries common parameters) from one of the components is cloned to assure common required parameters are copied to this object. This argument controls which of the three components passed through ts is used. Default is 0.

    : BasicTimeSeries(
          dynamic_cast<const BasicTimeSeries &>(ts[component_to_clone])),
      Metadata(dynamic_cast<const Metadata &>(ts[component_to_clone])), u() {
  const string base_error("CoreSeismogram constructor from 3 Time Series:  ");
  int i, j;
  /* This is needed in case nsamp does not match s.size(0) */
  int nstest = ts[component_to_clone].s.size();
  if (nsamp != nstest)
    this->nsamp = nstest;
  /* this method allocates u and sets the proper metadata for npts*/
  this->CoreSeismogram::set_npts(this->nsamp);
  /* beware irregular sample rates, but don' be too machevelian.
         Abort only if the mismatch is large defined as accumulated time
         over data range of this constructor is less than half a sample */
  if ((ts[0].dt() != ts[1].dt()) || (ts[1].dt() != ts[2].dt())) {
    double ddtmag1 = fabs(ts[0].dt() - ts[1].dt());
    double ddtmag2 = fabs(ts[1].dt() - ts[2].dt());
    double ddtmag;
    if (ddtmag1 > ddtmag1)
      ddtmag = ddtmag1;
    else
      ddtmag = ddtmag2;
    ddtmag1 = fabs(ts[0].dt() - ts[2].dt());
    if (ddtmag1 > ddtmag)
      ddtmag = ddtmag1;
    double ddtcum = ddtmag * ((double)ts[0].s.size());
    if (ddtcum > (ts[0].dt()) / 2.0) {
      stringstream ss;
      ss << base_error << "Sample intervals of components are not consistent"
         << endl;
      for (int ie = 0; ie < 3; ++ie)
        ss << "Component " << ie << " dt=" << ts[ie].dt() << " ";
      ss << endl;
      throw MsPASSError(ss.str(), ErrorSeverity::Invalid);
    }
  }
  // temporaries to hold component values
  double t0_component[3];
  double hang[3];
  double vang[3];
  // Load up these temporary arrays inside this try block and arrange to
  // throw an exception if required metadata are missing
  try {
    /* WARNING hang and vang attributes in Metadata
    are always assumed to have been read from a database where they
    were stored in degrees.  We convert these to radians below to
    compute the transformation matrix.
 
    Feb 2021:  converted to use keywords.h definitions assuming these
    came from the channel collection in MongoDB - Can be changed in
    kewords.h for application outside mspass */
    hang[0] = ts[0].get_double(SEISMICMD_hang);
    hang[1] = ts[1].get_double(SEISMICMD_hang);
    hang[2] = ts[2].get_double(SEISMICMD_hang);
    vang[0] = ts[0].get_double(SEISMICMD_vang);
    vang[1] = ts[1].get_double(SEISMICMD_vang);
    vang[2] = ts[2].get_double(SEISMICMD_vang);
  } catch (MetadataGetError &mde) {
    stringstream ss;
    ss << base_error
       << "missing hang or vang variable in component TimeSeries objects "
          "received"
       << endl;
    ss << "Message posted by Metadata::get_double:  " << mde.what() << endl;
    throw MsPASSError(ss.str(), ErrorSeverity::Invalid);
  }
  /* We couldn't get here if hang and vang were not set on comp 0 so
     we don't test for that condition.  We do need to clear hang and vang
     from result here, however, as both attributes are meaningless
     for a 3C seismogram */
  this->erase(SEISMICMD_hang);
  this->erase(SEISMICMD_vang);
  // These are loaded just for convenience
  t0_component[0] = ts[0].t0();
  t0_component[1] = ts[1].t0();
  t0_component[2] = ts[2].t0();
 
  // Treat the normal case specially and avoid a bunch of work unless
  // it is required
  if ((ts[0].s.size() == ts[1].s.size()) &&
      (ts[1].s.size() == ts[2].s.size()) &&
      (fabs((t0_component[0] - t0_component[1]) / dt()) < 1.0) &&
      (fabs((t0_component[1] - t0_component[2]) / dt()) < 1.0)) {
    /* Older code had this.   No longer needed with logic above that
    calls set_npts.  that method creates and initialized the u dmatrix*/
    // this->u=dmatrix(3,nsamp);
    //  Load data by a simple copy operation
    /* This is a simple loop version
    for(j=0;j<nsamp;++nsamp)
    {
            this->u(0,j)=ts[0].s[j];
            this->u(1,j)=ts[1].s[j];
            this->u(2,j)=ts[2].s[j];
    }
    */
    // This is a vector version that I'll use because it will
    // be faster albeit infinitely more obscure and
    // intrinsically more dangerous
    dcopy(nsamp, &(ts[0].s[0]), 1, u.get_address(0, 0), 3);
    dcopy(nsamp, &(ts[1].s[0]), 1, u.get_address(1, 0), 3);
    dcopy(nsamp, &(ts[2].s[0]), 1, u.get_address(2, 0), 3);
  } else {
    /*Land here if the start time or number of samples
    is irregular.  We cut the output to latest start time to earliest end time*/
    /* WARNING - debugging may be needed for this block. SEISPP versio of this
    used gaps.  Here we cut the output to match an irregularities. */
    double tsmax, temin;
    tsmax = max(t0_component[0], t0_component[1]);
    tsmax = max(tsmax, t0_component[2]);
    temin = min(ts[0].endtime(), ts[1].endtime());
    temin = min(temin, ts[2].endtime());
    nstest = (int)round((temin - tsmax) / mdt);
    if (nstest <= 0)
      throw MsPASSError(
          base_error + "Irregular time windows of components have no overlap",
          ErrorSeverity::Invalid);
    else
      this->CoreSeismogram::set_npts(nstest);
    // Now load the data.  Use the time and sample number methods
    // to simplify process
    double t;
    t = tsmax;
    this->set_t0(t);
    double delta = this->dt();
    for (int ic = 0; ic < 3; ++ic) {
      t = this->t0();
      for (j = 0; j < ts[ic].s.size(); ++j) {
        i = ts[ic].sample_number(t);
        // silently do nothing if outside bounds.  This
        // perhaps should be an error as it shouldn't really
        // happen with the above algorithm, but safety is good
        if ((i >= 0) && (i < nsamp))
          this->u(ic, j) = ts[ic].s[i];
        t += delta;
      }
    }
  }
  /* Finally we need to set the transformation matrix.
   This is a direct application of conversion of routines
  in spherical coordinate procedures.  They are procedural
  routines, not objects so the code is procedural.
  */
  SphericalCoordinate scor;
  double *nu;
  // convert all the hang values to spherical coordinate phi
  // (angle from postive east) from input assumed in degrees
  // azimuth from north.  At the same time convert vang to radians.
  for (i = 0; i < 3; ++i) {
    hang[i] = mspass::utility::rad(90.0 - hang[i]);
    vang[i] = mspass::utility::rad(vang[i]);
  }
  for (i = 0; i < 3; ++i) {
    scor.phi = hang[i];
    scor.theta = vang[i];
    nu = SphericalToUnitVector(scor);
    for (j = 0; j < 3; ++j)
      tmatrix[i][j] = nu[j];
    delete[] nu;
  }
  components_are_cardinal = this->tmatrix_is_cardinal();
  if (components_are_cardinal)
    components_are_orthogonal = true;
  else
    components_are_orthogonal = false;
  /* Last but not least set the datum live before returning */
  this->set_live();
}

References components_are_cardinal, components_are_orthogonal, mspass::seismic::BasicTimeSeries::dt(), endtime(), mspass::utility::Metadata::erase(), mspass::utility::Metadata::get(), mspass::utility::dmatrix::get_address(), mspass::seismic::BasicTimeSeries::mdt, mspass::seismic::BasicTimeSeries::nsamp, mspass::seismic::SEISMICMD_hang(), mspass::seismic::SEISMICMD_vang(), mspass::seismic::BasicTimeSeries::set_live(), set_npts(), set_t0(), mspass::seismic::BasicTimeSeries::t0(), tmatrix, and u.

CoreSeismogram() [4/5]

mspass::seismic::CoreSeismogram::CoreSeismogram	(	const mspass::utility::Metadata &	md,
		const bool	load_data = true
	)

Construct from Metadata definition that includes data path.

A Metadata object is sufficiently general that it can contain enough information to contruct an object from attributes contained in it. This constuctor uses that approach, with the actual loading of data being an option (on by default). In mspass this is constructor is used to load data with Metadata constructed from MongoDB and then using the path created from two parameters (dir and dfile used as in css3.0 wfdisc) to read data. The API is general but the implementation in mspass is very rigid. It blindly assumes the data being read are binary doubles in the right byte order and ordered in the native order for dmatrix (Fortran order). i.e. the constuctor does a raw fread of ns*3 doubles into the internal array used in the dmatrix implementation.

A second element of the Metadata that is special for MsPASS is the handling of the transformation matrix by this constructor. In MsPASS the transformation matrix is stored as a python object in MongoDB. This constructor aims to fetch that entity with the key 'tmatrix'. To be more robust and simpler to use with data not loaded from mongodb we default tmatrix to assume the data are in standard coordinates. That is, if the key tmatrix is not defined in Metadata passed as arg0, the constructor assumes it should set the transformation matrix to an identity. Use set_transformation_matrix if that assumption is wrong for your data.

Parameters

md	is the Metadata used for the construction.
load_data	if true (default) a file name is constructed from dir+"/"+dfile, the file is openned, fseek is called to foff, data are read with fread, and the file is closed. If false a dmatrix for u is still created of size 3xns, but the matrix is only initialized to all zeros.

Exceptions

Will	throw a MsPASSError if required metadata are missing.

    : Metadata(md) {
  string dfile, dir;
  long foff;
  FILE *fp;
  double *inbuffer;
 
  components_are_orthogonal = true;
  mlive = false;
  try {
    /* Names used are from mspass defintions as of Jan 2020.
    We don't need to call the set methods for these attributes as they
    would add the overhead of setting delta, startime, and npts to the
    same value passed. */
    this->mdt = this->get_double(SEISMICMD_dt);
    this->mt0 = this->get_double(SEISMICMD_t0);
    if (this->is_defined(SEISMICMD_time_standard)) {
      if (this->get_string(SEISMICMD_time_standard) == "UTC")
        this->set_tref(TimeReferenceType::UTC);
      else {
        this->set_tref(TimeReferenceType::Relative);
        /* For now we can't post an error because this is CoreSeismogram
        so elog is not defined.   For now let this error be silent as it
        is harmless */
        /*
        this->elog.log_error("CoreSeismogram Metadata constructor",
          SEISMICMD_time_standard+" attribute is not defined - set to Relative",
          ErrorSeverity::Complaint);
          */
      }
    }
    if (this->time_is_relative()) {
      /* It is not an error if a t0 shift is not defined and we are
      in relative time. That is the norm for active source data. */
      if (this->is_defined(SEISMICMD_t0_shift)) {
        double t0shift = this->get_double(SEISMICMD_t0_shift);
        this->force_t0_shift(t0shift);
      }
    }
    /* This section is done specially to handle interaction with MongoDB.
    We store tmatrix there as a python object so we use a get_any to fetch
    it.   Oct 22, 2021 added a bug fix to handle tmatrix not defined.
    We will take a null (undefined) tmatrix stored in the database to
    imply the data are cardinal and orthogonal (i.e. stardard geographic
    coordinates in e,n,z order.)*/
    if (this->is_defined("tmatrix")) {
      this->set_transformation_matrix(
          boost::any_cast<py::object>(this->get_any("tmatrix")));
      components_are_cardinal = this->tmatrix_is_cardinal();
      if (components_are_cardinal)
        components_are_orthogonal = true;
      else
        components_are_orthogonal = false; // May be wrong but cost is tiny
    } else {
      /* this might not be needed but best be explicit*/
      for (auto i = 0; i < 3; ++i)
        for (auto j = 0; j < 3; ++j)
          this->tmatrix[i][j] = 0.0;
      for (auto i = 0; i < 3; ++i)
        this->tmatrix[i][i] = 1.0;
      components_are_orthogonal = true;
      components_are_cardinal = true;
    }
    /* We have to handle nsamp specially in the case when load_data
    is false.  To be consistent with TimeSeries we use a feature that
    if the Metadata container does not define npts we default it.
    In this case that means the default constructor for u and set
    nsamp to 0 (via set_npts). */
    if (md.is_defined(SEISMICMD_npts)) {
      long int ns = md.get_long(SEISMICMD_npts);
      this->set_npts(ns); /* note this is assumed to initialize u*/
    } else {
      this->set_npts(0);
    }
    /* Note previous code had an else clause to to with the
    following conditional.  It used to zero the u matrix.
    The call to set_npts above will always do that so that would
    have been redundant and was removed June 2022*/
    if (load_data) {
      dir = this->get_string(SEISMICMD_dir);
      dfile = this->get_string(SEISMICMD_dfile);
      foff = this->get_long(SEISMICMD_foff);
      string fname = dir + "/" + dfile;
      if ((fp = fopen(fname.c_str(), "r")) == NULL)
        throw(MsPASSError(string("Open failure for file ") + fname,
                          ErrorSeverity::Invalid));
      if (foff > 0)
        fseek(fp, foff, SEEK_SET);
      /* The older seispp code allowed byte swapping here.   For
      efficiency we don't support that here and assume can do a
      raw fread from the file and get valid data.  If support for
      other types is needed this will need to be extended.  Here
      we just point fread at the internal u array. */
      inbuffer = this->u.get_address(0, 0);
      unsigned int nt = 3 * this->nsamp;
      if (fread((void *)(inbuffer), sizeof(double), nt, fp) != nt) {
        fclose(fp);
        throw(MsPASSError(
            string("CoreSeismogram constructor:  fread error on file ") + fname,
            ErrorSeverity::Invalid));
      }
      fclose(fp);
      mlive = true;
    }
  } catch (MsPASSError &mpe) {
    throw(mpe);
  } catch (boost::bad_any_cast &be) {
    throw(MsPASSError(
        string("CoreSeismogram constructor:  tmatrix type is not recognized"),
        ErrorSeverity::Invalid));
  } catch (...) {
    throw;
  };
}

References components_are_cardinal, components_are_orthogonal, mspass::seismic::BasicTimeSeries::force_t0_shift(), mspass::utility::Metadata::get(), mspass::utility::dmatrix::get_address(), mspass::utility::Metadata::get_any(), mspass::utility::Metadata::get_double(), mspass::utility::Metadata::get_long(), mspass::utility::Metadata::get_string(), mspass::utility::Metadata::is_defined(), mspass::utility::Metadata::md, mspass::seismic::BasicTimeSeries::mdt, mspass::seismic::BasicTimeSeries::mlive, mspass::seismic::BasicTimeSeries::mt0, mspass::seismic::BasicTimeSeries::nsamp, mspass::seismic::Relative, mspass::seismic::SEISMICMD_dfile(), mspass::seismic::SEISMICMD_dir(), mspass::seismic::SEISMICMD_dt(), mspass::seismic::SEISMICMD_foff(), mspass::seismic::SEISMICMD_npts(), mspass::seismic::SEISMICMD_t0(), mspass::seismic::SEISMICMD_t0_shift(), mspass::seismic::SEISMICMD_time_standard(), set_npts(), set_transformation_matrix(), mspass::seismic::BasicTimeSeries::set_tref(), mspass::seismic::BasicTimeSeries::t0shift, mspass::seismic::BasicTimeSeries::time_is_relative(), tmatrix, u, and mspass::seismic::UTC.

CoreSeismogram() [5/5]

mspass::seismic::CoreSeismogram::CoreSeismogram

(

const CoreSeismogram &

t3c

)

Standard copy constructor.

    : BasicTimeSeries(dynamic_cast<const BasicTimeSeries &>(t3c)),
      Metadata(dynamic_cast<const Metadata &>(t3c)), u(t3c.u) {
  int i, j;
  components_are_orthogonal = t3c.components_are_orthogonal;
  components_are_cardinal = t3c.components_are_cardinal;
  for (i = 0; i < 3; ++i)
    for (j = 0; j < 3; ++j)
      tmatrix[i][j] = t3c.tmatrix[i][j];
}

References components_are_cardinal, components_are_orthogonal, mspass::utility::Metadata::get(), and tmatrix.

~CoreSeismogram()

virtual mspass::seismic::CoreSeismogram::~CoreSeismogram ( )

inlinevirtual

Standard destructor.

{};

Member Function Documentation

cardinal()

bool mspass::seismic::CoreSeismogram::cardinal ( ) const

inline

Returns true of components are cardinal.

{ return components_are_cardinal; };

References components_are_cardinal.

endtime()

double mspass::seismic::CoreSeismogram::endtime ( ) const

inlinenoexcept

Returns the end time (time associated with last data sample) of this data object.

                                  {
    return (mt0 + mdt * static_cast<double>(u.columns() - 1));
  };

References mspass::utility::dmatrix::columns(), mspass::seismic::BasicTimeSeries::mdt, mspass::seismic::BasicTimeSeries::mt0, and u.

free_surface_transformation()

void mspass::seismic::CoreSeismogram::free_surface_transformation	(	const mspass::seismic::SlownessVector	u,
		const double	vp0,
		const double	vs0
	)

Computes and applies the Kennett [1991] free surface transformation matrix.

Kennett [1991] gives the form for a free surface transformation operator that reduces to a nonorthogonal transformation matrix when the wavefield is not evanescent. On output x1 will be transverse, x2 will be SV (radial), and x3 will be longitudinal.

Parameters

u	slowness vector off the incident wavefield
vp0	Surface P wave velocity
vs0	Surface S wave velocity.

                                                            {
  if ((u.size()[1] <= 0) || dead())
    return; // do nothing in these situations
  double a02, b02, pslow, p2;
  double qa, qb, vpz, vpr, vsr, vsz;
  pslow = uvec.mag();
  // silently do nothing if magnitude of the slowness vector is 0
  // (vertical incidence)
  if (pslow < DBL_EPSILON)
    return;
  // Can't handle evanescent waves with this operator
  double vapparent = 1.0 / pslow;
  if (vapparent < a0 || vapparent < b0) {
    stringstream ss;
    ss << "CoreSeismogram::free_surface_transformation method:  illegal input"
       << endl
       << "Apparent velocity defined by input slowness vector=" << vapparent
       << endl
       << "Smaller than specified surface P velocity=" << a0
       << " or S velocity=" << b0 << endl
       << "That implies evanescent waves that violate the assumption of this "
          "operator"
       << endl;
    throw MsPASSError(ss.str(), ErrorSeverity::Invalid);
  }
 
  // First the horizonal rotation
  SphericalCoordinate scor;
  // rotation angle is - azimuth to put x2 (north in standard coord)
  // in radial direction
  scor.phi = atan2(uvec.uy, uvec.ux);
  scor.theta = 0.0;
  scor.radius = 1.0;
  // after this transformation x1=transverse horizontal
  // x2=radial horizonal, and x3 is still vertical
  this->rotate(scor);
 
  a02 = a0 * a0;
  b02 = b0 * b0;
  p2 = pslow * pslow;
  qa = sqrt((1.0 / a02) - p2);
  qb = sqrt((1.0 / b02) - p2);
  vpz = -(1.0 - 2.0 * b02 * p2) / (2.0 * a0 * qa);
  vpr = pslow * b02 / a0;
  vsr = (1.0 - 2.0 * b02 * p2) / (2.0 * b0 * qb);
  vsz = pslow * b0;
  /* Now construct the transformation matrix
   This is different from Bostock's original code
   in sign and order.  Also note this transformation
       is not scaled to have a unit matrix norm so amplitudes
       after the transformation are distorted.  rotate_to_standard,
       however, should still restore original data within roundoff
       error if called on the result. */
  double fstran[3][3];
  fstran[0][0] = 0.5;
  fstran[0][1] = 0.0;
  fstran[0][2] = 0.0;
  fstran[1][0] = 0.0;
  fstran[1][1] = vsr;
  fstran[1][2] = vpr;
  fstran[2][0] = 0.0;
  fstran[2][1] = -vsz;
  fstran[2][2] = -vpz;
  this->transform(fstran);
 
  components_are_cardinal = false;
  components_are_orthogonal = false;
}

References components_are_cardinal, components_are_orthogonal, mspass::seismic::BasicTimeSeries::dead(), mspass::utility::Metadata::get(), mspass::utility::SphericalCoordinate::phi, rotate(), mspass::utility::dmatrix::size(), transform(), and u.

get_transformation_matrix()

mspass::utility::dmatrix mspass::seismic::CoreSeismogram::get_transformation_matrix ( ) const

inline

Return current transformation matrix.

The transformation matrix is maintained internally in this object. Transformations like rotations and the transform method can change make this matrix not an identity matrix. It should always be an identity matrix when the coordinates are cardinal (i.e. ENZ).

Returns

3x3 transformation matrix.

                                                         {
    mspass::utility::dmatrix result(3, 3);
    for (int i = 0; i < 3; ++i)
      for (int j = 0; j < 3; ++j)
        result(i, j) = tmatrix[i][j];
    return result;
  };

References tmatrix.

operator*=()

CoreSeismogram & mspass::seismic::CoreSeismogram::operator*=

(

const double

scale

)

Multiply data by a scalar.

                                                             {
  /* do nothing to empty data or data marked dead*/
  if ((this->npts() == 0) || (this->dead()))
    return (*this);
  /* We can use this dscal blas function because dmatrix puts all the data
  in a continguous block. Beware if there is am implementation change for
  the matrix data*/
  double *ptr = this->u.get_address(0, 0);
  dscal(3 * this->npts(), scale, ptr, 1);
  return (*this);
}

References mspass::seismic::BasicTimeSeries::dead(), mspass::utility::dmatrix::get_address(), mspass::seismic::BasicTimeSeries::npts(), and u.

operator+()

const CoreSeismogram mspass::seismic::CoreSeismogram::operator+

(

const CoreSeismogram &

other

)

const

Addition operator.

This operator is implemented in a standard way utilizing operator+=. For data with irregular start and end times that has an important consequence; the operator is not communative. i.e given x an y z=x+y will not yield the same result as z=y+x.

                                                           {
  CoreSeismogram result(*this);
  result += other;
  return result;
}

References mspass::utility::Metadata::get().

operator+=()

CoreSeismogram & mspass::seismic::CoreSeismogram::operator+=

(

const CoreSeismogram &

d

)

Summation operator.

Summing data from signals of irregular length requires handling potential mismatches in size and overlap. This behaves the way a += operator should logically behave in that situation. That is, because the lhs is where the sum is being accumulated, the size is always controlled by the left hand side of the operator. Any portions of the right hand side that are outside the t0 to endtime() of the left hand side are silently discarded. If the start time of the right hand side is greater than t0 or the endtime is less than endtime of the lhs there will be discontinuties in the sum there the ends of the rhs are inside the range of the lhs.

Parameters

d	is other signal to add to this.

Exceptions

MsPASSError

if lhs and rhs have incompatible time standards, sample intervals, or start-time grids.

                                                                  {
  // Silently do nothing if d or lhs is marked dead
  if (d.dead() || (this->dead()))
    return (*this);
  const auto overlap =
      detail::arithmetic_overlap(*this, d, "CoreSeismogram::operator+=");
  for (std::size_t i = 0; i < overlap.count; ++i) {
    const std::size_t lhs_index = overlap.lhs_begin + i;
    const std::size_t rhs_index = overlap.rhs_begin + i;
    this->u(0, lhs_index) += d.u(0, rhs_index);
    this->u(1, lhs_index) += d.u(1, rhs_index);
    this->u(2, lhs_index) += d.u(2, rhs_index);
  }
  return (*this);
}

References mspass::seismic::BasicTimeSeries::dead(), mspass::utility::Metadata::get(), and u.

operator-()

const CoreSeismogram mspass::seismic::CoreSeismogram::operator-

(

const CoreSeismogram &

other

)

const

Subtraction operator.

This operator is implemented in a standard way utilizing operator-=. For data with irregular start and end times that has an important consequence; the operator is not communative. i.e given x an y z=x-y will not yield the same result as z=-(y-x).

                                                           {
  CoreSeismogram result(*this);
  result -= other;
  return result;
}

References mspass::utility::Metadata::get().

operator-=()

CoreSeismogram & mspass::seismic::CoreSeismogram::operator-=

(

const CoreSeismogram &

d

)

Subtraction operator.

Differencing data from signals of irregular length requires handling potential mismatches in size and overlap. This behaves the way a -= operator should logically behave in that situation. That is, because the lhs is where the sum is being accumulated, the size is always controlled by the left hand side of the operator. Any portions of the right hand side that are outside the t0 to endtime() of the left hand side are silently discarded. If the start time of the right hand side is greater than t0 or the endtime is less than endtime of the lhs there will be discontinuties in the sum there the ends of the rhs are inside the range of the lhs.

Parameters

d	is other signal to subract from this.

Exceptions

MsPASSError

if lhs and rhs have incompatible time standards, sample intervals, or start-time grids.

                                                                  {
  // Silently do nothing if d is marked dead
  if (d.dead())
    return (*this);
  const auto overlap =
      detail::arithmetic_overlap(*this, d, "CoreSeismogram::operator-=");
  for (std::size_t i = 0; i < overlap.count; ++i) {
    const std::size_t lhs_index = overlap.lhs_begin + i;
    const std::size_t rhs_index = overlap.rhs_begin + i;
    this->u(0, lhs_index) -= d.u(0, rhs_index);
    this->u(1, lhs_index) -= d.u(1, rhs_index);
    this->u(2, lhs_index) -= d.u(2, rhs_index);
  }
  return (*this);
}

References mspass::seismic::BasicTimeSeries::dead(), mspass::utility::Metadata::get(), and u.

operator=()

CoreSeismogram & mspass::seismic::CoreSeismogram::operator=

(

const CoreSeismogram &

seisin

)

Standard assignment operator.

                                                                      {
  if (this != &seisin) {
    this->BasicTimeSeries::operator=(seisin);
    this->Metadata::operator=(seisin);
    components_are_orthogonal = seisin.components_are_orthogonal;
    components_are_cardinal = seisin.components_are_cardinal;
    for (int i = 0; i < 3; ++i) {
      for (int j = 0; j < 3; ++j) {
        tmatrix[i][j] = seisin.tmatrix[i][j];
      }
    }
    u = seisin.u;
  }
  return (*this);
}

References components_are_cardinal, components_are_orthogonal, mspass::utility::Metadata::get(), mspass::seismic::BasicTimeSeries::operator=(), mspass::utility::Metadata::operator=(), tmatrix, and u.

operator[]() [1/2]

vector< double > mspass::seismic::CoreSeismogram::operator[]

(

const double

time

)

const

Overloaded version of operator[] for time.

Sometimes it is useful to ask for data at a specified time without worrying about the time conversion. This simplifies that process. It is still subject to an exception if the the time requested is outside the data range.

Parameters

time	is the time of the requested sample

Returns

3 vector of data samples at requested time

Exceptions

MsPASSError

will be thrown if the time is outside the data range.

                                                              {
  try {
    vector<double> result;
    int i = this->sample_number(t);
    /* could test for a negative i and i too large but we assume
    the dmatrix container will throw an exception if it resolves that way*/
    for (int k = 0; k < 3; ++k)
      result.push_back(this->u(k, i));
    return result;
  } catch (...) {
    throw;
  };
}

References mspass::utility::Metadata::get(), and mspass::seismic::BasicTimeSeries::sample_number().

operator[]() [2/2]

vector< double > mspass::seismic::CoreSeismogram::operator[]

(

const int

sample

)

const

Extract a sample from data vector.

A sample in this context means a three-vector at a requested sample index. Range checking is implicit because of the internal use of the dmatrix to store the samples of data. This operator is an alternative to extracting samples through indexing of the internal dmatrix u that holds the data.

Parameters

sample

is the sample number requested (must be in range or an exception will be thrown)

Exceptions

MsPASSError

if the requested sample is outside the range of the data. Note this includes an implicit "outside" defined when the contents are marked dead. Note the code does this by catching an error thrown by dmatrix in this situation, printing the error message from the dmatrix object, and then throwing a new SeisppError with a shorter message.

Returns

std::vector containing a 3 vector of the samples at requested sample number

                                                           {
  try {
    vector<double> result;
    result.reserve(3);
    for (int k = 0; k < 3; ++k)
      result.push_back(this->u(k, i));
    return result;
  } catch (...) {
    throw;
  };
}

References mspass::utility::Metadata::get().

orthogonal()

bool mspass::seismic::CoreSeismogram::orthogonal ( ) const

inline

Return true if the components are orthogonal.

{ return components_are_orthogonal; };

References components_are_orthogonal.

rotate() [1/3]

void mspass::seismic::CoreSeismogram::rotate

(

const double

nu[3]

)

Rotate data using a P wave type coordinate definition.

In seismology the longitudinal motion direction of a P wave defines a direction in space. This method rotates the data into a coordinate system defined by a direction passed through the argument. The data are rotated such that x1 becomes the transverse component, x2 becomes radial, and x3 becomes longitudinal. In the special case for a vector pointing in the x3 direction the data are not altered.

This method effectively turns nu into a SphericalCoordinate object and calles the related rotate method that has a SphericalCoordinate object as an argument. The potential confusion of orientation is not as extreme here. After the transformation x3prime will point in the direction of nu, x2 will be in the x3-x3prime plane (rotation by theta) and orthogonal to x3prime, and x1 will be horizontal and perpendicular to x2prime and x3prime.

A VERY IMPORTANT thing to recognize about this tranformation is it will always yield a result relative to cardinal coordinates. i.e. if the data had been previously rotated or were not originally in ENZ form they will be first transformed to ENZ before actually performing this transformation. Use the transform or horizontal rotation method to

Parameters

nu	defines direction of x3 direction (longitudinal) as a unit vector with three components.

                                              {
  if ((u.size()[1] <= 0) || this->dead())
    return; // do nothing in these situations
  SphericalCoordinate xsc = UnitVectorToSpherical(nu);
  this->rotate(xsc);
}

References mspass::utility::Metadata::get(), rotate(), mspass::utility::dmatrix::size(), and u.

rotate() [2/3]

void mspass::seismic::CoreSeismogram::rotate

(

const double

phi

)

Rotate horizontals by a simple angle in degrees.

A common transformation in 3C processing is a rotation of the
horizontal components by an angle.  This leaves the vertical
(assumed here x3) unaltered.   This routine rotates the horizontals
by angle phi using with positive phi counterclockwise as in
polar coordinates and the azimuth angle of spherical coordinates.

Note this transformation is cummulative.  i.e. this transformation
is cumulative.  The internal transformation matrix will be updated.
This is a useful feature for things like incremental horizontal
rotation in rotational angle grid searches.

\param phi rotation angle around x3 axis in counterclockwise
  direction (in radians).

                                      {
  if ((u.size()[1] <= 0) || dead())
    return; // do nothing in these situations
  int i, j, k;
  double a, b;
  a = cos(phi);
  b = sin(phi);
  double tmnew[3][3];
  tmnew[0][0] = a;
  tmnew[1][0] = -b;
  tmnew[2][0] = 0.0;
  tmnew[0][1] = b;
  tmnew[1][1] = a;
  tmnew[2][1] = 0.0;
  tmnew[0][2] = 0.0;
  tmnew[1][2] = 0.0;
  tmnew[2][2] = 1.0;
 
  /* Now multiply the data by this transformation matrix.
   Note trick in this i only goes to 2 because 3 component
   is an identity.*/
  double *work[2];
  for (i = 0; i < 2; ++i)
    work[i] = new double[nsamp];
  for (i = 0; i < 2; ++i) {
    dcopy(nsamp, u.get_address(0, 0), 3, work[i], 1);
    dscal(nsamp, tmnew[i][0], work[i], 1);
    daxpy(nsamp, tmnew[i][1], u.get_address(1, 0), 3, work[i], 1);
  }
  for (i = 0; i < 2; ++i)
    dcopy(nsamp, work[i], 1, u.get_address(i, 0), 3);
  double tm_tmp[3][3];
  double prod;
  for (i = 0; i < 3; ++i)
    for (j = 0; j < 3; ++j) {
      for (prod = 0.0, k = 0; k < 3; ++k)
        prod += tmnew[i][k] * tmatrix[k][j];
      tm_tmp[i][j] = prod;
    }
  for (i = 0; i < 3; ++i)
    for (j = 0; j < 3; ++j)
      tmatrix[i][j] = tm_tmp[i][j];
  components_are_cardinal = false;
  for (i = 0; i < 2; ++i)
    delete[] work[i];
}

References components_are_cardinal, mspass::seismic::BasicTimeSeries::dead(), mspass::utility::Metadata::get(), mspass::utility::dmatrix::get_address(), mspass::seismic::BasicTimeSeries::nsamp, mspass::utility::dmatrix::size(), tmatrix, and u.

rotate() [3/3]

void mspass::seismic::CoreSeismogram::rotate

(

mspass::utility::SphericalCoordinate &

sc

)

Rotate data using a P wave type coordinate definition.

In seismology the longitudinal motion direction of a P wave defines a direction in space. This method rotates the data into a coordinate system defined by a direction passed through the argument. The data are rotated such that x1 becomes the transverse component, x2 becomes radial, and x3 becomes longitudinal. In the special case for a vector pointing in the x3 direction the data are not altered. The transformation matrix is effectively the matrix product of two coordinate rotations: (1) rotation around x3 by angle phi and (2) rotation around x1 by theta.

The sense of this transformation is confusing because of a difference in convention between spherical coordinates and standard earth coordinates. In particular, orientation on the earth uses a convention with x2 being the x2 axis and bearings are relative to that with a standard azimuth measured clockwise from north. Spherical coordinate angle phi (used here) is measured counterclockwise relative to the x1 axis, which is east in standard earth coordinates. This transformation is computed using a phi angle. To use this then to compute a transformation to standard ray coordinates with x2 pointing in the direction of wavefront advance, phi should be set to pi/2-azimuth which gives the phi angle needed to rotate x2 to radial. This is extremely confusing because in spherical coordinates it would be more intuitive to rotate x1 to radial, but this is NOT the convention used here. In general to use this feature the best way to avoid this confusion is to use the PMHalfSpaceModel procedure to compute a SphericalCoordinate object consistent with given propagation direction defined by a slowness vector. Alternatively, use the free_surface_transformation method defined below.

A VERY IMPORTANT thing to recognize about this tranformation is it will always yield a result relative to cardinal coordinates. i.e. if the data had been previously rotated or were not originally in ENZ form they will be first transformed to ENZ before actually performing this transformation. Use the transform or horizontal rotation method to create cummulative transformations.

Parameters

sc	defines final x3 direction (longitudinal) in a spherical coordinate structure.

                                                    {
  if ((u.size()[1] <= 0) || dead())
    return; // do nothing in these situations
 
  // Earlier version had a reset of the nsamp variable here - we need to trust
  // that is correct here for efficiency.  We the new API it would be hard
  // to have that happen. without a serious blunder
  int i;
  double theta, phi; /* corrected angles after dealing with signs */
  double a, b, c, d;
 
  //
  // Undo any previous transformations
  //
  this->rotate_to_standard();
  if (xsc.theta == M_PI) {
    // This will be left handed
    tmatrix[2][2] = -1.0;
    dscal(nsamp, -1.0, u.get_address(2, 0), 3);
    components_are_cardinal = false;
    components_are_orthogonal = true;
    return;
  }
 
  if (xsc.theta < 0.0) {
    theta = -(xsc.theta);
    phi = xsc.phi + M_PI;
    if (phi > M_PI)
      phi -= (2.0 * M_PI);
  } else if (xsc.theta > M_PI) {
    theta = xsc.theta - M_PI;
    phi = xsc.phi + M_PI;
    if (phi > M_PI)
      phi -= (2.0 * M_PI);
  } else {
    theta = xsc.theta;
    phi = xsc.phi;
  }
  /* Am using a formula here for azimuth with is pi/2 - phi*/
  double azimuth = M_PI_2 - phi;
  a = cos(azimuth);
  b = sin(azimuth);
  c = cos(theta);
  d = sin(theta);
 
  tmatrix[0][0] = a;
  tmatrix[1][0] = b * c;
  tmatrix[2][0] = b * d;
  tmatrix[0][1] = -b;
  tmatrix[1][1] = a * c;
  tmatrix[2][1] = a * d;
  tmatrix[0][2] = 0.0;
  tmatrix[1][2] = -d;
  tmatrix[2][2] = c;
 
  /* Now multiply the data by this transformation matrix.  */
  double *work[3];
  for (i = 0; i < 3; ++i)
    work[i] = new double[nsamp];
  for (i = 0; i < 3; ++i) {
    dcopy(nsamp, u.get_address(0, 0), 3, work[i], 1);
    dscal(nsamp, tmatrix[i][0], work[i], 1);
    daxpy(nsamp, tmatrix[i][1], u.get_address(1, 0), 3, work[i], 1);
    daxpy(nsamp, tmatrix[i][2], u.get_address(2, 0), 3, work[i], 1);
  }
  for (i = 0; i < 3; ++i)
    dcopy(nsamp, work[i], 1, u.get_address(i, 0), 3);
  components_are_cardinal = false;
  for (i = 0; i < 3; ++i)
    delete[] work[i];
}

References components_are_cardinal, components_are_orthogonal, mspass::seismic::BasicTimeSeries::dead(), mspass::utility::Metadata::get(), mspass::utility::dmatrix::get_address(), mspass::seismic::BasicTimeSeries::nsamp, rotate_to_standard(), mspass::utility::dmatrix::size(), mspass::utility::SphericalCoordinate::theta, tmatrix, and u.

rotate_to_standard()

void mspass::seismic::CoreSeismogram::rotate_to_standard

(

)

Apply inverse transformation matrix to return data to cardinal direction components.

It is frequently necessary to make certain a set of three component data are oriented to the standard reference frame (EW, NS, Vertical). This function does this. For efficiency it checks the components_are_cardinal variable and does nothing if it is set true. Otherwise, it applies the inverse transformation and then sets this variable true. Note even if the current transformation matrix is not orthogonal it will be put back into cardinal coordinates.

Exceptions

SeisppError

thrown if the an inversion of the transformation matrix is required and that matrix is singular. This can happen if the transformation matrix is incorrectly defined or the actual data are coplanar.

                                        {
  if ((u.size()[1] <= 0) || this->dead())
    return; // do nothing in these situations
  double *work[3];
  int i, j;
  if (components_are_cardinal)
    return;
  /* We assume nsamp is the number of samples = number of columns in u - we
  don't check here for efficiency */
  for (j = 0; j < 3; ++j)
    work[j] = new double[nsamp];
  if (components_are_orthogonal) {
    //
    // Use a daxpy algorithm.  tmatrix stores the
    // forward transformation used to get current
    // Use the transpose to get back
    //
    for (i = 0; i < 3; ++i) {
      // x has a stride of 3 because we store in fortran order in x
      dcopy(nsamp, u.get_address(0, 0), 3, work[i], 1);
      dscal(nsamp, tmatrix[0][i], work[i], 1);
      daxpy(nsamp, tmatrix[1][i], u.get_address(1, 0), 3, work[i], 1);
      daxpy(nsamp, tmatrix[2][i], u.get_address(2, 0), 3, work[i], 1);
    }
    for (i = 0; i < 3; ++i)
      dcopy(nsamp, work[i], 1, u.get_address(i, 0), 3);
  } else {
    //
    // Enter here only when the transformation matrix is
    // not orthogonal.  We have to construct a fortran
    // order matrix a to use LINPACK routine in sunperf/perf
    // This could be done with the matrix template library
    // but the overhead ain't worth it
    //
    double a[9];
    int ipivot[3];
    int info;
    a[0] = tmatrix[0][0];
    a[1] = tmatrix[1][0];
    a[2] = tmatrix[2][0];
    a[3] = tmatrix[0][1];
    a[4] = tmatrix[1][1];
    a[5] = tmatrix[2][1];
    a[6] = tmatrix[0][2];
    a[7] = tmatrix[1][2];
    a[8] = tmatrix[2][2];
    // LAPACK routine with FORTRAN interface using pass by reference and
    // pointers
    int three(3);
    dgetrf(three, three, a, three, ipivot, info);
    if (info != 0) {
      for (i = 0; i < 3; ++i)
        delete[] work[i];
      throw(MsPASSError(string("rotate_to_standard:  LU factorization of "
                               "transformation matrix failed"),
                        ErrorSeverity::Invalid));
    }
    // inversion routine after factorization from lapack FORT$RAN interface
    double awork[10]; // Larger than required but safety value small cost
    int ldwork(10);
    dgetri(three, a, three, ipivot, awork, ldwork, info);
    // This is the openblas version
    // info=LAPACKE_dgetri(LAPACK_COL_MAJOR,3,a,3,ipivot);
    if (info != 0) {
      for (i = 0; i < 3; ++i)
        delete[] work[i];
      throw(MsPASSError(string("rotate_to_standard:  LU factorization "
                               "inversion of transformation matrix failed"),
                        ErrorSeverity::Invalid));
    }
 
    tmatrix[0][0] = a[0];
    tmatrix[1][0] = a[1];
    tmatrix[2][0] = a[2];
    tmatrix[0][1] = a[3];
    tmatrix[1][1] = a[4];
    tmatrix[2][1] = a[5];
    tmatrix[0][2] = a[6];
    tmatrix[1][2] = a[7];
    tmatrix[2][2] = a[8];
    /* The inverse is now in tmatrix so we reverse the
       rows and columms from above loop */
 
    for (i = 0; i < 3; ++i) {
      dcopy(nsamp, u.get_address(0, 0), 3, work[i], 1);
      dscal(nsamp, tmatrix[i][0], work[i], 1);
      daxpy(nsamp, tmatrix[i][1], u.get_address(1, 0), 3, work[i], 1);
      daxpy(nsamp, tmatrix[i][2], u.get_address(2, 0), 3, work[i], 1);
    }
    for (i = 0; i < 3; ++i)
      dcopy(nsamp, work[i], 1, u.get_address(i, 0), 3);
    components_are_orthogonal = true;
  }
  //
  // Have to set the transformation matrix to an identity now
  //
  for (i = 0; i < 3; ++i)
    for (j = 0; j < 3; ++j)
      if (i == j)
        tmatrix[i][i] = 1.0;
      else
        tmatrix[i][j] = 0.0;
 
  components_are_cardinal = true;
  for (i = 0; i < 3; ++i)
    delete[] work[i];
}

References components_are_cardinal, components_are_orthogonal, mspass::utility::Metadata::get(), mspass::utility::dmatrix::get_address(), mspass::seismic::BasicTimeSeries::nsamp, mspass::utility::dmatrix::size(), tmatrix, and u.

set_dt()

void mspass::seismic::CoreSeismogram::set_dt ( const double  sample_interval )

virtual

Set the sample interval.

This method is complicated by the need to sync the changed value with Metadata. That is further complicated by the need to support aliases for the keys used to defined dt in Metadata. That is handled by first setting the internal dt value and then going through a fixed list of valid alias keys for dt. Any that exist are changed. If none were previously defined the unique name (see documentation) is added to Metadata.

Parameters

sample_interval

is the new data sample interval to be used.

Reimplemented from mspass::seismic::BasicTimeSeries.

                                                        {
  this->BasicTimeSeries::set_dt(sample_interval);
  /* This is the unique name defined in the mspass schema - we always set it. */
  this->put(SEISMICMD_dt, sample_interval);
  /* these are hard coded aliases for sample_interval */
  std::set<string> aliases;
  std::set<string>::iterator aptr;
  /* Note these aren't set in keywords - aliases are flexible and this
  can allow another way to make these attribute names more flexible. */
  aliases.insert("dt");
  for (aptr = aliases.begin(); aptr != aliases.end(); ++aptr) {
    if (this->is_defined(*aptr)) {
      this->put(*aptr, sample_interval);
    }
  }
}

References mspass::utility::Metadata::get(), mspass::utility::Metadata::is_defined(), mspass::utility::Metadata::put(), mspass::seismic::SEISMICMD_dt(), and mspass::seismic::BasicTimeSeries::set_dt().

set_npts()

void mspass::seismic::CoreSeismogram::set_npts ( const size_t  npts )

virtual

Set the number of samples attribute for data.

This method is complicated by the need to sync the changed value with Metadata. That is further complicated by the need to support aliases for the keys used to defined npts in Metadata. That is handled by first setting the internal npts value (actually ns) and then going through a fixed list of valid alias keys for npts. Any that exist are changed. If none were previously defined the unique name (see documentation) is added to Metadata.

This attribute has an additional complication compared to other setter that are overrides from BasicTimeSeries. That is, the number of points define the data buffer size to hold the sample data. To guarantee the buffer size and the internal remain consistent this method clears any existing content of the dmatrix u and initializes the 3xnpts matrix to 0s. Note this means if one is using this to assemble a data object in pieces you MUST call this method before loading any data or it will be cleared and you will mysteriously find the data are all zeros.

Parameters

npts	is the new number of points to set.

Reimplemented from mspass::seismic::BasicTimeSeries.

                                               {
  this->BasicTimeSeries::set_npts(npts);
  /* This is the unique name - we always set it.  The weird
  cast is necessary to avoid type mismatch with unsigned.
  We use the name defined in keywords.h which we can always assume
  matches the schema for the unique name*/
  this->put(SEISMICMD_npts, (long int)npts);
  /* these are hard coded aliases for sample_interval */
  std::set<string> aliases;
  std::set<string>::iterator aptr;
  aliases.insert("nsamp");
  aliases.insert("wfdisc.nsamp");
  for (aptr = aliases.begin(); aptr != aliases.end(); ++aptr) {
    if (this->is_defined(*aptr)) {
      this->put(*aptr, (long int)npts);
    }
  }
  /* this method has the further complication that npts sets the size of the
  data matrix.   Here we resize the matrix and initialize it to 0s.*/
  if (npts == 0) {
    this->u = dmatrix();
  } else {
    this->u = dmatrix(3, npts);
    this->u.zero();
  }
}

References mspass::utility::Metadata::get(), mspass::utility::Metadata::is_defined(), mspass::seismic::BasicTimeSeries::npts(), mspass::utility::Metadata::put(), mspass::seismic::SEISMICMD_npts(), mspass::seismic::BasicTimeSeries::set_npts(), u, and mspass::utility::dmatrix::zero().

set_t0()

void mspass::seismic::CoreSeismogram::set_t0 ( const double  t0in )

virtual

Set the data start time.

This method is complicated by the need to sync the changed value with Metadata. That is further complicated by the need to support aliases for the keys used to defined npts in Metadata. That is handled by first setting the internal t0 value and then going through a fixed list of valid alias keys for it. Any that exist are changed. If none were previously defined the unique name (see documentation) is added to Metadata.

This is a dangerous method to use on real data as it can mess up the time if not handled correctly. It should be used only when that sharp knife is needed such as in assembling data outside of constructors in a test program.

Parameters

t0in	is the new data sample interval to be used.

Reimplemented from mspass::seismic::BasicTimeSeries.

                                             {
  this->BasicTimeSeries::set_t0(t0in);
  /* This is the unique name - we always set it.  Pulled from keywords.h
  which should match the schema.  aliases are hard coded not defined as
  keywords */
  this->put(SEISMICMD_t0, t0in);
  /* these are hard coded aliases for sample_interval */
  std::set<string> aliases;
  std::set<string>::iterator aptr;
  aliases.insert("t0");
  aliases.insert("time");
  for (aptr = aliases.begin(); aptr != aliases.end(); ++aptr) {
    if (this->is_defined(*aptr)) {
      this->put(*aptr, t0in);
    }
  }
}

References mspass::utility::Metadata::get(), mspass::utility::Metadata::is_defined(), mspass::utility::Metadata::put(), mspass::seismic::SEISMICMD_t0(), and mspass::seismic::BasicTimeSeries::set_t0().

set_transformation_matrix() [1/3]

bool mspass::seismic::CoreSeismogram::set_transformation_matrix

(

const double

a[3][3]

)

Define the transformaton matrix with a C style 3x3 matrix.

Parameters

a	is a C style 3x3 matrix.

Returns

true if the given transformation matrix is an identity meaning components_are_cardinal gets set true. false if the test for an identity matrix fails.

Exceptions

Will	throw a MsPASSError if the input matrix is not 3x3.

                                                                   {
  for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j)
      tmatrix[i][j] = a[i][j];
  py::list tmatrix_l;
  for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j)
      tmatrix_l.append(a[i][j]);
  this->put_object(SEISMICMD_tmatrix, tmatrix_l);
  bool cardinal;
  cardinal = this->tmatrix_is_cardinal();
  if (cardinal) {
    components_are_cardinal = true;
    components_are_orthogonal = true;
  } else {
    components_are_cardinal = false;
    /* Not necessarily true, but small overhead cost*/
    components_are_orthogonal = false;
  }
  return components_are_cardinal;
}

References cardinal(), components_are_cardinal, components_are_orthogonal, mspass::utility::Metadata::get(), mspass::utility::Metadata::put_object(), mspass::seismic::SEISMICMD_tmatrix(), and tmatrix.

set_transformation_matrix() [2/3]

bool mspass::seismic::CoreSeismogram::set_transformation_matrix

(

const mspass::utility::dmatrix &

A

)

Define the transformaton matrix.

Occasionally we need to set the transformation matrix manually. The type example is input with a format where the component directions are embedded. We use a dmatrix as it is more easily wrapped for python than the raw C 2D array which really doesn't translate well between the languages.

Parameters

A	is the 3X3 matrix copied to the internal transformation matrix array.

Returns

true if the given transformation matrix is an identity meaning components_are_cardinal gets set true. false if the test for an identity matrix fails.

Exceptions

Will	throw a MsPASSError if the input matrix is not 3x3.

                                                               {
  for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j)
      tmatrix[i][j] = A(i, j);
  py::list tmatrix_l;
  for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j)
      tmatrix_l.append(A(i, j));
  this->put_object(SEISMICMD_tmatrix, tmatrix_l);
  bool cardinal;
  cardinal = this->tmatrix_is_cardinal();
  if (cardinal) {
    components_are_cardinal = true;
    components_are_orthogonal = true;
  } else {
    components_are_cardinal = false;
    /* Not necessarily true, but small overhead cost*/
    components_are_orthogonal = false;
  }
  return components_are_cardinal;
}

References cardinal(), components_are_cardinal, components_are_orthogonal, mspass::utility::Metadata::get(), mspass::utility::Metadata::put_object(), mspass::seismic::SEISMICMD_tmatrix(), and tmatrix.

set_transformation_matrix() [3/3]

bool mspass::seismic::CoreSeismogram::set_transformation_matrix

(

pybind11::object

a

)

Define the transformaton matrix with a python object.

Parameters

a	is a python object of 9 elements with types of dmatrix, numpy array, or list.

Returns

true if the given transformation matrix is an identity meaning components_are_cardinal gets set true. false if the test for an identity matrix fails.

Exceptions

Will	throw a MsPASSError if the input type or dimension is not recognized.

                                                               {
  if (py::isinstance<py::array>(a)) {
    auto tmatrix_ary = a.cast<
        py::array_t<double, py::array::c_style | py::array::forcecast>>();
    py::buffer_info info = tmatrix_ary.request();
    if ((info.ndim == 2 && info.shape[0] * info.shape[1] == 9) ||
        (info.ndim == 1 && info.shape[0] == 9))
      return this->set_transformation_matrix(
          static_cast<double(*)[3]>(info.ptr));
    else
      throw(MsPASSError(
          string("set_transformation_matrix: tmatrix should be a 3x3 matrix"),
          ErrorSeverity::Invalid));
  } else if (py::isinstance<dmatrix>(a)) {
    auto tmatrix_ary = a.cast<dmatrix>();
    if (tmatrix_ary.rows() != 3 || tmatrix_ary.columns() != 3)
      throw(MsPASSError(
          string("set_transformation_matrix: tmatrix should be a 3x3 matrix"),
          ErrorSeverity::Invalid));
    return this->set_transformation_matrix(tmatrix_ary);
  } else if (py::isinstance<py::list>(a)) {
    dmatrix tmatrix_ary(3, 3);
    double *ptr = tmatrix_ary.get_address(0, 0);
    if (py::len(a) == 9) {
      int i = 0;
      for (auto item : a) {
        try {
          *(ptr + i) = item.cast<double>();
          i++;
        } catch (...) {
          throw(MsPASSError(string("set_transformation_matrix: the elements of "
                                   "tmatrix should be float"),
                            ErrorSeverity::Invalid));
        }
      }
    } else if (py::len(a) == 3) {
      int i = 0;
      for (auto items : a) {
        if (!py::isinstance<py::list>(items))
          throw(MsPASSError(string("set_transformation_matrix: tmatrix should "
                                   "be a 3x3 list of list"),
                            ErrorSeverity::Invalid));
        else if (py::len(items) != 3)
          throw(MsPASSError(string("set_transformation_matrix: tmatrix should "
                                   "be a 3x3 list of list"),
                            ErrorSeverity::Invalid));
        else {
          for (auto item : items) {
            try {
              *(ptr + i) = item.cast<double>();
              i++;
            } catch (...) {
              throw(MsPASSError(string("set_transformation_matrix: the "
                                       "elements of tmatrix should be float"),
                                ErrorSeverity::Invalid));
            }
          }
        }
      }
    } else {
      throw(MsPASSError(string("set_transformation_matrix: tmatrix should be a "
                               "list of 9 floats or a 3x3 list of list"),
                        ErrorSeverity::Invalid));
    }
    return this->set_transformation_matrix(tr(tmatrix_ary));
  } else {
    throw(MsPASSError(
        string("set_transformation_matrix: tmatrix's type is not recognized"),
        ErrorSeverity::Invalid));
  }
}

References mspass::utility::Metadata::get(), and set_transformation_matrix().

sync_npts()

void mspass::seismic::CoreSeismogram::sync_npts

(

)

Sync the number of samples attribute with actual data size.

This method syncs the npts attribute with the actual size of the dmatrix u. It also syncs aliases in the same way as the set_npts method.

                               {
  if (nsamp != this->u.columns()) {
    this->BasicTimeSeries::set_npts(this->u.columns());
    /* This is the unique name - we always set it.  The weird
    cast is necessary to avoid type mismatch with unsigned.
    As above converted to keywords.h const string to make
    this easier to maintain*/
    this->put(SEISMICMD_npts, (long int)nsamp);
    /* these are hard coded aliases for sample_interval */
    std::set<string> aliases;
    std::set<string>::iterator aptr;
    aliases.insert("nsamp");
    aliases.insert("wfdisc.nsamp");
    for (aptr = aliases.begin(); aptr != aliases.end(); ++aptr) {
      if (this->is_defined(*aptr)) {
        this->put(*aptr, (long int)nsamp);
      }
    }
  }
}

References mspass::utility::dmatrix::columns(), mspass::utility::Metadata::get(), mspass::utility::Metadata::is_defined(), mspass::seismic::BasicTimeSeries::nsamp, mspass::utility::Metadata::put(), mspass::seismic::SEISMICMD_npts(), mspass::seismic::BasicTimeSeries::set_npts(), and u.

transform()

void mspass::seismic::CoreSeismogram::transform

(

const double

a[3][3]

)

Applies an arbitrary transformation matrix to the data. i.e. after calling this method the data will have been multiplied by the matrix a and the transformation matrix will be updated. The later allows cascaded transformations to data.

Parameters

a	is a C style 3x3 matrix.

                                                   {
  if ((u.size()[1] <= 0) || dead())
    return; // do nothing in these situations
  /* Older version had this - we need to trust ns is already u.columns().  */
  // size_t ns = u.size()[1];
  size_t i, j, k;
  double *work[3];
  for (i = 0; i < 3; ++i)
    work[i] = new double[nsamp];
  for (i = 0; i < 3; ++i) {
    dcopy(nsamp, u.get_address(0, 0), 3, work[i], 1);
    dscal(nsamp, a[i][0], work[i], 1);
    daxpy(nsamp, a[i][1], u.get_address(1, 0), 3, work[i], 1);
    daxpy(nsamp, a[i][2], u.get_address(2, 0), 3, work[i], 1);
  }
  for (i = 0; i < 3; ++i)
    dcopy(nsamp, work[i], 1, u.get_address(i, 0), 3);
  for (i = 0; i < 3; ++i)
    delete[] work[i];
  /* Hand code this rather than use dmatrix or other library.
     Probably dumb, but this is just a 3x3 system.  This
     is simply a multiply of a*tmatrix with result replacing
     the internal tmatrix */
  double tmnew[3][3];
  double prod;
  for (i = 0; i < 3; ++i)
    for (j = 0; j < 3; ++j) {
      for (prod = 0.0, k = 0; k < 3; ++k)
        prod += a[i][k] * tmatrix[k][j];
      tmnew[i][j] = prod;
    }
  for (i = 0; i < 3; ++i)
    for (j = 0; j < 3; ++j)
      tmatrix[i][j] = tmnew[i][j];
  components_are_cardinal = this->tmatrix_is_cardinal();
  /* Assume this method does not yield cartesian coordinate directions.*/
  if (!components_are_cardinal)
    components_are_orthogonal = false;
}

References components_are_cardinal, components_are_orthogonal, mspass::seismic::BasicTimeSeries::dead(), mspass::utility::Metadata::get(), mspass::utility::dmatrix::get_address(), mspass::seismic::BasicTimeSeries::nsamp, mspass::utility::dmatrix::size(), tmatrix, and u.

Member Data Documentation

components_are_cardinal

bool mspass::seismic::CoreSeismogram::components_are_cardinal

protected

Defines of the contents of the object are in Earth cardinal coordinates.

Cardinal means the cardinal directions at a point on the earth. That is, x1 is positive east, x2 is positive north, and x3 is positive up. Like the components_are_orthogonal variable the purpose of this variable is to simplify common tests for properties of a given data series.

components_are_orthogonal

bool mspass::seismic::CoreSeismogram::components_are_orthogonal

protected

Defines if the contents of this object are components of an orthogonal basis.

Most raw 3c seismic data use orthogonal components, but this is not universal. Furthermore, some transformations (e.g. the free surface transformation operator) define transformations to basis sets that are not orthogonal. Because detecting orthogonality from a transformation is a nontrivial thing (rounding error is the complication) this is made a part of the object to simplify a number of algorithms.

tmatrix

double mspass::seismic::CoreSeismogram::tmatrix[3][3]

protected

Transformation matrix.

This is a 3x3 transformation that defines how the data in this object is produced from cardinal coordinates. That is, if u is the contents of this object the data in cardinal directions can be produced by tmatrix^-1 * u.

u

mspass::utility::dmatrix mspass::seismic::CoreSeismogram::u

Holds the actual data.

Matrix is 3xns. Thus the rows are the component number and columns define time position. Note there is a redundancy in these definitions that must be watched if you manipulate the contents of this matrix. That is, BasicTimeSeries defines ns, but the u matrix has it's own internal size definitions. Currently no tests are done to validate this consistency. All constructors handle this, but again because u is public be very careful in altering u.

---

The documentation for this class was generated from the following files:

• /home/runner/work/mspass/mspass/cxx/include/mspass/seismic/CoreSeismogram.h
• /home/runner/work/mspass/mspass/cxx/src/lib/seismic/CoreSeismogram.cc