GreatCircle.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 GreatCircle_H
#define GreatCircle_H
 
// **** _SYSTEM INCLUDES_ ******************************************************
 
#include <vector>
#include <map>
#include <string>
#include <iostream>
 
using namespace std;
 
// **** _LOCAL INCLUDES_ *******************************************************
 
#include "SLBMGlobals.h"
#include "Location.h"
#include "LayerProfile.h"
//#include "CrustalProfile.h"
 
// **** _BEGIN SLBM NAMESPACE_ **************************************************
 
namespace slbm {
 
class Grid;
class CrustalProfile;
 
// **** _CLASS DEFINITION_ *****************************************************
//
class SLBM_EXP_IMP GreatCircle
{
 
public:
 
    GreatCircle(
        const int& _phase,
        Grid& _grid,
        const double& latSource,
        const double& lonSource,
        const double& depthSource,
        const double& latReceiver,
        const double& lonReceiver,
        const double& depthReceiver);
 
    virtual ~GreatCircle();
 
    GreatCircle(const GreatCircle &other);
 
    GreatCircle& operator=(const GreatCircle& other);
 
    virtual bool operator == (const GreatCircle& other) const;
    virtual bool operator != (const GreatCircle& other) const { return !(*this == other); } ;
 
    Grid& getGrid() { return grid; };
 
    int getPhase() { return phase; };
 
    const int& getHeadWaveInterface()  { return headWaveInterface; };
 
    string getPhaseString();
 
    CrustalProfile* getSourceProfile() { return source; };
 
    CrustalProfile* getReceiverProfile() { return receiver; };
 
    double getSourceRayParameter() { return sourceRayParameter; };
 
    double getReceiverRayParameter() { return receiverRayParameter; };
 
    LayerProfile* getProfile(const int& i);
 
    int getNProfiles() { return (int)profiles.size(); };
 
    void getLayerProfileLocation(const int& i, Location& loc);
 
    double getDistance() { return distance; };
 
    double getEsaz();
 
    double getSourceDistance() { return xSource; };
 
    double getSourceDistanceX() { return xSource; };
 
    double getSourceDistanceZ() { return zSource; };
 
    double getSourceDistanceS() { return sSource; };
    
    double getReceiverDistance() { return xReceiver; };
 
    double getReceiverDistanceX() { return xReceiver; };
 
    double getReceiverDistanceZ() { return zReceiver; };
 
    double getReceiverDistanceS() { return sReceiver; };
    
    double getHeadwaveDistance() { return distance-xSource-xReceiver; };
 
    double getHeadwaveDistanceKm() { return xHorizontal; };
 
    double getActualPathIncrement() { return actual_path_increment; };
 
    double getTravelTime();
 
    void getTravelTime(double& tTotal, double& tSource, double& tReceiver,
            double& tMantle, double& tGradient);
 
    virtual string toString(const int& verbosity) = 0;
 
    void getData(
        int& phase,
        double& actual_path_increment,
        double sourceDepth[NLAYERS],
        double sourceVelocity[NLAYERS],
        double receiverDepth[NLAYERS],
        double receiverVelocity[NLAYERS],
        int& npoints,
        double headWaveVelocity[],
        double gradient[]
        );
 
    void getData(
        int& phase,
        double& actual_path_increment,
        vector<double>& sourceDepth,
        vector<double>& sourceVelocity,
        vector<double>& receiverDepth,
        vector<double>& receiverVelocity,
        vector<double>& headWaveVelocity,
        vector<double>& gradient
        );
 
    void getNodeInfo(
        int** neighbors,
        double** coefficients,
        const int& maxpoints,
        const int& maxnodes,
        int& npoints,
        int* nnodes
        );
 
    void getNodeInfo( vector<vector<int> >& neighbors, vector<vector<double> >& coefficients );
 
    void getGreatCircleLocation(const double& distance, Location& loc);
 
    void getWeights(vector<int>& nodeids, vector<double>& weights, bool headWaveWeights = true);
 
    void getWeights(int nodeids[], double weights[], int& nweights);
 
    double getFractionActive();
 
    double getRayParameter() { return rayParameter; };
 
    double getTurningRadius() { return turningRadius; };
 
    virtual void getZhaoParameters(double& Vm, double& Gm, double& H, double& C, double& Cm, int& udSign) = 0;
 
    virtual void getPgLgComponents(double& tT,
                                double& tTaup, double& tHeadwave,
                                double& pTaup, double& pHeadwave,
                                double& trTaup, double& trHeadwave)
    { tT = tTaup = tHeadwave = pTaup = pHeadwave = trTaup = trHeadwave = NA_VALUE; };
 
    virtual size_t memSize();
 
    static int getClassCount();
 
 
    static void setDelDistance(const double& del_distance) { DEL_DISTANCE = del_distance; }
 
    static double getDelDistance() { return DEL_DISTANCE; }
 
    static void setDelDepth(const double& del_depth) { DEL_DEPTH = del_depth; }
 
    static double getDelDepth() { return DEL_DEPTH; }
 
    static void setPathIncrement(const double& path_increment) { PATH_INCREMENT = path_increment; }
 
    static double getPathIncrement() { return PATH_INCREMENT; }
 
    void get_dtt_ddist(double& dtt_ddist);
 
    void get_dtt_dlat(double& dtt_dlat);
 
    void get_dtt_dlon(double& dtt_dlon);
 
    void get_dtt_ddepth(double& dtt_ddepth);
 
    //void get_dtt_dlat_fast(double& dtt_dlat);
 
    //void get_dtt_dlon_fast(double& dtt_dlon);
 
    //void get_dsh_ddist(double& dsh_ddist);
 
    //void get_dsh_dlat(double& dsh_dlat);
 
    //void get_dsh_dlon(double& dsh_dlon);
 
    //void get_dsh_ddepth(double& dsh_ddepth);
 
    static double MAX_DISTANCE;
    static double MAX_DEPTH;
 
    void setNAValues();
 
protected:
 
    static int greatCircleClassCount;
 
    Grid& grid;
 
    int phase;
 
    double latSource;
    double lonSource;
    double depthSource;
    double latReceiver;
    double lonReceiver;
    double depthReceiver;
 
    int headWaveInterface;
 
    CrustalProfile* source;
 
    CrustalProfile* receiver;
 
    Location location;
 
    double vtp[3];
 
    vector<LayerProfile*> profiles;
 
    string solutionMethod;
 
    double tTotal;
 
    double actual_path_increment;
 
    double distance;
 
    double esaz;
 
    // \brief The ray parameter, in sec/radian or sec/km.
    double rayParameter;
 
    double xSource;
 
    double zSource;
 
    double sSource;
    
    double tSource;
 
    double xReceiver;
 
    double zReceiver;
 
    double sReceiver;
    
    double tReceiver;
 
    double xHorizontal;
 
    double tHorizontal;
 
    double tGamma;
 
    double sourceRayParameter;
    double receiverRayParameter;
 
    double turningRadius;
 
    int sourceIndex;
 
    int receiverIndex;
 
    double getActualPathIncrement(const int& i);
 
    virtual void computeTravelTime() = 0;
 
    double sqr(const double& x) { return x*x; };
 
    double ttHminus;
 
    double ttHplus;
 
    //double ttZplus;
 
    double ttHZplus;
 
    // Used in calculation of derivatives wrt source position.
    double ttEast;
 
    // Used in calculation of derivatives wrt source position.
    double ttWest;
 
    // Used in calculation of derivatives wrt source position.
    double ttEastZ;
 
    // Used in calculation of derivatives wrt source position.
    double ttNorth;
 
    // Used in calculation of derivatives wrt source position.
    double ttSouth;
 
    // Used in calculation of derivatives wrt source position.
    double ttNorthZ;
 
    double get_ttHminus();
 
    double get_ttHplus();
 
    double get_ttZplus();
 
    double get_ttZminus();
 
    double get_ttHZplus();
 
    double get_ttEast();
 
    double get_ttWest();
 
    double get_ttEastZ();
 
    double get_ttNorth();
 
    double get_ttSouth();
 
    double get_ttNorthZ();
 
    // when computing horizontal slowness, the source is moved distance
    // del_distance (radians) away from the receiver and a new travel time
    // is computed.  The difference in travel time divided by del_distance
    // is the horizontal slowness.  Units are radians.
    static double DEL_DISTANCE;
 
    // when computing derivative of travel time wrt depth, the depth of the
    // source is increased by del_depth and a new travel time computed.
    // The difference in travel time divided by del_depth is the
    // derivative of travel time wrt depth.  Units of DEL_DEPTH is km.
    static double DEL_DEPTH;
 
    static double PATH_INCREMENT;
 
};
 
inline double GreatCircle::getTravelTime()
{
    return tTotal;
}
 
inline void GreatCircle::getTravelTime(double& tT, double& tS, double& tR, double& tM, double& tG)
{
    tT = tTotal;
    tS = tSource;
    tR = tReceiver;
    tM = tHorizontal;
    tG = tGamma;
}
 
inline size_t GreatCircle::memSize()
{
    //size_t n =
    //4 // sizeof(grid)
    //+ sizeof(phase)
    //+ sizeof(source)
    //+ sizeof(receiver)
    //+ sizeof(location)
    //+ sizeof(vtp)+3*sizeof(double)
    //+ sizeof(tTotal)
    //+ sizeof(actual_path_increment)
    //+ sizeof(distance)
    //+ sizeof(esaz)
    //+ sizeof(xSource)
    //+ sizeof(tSource)
    //+ sizeof(xReceiver)
    //+ sizeof(tReceiver)
    //+ sizeof(xHorizontal)
    //+ sizeof(tHorizontal)
    //+ sizeof(sourceIndex)
    //+ sizeof(receiverIndex);
    //for (int i=0; i<(int)profiles.size(); i++)
    //  if (profiles[i]) n += profiles[i]->memSize();
 
    size_t n =
    sizeof(&grid) // Grid&
    + sizeof(phase)
    + sizeof(headWaveInterface)
    + sizeof(source)  // CrustalProfile*
    + sizeof(receiver) // CrustalProfile*
    + sizeof(location) // Location
    + sizeof(vtp) + 3 * sizeof(double) // vtp[3]
    + sizeof(profiles) // vector<LayerProfile*> profiles)
    + sizeof(solutionMethod) // string
    + sizeof(tTotal)
    + sizeof(actual_path_increment)
    + sizeof(distance)
    + sizeof(esaz)
    + sizeof(rayParameter)
    + sizeof(xSource)
    + sizeof(zSource)
    + sizeof(sSource)
    + sizeof(tSource)
    + sizeof(xReceiver)
    + sizeof(zReceiver)
    + sizeof(sReceiver)
    + sizeof(tReceiver)
    + sizeof(xHorizontal)
    + sizeof(tHorizontal)
    + sizeof(tGamma)
    + sizeof(sourceRayParameter)
    + sizeof(receiverRayParameter)
    + sizeof(turningRadius)
    + sizeof(sourceIndex)
    + sizeof(receiverIndex)
    + sizeof(ttHminus)
    + sizeof(ttHplus)
    //+ sizeof(ttZplus)
    + sizeof(ttHZplus)
    + sizeof(ttEast)
    + sizeof(ttWest)
    + sizeof(ttEastZ)
    + sizeof(ttNorth)
    + sizeof(ttSouth)
    + sizeof(ttNorthZ)
    ;
 
    return n;
}
 
inline double GreatCircle::getActualPathIncrement(const int& i)
{
    // this horizontal segment lies completely in between the
    // source and receiver.  return actual_path_increment.
    if (i > sourceIndex && i < receiverIndex)
        return actual_path_increment;
 
    // if the source and receiver have the same index, it means the ray
    // hit the headwave interface in some horizontal segment, and left
    // the headwave interface before exiting the horizontal segment.
    // In this case, return the source-receiver separation minus the
    // horizontal offsets below the source and the receiver.
    if (i == sourceIndex && i == receiverIndex)
        return distance - xSource - xReceiver;
 
    // the ray hit the interface in this segment.  return only the
    // part of the horizontal segment that starts at the source pierce point
    // and ends at the profile just to the right.
    if (i == sourceIndex)
        return (sourceIndex+1)*actual_path_increment - xSource;
 
    // this is the horizontal segment such that the ray left the
    // interface in this segment.  return only the part of the
    // horizontal segment that starts at the profile that is the
    // left side of the horizontal segment and ends at the
    // receiver pierce point.
    if (i == receiverIndex)
        return distance - xReceiver - receiverIndex*actual_path_increment;
 
    return 0.;
}
 
inline string GreatCircle::getPhaseString()
{
    return (phase==Pn ? "Pn" : (phase==Sn ? "Sn" : (phase==Pg ? "Pg" : (phase==Lg ? "Lg"
            : "unknown phase"))));
}
 
inline double GreatCircle::getFractionActive()
{
    int nactive = 0;
    for (int i=0; i<(int)profiles.size(); ++i)
        if (getProfile(i)->isActiveProfile())
            ++nactive;
    return ((double)nactive)/((double)profiles.size());
}
 
inline void GreatCircle::get_dtt_ddist(double& dtt_ddist)
{
    if (getTravelTime() > -1.)
      dtt_ddist = (get_ttHplus() - getTravelTime())/DEL_DISTANCE;
    else
        dtt_ddist = NA_VALUE;
}
 
inline void GreatCircle::get_dtt_dlat(double& dtt_dlat)
{
    if (getTravelTime() > -1.)
        dtt_dlat = (get_ttNorth() - getTravelTime())/DEL_DISTANCE;
    else
        dtt_dlat = NA_VALUE;
}
 
inline void GreatCircle::get_dtt_dlon(double& dtt_dlon)
{
    if (getTravelTime() > -1.)
        dtt_dlon = (get_ttEast() - getTravelTime())/DEL_DISTANCE;
    else
        dtt_dlon = NA_VALUE;
}
 
//inline void GreatCircle::get_dtt_dlat_fast(double& dtt_dlat)
//{
//  dtt_dlat = (get_ttHplus() - getTravelTime())/DEL_DISTANCE*cos(getEsaz()+PI);
//}
//
//inline void GreatCircle::get_dtt_dlon_fast(double& dtt_dlon)
//{
//  dtt_dlon = (get_ttHplus() - getTravelTime())/DEL_DISTANCE*sin(getEsaz()+PI);
//}
//
 
//inline void GreatCircle::get_dsh_ddist(double& dsh_ddist)
//{
//  dsh_ddist = (get_ttHminus() - 2*getTravelTime() + get_ttHplus())
//          /DEL_DISTANCE/DEL_DISTANCE;
//}
//
//inline void GreatCircle::get_dsh_dlat(double& dsh_dlat)
//{
//  dsh_dlat = (get_ttNorth() - 2*getTravelTime() + get_ttSouth())
//          /DEL_DISTANCE/DEL_DISTANCE;
//}
//
//inline void GreatCircle::get_dsh_dlon(double& dsh_dlon)
//{
//  dsh_dlon = (get_ttEast() - 2*getTravelTime() + get_ttWest())
//          /DEL_DISTANCE/DEL_DISTANCE;
//}
//
//inline void GreatCircle::get_dsh_ddepth(double& dsh_ddepth)
//{
//  dsh_ddepth = (get_ttHZplus() - get_ttZplus()- get_ttHplus()
//      + getTravelTime()) /DEL_DEPTH/DEL_DISTANCE;
//}
 
} // end slbm namespace
 
#endif // GreatCircle.h