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

#include <CoreSeismogram.h>

Inheritance 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. More...

	CoreSeismogram (const CoreSeismogram &)

void	set_dt (const double sample_interval)
	Set the sample interval. More...

void	set_npts (const size_t npts)
	Set the number of samples attribute for data. More...

void	sync_npts ()
	Sync the number of samples attribute with actual data size. More...

void	set_t0 (const double t0in)
	Set the data start time. More...

CoreSeismogram &	operator= (const CoreSeismogram &)

CoreSeismogram &	operator+= (const CoreSeismogram &d)
	Summation operator. More...

const CoreSeismogram	operator+ (const CoreSeismogram &other) const

CoreSeismogram &	operator*= (const double)

CoreSeismogram &	operator-= (const CoreSeismogram &d)
	Subtraction operator. More...

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. More...

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. More...

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. More...

bool	set_transformation_matrix (const double a[3][3])
	Define the transformaton matrix with a C style 3x3 matrix. More...

bool	set_transformation_matrix (pybind11::object a)
	Define the transformaton matrix with a python object. More...

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. More...

double	time (const int i) const

int	sample_number (double t) const

double	endtime () const noexcept

bool	shifted () const

double	time_reference () const

void	force_t0_shift (const double t)
	Force a t0 shift value on data. More...

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

void	set_tref (const TimeReferenceType newtref)
	Force the time standard. More...

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. More...

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. More...

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. More...

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, 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, 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.

◆ 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;
     };
 }

◆ 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, and tmatrix.

◆ ~CoreSeismogram()

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

inlinevirtual

Standard destructor.

301 {};

Member Function Documentation

◆ cardinal()

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

inline

Returns true of components are cardinal.

479 {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::seismic::SlownessVector::mag(), mspass::utility::SphericalCoordinate::phi, mspass::utility::SphericalCoordinate::radius, rotate(), mspass::utility::dmatrix::size(), mspass::utility::SphericalCoordinate::theta, transform(), u, mspass::seismic::SlownessVector::ux, and mspass::seismic::SlownessVector::uy.

◆ 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::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;
 }

◆ 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 can be thrown if lhs and rhs do not have matching time standards.

 {
     int i,iend,jend;
     size_t j,i0,j0;
     // Silently do nothing if d or lhs is marked dead
     if( d.dead() || (this->dead()) ) return(*this);
     // Silently do nothing if d does not overlap with data to contain sum
     if( (d.endtime()<mt0)
             || (d.mt0>(this->endtime())) ) return(*this);
     if(d.tref!=(this->tref))
         throw MsPASSError("CoreSeismogram += operator cannot handle data with inconsistent time base\n",
                           ErrorSeverity::Invalid);
     /* this defines the range of left and right hand sides to be summed */
     i=d.sample_number(this->mt0);
     if(i<0)
     {
       j0=this->sample_number(d.t0());
       i0=0;
     }
     else
     {
       j0=0;
       i0=i;
     }
     iend=d.sample_number(this->endtime());
     jend=this->sample_number(d.endtime());
     if(iend>=(d.npts()))
     {
         iend=d.npts()-1;
     }
     if(jend>=this->npts())
     {
       jend=this->npts()-1;
     }
     for(i=i0,j=j0; i<=iend && j<=jend; ++i,++j)
     {
         this->u(0,j)+=d.u(0,i);
         this->u(1,j)+=d.u(1,i);
         this->u(2,j)+=d.u(2,i);
     }
     return(*this);
 }

References mspass::seismic::BasicTimeSeries::dead(), endtime(), mspass::seismic::BasicTimeSeries::mt0, mspass::seismic::BasicTimeSeries::npts(), mspass::seismic::BasicTimeSeries::sample_number(), mspass::seismic::BasicTimeSeries::t0(), mspass::seismic::BasicTimeSeries::tref, 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;
 }

◆ 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 can be thrown if lhs and rhs do not have matching time standards.

 {
     int i,iend,jend;
     size_t j,i0,j0;
     // Sun's compiler complains about const objects without this.
     CoreSeismogram& d=const_cast<CoreSeismogram&>(data);
     // Silently do nothing if d is marked dead
     if(!d.mlive) return(*this);
     // Silently do nothing if d does not overlap with data to contain sum
     if( (d.endtime()<mt0)
             || (d.mt0>(this->endtime())) ) return(*this);
     if(d.tref!=(this->tref))
         throw MsPASSError("CoreSeismogram += operator cannot handle data with inconsistent time base\n",
                           ErrorSeverity::Invalid);
     /* this defines the range of left and right hand sides to be summed */
     i=d.sample_number(this->mt0);
     if(i<0)
     {
       j0=this->sample_number(d.t0());
       i0=0;
     }
     else
     {
       j0=0;
       i0=i;
     }
     iend=d.sample_number(this->endtime());
     jend=this->sample_number(d.endtime());
     if(iend>=(d.npts()))
     {
         iend=d.npts()-1;
     }
     if(jend>=this->npts())
     {
       jend=this->npts()-1;
     }
     for(i=i0,j=j0; i<=iend && j<=jend; ++i,++j)
     {
         this->u(0,j)-=d.u(0,i);
         this->u(1,j)-=d.u(1,i);
         this->u(2,j)-=d.u(2,i);
     }
     return(*this);
 }

References endtime(), mspass::seismic::BasicTimeSeries::mlive, mspass::seismic::BasicTimeSeries::mt0, mspass::seismic::BasicTimeSeries::npts(), mspass::seismic::BasicTimeSeries::sample_number(), mspass::seismic::BasicTimeSeries::t0(), mspass::seismic::BasicTimeSeries::tref, 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::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::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

Parameters

sample is the integer sample number of data desired.

 {
     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;};
 }

◆ orthogonal()

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

inline

Return true if the components are orthogonal.

481 {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 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::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);
         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, mspass::seismic::BasicTimeSeries::dead(), mspass::seismic::BasicTimeSeries::nsamp, mspass::utility::SphericalCoordinate::phi, 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::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::is_defined(), 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::is_defined(), mspass::seismic::BasicTimeSeries::npts(), 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::is_defined(), 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::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::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.

◆ 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::is_defined(), mspass::seismic::BasicTimeSeries::nsamp, 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::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

Public Member Functions

Public Attributes

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ CoreSeismogram() [1/5]

◆ CoreSeismogram() [2/5]

◆ CoreSeismogram() [3/5]

◆ CoreSeismogram() [4/5]

◆ CoreSeismogram() [5/5]

◆ ~CoreSeismogram()

Member Function Documentation

◆ cardinal()

◆ endtime()

◆ free_surface_transformation()

◆ get_transformation_matrix()

◆ operator*=()

◆ operator+()

◆ operator+=()

◆ operator-()

◆ operator-=()

◆ operator=()

◆ operator[]() [1/2]

◆ operator[]() [2/2]

◆ orthogonal()

◆ rotate() [1/3]

◆ rotate() [2/3]

◆ rotate() [3/3]

◆ rotate_to_standard()

◆ set_dt()

◆ set_npts()

◆ set_t0()

◆ set_transformation_matrix() [1/3]

◆ set_transformation_matrix() [2/3]

◆ set_transformation_matrix() [3/3]

◆ sync_npts()

◆ transform()

Member Data Documentation

◆ components_are_cardinal

◆ components_are_orthogonal

◆ tmatrix

◆ u