TPVelocityModels.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 TPVELOCITYMODELS_H
#define TPVELOCITYMODELS_H
 
// **** _SYSTEM INCLUDES_ ******************************************************
#include <float.h>
#include <vector>
#include <cmath>
#include <iostream>
 
using namespace std;
// use standard library objects
 
// **** _LOCAL INCLUDES_ *******************************************************
 
#include "TauPGlobals.h"
#include "IntegrateFunction.h"
#include "DataBuffer.h"
#include "CPPUtils.h"
 
using util::DataBuffer;
 
using namespace geotess;
 
//template<class F> class util::IntegrateFunction;
 
// **** _BEGIN TAUP NAMESPACE_ *************************************************
 
namespace taup {
 
// **** _CLASS CONSTANTS_ ******************************************************
// **** _FORWARD REFERENCES_ ***************************************************
 
class TPVelocityLayer;
class VelocityConst;
class VelocityPower;
class VelocityLinear;
class VelocityQuadratic;
class VelocityCubic;
 
template<class V> class TPdDistdr;
template<class V> class TPdTaudr;
 
// *****************************************************************************
//
//
//
// *****************************************************************************
template<class V>
class TPdDistdr
{
  public:
 
    //
    TPdDistdr(V& v) : diV(v), diP(0.0) {};
 
    TPdDistdr(const TPdDistdr<V>& di) : diV(di.diV), diP(di.diP) {};
 
    virtual ~TPdDistdr() {};
 
    //
    TPdDistdr&     operator=(const TPdDistdr<V>& di)
    {
      diP = di.diP;
      return *this;
    };
 
    double         operator()(double r)
    {
      double pv = diP * diV(r);
      double d  = fabs(r - pv) * (r + pv);
      if (d == 0.0)
        return 1.0 / sqrt(DBL_EPSILON * (r + pv));
      else
        return pv / r / sqrt(d);
    };
 
    void           setP(double p) {diP = p;};
 
  private:
 
    double         diP;
 
    V&             diV;
};
 
// *****************************************************************************
//
//
//
// *****************************************************************************
template<class V>
class TPdTaudr
{
  public:
 
    //
    TPdTaudr(V& v) : tiV(v), tiP(0.0) {};
 
    TPdTaudr(const TPdTaudr& ti) : tiV(ti.tiV), tiP(ti.tiP) {};
 
    virtual ~TPdTaudr() {};
 
    //
    TPdTaudr&     operator=(const TPdTaudr& ti)
    {
      tiP = ti.tiP;
      return *this;
    };
 
    double         operator()(double r)
    {
      double vr = tiV(r);
      double pv = tiP * vr;
 
      return sqrt(fabs(r - pv) * (r + pv)) / r / vr;
    };
 
    void           setP(double p) {tiP = p;};
 
  private:
 
    double         tiP;
 
    V&             tiV;
};
 
// *****************************************************************************
//
//
//
// *****************************************************************************
class TAUP_EXP TPVelocityLayer
{
  public:
 
    // **** _PUBLIC LIFCYCLES_ *************************************************
 
    TPVelocityLayer() : vlRt(0.0), vlRb(0.0), vlVt(-1.0), vlVb(-1.0),
                        vlPt(-1.0), vlPb(-1.0), vlRLast(-1.0),
                        vlDistT(0.0), vlDistB(0.0), vlIRt(0.0), vlIRb(0.0),
                        vlPmin(-1.0), vlLayerType(-1),
                        vlPCrit(-1.0), vlRTurn(-1.0), vlDistCrit(-1.0),
                        vlSplitLayer(false), vlInvalidRay(false),
                        vlPassingRay(false), vlTurningRay(false),
                        vldDistdP_T(0.0), vldDistdP_B(0.0), vlVTurn(-1.0),
                        vlPhaseUpperIDef(false), vlPhaseLowerIDef(false),
                        vlLayerName("")
    {
      vlPhase = vlPhaseUpper = vlPhaseLower = vlPhaseType= "";
      vlPhaseIUpper = vlPhaseILower = "";
    };
 
    TPVelocityLayer(double rt, double rb, const string& layrnam) :
                    vlRt(rt), vlRb(rb), vlVt(-1.0), vlVb(-1.0),
                    vlPt(-1.0), vlPb(-1.0), vlRLast(-1.0),
                    vlDistT(0.0), vlDistB(0.0), vlIRt(0.0), vlIRb(0.0),
                    vldDistdP_T(0.0), vldDistdP_B(0.0), vlRTurn(-1.0),
                    vlPmin(-1.0), vlLayerType(-1), vlPCrit(-1.0),
                    vlDistCrit(-1.0), vlSplitLayer(false), vlVTurn(-1.0),
                    vlInvalidRay(false), vlPassingRay(false),
                    vlPhaseUpperIDef(false), vlPhaseLowerIDef(false),
                    vlTurningRay(false), vlLayerName(layrnam)
    {
      vlPhase = vlPhaseUpper = vlPhaseLower = vlPhaseType = "";
      vlPhaseIUpper = vlPhaseILower = "";
    };
 
    TPVelocityLayer(const TPVelocityLayer& vl) : vlRt(vl.vlRt), vlRb(vl.vlRb),
                    vlVt(vl.vlVt), vlVb(vl.vlVb), vlPt(vl.vlPt), vlPb(vl.vlPb),
                    vlRLast(vl.vlRLast), vlIRt(vl.vlIRt), vlIRb(vl.vlIRb),
                    vlDistT(vl.vlDistT), vlDistB(vl.vlDistB),
                    vldDistdP_T(vl.vldDistdP_T), vldDistdP_B(vl.vldDistdP_B),
                    vlInvalidRay(vl.vlInvalidRay), vlRTurn(vl.vlRTurn),
                    vlPmin(vl.vlPmin), vlLayerType(vl.vlLayerType),
                    vlPCrit(vl.vlPCrit), vlDistCrit(vl.vlDistCrit),
                    vlSplitLayer(vl.vlSplitLayer), vlVTurn(vl.vlVTurn),
                    vlPassingRay(vl.vlPassingRay),
                    vlTurningRay(vl.vlTurningRay), vlLayerName(vl.vlLayerName),
                    vlPhase(vl.vlPhase), vlPhaseUpper(vl.vlPhaseUpper),
                    vlPhaseLower(vl.vlPhaseLower), vlPhaseType(vl.vlPhaseType),
                    vlPhaseIUpper(vl.vlPhaseIUpper),
                    vlPhaseUpperIDef(vl.vlPhaseUpperIDef),
                    vlPhaseLowerIDef(vl.vlPhaseLowerIDef),
                    vlPhaseILower(vl.vlPhaseILower)
    {};
 
    virtual ~TPVelocityLayer() {};
 
    // **** _PUBLIC OPERATORS_ *************************************************
 
    TPVelocityLayer&       operator=(const TPVelocityLayer& vl)
    {
      vlRt          = vl.vlRt;
      vlRb          = vl.vlRb;
      vlVt          = vl.vlVt;
      vlVb          = vl.vlVb;
      vlPt          = vl.vlPt;
      vlPb          = vl.vlPb;
      vlRTurn       = vl.vlRTurn;
      vlVTurn       = vl.vlVTurn;
 
      vlLayerName      = vl.vlLayerName;
      vlPhaseType      = vl.vlPhaseType;
      vlPhase          = vl.vlPhase;
      vlPhaseUpper     = vl.vlPhaseUpper;
      vlPhaseLower     = vl.vlPhaseLower;
      vlPhaseIUpper    = vl.vlPhaseIUpper;
      vlPhaseILower    = vl.vlPhaseILower;
      vlPhaseUpperIDef = vl.vlPhaseUpperIDef;
      vlPhaseLowerIDef = vl.vlPhaseLowerIDef;
 
      vlIRt         = vl.vlIRt;
      vlIRb         = vl.vlIRb;
      vlRLast       = vl.vlRLast;
      vlDistT       = vl.vlDistT;
      vlDistB       = vl.vlDistB;
      vldDistdP_T   = vl.vldDistdP_T;
      vldDistdP_B   = vl.vldDistdP_B;
 
      vlPmin        = vl.vlPmin;
      vlLayerType   = vl.vlLayerType;
 
      vlPCrit       = vl.vlPCrit;
      vlDistCrit    = vl.vlDistCrit;
      vlSplitLayer  = vl.vlSplitLayer;
 
      vlPassingRay  = vl.vlPassingRay;
      vlTurningRay  = vl.vlTurningRay;
      vlInvalidRay  = vl.vlInvalidRay;
 
      return *this;
    };
 
    virtual double         operator()(double r) = ABSTRACT;
 
    // **** _PUBLIC METHODS_ ***************************************************
 
    void                   init()
    {
      vlVt = operator()(vlRt);
      vlVb = operator()(vlRb);
      vlPt = pAtR(vlRt);
      vlPb = pAtR(vlRb);
    };
 
    static TPVelocityLayer* newModelCopy(TPVelocityLayer* tpvl);
 
    static TPVelocityLayer* newModelCopy(const string& cnam,
                                         DataBuffer& buffer);
 
    virtual void           writeNormRadius(ostream& os) const {};
 
    virtual void           writeVelocity(ostream& os) const = ABSTRACT;
 
    virtual double         rAtP(double p) = ABSTRACT;
 
    double                 pAtR(double r) {return r / operator()(r);};
 
    double                 integDistance(double p, double r = -1.0)
    {
      double dist = 0.0;
 
      // set ray and integrate if valid
 
      setRay(p);
      if (!vlInvalidRay)
      {
        // use the turning depth (or bottom of the layer) vlRTurn if
        // r is not defined (-1.0) or is smaller than vlRTurn. Otherwise,
        // integrate from r to the top of the layer.
 
        double rb = vlRTurn;
        bool   tflg = vlTurningRay;
        if ((r != -1.0) && (r > vlRTurn))
        {
          rb = r;
          tflg = false;
        }
 
        // assign integration parameters and perform integration
 
        vlRLast = rb;
        dist = integrateDistance(p, rb, tflg);
      }
 
      // return distance
 
      return dist;
    };
 
    double                 integDistance(double p, double r1, double r2)
    {
      double dist = 0.0;
 
      // set updown ray flags and integrate if valid or turning
 
      setUpDownRay(p, r1, r2);
      if (!vlInvalidRay)
        dist = integrateDistance(p, vlIRb, vlIRt, vlTurningRay);
      else if (vlTurningRay)
        dist = integrateDistance(p, vlRTurn, vlIRt, true);
 
      // return distance
 
      return dist;
    };
 
    double                 integTime(double p, double r = -1.0)
    {
      double tim = 0.0;
 
      // get distance component and integrate Tau if valid
 
      double d = integDistance(p, r);
      if (!vlInvalidRay)
      {
        // ray is valid calculate time from tau or straight-away
 
        tim = integrateTime(p, vlRLast);
        if (isTimeIntegralTau()) tim += p * d;
      }
 
      return tim;
    };
 
    double                 integTime(double p, double r1, double r2)
    {
      double tim = 0.0;
 
      // set updown ray flags and integrate if valid or turning
 
      setUpDownRay(p, r1, r2);
      double r;
      if (!vlInvalidRay)
        r = vlIRb;
      else if (vlTurningRay)
        r = vlRTurn;
      else
        return tim;
 
      tim = integrateTime(p, r, vlIRt);
      if (isTimeIntegralTau()) tim += p * integrateDistance(p, r, vlIRt);
 
      return tim;
    };
 
    //
    virtual double         integrateDistance(double p, double ra,
                                             bool r_open = false) = ABSTRACT;
 
    //
    virtual double         integrateDistance(double p, double ra, double rb,
                                             bool r_open = false) = ABSTRACT;
 
    virtual double         integrateTime(double p, double ra) = ABSTRACT;
 
    virtual double         integrateTime(double p, double ra,
                                         double rb) = ABSTRACT;
 
    virtual bool           isTimeIntegralTau() const {return true;};
 
    double                 getRt() const {return vlRt;};
 
    double                 getRb() const {return vlRb;};
 
    double                 getVt() const {return vlVt;};
 
    double                 getVb() const {return vlVb;};
 
    double                 getPt() const {return vlPt;};
 
    double                 getPb() const {return vlPb;};
 
    double                 getTurningRadius() const {return vlRTurn;};
 
    double                 getTurningVelocity() const {return vlVTurn;};
 
    void                   setPmin(double pmin) {vlPmin = pmin;};
 
    double                 getPmin() const {return vlPmin;};
 
    void                   setDistT(double d) {vlDistT = d;};
 
    double                 getDistT() const {return vlDistT;};
 
    void                   setDistB(double d) {vlDistB = d;};
 
    double                 getDistB() const {return vlDistB;};
 
    void                   setdDistdPT(double dDdP) {vldDistdP_T = dDdP;};
 
    double                 getdDistdPT() const {return vldDistdP_T;};
 
    void                   setdDistdPB(double dDdP) {vldDistdP_B = dDdP;};
 
    double                 getdDistdPB() const {return vldDistdP_B;};
 
    void                   setLayerType(int lt) {vlLayerType = lt;};
 
    int                    getLayerType() const {return vlLayerType;};
 
    const string&          getLayerName() const {return vlLayerName;};
 
    void                   setPhaseTypeP() {vlPhaseType = "P";};
 
    void                   setPhaseTypeS() {vlPhaseType = "S";};
 
    void                   setPhaseType(const string& phtype)
    {vlPhaseType = phtype;};
 
    const string&          getPhaseType() const {return vlPhaseType;};
 
    void                   setPhaseName(const string& name) {vlPhase = name;};
 
    const string&          getPhaseName() const {return vlPhase;};
 
    void                   setPhaseNameUpper(const string& name)
    {vlPhaseUpper = name;};
 
    const string&          getPhaseNameUpper() const {return vlPhaseUpper;};
 
    void                   setPhaseNameLower(const string& name)
    {vlPhaseLower = name;};
 
    const string&          getPhaseNameLower() const {return vlPhaseLower;};
 
    void                   setPhaseNameDiff(const string& name)
    {vlPhaseIUpper = name;};
 
    const string&          getPhaseNameDiff() const {return vlPhaseIUpper;};
 
    void                   setPhaseNameDiffLower(const string& name)
    {vlPhaseILower = name;};
 
    const string&          getPhaseNameDiffLower() const {return vlPhaseILower;};
 
    void                   setPhaseDiffDef(bool def)
    {vlPhaseUpperIDef = def;};
 
    void                   setPhaseDiffLowerDef(bool def)
    {vlPhaseLowerIDef = def;};
 
    bool                   isPhaseDiffDefined() const
    {return vlPhaseUpperIDef;};
 
    bool                   isPhaseDiffLowerDefined() const
    {return vlPhaseLowerIDef;};
 
    double                 getPCrit() const {return vlPCrit;};
 
    double                 getDistCrit() const {return vlDistCrit;};
 
    bool                   isSplitLayer() const {return vlSplitLayer;};
 
    void                   setSplitLayer(double pcrit, double dcrit)
    {
      vlPCrit    = pcrit;
      vlDistCrit = dcrit;
      vlSplitLayer = true;
    };
 
    bool                   invalidRay() const {return vlInvalidRay;};
 
    bool                   passingRay() const {return vlPassingRay;};
 
    bool                   turningRay() const {return vlTurningRay;};
 
    string                 toString() const;
 
    virtual void           toStream(ostream& os, string indent) const;
 
    static  string         class_name() {return "TPVelocityLayer";};
 
    virtual string         get_class_name() const {return class_name();};
 
    virtual int            classSize() const
    {return (int) sizeof(TPVelocityLayer);};
 
    virtual bool           isVelocityConstant() const {return false;};
    virtual bool           isVelocityPowerLaw() const {return false;};
    virtual bool           isVelocityLinear() const {return false;};
    virtual bool           isVelocityQuadratic() const {return false;};
    virtual bool           isVelocityCubic() const {return false;};
 
    virtual int            bufferSize() const
    {
      int bs = (int) vlLayerName.size() + (int) vlPhaseType.size() +
               (int) vlPhase.size() + (int) vlPhaseUpper.size() +
               (int) vlPhaseLower.size() + (int) vlPhaseIUpper.size() +
               (int) vlPhaseILower.size() + 14 * sizeof(int);
      bs += 2 * sizeof(int) + 3 * sizeof(char) + 16 * sizeof(double);
 
      return bs;
    };
 
    virtual void           serialize(DataBuffer& buffer)
    {
      buffer.writeString(vlLayerName);
      buffer.writeString(vlPhaseType);
      buffer.writeString(vlPhase);
      buffer.writeString(vlPhaseUpper);
      buffer.writeString(vlPhaseLower);
      buffer.writeString(vlPhaseIUpper);
      buffer.writeString(vlPhaseILower);
 
      buffer.writeInt32(vlLayerType);
 
      buffer.writeByte(vlSplitLayer);
      buffer.writeByte(vlPhaseUpperIDef);
      buffer.writeByte(vlPhaseLowerIDef);
 
      buffer.writeDouble(vlRt);
      buffer.writeRawDouble(vlRb);
      buffer.writeRawDouble(vlVt);
      buffer.writeRawDouble(vlVb);
      buffer.writeRawDouble(vlPt);
      buffer.writeRawDouble(vlPb);
      buffer.writeRawDouble(vlRTurn);
      buffer.writeRawDouble(vlVTurn);
      buffer.writeRawDouble(vlDistT);
      buffer.writeRawDouble(vlDistB);
      buffer.writeRawDouble(vldDistdP_T);
      buffer.writeRawDouble(vldDistdP_B);
      buffer.writeRawDouble(vlPmin);
      buffer.writeRawDouble(vlPCrit);
      buffer.writeRawDouble(vlDistCrit);
    };
 
    virtual void           deserialize(DataBuffer& buffer)
    {
      vlLayerName   = buffer.readString();
      vlPhaseType   = buffer.readString();
      vlPhase       = buffer.readString();
      vlPhaseUpper  = buffer.readString();
      vlPhaseLower  = buffer.readString();
      vlPhaseIUpper = buffer.readString();
      vlPhaseILower = buffer.readString();
 
      vlLayerType   = buffer.readInt32();
 
      vlSplitLayer     = buffer.readByte();
      vlPhaseUpperIDef = buffer.readByte();
      vlPhaseLowerIDef = buffer.readByte();
 
      vlRt        = buffer.readDouble();
      vlRb        = buffer.readRawDouble();
      vlVt        = buffer.readRawDouble();
      vlVb        = buffer.readRawDouble();
      vlPt        = buffer.readRawDouble();
      vlPb        = buffer.readRawDouble();
      vlRTurn     = buffer.readRawDouble();
      vlVTurn     = buffer.readRawDouble();
      vlDistT     = buffer.readRawDouble();
      vlDistB     = buffer.readRawDouble();
      vldDistdP_T = buffer.readRawDouble();
      vldDistdP_B = buffer.readRawDouble();
      vlPmin      = buffer.readRawDouble();
      vlPCrit     = buffer.readRawDouble();
      vlDistCrit  = buffer.readRawDouble();
    };
 
  protected:
 
    void                   setRay(double p)
    {
      vlInvalidRay = vlPassingRay = vlTurningRay = false;
 
      // determine limits of integration
 
      if (p < vlPt)
      {
        // valid integral determine bottom limit
 
        if (p < vlPb)
        {
          // ray passes through layer
 
          vlRTurn = vlRb;
          vlVTurn = vlVb;
          vlPassingRay = true;
        }
        else if (p == vlPb)
        {
          // ray turns in layer at bottom
 
          vlRTurn = vlRb;
          vlVTurn = vlVb;
          vlTurningRay = true;
        }
        else
        {
          // ray turns in layer
 
          vlRTurn = rAtP(p);
          vlVTurn = operator()(vlRTurn);
          vlTurningRay = true;
        }
      }
      else
      {
        // invalid ray
 
        vlInvalidRay = true;
      }
    };
 
    void                 setUpDownRay(double p, double r1, double r2)
    {
      vlPassingRay = vlTurningRay = false;
      vlInvalidRay = true;
 
      // ensure validity of radial bounds r1 and r2
 
      if ((r2 < vlRt) && (r1 > vlRb))
      {
        // set top radius to the smaller of vlRt and r1
 
        double rtt = vlRt;
        if (r1 < vlRt) rtt = r1;
 
        // set bottom radius to the larger of vlRb and r2
 
        double rbb = vlRb;
        if (r2 > vlRb) rbb = r2;
 
        // find minimum and maximum ray paramameter at rtt and rbb
 
        double ptt = pAtR(rtt);
        double pbb = pAtR(rbb);
        double pmin = pbb;
        double pmax = ptt;
        if (ptt < pbb)
        {
          pmin = ptt;
          pmax = pbb;
        }
 
        // if p > pmax ray is invalid so just return 0 ... otherwise
 
        if (p < pmin)
        {
          //valid up/down
 
          vlRLast = rbb;
          vlIRt   = rtt;
          vlIRb   = rbb;
          vlPassingRay = true;
          vlInvalidRay = false;
        }
        else if (p < pmax)
        {
          //valid turning ray
 
          vlRTurn = rAtP(p);
          vlVTurn = operator()(vlRTurn);
          vlRLast = vlRTurn;
          vlIRt   = rtt;
          vlIRb   = vlRTurn;
          vlTurningRay = true;
          if (p == pmin) vlInvalidRay = false;
        }
      }
    };
 
    string                 vlLayerName;
 
    string                 vlPhaseType;
 
    string                 vlPhase;
 
    string                 vlPhaseUpper;
 
    string                 vlPhaseLower;
 
    string                 vlPhaseIUpper;
 
    string                 vlPhaseILower;
 
    //
    int                    vlLayerType;
 
    bool                   vlSplitLayer;
 
    bool                   vlInvalidRay;
 
    bool                   vlPassingRay;
 
    bool                   vlTurningRay;
 
    bool                   vlPhaseUpperIDef;
 
    bool                   vlPhaseLowerIDef;
 
    double                 vlRt;
 
    double                 vlRb;
 
    double                 vlIRt;
 
    double                 vlIRb;
 
    double                 vlVt;
 
    double                 vlVb;
 
    double                 vlPt;
 
    double                 vlPb;
 
    // Turning radius or bottom of layer if this is a passing ray
    double                 vlRTurn;
 
    // Turning velocity or layer bottom velocity if this is a passing ray
    double                 vlVTurn;
 
    double                 vlRLast;
 
    double                 vlDistT;
 
    double                 vlDistB;
 
    double                 vldDistdP_T;
 
    double                 vldDistdP_B;
 
    double                 vlPmin;
 
    double                 vlPCrit;
 
    double                 vlDistCrit;
 
  private:
};
 
// *****************************************************************************
//
//
//
// *****************************************************************************
template<class V>
class TAUP_EXP VelocityIntegrate : public TPVelocityLayer
{
  public:
 
    // **** _PUBLIC LIFCYCLES_ *************************************************
 
    VelocityIntegrate() : vmDist(NULL), vmTau(NULL),
                          vmDistNI(NULL), vmTauNI(NULL), TPVelocityLayer() {};
 
    VelocityIntegrate(double rt, double rb, const string& layrnam) :
                      vmDist(NULL), vmTau(NULL),
                      vmDistNI(NULL), vmTauNI(NULL),
                      TPVelocityLayer(rt, rb, layrnam) {};
 
    VelocityIntegrate(const VelocityIntegrate& vi) :
                      vmDist(NULL), vmTau(NULL),
                      vmDistNI(NULL), vmTauNI(NULL), TPVelocityLayer(vi) {};
 
    virtual ~VelocityIntegrate()
    {
      // Delete integrand and integration function objects if they were
      // defined.
 
      if (vmDist)
      {
        delete vmDist;
        delete vmTau;
        delete vmDistNI;
        delete vmTauNI;
      }
    };
 
    // **** _PUBLIC OPERATORS_ *************************************************
 
    VelocityIntegrate&     operator=(const VelocityIntegrate& vi)
    {
      // delete old definitions and create new ones if they were defined
 
      if (vmDist)
      {
        delete vmDist;
        vmDist = (TPdDistdr<V>*) NULL;
        delete vmTau;
        vmTau  = (TPdTaudr<V>*) NULL;
        delete vmDistNI;
        vmDistNI = (util::IntegrateFunction<TPdDistdr<V> >*) NULL;
        delete vmTauNI;
        vmTauNI  = (util::IntegrateFunction<TPdTaudr<V> >*) NULL;
      }
 
      TPVelocityLayer::operator=(vi);
 
      return *this;
    };
 
    // **** _PUBLIC METHODS_ ***************************************************
 
    //
    virtual double         integrateDistance(double p, double ra, double rb,
                                             bool r_open = false)
    {
      // create the numerical integration objects if they do not yet exist
 
      if (!vmDist) createNumericObjects(*((V*) this));
 
      // set p and integrate distance
 
      vmDist->setP(p);
      if (r_open)
        return vmDistNI->integrateAOpenS(ra, rb);
      else
        return vmDistNI->integrateClosed(ra, rb);
    };
 
    //
    virtual double         integrateDistance(double p, double ra,
                                             bool r_open = false)
    {
      return integrateDistance(p, ra, vlRt, r_open);
    };
 
    virtual double         integrateTime(double p, double ra, double rb)
    {
      // create the numerical integration objects if they do not yet exist
 
      if (!vmTau) createNumericObjects(*((V*) this));
 
      // set p and integrate time
 
      vmTau->setP(p);
      if (ra == 0.0)
        return vmTauNI->integrateAOpenS(ra, rb);
      else
        return vmTauNI->integrateClosed(ra, rb);
    };
 
    virtual double         integrateTime(double p, double ra)
    {
      return integrateTime(p, ra, vlRt);
    };
 
    void                   setIntegTolerance(double tol)
    {
      // create the numerical integration objects if they do not yet exist
 
      if (!vmDist) createNumericObjects(*((V*) this));
 
      // set tolerance
 
      vmDistNI->setTolerance(tol);
      vmTauNI->setTolerance(tol);
    };
 
  protected:
 
    void                   createNumericObjects(V& v);
 
    static const double    vmIntegTol;
 
    TPdDistdr<V>*          vmDist;
 
    TPdTaudr<V>*           vmTau;
 
    util::IntegrateFunction<TPdDistdr<V> >* vmDistNI;
 
    util::IntegrateFunction<TPdTaudr<V> >*  vmTauNI;
};
 
// *****************************************************************************
//
//
// *****************************************************************************
class TAUP_EXP VelocityConst : public VelocityIntegrate<VelocityConst>
{
  public:
 
    // **** _PUBLIC LIFCYCLES_ *************************************************
 
    VelocityConst() : vc(0.0), VelocityIntegrate<VelocityConst>() {};
 
    VelocityConst(double c, double rt, double rb, const string& layrnam = "") :
                  vc(c), VelocityIntegrate<VelocityConst>(rt, rb, layrnam) {init();};
 
    VelocityConst(const VelocityConst& vcst) : vc(vcst.vc),
                  VelocityIntegrate<VelocityConst>(vcst) {};
 
    VelocityConst(DataBuffer& buffer) {deserialize(buffer);};
 
    virtual ~VelocityConst() {};
 
    // **** _PUBLIC OPERATORS_ *************************************************
 
    VelocityConst&         operator=(const VelocityConst& vcst)
    {
      VelocityIntegrate<VelocityConst>::operator=(vcst);
 
      vc = vcst.vc;
      return *this;
    };
 
    virtual double         operator()(double r) {return vc;};
 
    // **** _PUBLIC METHODS_ ***************************************************
 
    virtual double         rAtP(double p) {return p * vc;};
 
    //
    virtual double         integrateDistance(double p, double ra,
                                             bool r_open = false)
    {
      return integrateDistance(p, ra, vlRt, r_open);
    };
 
    //
    virtual double         integrateDistance(double p, double ra, double rb,
                                             bool r_open = false)
    {
      if (vc == 0.0) return 0.0;
      double pv = p * vc;
 
      double prt = pv / rb;
      double prb = pv / ra;
 
      return asin(min(prb, 1.0)) - asin(min(prt, 1.0));
    };
 
    //
    virtual double         integrateTime(double p, double ra)
    {
      return integrateTime(p, ra, vlRt);
    };
 
    //
    virtual double         integrateTime(double p, double ra, double rb)
    {
      if (vc == 0.0) return 0.0;
      double pt = rb / vc;
      double pr = ra / vc;
      return sqrt(fabs(pt - p) * (pt + p)) - sqrt(fabs(pr - p) * (pr + p));
    };
 
    virtual bool           isTimeIntegralTau() const {return false;};
 
    virtual void           writeVelocity(ostream& os) const;
 
    virtual void           toStream(ostream& os, string indent) const;
 
    static  string         class_name() {return "VelocityConst";};
 
    virtual string         get_class_name() const {return class_name();};
 
    virtual int            classSize() const
    {return (int) sizeof(VelocityConst);};
 
    virtual bool           isVelocityConstant() const {return true;};
 
    virtual int            bufferSize() const
    {
      return 2 * sizeof(double) + TPVelocityLayer::bufferSize();
    };
 
    virtual void           serialize(DataBuffer& buffer)
    {
      buffer.writeDouble(vc);
      TPVelocityLayer::serialize(buffer);
    };
 
    virtual void           deserialize(DataBuffer& buffer)
    {
      vc = buffer.readDouble();
      TPVelocityLayer::deserialize(buffer);
    };
 
  private:
 
    // **** _PRIVATE DATA_ *****************************************************
 
    double                 vc;
};
 
// *****************************************************************************
//
//
// *****************************************************************************
class TAUP_EXP VelocityPower : public VelocityIntegrate<VelocityPower>
{
  public:
 
    // **** _PUBLIC LIFCYCLES_ *************************************************
 
    VelocityPower() : vp0(0.0), vp1(0.0), vpB(0.0),
                      VelocityIntegrate<VelocityPower>() {};
 
    VelocityPower(double vt, double vb, double rt, double rb,
                  const string& layrnam = "") :
                  vp0(vt), vp1(vb),
                  VelocityIntegrate<VelocityPower>(rt, rb, layrnam)
    {
      vpB = log(vb / vt) / log(vlRb / vlRt);
      vp1_B = 1.0 - vpB;
      init();
    };
 
    VelocityPower(const VelocityPower& vp) :
                  vp0(vp.vp0), vp1(vp.vp1), vp1_B(vp.vp1_B),
                  vpB(vp.vpB), VelocityIntegrate<VelocityPower>(vp) {};
 
    VelocityPower(DataBuffer& buffer) {deserialize(buffer);};
 
    virtual ~VelocityPower() {};
 
    // **** _PUBLIC OPERATORS_ *************************************************
 
    VelocityPower& operator=(const VelocityPower& vp)
    {
      VelocityIntegrate<VelocityPower>::operator=(vp);
 
      vp0   = vp.vp0;
      vp1   = vp.vp1;
      vpB   = vp.vpB;
      vp1_B = vp.vp1_B;
 
      return *this;
    };
 
    virtual double         operator()(double r)
    {
      return vp0 * pow(r / vlRt, vpB);
    };
 
    // **** _PUBLIC METHODS_ ***************************************************
 
    virtual double         rAtP(double p)
    {
      return pow(p * vp0 * pow(1.0 / vlRt, vpB), 1.0 / vp1_B);
    };
 
    double                 B() const {return vpB;};
 
    //
    virtual double         integrateDistance(double p, double ra,
                                             bool r_open = false)
    {
      double pva = p * operator()(ra);
      double pvb = p * vlVt;
      return (asin(min(pva / ra, 1.0)) - asin(min(pvb / vlRt, 1.0))) / vp1_B;
    };
 
    //
    virtual double         integrateDistance(double p, double ra, double rb,
                                             bool r_open = false)
    {
      double pva = p * operator()(ra);
      double pvb = p * operator()(rb);
      return (asin(min(pva / ra, 1.0)) - asin(min(pvb / rb, 1.0))) / vp1_B;
    };
 
    //
    virtual double         integrateTime(double p, double ra)
    {
      double prb = vlRt / vlVt;
      double pra = ra / operator()(ra);
      return (sqrt(fabs(prb - p) * (prb + p)) -
              sqrt(fabs(pra - p) * (pra + p))) / vp1_B;
    };
 
    //
    virtual double         integrateTime(double p, double ra, double rb)
    {
      double pra = ra / operator()(ra);
      double prb = rb / operator()(rb);
      //double p2 = p * p;
      return (sqrt(fabs(pra - p) * (pra + p)) -
              sqrt(fabs(prb - p) * (prb + p))) / vp1_B;
    };
 
    virtual bool           isTimeIntegralTau() const {return false;};
 
    virtual void           writeVelocity(ostream& os) const;
 
    virtual void           toStream(ostream& os, string indent) const;
 
    static  string         class_name() {return "VelocityPower";};
 
    virtual string         get_class_name() const {return class_name();};
 
    virtual int            classSize() const
    {return (int) sizeof(VelocityPower);};
 
    virtual bool           isVelocityPowerLaw() const {return true;};
 
    virtual int            bufferSize() const
    {
      return 5 * sizeof(double) + TPVelocityLayer::bufferSize();
    };
 
    virtual void           serialize(DataBuffer& buffer)
    {
      buffer.writeDouble(vp0);
      buffer.writeRawDouble(vp1);
      buffer.writeRawDouble(vpB);
      buffer.writeRawDouble(vp1_B);
      TPVelocityLayer::serialize(buffer);
    };
 
    virtual void           deserialize(DataBuffer& buffer)
    {
      vp0 = buffer.readDouble();
      vp1 = buffer.readRawDouble();
      vpB = buffer.readRawDouble();
      vp1_B = buffer.readRawDouble();
      TPVelocityLayer::deserialize(buffer);
    };
 
  private:
 
    // **** _PRIVATE DATA_ *****************************************************
 
    double                 vp0;
 
    double                 vp1;
 
    double                 vpB;
 
    double                 vp1_B;
};
 
// *****************************************************************************
//
//
// *****************************************************************************
class TAUP_EXP VelocityLinear : public VelocityIntegrate<VelocityLinear>
{
  public:
 
    // **** _PUBLIC LIFCYCLES_ *************************************************
 
    VelocityLinear() : va0(0.0), va1(0.0), VelocityIntegrate<VelocityLinear>(){};
 
    VelocityLinear(double a0, double a1, double rt, double rb,
                   const string& layrnam = "", double normradius = 1.0) :
                   va0(a0), va1(a1), vNormRadius(normradius),
                   VelocityIntegrate<VelocityLinear>(rt, rb, layrnam) {init();};
 
    VelocityLinear(const VelocityLinear& vl) : va0(vl.va0), va1(vl.va1),
                   vNormRadius(vl.vNormRadius), VelocityIntegrate<VelocityLinear>(vl) {};
 
    VelocityLinear(DataBuffer& buffer) {deserialize(buffer);};
 
    virtual ~VelocityLinear() {};
 
    // **** _PUBLIC OPERATORS_ *************************************************
 
    VelocityLinear&        operator=(const VelocityLinear& vl)
    {
      VelocityIntegrate<VelocityLinear>::operator=(vl);
 
      va0 = vl.va0;
      va1 = vl.va1;
      vNormRadius = vl.vNormRadius;
 
      return *this;
    };
 
    virtual double         operator()(double r)
    {
      double rn = r / vNormRadius;
      return va1 * rn + va0;
    };
 
    // **** _PUBLIC METHODS_ ***************************************************
 
    virtual double         rAtP(double p)
    {
      return (p * va0) / (1 - p * va1 / vNormRadius);
    };
 
    //
    virtual double         integrateDistance(double p, double r,
                                             bool r_open = false)
    {
      double v1, pv, pv1, pv0, rslt, a, b, c;
 
      // first set result to analytic solution from the constant coefficient
      // component
 
      v1 = va1 / vNormRadius;
      pv = p * (va0 + v1 * r) / r;
      if (pv >= 1.0)
        rslt = PI / 2.0;
      else
        rslt = asin(pv);
      pv = p * vlVt / vlRt;
      rslt -= asin(pv);
 
      // next add to the result the analytic solution from the gradient
      // component ... note: that the additive component depends on the
      // value of c
 
      pv0 = p * va0;
      pv1 = p * v1;
      pv = pv1 * pv1;
      c = 1.0 - pv;
      if (c < 0.0)
      {
        b = - pv0 * pv1;
        double arg = (c * r + b) / pv0;
        if (arg >= 1.0)
          rslt += pv1 * (PI / 2.0  - asin((c * vlRt + b) / pv0)) / sqrt(-c);
        else
          rslt += pv1 * (asin(arg) - asin((c * vlRt + b) / pv0)) / sqrt(-c);
      }
      else if (c == 0.0)
      {
        rslt += sqrt(-2.0 * v1 * r / va0 - 1.0) -
                sqrt(-2.0 * v1 * vlRt / va0 - 1.0);
      }
      else // (c > 0.0)
      {
        b = - 2.0 * pv0 * pv1;
        a = - pv0 * pv0;
 
        double sc  = sqrt(c);
        double Rrt = a + vlRt * (b + c * vlRt);
        double Rr  = a + r * (b + c * r);
        rslt += pv1 * (log(2.0 * sc * sqrt(fabs(Rrt)) + 2.0 * c * vlRt + b) -
                       log(2.0 * sc * sqrt(fabs(Rr))  + 2.0 * c * r + b)) / sc;
      }
 
      // return the result
 
      return rslt;
    };
 
    //
    virtual double         integrateDistance(double p, double ra, double rb,
                                             bool r_open = false)
    {
      double v1, pv, pv1, pv0, rslt, a, b, c;
 
      // first set result to analytic solution from the constant coefficient
      // component
 
      v1 = va1 / vNormRadius;
      pv = p * (va0 + v1 * ra) / ra;
      if (pv >= 1.0)
        rslt = PI / 2.0;
      else
        rslt = asin(pv);
 
      pv = p * (va0 + v1 * rb) / rb;
      if (pv >= 1.0)
        rslt -= PI / 2.0;
      else
        rslt -= asin(pv);
 
      // next add to the result the analytic solution from the gradient
      // component ... note: that the additive component depends on the
      // value of c
 
      pv0 = p * va0;
      pv1 = p * v1;
      pv = pv1 * pv1;
      c = 1.0 - pv;
      if (c < 0.0)
      {
        b = - pv0 * pv1;
        double arg = (c * ra + b) / pv0;
        if (arg >= 1.0)
          rslt += pv1 * (PI / 2.0  - asin((c * rb + b) / pv0)) / sqrt(-c);
        else
          rslt += pv1 * (asin(arg) - asin((c * rb + b) / pv0)) / sqrt(-c);
      }
      else if (c == 0.0)
      {
        rslt += sqrt(-2.0 * v1 * ra / va0 - 1.0) -
                sqrt(-2.0 * v1 * rb / va0 - 1.0);
      }
      else // (c > 0.0)
      {
        b = - 2.0 * pv0 * pv1;
        a = - pv0 * pv0;
 
        double sc  = sqrt(c);
        double Rrb = a + rb * (b + c * rb);
        double Rra = a + ra * (b + c * ra);
        rslt += pv1 * (log(2.0 * sc * sqrt(fabs(Rrb)) + 2.0 * c * rb + b) -
                       log(2.0 * sc * sqrt(fabs(Rra)) + 2.0 * c * ra + b)) / sc;
      }
 
      // return the result
 
      return rslt;
    };
 
    virtual void           writeNormRadius(ostream& os) const;
 
    virtual void           writeVelocity(ostream& os) const;
 
    virtual void           toStream(ostream& os, string indent) const;
 
    static  string         class_name() {return "VelocityLinear";};
 
    virtual string         get_class_name() const {return class_name();};
 
    virtual int            classSize() const
    {return (int) sizeof(VelocityLinear);};
 
    virtual bool           isVelocityLinear() const {return true;};
 
    virtual int            bufferSize() const
    {
      return 4 * sizeof(double) + TPVelocityLayer::bufferSize();
    };
 
    virtual void           serialize(DataBuffer& buffer)
    {
      buffer.writeDouble(va0);
      buffer.writeRawDouble(va1);
      buffer.writeRawDouble(vNormRadius);
      TPVelocityLayer::serialize(buffer);
    };
 
    virtual void           deserialize(DataBuffer& buffer)
    {
      va0 = buffer.readDouble();
      va1 = buffer.readRawDouble();
      vNormRadius = buffer.readRawDouble();
      TPVelocityLayer::deserialize(buffer);
    };
  private:
 
    // **** _PRIVATE DATA_ *****************************************************
 
    double                 vNormRadius;
 
    double                 va0;
 
    double                 va1;
};
 
// *****************************************************************************
//
//
// *****************************************************************************
class TAUP_EXP VelocityQuadratic : public VelocityIntegrate<VelocityQuadratic>
{
  public:
 
    // **** _PUBLIC LIFCYCLES_ *************************************************
 
    VelocityQuadratic() : va0(0.0), va1(0.0), va2(0.0),
                          VelocityIntegrate<VelocityQuadratic>() {};
 
    VelocityQuadratic(double a0, double a1, double a2, double rt, double rb,
                      const string& layrnam = "", double normradius = 1.0) :
                      va0(a0), va1(a1), va2(a2), vNormRadius(normradius),
                      VelocityIntegrate<VelocityQuadratic>(rt, rb, layrnam) {init();};
 
    VelocityQuadratic(const VelocityQuadratic& vq) :
                      va0(vq.va0), va1(vq.va1), va2(vq.va2),
                      vNormRadius(vq.vNormRadius), VelocityIntegrate<VelocityQuadratic>(vq) {};
 
    VelocityQuadratic(DataBuffer& buffer) {deserialize(buffer);};
 
    virtual ~VelocityQuadratic() {};
 
    // **** _PUBLIC OPERATORS_ *************************************************
 
    VelocityQuadratic&     operator=(const VelocityQuadratic& vq)
    {
      VelocityIntegrate<VelocityQuadratic>::operator=(vq);
 
      va0 = vq.va0;
      va1 = vq.va1;
      va2 = vq.va2;
      vNormRadius = vq.vNormRadius;
 
      return *this;
    };
 
    virtual double         operator()(double r)
    {
      double rn = r / vNormRadius;
      return va0 + rn * (va1 + rn * va2);
    };
 
    // **** _PUBLIC METHODS_ ***************************************************
 
    virtual double         rAtP(double p)
    {
      // return 0 if p = 0
 
      if (p == 0.0) return 0.0;
 
      // calculate quadratic coefficient (a +- sqrt(a*a - c)) / b
      // calculate quadratic solution from ar^2 + br + c = 0 where
 
      double a = p * va2 / vNormRadius / vNormRadius;
      double b = p * va1 / vNormRadius - 1.0;
      double c = p * va0;
 
      // solution is r = (-b +- sqrt(b^2 - 4ac)) / 2a
      // first check for error
 
      double sqarg = b * b - 4.0 * a * c;
      if (sqarg < 0.0)
      {
        // error: p is not a valid ray parameter to calculate r
        return 0.0;
      }
      else if (sqarg == 0.0)
      {
        // one root ... return result
 
        return -b / 2.0 / a;
      }
      else
      {
        // two roots ... get both parts
 
        double cc = 2.0 * a;
        double aa = -b / cc;
        double sq = sqrt(sqarg) / cc;
        double r1 = aa - sq;
        if ((r1 <= vlRt) && (r1 >= vlRb))
          return r1;
        else
          return aa + sq;
      }
    };
 
    virtual void           writeNormRadius(ostream& os) const;
 
    virtual void           writeVelocity(ostream& os) const;
 
    virtual void           toStream(ostream& os, string indent) const;
 
    static  string         class_name() {return "VelocityQuadratic";};
 
    virtual string         get_class_name() const {return class_name();};
 
    virtual int            classSize() const
    {return (int) sizeof(VelocityQuadratic);};
 
    virtual bool           isVelocityQuadratic() const {return true;};
 
    virtual int            bufferSize() const
    {
      return 5 * sizeof(double) + TPVelocityLayer::bufferSize();
    };
 
    virtual void           serialize(DataBuffer& buffer)
    {
      buffer.writeDouble(va0);
      buffer.writeRawDouble(va1);
      buffer.writeRawDouble(va2);
      buffer.writeRawDouble(vNormRadius);
      TPVelocityLayer::serialize(buffer);
    };
 
    virtual void           deserialize(DataBuffer& buffer)
    {
      va0 = buffer.readDouble();
      va1 = buffer.readRawDouble();
      va2 = buffer.readRawDouble();
      vNormRadius = buffer.readRawDouble();
      TPVelocityLayer::deserialize(buffer);
    };
 
  private:
 
    // **** _PRIVATE DATA_ *****************************************************
 
    double                 vNormRadius;
 
    double                 va0;
 
    double                 va1;
 
    double                 va2;
};
 
// *****************************************************************************
//
//
// *****************************************************************************
class TAUP_EXP VelocityCubic : public VelocityIntegrate<VelocityCubic>
{
  public:
 
    // **** _PUBLIC LIFCYCLES_ *************************************************
 
    VelocityCubic() : va0(0.0), va1(0.0), va2(0.0), va3(0.0), vNormRadius(1.0),
                      VelocityIntegrate<VelocityCubic>() {};
 
    VelocityCubic(double a0, double a1, double a2, double a3,
                  double rt, double rb, const string& layrnam = "",
                  double normradius = 1.0) :
                  va0(a0), va1(a1), va2(a2), va3(a3), vNormRadius(normradius),
                  VelocityIntegrate<VelocityCubic>(rt, rb, layrnam) {init();};
 
    VelocityCubic(const VelocityCubic& vc) :
                  va0(vc.va0), va1(vc.va1), va2(vc.va2), va3(vc.va3),
                  vNormRadius(vc.vNormRadius), VelocityIntegrate<VelocityCubic>(vc) {};
 
    VelocityCubic(DataBuffer& buffer) {deserialize(buffer);};
 
    virtual ~VelocityCubic() {};
 
    // **** _PUBLIC OPERATORS_ *************************************************
 
    VelocityCubic&         operator=(const VelocityCubic& vc)
    {
      VelocityIntegrate<VelocityCubic>::operator=(vc);
 
      va0 = vc.va0;
      va1 = vc.va1;
      va2 = vc.va2;
      va3 = vc.va3;
      vNormRadius = vc.vNormRadius;
 
      return *this;
    };
 
    virtual double         operator()(double r)
    {
      double rn = r / vNormRadius;
      return va0 + rn * (va1 + rn * (va2 + rn * va3));
    };
 
    // **** _PUBLIC METHODS_ ***************************************************
 
    virtual double         rAtP(double p);
 
    virtual void           writeNormRadius(ostream& os) const;
 
    virtual void           writeVelocity(ostream& os) const;
 
    virtual void           toStream(ostream& os, string indent) const;
 
    static  string         class_name() {return "VelocityCubic";};
 
    virtual string         get_class_name() const {return class_name();};
 
    virtual int            classSize() const
    {return (int) sizeof(VelocityCubic);};
 
    virtual bool           isVelocityCubic() const {return true;};
 
    virtual int            bufferSize() const
    {
      return 6 * sizeof(double) + TPVelocityLayer::bufferSize();
    };
 
    virtual void           serialize(DataBuffer& buffer)
    {
      buffer.writeDouble(va0);
      buffer.writeRawDouble(va1);
      buffer.writeRawDouble(va2);
      buffer.writeRawDouble(va3);
      buffer.writeRawDouble(vNormRadius);
      TPVelocityLayer::serialize(buffer);
    };
 
    virtual void           deserialize(DataBuffer& buffer)
    {
      va0 = buffer.readDouble();
      va1 = buffer.readRawDouble();
      va2 = buffer.readRawDouble();
      va3 = buffer.readRawDouble();
      vNormRadius = buffer.readRawDouble();
      TPVelocityLayer::deserialize(buffer);
    };
 
  private:
 
    // **** _PRIVATE DATA_ *****************************************************
 
    double                 vNormRadius;
 
    double                 va0;
 
    double                 va1;
 
    double                 va2;
 
    double                 va3;
};
 
// *****************************************************************************
//
//
// *****************************************************************************
template<class V>
class VZero
{
  public:
 
    // **** _PUBLIC LIFCYCLES_ *************************************************
 
    //
    VZero(double p, V& v) : vzP(p), vzV(v) {};
 
    VZero(const VZero& vz) : vzP(vz.vzP), vzV(vz.vzV) {};
 
    virtual ~VZero() {};
 
    // **** _PUBLIC OPERATORS_ *************************************************
 
    VZero&                 operator=(const VZero& vz)
    {
      vzP = vz.vzP;
      return *this;
    };
 
    double                 operator()(double r) {return r - vzP * vzV(r);};
 
  private:
 
    double                 vzP;
 
    V&                     vzV;
};
 
} // end namespace taup
 
#endif // TPVelocityModels.h