TauPSite.h

Go to the documentation of this file.

//- ****************************************************************************
//-
//- Copyright 2009 National Technology & Engineering Solutions of Sandia, LLC
//- (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
//- Government retains certain rights in this software.
//-
//- BSD Open Source License
//- All rights reserved.
//-
//- Redistribution and use in source and binary forms, with or without
//- modification, are permitted provided that the following conditions are met:
//-
//-   1. Redistributions of source code must retain the above copyright notice,
//-      this list of conditions and the following disclaimer.
//-
//-   2. Redistributions in binary form must reproduce the above copyright
//-      notice, this list of conditions and the following disclaimer in the
//-      documentation and/or other materials provided with the distribution.
//-
//-   3. Neither the name of the copyright holder nor the names of its
//-      contributors may be used to endorse or promote products derived from
//-      this software without specific prior written permission.
//-
//- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
//- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
//- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
//- ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
//- LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
//- CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
//- SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
//- INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
//- CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
//- ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//- POSSIBILITY OF SUCH DAMAGE.
//-
//- ****************************************************************************
 
#ifndef TAUPSITE_H
#define TAUPSITE_H
 
// **** _SYSTEM INCLUDES_ ******************************************************
 
#include <vector>
#include <map>
#include <set>
 
using namespace std;
// use standard library objects
 
// **** _LOCAL INCLUDES_ *******************************************************
 
#include "TauPGlobals.h"
 
#include "Brents.h"
#include "TauPSiteFunctionals.h"
 
// **** _BEGIN TAUP NAMESPACE_ *************************************************
 
namespace taup {
 
// **** _FORWARD REFERENCES_ ***************************************************
 
class TPVelocityLayer;
class TravelTimeResult;
 
// **** _CLASS CONSTANTS_ ******************************************************
 
 
// *****************************************************************************
//
//
//
// *****************************************************************************
class TAUP_EXP TauPSite
{
  public:
 
    // **** _PUBLIC LIFECYCLES_ ************************************************
 
    TauPSite();
 
    TauPSite(const string& staname, const string& phase);
 
    TauPSite(const TauPSite& tps);
 
    virtual ~TauPSite();
 
    // **** _PUBLIC OPERATORS_ *************************************************
 
    TauPSite&         operator=(const TauPSite& tps);
 
    // **** _PUBLIC METHODS_ ***************************************************
 
    void              setSiteDepth(double depth);
 
    void              setVelocityModels(const vector<TPVelocityLayer*>& vm);
 
    void              appendVelocityModel(TPVelocityLayer* vm);
 
    void              appendConstVelocityModel(double c, double rt, double rb,
                                               const string& layrnam = "");
 
    void              appendLinearVelocityModel(double a0, double a1,
                                                double rt, double rb,
                                                const string& layrnam = "",
                                                double normradius = 1.0);
 
    void              appendQuadraticVelocityModel(double a0, double a1,
                                                   double a2,
                                                   double rt, double rb,
                                                   const string& layrnam = "",
                                                   double normradius = 1.0);
 
    void              appendCubicVelocityModel(double a0, double a1,
                                               double a2, double a3,
                                               double rt, double rb,
                                               const string& layrnam = "",
                                               double normradius = 1.0);
 
    void              clearVelocityModels();
 
    void              calculateTravelTimes(double srcdist, double srcdepth,
                                           bool evalderivs = true);
 
    TravelTimeResult* getFirstTravelTimeResult();
 
    double            getFirstTravelTime() const;
 
    double            getFirstNonDiffractedTravelTime() const;
 
    TravelTimeResult* getTravelTimeResult(const string& phase,
                                          bool matchPhase);
    //- Returns the travel time result with the input matching phase.
    //- If the matching phase was not found and matchPhase is true a NULL
    //- pointer is returned. If the matching phase was not found and matchPhase
    //- is false then the first arrival travel time result is returned.
 
    double            getTravelTime(const string& phase,
                                    bool matchPhase) const;
    //- Returns the travel time with the input matching phase.
    //- If the matching phase was not found and matchPhase is true a travel
    //- time result of -1 is returned. If the matching phase was not found and
    //- matchPhase is false then the first arrival travel time is returned.
 
    const map<double, TravelTimeResult*>& getAllTravelTimes() const;
 
    //
    bool              integrateDistance(double p, double& d,
                                        bool bottom_pass = false);
 
    bool              integrateDistance(double p, double r, double& d);
 
    bool              integrateDistance(double p, double r1, double r2, double& d);
 
    //
    bool              integrateTime(double p, double& d,
                                    bool bottom_pass = false);
 
    bool              integrateTime(double p, double r, double& d);
 
    bool              integrateTime(double p, double r1, double r2, double& t);
 
    virtual void      writeData(const string& filename) const;
 
    //
    void              writeLayerData(int i, ostream& os, int n,
                                     double f0 = 0.0, double f1 = 1.0);
 
    void              dumpLocalSrcRcvrLayers(ostream& os);
 
    void              dumpLayerInfo(ostream& os);
 
    string            toString() const;
 
    virtual void      toStream(ostream& os, string indent) const;
 
    // **** _PUBLIC PROPERTIES_ ************************************************
 
    static  string    class_name();
 
    virtual string    get_class_name() const;
 
    virtual int       class_size() const;
 
    static  string    commonName();
 
    virtual string    getCommonName() const;
 
    const string&     getName() const;
 
    const string&     getPhase() const;
 
    void              setModelName(const string& modname);
 
    const string&     getModelName() const;
 
    void              setEarthRadius(double r);
 
    double            getEarthRadius() const;
 
    void              setBrentsTolerance(double tol);
 
    double            getBrentsTolerance() const;
 
    const vector<TPVelocityLayer*>& getVelocityModels() const;
 
    virtual int64     get_memory() const;
 
    virtual int64     get_alloc_memory() const;
 
    // **** _PUBLIC INQUIRIES_ *************************************************
 
  private:
 
    // **** _PRIVATE METHODS_ **************************************************
 
    virtual void      writeData(ostream& os) const;
 
    void              clearTTRMap();
 
    TravelTimeResult* saveResultI(double pB, double pI, double dist,
                                  int i, bool isUpper, bool isLower);
 
    TravelTimeResult* saveResult(double p, int i,
                                 bool isintrfcupper, bool isintrfclower,
                                 bool isspltlower, bool isspltupper,
                                 bool isturningzero);
 
    void              findLimits();
 
    void              debugOutZeros(int n, double p0, double p1);
 
    //**************************************************************************
    //**** Begin: Derivative Evaluation Functions ******************************
    //**************************************************************************
 
    // Derivative calculation occurs in the function evaluateDerivatives(...).
    // The initial call is made by calculateTravelTimes(...). The arguments
    // assume that a valid travel time has been evaluated for some layer whose
    // top and bottom ray parameter is given by pT and pB, respectively. The
    // valid travel time is passed in as T00. The travel time was evaluated at
    // dist and depth. The derivatives will be determined numerically
    // by calculating travel times at steps away from dist and depth using the
    // step sizes edist and edepth. The derivatives are returned in the array
    // derivs[] and include derivs[0] = dT/dDist, derivs[1] = dT/dDpth,
    // derivs[2] = d^2T/dDist/dDpth, and derivs[3] = d^2T/dDist^2. The travel
    // time is assumed to correspond to a turning ray unless the diffracted
    // ray parameter (pI) is input as a non-negative value (defaults to -1.0).
    //
    // On entry the function evaluateDerivatives(...) assumes that the
    // TPZeroFunctional object has had its top ray parameter set to pT.
    // It immediately calls evalDerivsPrimary(...) which calls the function
    // layeredTravelTime(...) at the four primary stencil locations. If
    // successful it immediately returns with the derivatives stored in derivs[].
    //
    // If one or more of the stencil locations cannot be evaluated for the input
    // layer then one or two alternate stencils are tried, These include B and A,
    // B and C, E and F, and E and D. If any are successful the function returns
    // with the evaluated derivatives. Each pair is selected based on the point
    // that failed during the primary evaluation. If both alternates also fail
    // then the function reEvaluateDerivatives(...) is called which reduces the
    // step sizes (edist and edepth) by 1/2 and recalls evaluateDerivatives(...)
    // recursively. Each recursive call increments a counter which, if it exceeds
    // 5, will trip an error condition which causes a return without evaluting
    // the derivatives. The derivative function calling sequence is shown in the
    // diagram below.
    //
    //                         initial call
    //                              |
    //  |--> evaluateDerivatives <--|
    //  |    {
    //  |      evalDerivsPrimary ------->|
    //  |      evalDerivsAlternateA ---->|
    //  |      evalDerivsAlternateB ---->|
    //  |      evalDerivsAlternateC ---->|----> layeredTravelTime
    //  |      evalDerivsAlternateD ---->|
    //  |      evalDerivsAlternateE ---->|
    //  |      evalDerivsAlternateF ---->|
    //  |
    //  |<--   reEvaluateDerivatives
    //       }
    //
 
    //
    void              evaluateDerivatives(double pT, double pB, double T00,
                                          double dist, double depth,
                                          double edist, double edepth,
                                          double* derivs, double pI = -1.0,
                                          int derivcnt = 0);
 
    void              reEvaluateDerivatives(double pT, double pB, double T00,
                                            double dist, double depth,
                                            double edist, double edepth,
                                            double* derivs, double pI,
                                            int derivcnt);
 
    //
    int               evalDerivsPrimary(double pT, double pB, double T00,
                                        double dist, double depth,
                                        double edist, double edepth,
                                        double* derivs, double pI);
 
    //
    int               evalDerivsAlternateA(double pT, double pB, double T00,
                                           double dist, double depth,
                                           double edist, double edepth,
                                           double* derivs, double pI);
 
    //
    int               evalDerivsAlternateB(double pT, double pB, double T00,
                                           double dist, double depth,
                                           double edist, double edepth,
                                           double* derivs, double pI);
 
    //
    int               evalDerivsAlternateC(double pT, double pB, double T00,
                                           double dist, double depth,
                                           double edist, double edepth,
                                           double* derivs, double pI);
 
    //
    int               evalDerivsAlternateD(double pT, double pB, double T00,
                                           double dist, double depth,
                                           double edist, double edepth,
                                           double* derivs, double pI);
 
    //
    int               evalDerivsAlternateE(double pT, double pB, double T00,
                                           double dist, double depth,
                                           double edist, double edepth,
                                           double* derivs, double pI);
 
    //
    int               evalDerivsAlternateF(double pT, double pB, double T00,
                                           double dist, double depth,
                                           double edist, double edepth,
                                           double* derivs, double pI);
 
    double            layeredTravelTime(double pT, double pB, double dist,
                                        double depth, double& T, double pI);
 
    //**************************************************************************
    //**** End: Derivative Evaluation Functions ********************************
    //**************************************************************************
 
    // **** _PRIVATE DATA_ *****************************************************
 
    static int        tpsClassCount;
 
    static const double tpsDerivStepSize;
 
    static const double tpsBrentsZeroInTol;
 
    static vector<TravelTimeResult*> tpsReuseTTR;
 
    vector<TPVelocityLayer*> tpsVLayer;
 
    bool              tpsIsVelModlOwned;
 
    string            tpsRcvrName;
 
    string            tpsRcvrPhase;
 
    string            tpsModelName;
 
    double            tpsRcvrRad;
 
    double            tpsRcvrDepth;
 
    double            tpsSrcRad;
 
    double            tpsSrcDepth;
 
    double            tpsDist;
 
    double            tpsLastDist;
 
    double            tpsLastDepth;
 
    TPZeroFunctional  tpsZeroF;
 
    util::Brents<TPZeroFunctional> tpsZeroIn;
 
    map<double, TravelTimeResult*> tpsTTR;
 
}; // TauPSite End Definition
 
// *****************************************************************************
// **** TravelTimeResult Definition ********************************************
// *****************************************************************************
 
// *****************************************************************************
//
//
// *****************************************************************************
class TAUP_EXP TravelTimeResult
{
  public:
 
    TravelTimeResult() {};
 
    string           toString(string indent = "") const;
 
    void             toStream(ostream& os, string indent = "") const;
 
    TPVelocityLayer* ttrVelLayer;
 
    string           ttrPhaseName;
 
    string           ttrPhaseEval;
 
    string           ttrRayType;
 
    int              ttrLayerIndex;             // Xn
 
    bool             ttrIsInterfaceUpper;       // Xni+
 
    bool             ttrIsInterfaceLower;       // Xni-
 
    bool             ttrIsSplitUpper;           // Xnb+
 
    bool             ttrIsSplitLower;           // Xnb- ... else Xnb
 
    double           ttrP;
 
    double           ttrR;
 
    double           ttrV;
 
    double           ttrT;
 
    double           ttrDRay;
 
    double           ttrTRay;
 
    double           ttrDSrc;
 
    double           ttrTSrc;
 
    double           ttrDRcvr;
 
    double           ttrTRcvr;
 
    double           ttrDIntrfc;
 
    double           ttrTIntrfc;
 
    double           ttrDerivs[4];
 
  private:
};
 
// *****************************************************************************
// **** _INLINE FUNCTION IMPLEMENTATIONS_ **************************************
// *****************************************************************************
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
//
// *****************************************************************************
inline string TauPSite::class_name()
{
  return "TauPSite";
}
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
//
// *****************************************************************************
inline string TauPSite::get_class_name() const
{
  return class_name();
}
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
//
// *****************************************************************************
inline int TauPSite::class_size() const
{
  return (int) sizeof(TauPSite);
}
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
//
// *****************************************************************************
inline string TauPSite::commonName()
{
  return ("1-D Radial TauP Travel Time Site Model");
}
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
//
// *****************************************************************************
inline string TauPSite::getCommonName() const
{
  return commonName();
}
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
//
// *****************************************************************************
inline int64 TauPSite::get_memory() const
{
  return (sizeof(TauPSite) + TauPSite::get_alloc_memory());
}
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
//
// *****************************************************************************
inline const string& TauPSite::getName() const
{
  return tpsRcvrName;
}
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
//
// *****************************************************************************
inline const string& TauPSite::getPhase() const
{
  return tpsRcvrPhase;
}
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
//
// *****************************************************************************
inline void TauPSite::setModelName(const string& modname)
{
  tpsModelName = modname;
}
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
//
// *****************************************************************************
inline const string& TauPSite::getModelName() const
{
  return tpsModelName;
}
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
//
// *****************************************************************************
inline const vector<TPVelocityLayer*>& TauPSite::getVelocityModels() const
{
  return tpsVLayer;
}
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
//
// *****************************************************************************
inline void TauPSite::setEarthRadius(double r)
{
  tpsZeroF.setPlanetRadius(r);
}
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
//
// *****************************************************************************
inline double TauPSite::getEarthRadius() const
{
  return tpsZeroF.getPlanetRadius();
}
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
//
// *****************************************************************************
inline void TauPSite::setBrentsTolerance(double tol)
{
  tpsZeroIn.setTolerance(tol);
}
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
//
// *****************************************************************************
inline double TauPSite::getBrentsTolerance() const
{
  return tpsZeroIn.getTolerance();
}
 
// **** _FUNCTION DESCRIPTION_ *************************************************
//
// *****************************************************************************
inline const map<double, TravelTimeResult*>& TauPSite::getAllTravelTimes() const
{
  return tpsTTR;
}
 
} // end namespace taup
 
#endif // TAUPSITE_H