SlbmInterface.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 SlbmInterface_H
#define SlbmInterface_H
 
// **** _SYSTEM INCLUDES_ ******************************************************
#include <map>
#include <cmath>
 
using namespace std;
 
// **** _LOCAL INCLUDES_ *******************************************************
 
#include "GreatCircle.h"
#include "GreatCircleFactory.h"
#include "GridGeoTess.h"
#include "UncertaintyPIU.h"
#include "UncertaintyPDU.h"
#include "SLBMException.h"
 
// **** _BEGIN SLBM NAMESPACE_ **************************************************
 
namespace slbm {
 
class SLBM_EXP SlbmInterface
{
 
public:
 
    SlbmInterface();
 
    SlbmInterface(const double& earthRadius);
 
    virtual ~SlbmInterface();
 
    bool isEqual(SlbmInterface *other);
 
    bool operator == (SlbmInterface &other) { return isEqual(&other); };
    bool operator != (SlbmInterface &other) { return !(*this == other); };
 
 
    static bool modelsEqual(const string modelPath1, const string modelPath2);
 
    string getVersion() { return SlbmVersion; };
 
    void loadVelocityModel(const string& modelPath);
 
    void saveVelocityModel(const string& modelFileName, const int& format=4);
 
 
    void loadVelocityModelBinary(const string& modelPath)
    { loadVelocityModel(modelPath); };
 
    void specifyOutputDirectory(const string& directoryName);
 
    void saveVelocityModelBinary();
 
    void loadVelocityModelBinary(util::DataBuffer& buffer);
 
    void saveVelocityModelBinary(util::DataBuffer& buffer);
 
 
    int  getBufferSize() const;
 
    void setInterpolatorType(const string& interpolatorType);
 
    string getInterpolatorType();
 
    void createGreatCircle(const string& phase,
                    const double& sourceLat,
                    const double& sourceLon,
                    const double& sourceDepth,
                    const double& receiverLat,
                    const double& receiverLon,
                    const double& receiverDepth);
 
    void createGreatCircle(const int& phase,
                    const double& sourceLat,
                    const double& sourceLon,
                    const double& sourceDepth,
                    const double& receiverLat,
                    const double& receiverLon,
                    const double& receiverDepth);
 
    void clear();
 
    bool isValid() { return valid; };
 
    string getPhase() { return greatCircle->getPhaseString(); };
 
    void getDistance(double& distance) { distance = getDistance(); };
 
    double getDistance();
 
    void getSourceDistance(double& dist);
    double getSourceDistance() { double x; getSourceDistance(x); return x; };
 
    void getReceiverDistance(double& dist);
    double getReceiverDistance() { double x; getReceiverDistance(x); return x; };
 
    void getHeadwaveDistance(double& dist);
    double getHeadwaveDistance() { double x; getHeadwaveDistance(x); return x; };
 
    void getHeadwaveDistanceKm(double& dist);
    double getHeadwaveDistanceKm() { double x; getHeadwaveDistanceKm(x); return x; };
 
    void getTravelTime(double& travelTime);
    double getTravelTime() { double x; getTravelTime(x); return x; };
 
    void getTravelTimeComponents(double& tTotal, double& tSource, double& tReceiver,
            double& tHeadwave, double& tGradient);
 
    void getSlowness(double& slowness);
    double getSlowness() { double x; getSlowness(x); return x; };
 
    void get_dtt_ddist(double& slowness) { getSlowness(slowness); };
    double get_dtt_ddist() { double x; getSlowness(x); return x; };
 
    void get_dtt_dlat(double& dtt_dlat);
    double get_dtt_dlat() { double x; get_dtt_dlat(x); return x; };
 
    void get_dtt_dlon(double& dtt_dlon);
    double get_dtt_dlon() { double x; get_dtt_dlon(x); return x; };
 
    void get_dtt_ddepth(double& dtt_ddepth);
    double get_dtt_ddepth() { double x; get_dtt_ddepth(x); return x; };
 
    //void get_dtt_dlat_fast(double& dtt_dlat);
 
    //void get_dtt_dlon_fast(double& dtt_dlon);
 
//  //! \brief Retrieve the derivative of horizontal slowness wrt to source-receiver distance,
//  //! in seconds/radian^2.
//  //!
//  //! Retrieve the derivative of horizontal slowness wrt to source-receiver distance,
//  //! in seconds/radian^2.
//  void get_dsh_ddist(double& dsh_ddist);
//  double get_dsh_ddist() { double x; get_dsh_ddist(x); return x; };
//
//  //! \brief Retrieve the derivative of horizontal slowness wrt to source latitude,
//  //! in seconds/radian^2.
//  //!
//  //! Retrieve the derivative of horizontal slowness wrt to source latitude,
//  //! in seconds/radian^2.
//  //! @param dsh_dlat the derivative of horizontal slowness wrt to source latitude.
//  void get_dsh_dlat(double& dsh_dlat);
//
//  //! \brief Retrieve the derivative of horizontal slowness wrt to source longitude,
//  //! in seconds/radian^2.
//  //!
//  //! Retrieve the derivative of horizontal slowness wrt to source longitude,
//  //! in seconds/radian^2.
//  //! @param dsh_dlon the derivative of horizontal slowness wrt to source longitude.
//  void get_dsh_dlon(double& dsh_dlon);
//
//  //! \brief Retrieve the derivative of horizontal slowness wrt to source depth,
//  //! in seconds/radian-km.
//  //!
//  //! Retrieve the derivative of horizontal slowness wrt to source depth,
//  //! in seconds/radian-km.
//  //! @param dsh_ddepth the derivative of horizontal slowness wrt to source depth.
//  void get_dsh_ddepth(double& dsh_ddepth);
//  double get_dsh_ddepth() { double x; get_dsh_ddepth(x); return x; }
 
    void getWeights(int nodeId[], double weight[], int& nWeights);
 
    void getActiveNodeWeights(int nodeId[], double weight[], int& nWeights);
 
    void getWeights(vector<int>& nodeId, vector<double>& weight);
 
    void getActiveNodeWeights(vector<int>& nodeId, vector<double>& weight);
 
    void getWeightsSource(int nodeids[], double weights[], int& nWeights);
    void getWeightsSource(vector<int> &nodeids, vector<double> &weights);
 
    void getActiveNodeWeightsSource(int nodeids[], double weights[], int& nWeights);
    void getActiveNodeWeightsSource(vector<int> &nodeids, vector<double> &weights);
 
    void getWeightsReceiver(int nodeids[], double weights[], int& nWeights);
    void getWeightsReceiver(vector<int> &nodeids, vector<double> &weights);
 
    void getActiveNodeWeightsReceiver(int nodeids[], double weights[], int& nWeights);
    void getActiveNodeWeightsReceiver(vector<int> &nodeids, vector<double> &weights);
 
    string toString(const int& verbosity);
 
    void getNGridNodes(int& n);
 
    int getNGridNodes();
 
    void getNHeadWavePoints(int& nHeadWavePoints);
    int getNHeadWavePoints() { int x; getNHeadWavePoints(x); return x; };
 
    void getGridData(
        const int& nodeId,
        double& latitude,
        double& longitude,
        double depth[NLAYERS],
        double pvelocity[NLAYERS],
        double svelocity[NLAYERS],
        double gradient[2]);
 
    void getActiveNodeData(
        const int& nodeId,
        double& latitude,
        double& longitude,
        double depth[NLAYERS],
        double pvelocity[NLAYERS],
        double svelocity[NLAYERS],
        double gradient[2]);
 
    void setGridData(
        const int& nodeId,
        double depths[NLAYERS],
        double pvelocity[NLAYERS],
        double svelocity[NLAYERS],
        double gradient[2]);
 
    void setActiveNodeData(
        const int& nodeId,
        double depths[NLAYERS],
        double pvelocity[NLAYERS],
        double svelocity[NLAYERS],
        double gradient[2]);
 
    void getGreatCircleData(
        string& phase,
        double& actual_path_increment,
        double sourceDepth[NLAYERS],
        double sourceVelocity[NLAYERS],
        double receiverDepth[NLAYERS],
        double receiverVelocity[NLAYERS],
        int& npoints,
        double headWaveVelocity[],
        double gradient[]
        );
    void getGreatCircleData(
        string &phase,
        double &actual_path_increment,
        vector<double> &sourceDepth,
        vector<double> &sourceVelocity,
        vector<double> &receiverDepth,
        vector<double> &receiverVelocity,
        int &npoints,
        vector<double> &headWaveVelocity,
        vector<double> &gradient
        );
 
    void getGreatCircleLocations(double lat[], double lon[], double depth[], int& npoints);
    void getGreatCircleLocations(vector<double> &lat, vector<double> &lon, vector<double> &depth);
 
    void getGreatCircleNodeInfo(
        int** neighbors,
        double** coefficients,
        const int& maxpoints,
        const int& maxnodes,
        int& npoints,
        int* nnodes
        );
    void getGreatCircleNodeInfo(vector<vector<int> >& neighbors, vector<vector<double> >& coefficients);
 
    void getInterpolatedPoint(
        const double& lat,
        const double& lon,
        int* nodeIds,
        double* coefficients,
        int& nWeights,
        double depth[NLAYERS],
        double pvelocity[NLAYERS],
        double svelocity[NLAYERS],
        double& pgradient,
        double& sgradient
        );
 
    void getInterpolatedPoint(
        const double &lat,
        const double &lon,
        vector<int> &nodeId,
        vector<double> &coefficients,
        vector<double> &depth,
        vector<double> &pvelocity,
        vector<double> &svelocity,
        double &pgradient,
        double &sgradient
        );
 
    void getInterpolatedTransect(
        double lat[],
        double lon[],
        const int& nLatLon,
        int** neighbors,
        double** coefficients,
        int* nNeighbors,
        double depth[][NLAYERS],
        double pvelocity[][NLAYERS],
        double svelocity[][NLAYERS],
        double pgradient[NLAYERS],
        double sgradient[NLAYERS],
        int& nInvalid
        );
 
    void getInterpolatedTransect(
        vector<double> lat,
        vector<double> lon,
        vector<vector<int> > &nodeId,
        vector<vector<double> > &coefficients,
        vector<vector<double> > &depth,
        vector<vector<double> > &pvelocity,
        vector<vector<double> > &svelocity,
        vector<double> &pgradient,
        vector<double> &sgradient,
        int& nInvalid
        );
 
    void initializeActiveNodes(const double& latmin, const double& lonmin,
        const double& latmax, const double& lonmax);
 
    void initializeActiveNodes(const string& polygonFileName);
 
    void initializeActiveNodes(GeoTessPolygon* polygon);
 
    void initializeActiveNodes(double* lat, double* lon, const int& npoints, const bool& inDegrees=true);
 
    void initializeActiveNodes(const vector<double> lat, const vector<double> lon, const bool& inDegrees=true);
 
    void initializeActiveNodes(vector<double*>& unitVectors);
 
    int getNActiveNodes();
 
    void clearActiveNodes();
 
    int getGridNodeId(int activeNodeId);
 
    int getActiveNodeId(int gridNodeId);
 
    void getNodeHitCount(const int& nodeId, int& hitCount);
    int getNodeHitCount(const int &nodeId) { int x; getNodeHitCount(nodeId, x); return x; };
 
 
    void clearNodeHitCount();
 
    void getNodeNeighbors(const int& nid, int neighbors[], int& nNeighbors);
 
    void getActiveNodeNeighbors(const int& nid, int neighbors[], int& nNeighbors);
 
    void getNodeNeighbors(const int& nid, vector<int>& neighbors);
 
    void getActiveNodeNeighbors(const int& nid, vector<int>& neighbors);
 
    void getNodeNeighborInfo(const int& nid, int neighbors[], double distance[],
        double azimuth[], int& nNeighbors);
    void getNodeNeighborInfo(const int nid,
        vector<int> &neighbors, vector<double> &distance, vector<double> &azimuth);
 
    void getActiveNodeNeighborInfo(const int& nid, int neighbors[], double distance[],
        double azimuth[], int& nNeighbors);
 
    void getActiveNodeNeighborInfo(const int& nid,
            vector<int>& neighbors, vector<double>& distance, vector<double>& azimuth);
 
    void getNodeSeparation(const int& node1, const int& node2, double& distance);
 
    void getNodeAzimuth(const int& node1, const int& node2, double& azimuth);
 
    GreatCircle* getGreatCircleObject();
 
    Grid* getGridObject() { return grid; };
 
    GeoTessModelSLBM* getModelObject() { return grid->getModel(); }
 
    void getTravelTimeUncertainty( const int& phase, const double& distance, double& uncert );
 
    void getTravelTimeUncertainty(double& travelTimeUncertainty, bool calcRandomError = false);
 
    void getTravelTimeUncertainty1D(double& travelTimeUncertainty);
 
    void getSlownessUncertainty( const int& phase, const double& distance, double& uncert );
 
    void getSlownessUncertainty(double& slownessUncertainty);
 
    string getUncertaintyTable(const int& phase, const int& attribute);
 
    string getUncertaintyFileFormat(const int& phase, const int& attribute);
 
    void getZhaoParameters(double& Vm, double& Gm, double& H, double& C, double& Cm, int& udSign);
 
    void getPgLgComponents(double& tTotal, double& tTaup, double& tHeadwave,
        double& pTaup, double& pHeadwave, double& trTaup, double& trHeadwave);
 
    void static setCHMax(const double& chMax);
 
    void static getCHMax(double& chMax);
    double static getCHMax() { double x; getCHMax(x); return x; }
 
    void getAverageMantleVelocity(const int& type, double& velocity);
    double getAverageMantleVelocity(const int &type);
 
    void setAverageMantleVelocity(const int& type, const double& velocity);
 
    void getTessId(string& tessId);
    string getTessId() { string s; getTessId(s); return s; };
 
    void getFractionActive(double& fractionActive);
    double getFractionActive() { double x; getFractionActive(x); return x; }
 
    void static setMaxDistance(const double& maxDistance);
 
    void static getMaxDistance(double& maxDistance);
    double static getMaxDistance() { double x; getMaxDistance(x); return x; }
 
    void static setMaxDepth(const double& maxDepth);
 
    void static getMaxDepth(double& maxDepth);
    double static getMaxDepth() { double x; getMaxDepth(x); return x; }
 
    string getClassCount();
 
    const string& getModelPath() const;
 
    void getDistAz(const double& aLat, const double& aLon,
        const double& bLat, const double& bLon,
        double& distance, double& azimuth, const double& naValue);
 
    void movePoint(const double& aLat, const double& aLon,
        const double& distance, const double& azimuth,
        double& bLat, double& bLon);
 
    void getPiercePointSource(double& lat, double& lon, double& depth);
 
    void getPiercePointReceiver(double& lat, double& lon, double& depth);
 
    void getGreatCirclePoints(
        const double& aLat,
        const double& aLon,
        const double& bLat,
        const double& bLon,
        const int& npoints,
        double latitude[],
        double longitude[]);
    void getGreatCirclePoints(
        const double& aLat,
        const double& aLon,
        const double& bLat,
        const double& bLon,
        const int& npoints,
        vector<double> &latitude,
        vector<double> &longitude);
 
    void getGreatCirclePointsOnCenters(
        const double& aLat,
        const double& aLon,
        const double& bLat,
        const double& bLon,
        const int& npoints,
        double latitude[],
        double longitude[]);
    void getGreatCirclePointsOnCenters(
        const double& aLat,
        const double& aLon,
        const double& bLat,
        const double& bLon,
        const int& npoints,
        vector<double> &latitude,
        vector<double> &longitude);
 
    void setDelDistance(const double& del_distance);
 
    void getDelDistance(double& del_distance);
    double getDelDistance() { double x; getDelDistance(x); return x; };
 
    void setDelDepth(const double& del_depth);
 
    void getDelDepth(double& del_depth);
    double getDelDepth() { double x; getDelDepth(x); return x; };
 
    void getRayParameter(double& ray_parameter);
    double getRayParameter() { double x; getRayParameter(x); return x; };
 
    void getTurningRadius(double& turning_radius);
    double getTurningRadius() { double x; getTurningRadius(x); return x; };
 
    void setPathIncrement(const double& pathIncrement);
 
    void getPathIncrement(double& pathIncrement);
 
    double getPathIncrement();
 
    string getModelString() { return grid->toString(); }
 
    double getSrcLat() { return srcLat; }
 
    double getSrcLon() { return srcLon; }
 
    double getSrcDep() { return srcDep; }
 
    double getRcvLat() { return rcvLat; }
 
    double getRcvLon() { return rcvLon; }
 
    double getRcvDep() { return rcvDep; }
 
 
protected:
 
    Grid* grid;
 
    GreatCircle* greatCircle;
 
    static double CH_MAX;
 
    bool valid;
 
    void clearGreatCircles();
 
private:
 
    // copies of the source and receiver locations that were
    // passed to createGreatCircle the last time it was called.
    // Used to compute distance in instances when createGreateCircle
    // fails. Units are radians.
    double srcLat, srcLon, srcDep, rcvLat, rcvLon, rcvDep;
};
 
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
//  INLINE FUNCTIONS
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
inline void SlbmInterface::createGreatCircle(
                    const string& p,
                    const double& sourceLat,
                    const double& sourceLon,
                    const double& sourceDepth,
                    const double& receiverLat,
                    const double& receiverLon,
                    const double& receiverDepth)
{
    int iphase = (p=="Pn" ? Pn : (p=="Sn" ? Sn : (p=="Pg" ? Pg : (p=="Lg" ? Lg : -1))));
    if (iphase == -1)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
         os << endl << "ERROR in SlbmInterface::createGreatCircle" << endl
            << p << " is not a recognized phase.  Must be one of Pn, Sn, Pg, Lg" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),112);
    }
 
    createGreatCircle(
    iphase,
    sourceLat,
    sourceLon,
    sourceDepth,
    receiverLat,
    receiverLon,
    receiverDepth
    );
}
 
inline GreatCircle* SlbmInterface::getGreatCircleObject()
{
    return greatCircle;
}
 
inline void SlbmInterface::get_dtt_dlat(double& value)
{
    if (!isValid())
    {
        value = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::get_dtt_dlat" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    greatCircle->get_dtt_dlat(value);
}
 
inline void SlbmInterface::get_dtt_dlon(double& value)
{
    if (!isValid())
    {
        value = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::get_dtt_dlon" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    greatCircle->get_dtt_dlon(value);
}
 
inline void SlbmInterface::get_dtt_ddepth(double& value)
{
    if (!isValid())
    {
        value = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::get_dtt_ddepth" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    greatCircle->get_dtt_ddepth(value);
}
 
//inline void SlbmInterface::get_dtt_dlat_fast(double& value)
//{
//  if (!isValid())
//  {
//      value = NA_VALUE;
//      ostringstream os;
//      os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
//        os << endl << "ERROR in SlbmInterface::get_dtt_dlat" << endl
//          << "GreatCircle is invalid." << endl
//            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
//      throw SLBMException(os.str(),113);
//  }
//  greatCircle->get_dtt_dlat_fast(value);
//}
//
//inline void SlbmInterface::get_dtt_dlon_fast(double& value)
//{
//  if (!isValid())
//  {
//      value = NA_VALUE;
//      ostringstream os;
//      os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
//        os << endl << "ERROR in SlbmInterface::get_dtt_dlon" << endl
//          << "GreatCircle is invalid." << endl
//            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
//      throw SLBMException(os.str(),113);
//  }
//  greatCircle->get_dtt_dlon_fast(value);
//}
//
//inline void SlbmInterface::get_dsh_ddist(double& value)
//{
//  if (!isValid())
//  {
//      value = NA_VALUE;
//      ostringstream os;
//      os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
//        os << endl << "ERROR in SlbmInterface::get_dsh_ddist" << endl
//          << "GreatCircle is invalid." << endl
//            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
//      throw SLBMException(os.str(),113);
//  }
//  greatCircle->get_dsh_ddist(value);
//}
//
//inline void SlbmInterface::get_dsh_dlat(double& value)
//{
//  if (!isValid())
//  {
//      value = NA_VALUE;
//      ostringstream os;
//      os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
//        os << endl << "ERROR in SlbmInterface::get_dsh_dlat" << endl
//          << "GreatCircle is invalid." << endl
//            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
//      throw SLBMException(os.str(),113);
//  }
//  greatCircle->get_dsh_dlat(value);
//}
//
//inline void SlbmInterface::get_dsh_dlon(double& value)
//{
//  if (!isValid())
//  {
//      value = NA_VALUE;
//      ostringstream os;
//      os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
//        os << endl << "ERROR in SlbmInterface::get_dsh_dlon" << endl
//          << "GreatCircle is invalid." << endl
//            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
//      throw SLBMException(os.str(),113);
//  }
//  greatCircle->get_dsh_dlon(value);
//}
//
//inline void SlbmInterface::get_dsh_ddepth(double& value)
//{
//  if (!isValid())
//  {
//      value = NA_VALUE;
//      ostringstream os;
//      os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
//        os << endl << "ERROR in SlbmInterface::get_dsh_ddepth" << endl
//          << "GreatCircle is invalid." << endl
//            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
//      throw SLBMException(os.str(),113);
//  }
//  greatCircle->get_dsh_ddepth(value);
//}
 
inline void SlbmInterface::getSlowness(double& value)
{
    if (!isValid())
    {
        value = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getSlowness" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    greatCircle->get_dtt_ddist(value);
}
 
inline double SlbmInterface::getDistance()
{
    if (greatCircle == NULL)
    {
        Location s(srcLat, srcLon, 0.);
        Location r(rcvLat, rcvLon, 0.);
        return s.distance(r);
    }
    else
        return greatCircle->getDistance();
}
 
inline void SlbmInterface::getSourceDistance(double& distance)
{
    if (!isValid())
    {
        distance = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getSourceDistance" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    distance = greatCircle->getSourceDistance();
}
 
inline void SlbmInterface::getReceiverDistance(double& distance)
{
    if (!isValid())
    {
        distance = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getReceiverDistance" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    distance = greatCircle->getReceiverDistance();
}
 
inline void SlbmInterface::getHeadwaveDistance(double& distance)
{
    if (!isValid())
    {
        distance = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getHeadwaveDistance" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    distance = greatCircle->getHeadwaveDistance();
}
 
inline void SlbmInterface::getHeadwaveDistanceKm(double& distance)
{
    if (!isValid())
    {
        distance = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getHeadwaveDistanceKm" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    distance = greatCircle->getHeadwaveDistanceKm();
}
 
inline void SlbmInterface::getTravelTime(double& tTotal)
{
    if (!isValid())
    {
        tTotal = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getTravelTime" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    tTotal=greatCircle->getTravelTime();
}
 
inline void SlbmInterface::getTravelTimeComponents(
            double& tTotal, double& tSource, double& tReceiver,
            double& tHeadwave, double& tGradient)
{
    if (!isValid())
    {
        tTotal = NA_VALUE;
        tSource = NA_VALUE;
        tReceiver = NA_VALUE;
        tHeadwave = NA_VALUE;
        tGradient = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getTravelTimeComponents" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    greatCircle->getTravelTime(tTotal, tSource, tReceiver,
        tHeadwave, tGradient);
}
 
inline void SlbmInterface::getWeights(int nodeId[], double weight[], int& nWeights)
{
    if (!isValid())
    {
        nWeights=-1;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getWeights" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    greatCircle->getWeights(nodeId, weight, nWeights);
}
 
inline void SlbmInterface::getActiveNodeWeights(int nodeId[], double weight[], int& nWeights)
{
    getWeights(nodeId, weight, nWeights);
    for (int i=0; i<nWeights; ++i)
        nodeId[i] = grid->getActiveNodeId(nodeId[i]);
}
 
inline void SlbmInterface::getWeights(vector<int>& nodeId, vector<double>& weight)
{
    if (!isValid())
    {
        nodeId.clear();
        weight.clear();
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getWeights" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    greatCircle->getWeights(nodeId, weight);
}
 
inline void SlbmInterface::getActiveNodeWeights(vector<int>& nodeId, vector<double>& weight)
{
    getWeights(nodeId, weight);
    for (int i=0; i<(int)nodeId.size(); ++i)
        nodeId[i] = grid->getActiveNodeId(nodeId[i]);
}
 
inline void SlbmInterface::getWeightsSource(int nodeids[], double weights[], int& nWeights)
{
    if (!isValid())
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getWeightsSource" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    greatCircle->getSourceProfile()->getWeights(nodeids, weights, nWeights);
}
 
inline void SlbmInterface::getWeightsSource(vector<int> &nodeids, vector<double> &weights)
{
    // get the number of weights
    int nWeights = greatCircle->getSourceProfile()->getNCoefficients();
 
    // resize arrays to proper size
    nodeids.resize(nWeights);
    weights.resize(nWeights);
 
    // call getWeightsSource (v.data() sends a pointer similar to a c++ array)
    getWeightsSource(nodeids.data(), weights.data(), nWeights);
}
 
inline void SlbmInterface::getActiveNodeWeightsSource(int nodeids[], double weights[], int& nWeights)
{
    getWeightsSource(nodeids, weights, nWeights);
    for (int i=0; i<nWeights; ++i)
        nodeids[i] = grid->getActiveNodeId(nodeids[i]);
}
 
inline void SlbmInterface::getActiveNodeWeightsSource(vector<int> &nodeids, vector<double> &weights)
{
    // get the number of weights
    int nWeights = greatCircle->getSourceProfile()->getNCoefficients();
 
    // resize arrays to proper size
    nodeids.resize(nWeights);
    weights.resize(nWeights);
 
    // call getActiveNodeWeightsSource (v.data() sends a pointer similar to a c++ array)
    getActiveNodeWeightsSource(nodeids.data(), weights.data(), nWeights);
}
 
inline void SlbmInterface::getWeightsReceiver(int nodeids[], double weights[], int& nWeights)
{
    if (!isValid())
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getWeightsReceiver" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    greatCircle->getReceiverProfile()->getWeights(nodeids, weights, nWeights);
}
 
inline void SlbmInterface::getWeightsReceiver(vector<int> &nodeids, vector<double> &weights)
{
    // get the number of weights
    int nWeights = greatCircle->getReceiverProfile()->getNCoefficients();
 
    // resize arrays to proper size
    nodeids.resize(nWeights);
    weights.resize(nWeights);
 
    // call getActiveNodeWeightsSource (v.data() sends a pointer similar to a c++ array)
    getWeightsReceiver(nodeids.data(), weights.data(), nWeights);
}
 
inline void SlbmInterface::getActiveNodeWeightsReceiver(int nodeids[], double weights[], int& nWeights)
{
    getWeightsReceiver(nodeids, weights, nWeights);
    for (int i=0; i<nWeights; ++i)
        nodeids[i] = grid->getActiveNodeId(nodeids[i]);
}
 
inline void SlbmInterface::getActiveNodeWeightsReceiver(vector<int> &nodeids, vector<double> &weights)
{
    // get the number of weights
    int nWeights = greatCircle->getReceiverProfile()->getNCoefficients();
 
    // resize arrays to proper size
    nodeids.resize(nWeights);
    weights.resize(nWeights);
 
    // call getActiveNodeWeightsSource (v.data() sends a pointer similar to a c++ array)
    getActiveNodeWeightsReceiver(nodeids.data(), weights.data(), nWeights);
}
 
inline string SlbmInterface::toString(const int& verbosity)
{
    if (!isValid())
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::toString" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
 
    ostringstream os;
    if (verbosity > 0)
    {
        os << endl
            << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl
            << "Great Circle " << endl << endl
            << greatCircle->toString(verbosity);
 
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint)
            << setprecision(4);
        double slowness;
        getSlowness(slowness);
        os << endl << "Horizontal slowness = "
            << slowness << " sec/radian"
            << endl << endl;
    }
 
    return os.str();
}
 
inline void SlbmInterface::getGreatCircleData(
        string& phase,
        double& actual_path_increment,
        double sourceDepth[NLAYERS],
        double sourceVelocity[NLAYERS],
        double receiverDepth[NLAYERS],
        double receiverVelocity[NLAYERS],
        int& npoints,
        double headWaveVelocity[],
        double gradient[]
        )
{
    if (!isValid())
    {
        phase = "";
        actual_path_increment = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getGreatCircleData" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
 
    int p;
    greatCircle->getData(p, actual_path_increment,
        sourceDepth, sourceVelocity, receiverDepth, receiverVelocity,
        npoints, headWaveVelocity, gradient);
    phase = (p==Pn ? "Pn" : (p==Sn ? "Sn" : (p==Pg ? "Pg" : (p==Lg ? "Lg"
        : "unknown phase"))));
}
 
inline void SlbmInterface::getGreatCircleData(
        string &phase,
        double &actual_path_increment,
        vector<double> &sourceDepth,
        vector<double> &sourceVelocity,
        vector<double> &receiverDepth,
        vector<double> &receiverVelocity,
        int &npoints,
        vector<double> &headWaveVelocity,
        vector<double> &gradient
        )
{
    if (!isValid())
    {
        phase = "";
        actual_path_increment = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getGreatCircleData" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
 
    // figure out how many points there are
    getNHeadWavePoints(npoints);
 
    // make some room in memory
    sourceDepth.resize(NLAYERS);
    sourceVelocity.resize(NLAYERS);
    receiverDepth.resize(NLAYERS);
    receiverVelocity.resize(NLAYERS);
    headWaveVelocity.resize(npoints);
    gradient.resize(npoints);
 
    int p;
    greatCircle->getData(p, actual_path_increment,
        sourceDepth.data(), sourceVelocity.data(), receiverDepth.data(),
        receiverVelocity.data(), npoints, headWaveVelocity.data(), gradient.data());
    phase = (p==Pn ? "Pn" : (p==Sn ? "Sn" : (p==Pg ? "Pg" : (p==Lg ? "Lg"
        : "unknown phase"))));
}
 
inline void SlbmInterface::getGreatCircleNodeInfo(
    int** neighbors,
    double** coefficients,
    const int& maxpoints,
    const int& maxnodes,
    int& npoints,
    int* nnodes
    )
{
    if (!isValid())
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getGreatCircleNodeInfo" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
 
    greatCircle->getNodeInfo(neighbors, coefficients, maxpoints, maxnodes, npoints, nnodes);
}
inline void SlbmInterface::getGreatCircleNodeInfo(
    vector<vector<int> >& neighbors, vector<vector<double> >& coefficients)
    {
        if (!isValid())
        {
            ostringstream os;
            os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
              os << endl << "ERROR in SlbmInterface::getGreatCircleNodeInfo" << endl
                << "GreatCircle is invalid." << endl
                  << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
            throw SLBMException(os.str(),113);
        }
 
        greatCircle->getNodeInfo(neighbors, coefficients);
    }
 
inline void SlbmInterface::getGreatCircleLocations(
        double lat[],
        double lon[],
        double depth[],
        int& npoints
        )
{
    if (!isValid())
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getGreatCircleData" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
 
    npoints = greatCircle->getNProfiles();
    Location loc;
    for (int i=0; i<greatCircle->getNProfiles(); i++)
    {
        // get the location of this layer profile
        greatCircle->getLayerProfileLocation(i, loc);
        lat[i] = loc.getLat();
        lon[i] = loc.getLon();
        depth[i] = loc.getDepth();
    }
}
 
inline void SlbmInterface::getGreatCircleLocations(
        vector<double> &lat,
        vector<double> &lon,
        vector<double> &depth
        )
{
    if (!isValid())
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getGreatCircleData" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    int npoints = getNHeadWavePoints();
    lat.resize(npoints);
    lon.resize(npoints);
    depth.resize(npoints);
    getGreatCircleLocations(lat.data(), lon.data(), depth.data(), npoints);
}
 
 
inline void SlbmInterface::getInterpolatedPoint(
        const double& lat,
        const double& lon,
        int* nodeIds,
        double* coefficients,
        int& nWeights,
        double depth[NLAYERS],
        double pvelocity[NLAYERS],
        double svelocity[NLAYERS],
        double& pgradient,
        double& sgradient
        )
{
    Location location(lat, lon, 0.);
    QueryProfile* profile = grid->getQueryProfile(location);
    profile->getData(nodeIds, coefficients, nWeights, depth, pvelocity, svelocity, pgradient, sgradient);
    delete profile;
}
 
inline void SlbmInterface::getInterpolatedPoint(
        const double &lat,
        const double &lon,
        vector<int> &nodeIds,
        vector<double> &coefficients,
        vector<double> &depth,
        vector<double> &pvelocity,
        vector<double> &svelocity,
        double &pgradient,
        double &sgradient
        )
{
    // create a new location at the lat/lon
    Location location(lat, lon, 0.);
 
    // interpolate a profile at that location
    QueryProfile* profile = grid->getQueryProfile(location);
 
    // get the number of points used for that interpolation
    int nWeights = profile->getNCoefficients();
 
    // allocate some memory
    nodeIds.resize(nWeights);
    coefficients.resize(nWeights);
    depth.resize(NLAYERS);
    pvelocity.resize(NLAYERS);
    svelocity.resize(NLAYERS);
 
    // call the method
    profile->getData(nodeIds.data(), coefficients.data(), nWeights,
        depth.data(), pvelocity.data(), svelocity.data(), pgradient, sgradient);
    delete profile;
}
 
inline void SlbmInterface::getInterpolatedTransect(
        double lat[],
        double lon[],
        const int& nLatLon,
        int** nodeId,
        double** coefficients,
        int* nNeighbors,
        double depth[][NLAYERS],
        double pvelocity[][NLAYERS],
        double svelocity[][NLAYERS],
        double pgradient[NLAYERS],
        double sgradient[NLAYERS],
        int& nInvalid
        )
{
    nInvalid = 0;
    for (int i=0; i<nLatLon; i++)
    {
        try
        {
            getInterpolatedPoint(lat[i], lon[i], nodeId[i],
                coefficients[i], nNeighbors[i], depth[i], pvelocity[i], svelocity[i],
                pgradient[i], sgradient[i]);
        }
        catch (SLBMException ex)
        {
            for (int j=0; j<nNeighbors[i]; j++)
            {
                nodeId[i][j] = -1;
                coefficients[i][j] = NA_VALUE;
            }
 
            for (int j=0; j<NLAYERS; j++)
            {
                depth[i][j] = NA_VALUE;
                pvelocity[i][j] = NA_VALUE;
                svelocity[i][j] = NA_VALUE;
            }
 
            pgradient[i] = NA_VALUE;
            sgradient[i] = NA_VALUE;
 
            ++nInvalid;
        }
    }
}
inline void SlbmInterface::getInterpolatedTransect(
        vector<double> lat,
        vector<double> lon,
        vector<vector<int> > &nodeId,
        vector<vector<double> > &coefficients,
        vector<vector<double> > &depth,
        vector<vector<double> > &pvelocity,
        vector<vector<double> > &svelocity,
        vector<double> &pgradient,
        vector<double> &sgradient,
        int &nInvalid
        )
{
    nInvalid = 0;
 
    // check lat/lon sizes match
    int nLatLon;
    if (lat.size() == lon.size())
    {
        nLatLon = lat.size();
    }
    else
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getInterpolatedTransect" << endl
            << "Size of lat and lon vectors do not match." << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
 
    // resize vectors
    nodeId.resize(nLatLon);
    coefficients.resize(nLatLon);
    depth.resize(nLatLon);
    pvelocity.resize(nLatLon);
    svelocity.resize(nLatLon);
    pgradient.resize(nLatLon);
    sgradient.resize(nLatLon);
    for (int i=0; i<nLatLon; i++)
    {
        depth[i].resize(NLAYERS);
        pvelocity[i].resize(NLAYERS);
        svelocity[i].resize(NLAYERS);
    }
 
 
    for (int i=0; i<nLatLon; i++)
    {
        try
        {
            // create a new Location at the input lat/lon
            Location location(lat[i], lon[i], 0.);
 
            // interpolate a profile at that location
            QueryProfile profile(*grid, location);  // create a query profile at that location
 
            // get the number of points used for that interpolation
            int nWeights = profile.getNCoefficients();
 
            // resize vectors
            nodeId[i].resize(nWeights);
            coefficients[i].resize(nWeights);
 
            getInterpolatedPoint(lat[i], lon[i], nodeId[i].data(), coefficients[i].data(), nWeights,
                depth[i].data(), pvelocity[i].data(), svelocity[i].data(), pgradient[i], sgradient[i]);
 
        }
        catch (SLBMException ex)
        {
            nodeId[i].push_back(-1);
            coefficients[i].push_back(NA_VALUE);
            pgradient[i] = NA_VALUE;
            sgradient[i] = NA_VALUE;
 
            for (int j=0; j<NLAYERS; j++)
            {
                depth[i][j]     = NA_VALUE;
                pvelocity[i][j] = NA_VALUE;
                svelocity[i][j] = NA_VALUE;
            }
 
            ++nInvalid;
        }
    }
}
 
inline void SlbmInterface::getNGridNodes(int& n)
{
    if (!grid)
    {
        n = -1;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getNGridNodes" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    n = grid->getNNodes();
}
 
inline int SlbmInterface::getNGridNodes()
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getNGridNodes" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    return grid->getNNodes();
}
 
inline void SlbmInterface::getNHeadWavePoints(int& nHeadWavePoints)
{
    if (!greatCircle)
    {
        nHeadWavePoints = -1;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getNHeadWavePoints" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    nHeadWavePoints = greatCircle->getNProfiles();
}
 
//inline void SlbmInterface::getMaxInterpotionNodes(int& maxNodes)
//{
//  if (!greatCircle)
//  {
//      nHeadWavePoints = -1;
//      ostringstream os;
//      os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
//        os << endl << "ERROR in SlbmInterface::getMaxInterpotionNodes" << endl
//          << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
//          << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
//      throw SLBMException(os.str(),113);
//  }
//
//  maxNodes = 1;
//  for (int i=0; i<greatCircle->getNProfiles(); i++)
//      if (greatCircle->profiles[i]->getNodes().size() > maxNodes)
//          maxNodes = greatCircle->profiles[i]->getNodes().size();
//}
 
inline void SlbmInterface::getSlownessUncertainty( const int& phase, const double& distance, double& uncert )
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getSlownessUncertainty" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
 
    if (grid->getUncertainty()[phase][SH] == NULL)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
        os << endl << "ERROR in SlbmInterface::getSlownessUncertainty" << endl
            << "Uncertainty object is invalid.." << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),602);
    }
    else
    {
        // retrieve slowness uncertainty in sec/radian
        uncert =  grid->getUncertainty()[phase][SH]->getUncertainty( distance );
    }
}
 
inline void SlbmInterface::getSlownessUncertainty(double& slownessUncertainty)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getSlownessUncertainty" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
 
    getSlownessUncertainty(greatCircle->getPhase(),
        greatCircle->getDistance(), slownessUncertainty);
}
 
inline  string SlbmInterface::getUncertaintyTable(const int& phase, const int& attribute)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getUncertaintyTable" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
 
    if (grid->getUncertainty()[phase][attribute] == NULL)
    {
        ostringstream os;
        os << "No uncertainty information is available for phase " <<
                UncertaintyPIU::getPhase(phase) << " attribute " <<
                UncertaintyPIU::getAttribute(attribute) << endl;
        return os.str();
    }
    return grid->getUncertainty()[phase][attribute]->toStringTable();
}
 
inline  string SlbmInterface::getUncertaintyFileFormat(const int& phase, const int& attribute)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getUncertaintyFileFormat" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
 
    if (grid->getUncertainty()[phase][attribute] == NULL)
    {
        ostringstream os;
        os << "No uncertainty information is available for phase " << phase << " attribute " << attribute << endl;
        return os.str();
    }
    return grid->getUncertainty()[phase][attribute]->toStringFile();
}
 
inline  void SlbmInterface::getZhaoParameters(double& Vm, double& Gm, double& H, double& C, double& Cm, int& udSign)
{
    if (!greatCircle)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getZhaoParameters" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    greatCircle->getZhaoParameters(Vm, Gm, H, C, Cm, udSign);
}
 
inline  void SlbmInterface::getPgLgComponents(double& tTotal,
                                              double& tTaup, double& tHeadwave,
                                              double& pTaup, double& pHeadwave,
                                              double& trTaup, double& trHeadwave)
{
    if (!greatCircle)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getPgLgComponents" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    greatCircle->getPgLgComponents(tTotal, tTaup, tHeadwave,
        pTaup, pHeadwave, trTaup, trHeadwave);
}
 
 
inline  void SlbmInterface::getNodeNeighbors(const int& nid, int neighbors[], int& nNeighbors)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getNodeNeighbors" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    grid->getNodeNeighbors(nid, neighbors, nNeighbors);
}
 
inline  void SlbmInterface::getActiveNodeNeighbors(const int& nid, int neighbors[], int& nNeighbors)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getActiveNodeNeighbors" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    grid->getActiveNodeNeighbors(nid, neighbors, nNeighbors);
}
 
inline  void SlbmInterface::getNodeNeighbors(const int& nid, vector<int>& neighbors)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getNodeNeighbors" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    grid->getNodeNeighbors(nid, neighbors);
}
 
inline  void SlbmInterface::getActiveNodeNeighbors(const int& nid, vector<int>& neighbors)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getActiveNodeNeighbors" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    grid->getActiveNodeNeighbors(nid, neighbors);
}
 
inline  void SlbmInterface::getNodeNeighborInfo(const int& nid, int neighbors[],
            double distance[], double azimuth[], int& nNeighbors)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getNodeNeighborInfo" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    grid->getNodeNeighborInfo(nid, neighbors, distance, azimuth, nNeighbors);
}
 
inline  void SlbmInterface::getNodeNeighborInfo(const int nid,
        vector<int> &neighbors, vector<double> &distance, vector<double> &azimuth)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getNodeNeighborInfo" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    grid->getNodeNeighborInfo(nid, neighbors, distance, azimuth);
}
 
inline  void SlbmInterface::getActiveNodeNeighborInfo(const int& nid, int neighbors[],
            double distance[], double azimuth[], int& nNeighbors)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getActiveNodeNeighborInfo" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    grid->getActiveNodeNeighborInfo(nid, neighbors, distance, azimuth, nNeighbors);
}
 
inline  void SlbmInterface::getActiveNodeNeighborInfo(const int& nid,
            vector<int>& neighbors, vector<double>& distance, vector<double>& azimuth)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getActiveNodeNeighborInfo" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    grid->getActiveNodeNeighborInfo(nid, neighbors, distance, azimuth);
}
 
inline  void SlbmInterface::getNodeSeparation(const int& node1, const int& node2, double& distance)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getNodeSeparation" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    grid->getNodeSeparation(node1, node2, distance);
}
 
inline  void SlbmInterface::getNodeAzimuth(const int& node1, const int& node2, double& azimuth)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getNodeAzimuth" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    grid->getNodeAzimuth(node1, node2, azimuth);
}
 
inline void SlbmInterface::initializeActiveNodes(const double& latmin, const double& lonmin,
        const double& latmax, const double& lonmax)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::initializeActiveNodes" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    grid->initializeActiveNodes(latmin, lonmin, latmax, lonmax);
}
 
inline void SlbmInterface::initializeActiveNodes(double* lat, double* lon, const int& npoints, const bool& inDegrees)
{
    vector<double*> v;
    v.reserve(npoints);
    if (inDegrees)
        for (int i=0; i<npoints; ++i)
            v.push_back(GeoTessUtils::getVectorDegrees(lat[i], lon[i]));
    else
        for (int i=0; i<npoints; ++i)
            v.push_back(GeoTessUtils::getVector(lat[i], lon[i]));
    initializeActiveNodes(v);
}
inline void SlbmInterface::initializeActiveNodes(const vector<double> lat, const vector<double> lon, const bool& inDegrees)
{
    int npoints;
    if (lat.size() == lon.size())
    {
        npoints = lat.size();
    }
    else
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getInterpolatedTransect" << endl
            << "Size of lat and lon vectors do not match." << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    vector<double*> v;
    v.reserve(npoints);
    if (inDegrees)
        for (int i=0; i<npoints; ++i)
            v.push_back(GeoTessUtils::getVectorDegrees(lat[i], lon[i]));
    else
        for (int i=0; i<npoints; ++i)
            v.push_back(GeoTessUtils::getVector(lat[i], lon[i]));
    initializeActiveNodes(v);
}
 
inline void SlbmInterface::initializeActiveNodes(vector<double*>& unitVectors)
{
    initializeActiveNodes(new GeoTessPolygon(unitVectors));
}
 
inline void SlbmInterface::initializeActiveNodes(GeoTessPolygon* polygon)
{
    grid->initializeActiveNodes(polygon);
}
 
inline void SlbmInterface::initializeActiveNodes(const string& fileName)
{
    grid->initializeActiveNodes(new GeoTessPolygon(fileName));
}
 
inline void SlbmInterface::clearActiveNodes()
{
    grid->clearActiveNodes();
}
 
 
inline int SlbmInterface::getNActiveNodes()
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::nextActiveNode" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    return grid->getNActiveNodes();
}
 
inline int SlbmInterface::getGridNodeId(int activeNodeId)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getGridNodeId" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    return grid->getGridNodeId(activeNodeId);
}
 
inline int SlbmInterface::getActiveNodeId(int gridNodeId)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getActiveNodeId" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    return grid->getActiveNodeId(gridNodeId);
}
 
inline void SlbmInterface::getNodeHitCount(const int& nodeId, int& hitCount)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getNodeHitCount" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    grid->getNodeHitCount(nodeId, hitCount);
}
 
inline void SlbmInterface::clearNodeHitCount()
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::clearNodeHitCount" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    grid->clearNodeHitCount();
}
 
inline void SlbmInterface::getAverageMantleVelocity(const int& type, double& velocity)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::setAverageMantleVelocity" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    velocity = grid->getAverageMantleVelocity(type);
}
 
inline double SlbmInterface::getAverageMantleVelocity(const int &type)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::setAverageMantleVelocity" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    return grid->getAverageMantleVelocity(type);
}
 
inline void SlbmInterface::setAverageMantleVelocity(const int& type,
                                                    const double& velocity)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::setAverageMantleVelocity" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    grid->setAverageMantleVelocity(type, velocity);
}
 
inline void SlbmInterface::getTessId(string& tessId)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getTessId" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    tessId = grid->getTessId();
}
 
inline void SlbmInterface::getFractionActive(double& fractionActive)
{
    if (!isValid())
    {
        fractionActive = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getFractionActive" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
    fractionActive = greatCircle->getFractionActive();
}
 
inline string SlbmInterface::getClassCount()
{
    ostringstream os;
    os << "Class counts:" << endl;
    os << "GreatCircle         = " << GreatCircle::getClassCount()  << endl;
    os << "GridProfile         = " << GridProfile::getClassCount()  << endl;
    os << "GeoStack            = " << GeoStack::getClassCount() << endl;
    os << "InterpolatedProfile = " << InterpolatedProfile::getClassCount()  << endl;
    os << "CrustalProfile      = " << CrustalProfile::getClassCount()   << endl;
    os << "LayerProfile        = " << LayerProfile::getClassCount() << endl;
    os << "QueryProfile        = " << QueryProfile::getClassCount() << endl;
    os << "Location            = " << Location::getClassCount() << endl;
    return os.str();
}
 
inline const string& SlbmInterface::getModelPath() const { return grid->getModelPath(); }
 
inline void SlbmInterface::getPiercePointSource(double& lat, double& lon, double& depth)
{
    if (!isValid())
    {
        lat = lon = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getPiercePointSource" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
 
 
    if ((greatCircle->getPhase() == Pg || greatCircle->getPhase() == Lg) && !greatCircle->getSourceProfile()->isInCrust())
    {
        lat = lon = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getPiercePointSource" << endl
            << "Cannot compute moho pierce points for crustal phase (" << getPhase() << ") when source is in the mantle." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
 
    Location loc;
    greatCircle->getGreatCircleLocation(greatCircle->getSourceDistance(), loc);
 
    lat = loc.getLat();
    lon = loc.getLon();
 
    QueryProfile* profile = NULL;
    profile = grid->getQueryProfile(loc);
    if (greatCircle->getPhase() == Pg || greatCircle->getPhase() == Lg)
        depth = profile->getDepth()[MIDDLE_CRUST_G];
    else
        depth = profile->getDepth()[MANTLE];
    delete profile;
}
 
inline void SlbmInterface::getPiercePointReceiver(double& lat, double& lon, double& depth)
{
    if (!isValid())
    {
        lat = lon = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getPiercePointReceiver" << endl
            << "GreatCircle is invalid." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
 
 
    if ((greatCircle->getPhase() == Pg || greatCircle->getPhase() == Lg) && !greatCircle->getReceiverProfile()->isInCrust())
    {
        lat = lon = NA_VALUE;
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getPiercePointReceiver" << endl
            << "Cannot compute moho pierce points for crustal phase (" << getPhase() << ") when receiver is in the mantle." << endl
              << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),113);
    }
 
    Location loc;
    greatCircle->getGreatCircleLocation(
        greatCircle->getDistance()-greatCircle->getReceiverDistance(), loc);
 
    lat = loc.getLat();
    lon = loc.getLon();
 
    QueryProfile* profile = NULL;
    profile = grid->getQueryProfile(loc);
    if (greatCircle->getPhase() == Pg || greatCircle->getPhase() == Lg)
        depth = profile->getDepth()[MIDDLE_CRUST_G];
    else
        depth = profile->getDepth()[MANTLE];
    delete profile;
}
 
inline void SlbmInterface::getGreatCirclePoints(
        const double& aLat,
        const double& aLon,
        const double& bLat,
        const double& bLon,
        const int& npoints,
        double latitude[],
        double longitude[])
{
    Location a(aLat, aLon, 0.);
    Location b(bLat, bLon, 0.);
    double dx = a.distance(b)/(npoints-1);
 
    double moveDirection[3];
    a.vectorTripleProduct(b, moveDirection);
    for (int i=0; i<npoints; i++)
    {
        a.move(moveDirection, i*dx, b);
        latitude[i] = b.getLat();
        longitude[i] = b.getLon();
    }
}
 
inline void SlbmInterface::getGreatCirclePoints(
        const double& aLat,
        const double& aLon,
        const double& bLat,
        const double& bLon,
        const int& npoints,
        vector<double> &latitude,
        vector<double> &longitude)
{
    latitude.resize(npoints);
    longitude.resize(npoints);
 
    Location a(aLat, aLon, 0.);
    Location b(bLat, bLon, 0.);
    double dx = a.distance(b)/(npoints-1);
 
    double moveDirection[3];
    a.vectorTripleProduct(b, moveDirection);
    for (int i=0; i<npoints; i++)
    {
        a.move(moveDirection, i*dx, b);
        latitude[i] = b.getLat();
        longitude[i] = b.getLon();
    }
}
 
inline void SlbmInterface::getGreatCirclePointsOnCenters(
        const double& aLat,
        const double& aLon,
        const double& bLat,
        const double& bLon,
        const int& npoints,
        double latitude[],
        double longitude[])
{
    Location a(aLat, aLon, 0.);
    Location b(bLat, bLon, 0.);
    double dx = a.distance(b)/npoints;
 
    double moveDirection[3];
    a.vectorTripleProduct(b, moveDirection);
    for (int i=0; i<npoints; i++)
    {
        a.move(moveDirection, (i+0.5)*dx, b);
        latitude[i] = b.getLat();
        longitude[i] = b.getLon();
    }
}
 
inline void SlbmInterface::getGreatCirclePointsOnCenters(
        const double& aLat,
        const double& aLon,
        const double& bLat,
        const double& bLon,
        const int& npoints,
        vector<double> &latitude,
        vector<double> &longitude)
{
    latitude.resize(npoints);
    longitude.resize(npoints);
 
    Location a(aLat, aLon, 0.);
    Location b(bLat, bLon, 0.);
    double dx = a.distance(b)/npoints;
 
    double moveDirection[3];
    a.vectorTripleProduct(b, moveDirection);
    for (int i=0; i<npoints; i++)
    {
        a.move(moveDirection, (i+0.5)*dx, b);
        latitude[i] = b.getLat();
        longitude[i] = b.getLon();
    }
}
 
inline void SlbmInterface::getDistAz(const double& aLat, const double& aLon,
    const double& bLat, const double& bLon,
    double& distance, double& azimuth, const double& naValue)
{
    Location ptA(aLat, aLon);
    Location ptB(bLat, bLon);
    distance = ptA.distance(ptB);
    azimuth = ptA.azimuth(ptB, naValue);
}
 
inline void SlbmInterface::movePoint(const double& aLat, const double& aLon,
        const double& distance, const double& azimuth,
        double& bLat, double& bLon)
{
    Location ptA(aLat, aLon);
    Location ptB;
    ptA.move(azimuth, distance, ptB);
    bLat = ptB.getLat();
    bLon = ptB.getLon();
}
 
inline void SlbmInterface::setInterpolatorType(const string& interpolatorType)
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::setInterpolatorType" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    grid->setInterpolatorType(interpolatorType);
}
 
inline string SlbmInterface::getInterpolatorType()
{
    if (!grid)
    {
        ostringstream os;
        os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(9);
          os << endl << "ERROR in SlbmInterface::getInterpolatorType" << endl
            << "Grid is invalid.  Has the earth model been loaded with call to loadVelocityModel()?" << endl
            << "Version " << SlbmVersion << "  File " << __FILE__ << " line " << __LINE__ << endl << endl;
        throw SLBMException(os.str(),114);
    }
    return grid->getInterpolatorType();
}
 
 
 
} // end slbm namespace
 
#endif // SlbmInterface_H