RSTT  3.2.0
Regional Seismic Travel Time
All Classes Namespaces Files Functions Variables Typedefs Friends Macros
TauPSiteFunctionals.h
Go to the documentation of this file.
1 //- ****************************************************************************
2 //-
3 //- Copyright 2009 National Technology & Engineering Solutions of Sandia, LLC
4 //- (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
5 //- Government retains certain rights in this software.
6 //-
7 //- BSD Open Source License
8 //- All rights reserved.
9 //-
10 //- Redistribution and use in source and binary forms, with or without
11 //- modification, are permitted provided that the following conditions are met:
12 //-
13 //- 1. Redistributions of source code must retain the above copyright notice,
14 //- this list of conditions and the following disclaimer.
15 //-
16 //- 2. Redistributions in binary form must reproduce the above copyright
17 //- notice, this list of conditions and the following disclaimer in the
18 //- documentation and/or other materials provided with the distribution.
19 //-
20 //- 3. Neither the name of the copyright holder nor the names of its
21 //- contributors may be used to endorse or promote products derived from
22 //- this software without specific prior written permission.
23 //-
24 //- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
25 //- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
26 //- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
27 //- ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
28 //- LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
29 //- CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
30 //- SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
31 //- INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
32 //- CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
33 //- ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
34 //- POSSIBILITY OF SUCH DAMAGE.
35 //-
36 //- ****************************************************************************
37 
38 #ifndef TAUPSITEFNCTNLS_H
39 #define TAUPSITEFNCTNLS_H
40 
41 // **** _SYSTEM INCLUDES_ ******************************************************
42 
43 #include <vector>
44 #include <map>
45 #include <set>
46 #include <cmath>
47 
48 using namespace std;
49 // use standard library objects
50 
51 // **** _LOCAL INCLUDES_ *******************************************************
52 
53 // **** _BEGIN TAUP NAMESPACE_ *************************************************
54 
55 namespace taup {
56 
57 // **** _FORWARD REFERENCES_ ***************************************************
58 
59 class TauPSite;
60 
61 // **** _CLASS CONSTANTS_ ******************************************************
62 
63 // *****************************************************************************
64 // **** SplitDistance Definition ***********************************************
65 // *****************************************************************************
66 
67 // *****************************************************************************
68 //
74 //
75 // *****************************************************************************
77 {
78  public:
79 
82  SplitDistance(TauPSite* tps) : sdTPS(tps) {};
83 
85  SplitDistance(const SplitDistance& sd) : sdTPS(sd.sdTPS) {};
86 
88  virtual ~SplitDistance() {};
89 
91  SplitDistance& operator=(const SplitDistance& sd) {return *this;};
92 
96  double operator()(double p);
97 
98  private:
99 
102  TauPSite* sdTPS;
103 
104 }; // SplitDistance End Definition
105 
106 // *****************************************************************************
107 // **** TPZeroFunctional Definition ********************************************
108 // *****************************************************************************
109 
110 // *****************************************************************************
111 //
115 //
148 //
149 // *****************************************************************************
151 {
152  public:
153 
154  // **** _PUBLIC LIFECYCLES_ ************************************************
155 
157  TPZeroFunctional() : tpzTPS(NULL), tpzRSrc(0.0), tpzRSrcSgn(1.0),
158  tpzRRcvr(0.0), tpzRRcvrSgn(1.0), tpzPLast(-1.0),
159  tpzDLast(-1.0), tpzD(0.0), tpzPT(0.0),
160  tpzRayLegDist(0.0), tpzRayLegTime(0.0),
161  tpzSrcLegDist(0.0), tpzSrcLegTime(0.0),
162  tpzRcvrLegDist(0.0), tpzRcvrLegTime(0.0),
163  tpzIsTurningZero(true), tpzIsRayLegValid(false),
164  tpzIsRcvrLegValid(false), tpzIsSrcLegValid(false),
165  tpzRadius(EARTH_RAD) {};
166 
170  tpzTPS(tps), tpzRSrc(0.0), tpzRSrcSgn(1.0),
171  tpzRRcvr(0.0), tpzRRcvrSgn(1.0), tpzPLast(-1.0),
172  tpzDLast(-1.0), tpzD(0.0), tpzPT(0.0),
173  tpzRayLegDist(0.0), tpzRayLegTime(0.0),
174  tpzSrcLegDist(0.0), tpzSrcLegTime(0.0),
175  tpzRcvrLegDist(0.0), tpzRcvrLegTime(0.0),
176  tpzIsTurningZero(true), tpzIsRayLegValid(false),
177  tpzIsRcvrLegValid(false), tpzIsSrcLegValid(false),
178  tpzRadius(EARTH_RAD) {};
179 
182  tpzTPS(tpzf.tpzTPS),
183  tpzRSrc(tpzf.tpzRSrc), tpzRSrcSgn(tpzf.tpzRSrcSgn),
184  tpzRRcvr(tpzf.tpzRRcvr), tpzRRcvrSgn(tpzf.tpzRRcvrSgn),
185  tpzPLast(tpzf.tpzPLast), tpzDLast(tpzf.tpzDLast),
186  tpzD(tpzf.tpzD), tpzPT(tpzf.tpzPT),
187  tpzRayLegDist(tpzf.tpzRayLegDist),
188  tpzRayLegTime(tpzf.tpzRayLegTime),
189  tpzSrcLegDist(tpzf.tpzSrcLegDist),
190  tpzSrcLegTime(tpzf.tpzSrcLegTime),
191  tpzRcvrLegDist(tpzf.tpzRcvrLegDist),
192  tpzRcvrLegTime(tpzf.tpzRcvrLegTime),
193  tpzIsTurningZero(tpzf.tpzIsTurningZero),
194  tpzIsRayLegValid(tpzf.tpzIsRayLegValid),
195  tpzIsRcvrLegValid(tpzf.tpzIsRcvrLegValid),
196  tpzIsSrcLegValid(tpzf.tpzIsSrcLegValid),
197  tpzRadius(tpzf.tpzRadius) {};
198 
200  virtual ~TPZeroFunctional() {};
201 
202  // **** _PUBLIC OPERATORS_ *************************************************
203 
206  {
207  // set all values from tpzf into this TPZeroFunctional
208 
209  tpzD = tpzf.tpzD;
210  tpzPT = tpzf.tpzPT;
211  tpzRSrc = tpzf.tpzRSrc;
212  tpzRSrcSgn = tpzf.tpzRSrcSgn;
213  tpzRRcvr = tpzf.tpzRRcvr;
214  tpzRRcvrSgn = tpzf.tpzRRcvrSgn;
215  tpzPLast = tpzf.tpzPLast;
216  tpzDLast = tpzf.tpzDLast;
217  tpzRayLegDist = tpzf.tpzRayLegDist;
218  tpzRayLegTime = tpzf.tpzRayLegTime;
219  tpzSrcLegDist = tpzf.tpzSrcLegDist;
220  tpzSrcLegTime = tpzf.tpzSrcLegTime;
221  tpzRcvrLegDist = tpzf.tpzRcvrLegDist;
222  tpzRcvrLegTime = tpzf.tpzRcvrLegTime;
223  tpzRadius = tpzf.tpzRadius;
224  tpzTPS = tpzf.tpzTPS;
225  tpzIsTurningZero = tpzf.tpzIsTurningZero;
226  tpzIsRayLegValid = tpzf.tpzIsRayLegValid;
227  tpzIsRcvrLegValid = tpzf.tpzIsRcvrLegValid;
228  tpzIsSrcLegValid = tpzf.tpzIsSrcLegValid;
229  tpzPLast = -1.0;
230 
231  // exit
232 
233  return *this;
234  };
235 
240  double operator()(double p)
241  {
242  // calculate new zero in condition if p is different than last
243  // call
244 
245  if (p != tpzPLast)
246  {
247  // save p and new zero in condition
248 
249  tpzPLast = p;
250  distance(p);
251  if (tpzIsTurningZero)
252  tpzDLast = getTurningZero();
253  else
254  tpzDLast = getUpGoingZero();
255  }
256 
257  // return zero in condition at p
258 
259  return tpzDLast;
260  };
261 
262  // **** _PUBLIC METHODS_ ***************************************************
263 
267  {
268  tpzIsTurningZero = true;
269  };
270 
274  {
275  return tpzIsTurningZero;
276  };
277 
280  double getTurningZero()
281  {
282  return (tpzD - (2.0 * tpzRayLegDist - tpzRSrcSgn * tpzSrcLegDist -
283  tpzRRcvrSgn * tpzRcvrLegDist));
284  };
285 
289  {
290  return (tpzIsRayLegValid && tpzIsRcvrLegValid && tpzIsSrcLegValid);
291  };
292 
296  {
297  tpzIsTurningZero = false;
298  };
299 
303  {
304  return !tpzIsTurningZero;
305  };
306 
309  double getUpGoingZero()
310  {
311  //return (tpzD - fabs(tpzRSrcSgn * tpzSrcLegDist -
312  // tpzRRcvrSgn * tpzRcvrLegDist));
313  return (tpzD - tpzSrcLegDist);
314  };
315 
319  {
320  return (tpzIsRcvrLegValid && tpzIsSrcLegValid);
321  };
322 
325  double getMinP();
326 
329  {
330  tpzTPS = tps;
331  tpzPLast = -1.0;
332  };
333 
336  {
337  return *tpzTPS;
338  };
339 
341  void setPlanetRadius(double pr)
342  {
343  tpzRadius = pr;
344  tpzPLast = -1.0;
345  };
346 
348  double getPlanetRadius() const
349  {
350  return tpzRadius;
351  };
352 
354  void setSourceRadius(double r)
355  {
356  tpzRSrcSgn = setRadius(r);
357  tpzRSrc = r;
358  tpzPLast = -1.0;
359  };
360 
362  double getSourceRadius() const
363  {
364  return tpzRSrc;
365  }
366 
368  void setSourceDepth(double d)
369  {
370  setSourceRadius(tpzRadius - d);
371  };
372 
374  double getSourceDepth() const
375  {
376  return tpzRadius - tpzRSrc;
377  }
378 
380  void setReceiverRadius(double r)
381  {
382  tpzRRcvrSgn = setRadius(r);
383  tpzRRcvr = r;
384  tpzPLast = -1.0;
385  };
386 
388  double getReceiverRadius() const
389  {
390  return tpzRRcvr;
391  }
392 
394  void setReceiverDepth(double d)
395  {
396  setReceiverRadius(tpzRadius - d);
397  };
398 
400  double getReceiverDepth() const
401  {
402  return tpzRadius - tpzRRcvr;
403  }
404 
406  void setDist(double d)
407  {
408  tpzD = d;
409  tpzPLast = -1.0;
410  };
411 
413  double getDist() const
414  {
415  return tpzD;
416  };
417 
419  void setPTop(double p)
420  {
421  tpzPT = p;
422  tpzPLast = -1.0;
423  };
424 
426  double getPTop() const
427  {
428  return tpzPT;
429  };
430 
436  void distance(double p);
437 
443  double time(double p);
444 
446  double getRayDistance() const {return tpzRayLegDist;};
447 
449  double getRayTime() const {return tpzRayLegTime;};
450 
452  double getSourceLegDistance() const {return tpzRSrcSgn * tpzSrcLegDist;};
453 
455  double getSourceLegTime() const {return tpzRSrcSgn * tpzSrcLegTime;};
456 
458  double getReceiverLegDistance() const
459  {
460  return tpzRRcvrSgn * tpzRcvrLegDist;
461  };
462 
464  double getReceiverLegTime() const {return tpzRRcvrSgn * tpzRcvrLegTime;};
465 
467  int getSourceLayerId() const {return getRadiusLayerId(tpzRSrc);};
468 
470  int getReceiverLayerId() const {return getRadiusLayerId(tpzRRcvr);};
471 
473  int getRadiusLayerId(double r) const;
474 
475  private:
476 
477  // **** _PRIVATE METHODS_ **************************************************
478 
484  double setRadius(double& r)
485  {
486  double rsgn = 1.0;
487 
488  // if the input radius exceeds the maximum planet radius then
489  // this is an elevated source ... set r to depth and
490  // reverse the sign
491 
492  if (r > tpzRadius)
493  {
494  r = 2.0 * tpzRadius - r;
495  rsgn = -1.0;
496  }
497  return rsgn;
498  };
499 
500  // **** _PRIVATE DATA_ *****************************************************
501 
503  TauPSite* tpzTPS;
504 
506  double tpzRSrc;
507 
510  double tpzRSrcSgn;
511 
513  double tpzRRcvr;
514 
517  double tpzRRcvrSgn;
518 
520  double tpzPLast;
521 
523  double tpzDLast;
524 
526  double tpzD;
527 
530  double tpzPT;
531 
533  double tpzRayLegDist;
534 
536  double tpzRayLegTime;
537 
539  double tpzSrcLegDist;
540 
542  double tpzSrcLegTime;
543 
545  double tpzRcvrLegDist;
546 
548  double tpzRcvrLegTime;
549 
552  bool tpzIsTurningZero;
553 
556  bool tpzIsRayLegValid;
557 
560  bool tpzIsRcvrLegValid;
561 
564  bool tpzIsSrcLegValid;
565 
567  double tpzRadius;
568 
569 }; // end TPZeroFunctional Definition
570 
571 } // end namespace taup
572 
573 #endif // TAUPSITEFNCTNLS_H
Function object used by the Brents minF function to find the minimum (or maximum) of a retrograde lay...
SplitDistance(TauPSite *tps)
Standard constructor sets the internal TauPModel and a set of velocity layers.
double operator()(double p)
Function object operator() definition. This operator returns the ray distance as a function of the in...
SplitDistance & operator=(const SplitDistance &sd)
Assignment operator.
SplitDistance(const SplitDistance &sd)
Copy constructor.
virtual ~SplitDistance()
Destructor.
The primary layer search functional used by Brents zeroIn(...) function to find layers that contain a...
void setTauPSite(TauPSite *tps)
Sets the TauPSite.
TPZeroFunctional & operator=(const TPZeroFunctional &tpzf)
Assignment operator.
double getSourceDepth() const
Returns the source depth.
void setSourceRadius(double r)
Sets the source radius and sign.
void distance(double p)
The primary function of this object which calculates the ray travel distance between the source and t...
void setPTop(double p)
Sets the layer top ray parameter p for the current search layer.
void setDist(double d)
Sets the search distance between the source and receiver.
TPZeroFunctional(TauPSite *tps)
Standard constructor. Assigns the TauPSite for this TPZeroFunctional.
void setReceiverDepth(double d)
Sets the receiver radius and sign from the input depth.
double operator()(double p)
The function objects operator() definition which is used by a Brents::zeroIn(...) function to find th...
void setReceiverRadius(double r)
Sets the receiver radius and sign.
double getPTop() const
Gets the layer top ray parameter p for the current search layer.
double getReceiverLegTime() const
Returns the surface-to-receiver ray time.
double getRayDistance() const
Returns the surface-to-surface ray distance.
double getDist() const
Sets the search distance between the source and receiver.
double getTurningZero()
Returns the turning leg zero evaluation from the last distance() function evaluation.
double getReceiverRadius() const
Returns the receiver radius.
double getReceiverDepth() const
Returns the receiver depth.
TPZeroFunctional()
Default constructor.
bool isTurningZero()
Returns true if the turning zero is set for return by the operator() function (tpzIsTurningZero is tr...
double getSourceLegTime() const
Returns the surface-to-source ray time.
double getPlanetRadius() const
Returns the planet radius (default to Earth = 6371.0 km)
double getSourceRadius() const
Returns the source radius.
TPZeroFunctional(const TPZeroFunctional &tpzf)
Copy constructor.
double getSourceLegDistance() const
Returns the surface-to-source ray distance.
bool isUpGoingRayValid()
Returns true if the upgoing receiver and source legs are valid in the last distance() function call e...
double getReceiverLegDistance() const
Returns the surface-to-receiver ray distance.
bool isTurningRayValid()
Returns true if the turning ray, receiver, and source legs are valid in the last distance() function ...
int getReceiverLayerId() const
Return the layer containing the receiver position.
double getMinP()
Returns the minimum allowable ray parameter for a ray to transfer between the source and receiver dep...
void setUpGoingZero()
Sets the upgoing zero for the operator() function (tpzIsTurningZero is set to false).
TauPSite & getTauPSite()
Gets the TauPSite.
int getRadiusLayerId(double r) const
Return the layer containing the input radius r.
double getRayTime() const
Returns the surface-to-surface ray time.
double getUpGoingZero()
Returns the upgoing leg zero evaluation from the last distance() function evaluation.
bool isUpGoingZero()
Returns true if the upgoing zero is set for return by the operator() function (tpzIsTurningZero is fa...
void setSourceDepth(double d)
Sets the source radius and sign from the input depth.
virtual ~TPZeroFunctional()
Destructor.
void setPlanetRadius(double pr)
Sets the planet radius (default to Earth = 6371.0 km)
double time(double p)
Calculates the travel time between the source and receiver for the current layer at the zero in ray p...
int getSourceLayerId() const
Return the layer containing the source position.
void setTurningZero()
Sets the turning zero for the operator() function (tpzIsTurningZero is set to true).