RSTT 3.2.1 Regional Seismic Travel Time

Loading...
Searching...
No Matches
Location.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 Location_H
39#define Location_H
40
41// **** _SYSTEM INCLUDES_ ******************************************************
42#include <vector>
43#include <cmath>
44#include <sstream>
45#include <iomanip>
46
47#include "SLBMGlobals.h"
48
49#include "CPPUtils.h"
50
51 using namespace std;
52
53 using namespace geotess;
54
55// **** _BEGIN SLBM NAMESPACE_ **************************************************
56
57 namespace slbm {
58
59 class SLBM_EXP_IMP Location
60// **** _CLASS DEFINITION_ *****************************************************
61//
64//
76//
77// *****************************************************************************
78{
79 public:
80
81 // Default constructor
82 Location();
83
84 //destructor
85 virtual ~Location();
86
87 // copy constructor
88 Location(const Location& location);
89
90 // constructor specifying latitude, longitude (in radians), and depth in km
91 // Depth defaults to zero if not specified.
92 Location(const double &lat, const double &lon, const double &depth=0);
93
94 // constructor specifying 3D unit vector with origin at center of earth.
95 Location(const double v[], const double &radius);
96
101 Location(const Location& loc1, const Location& loc2);
102
103 // equal operator
104 Location& operator=(const Location& other);
105
106 // equality operator
107 bool operator==(const Location& other) const;
108
113 bool operator!=(const Location& other) {return ! (*this == other);};
114
118 string toString() const;
119
126 void setLocation(const double& lat, const double& lon, const double& depth);
127
133 void setLocation(const double* u, const double& r)
134 { v[0]=u[0]; v[1]=u[1]; v[2]=u[2]; radius=r; }
135
140 void setRadius(const double &r);
141
146 double getRadius() const;
147
153 double getEarthRadius() const;
154
159 double getDepth() const;
160
165 void setDepth(const double &depth);
166
173 double distance(const Location& other) const;
174
175 // returns the distance in km between this Location and some other Location,
176 // measured at the surface of the ellipsoid.
177 // The product of distance * getEarthRadius() is integrated along the great
178 // circle path from source to receiver (in approximately 1 degree increments).
179 double distanceKm(Location& other) const;
180
186 double distanceDegrees(const Location& other) const;
187
195 double azimuth(const Location& other) const;
196
204 double azimuthDegrees(const Location& other) const;
205
221 double azimuth(const Location& other, const double& errorValue) const;
222
238 double azimuthDegrees(const Location& other, const double& errorValue) const;
239
244 double getLat() const;
245
250 double getGeocentricLat() const
251 { return asin(v[2]); }
252
257 double getGeocentricLatDegrees() const
258 { return asin(v[2]) * RAD_TO_DEG; }
259
266 double getLon() const;
267
272 double getLatDegrees() const;
273
280 double getLonDegrees() const;
281
282 const double* getUnitVector() { return v; };
283
284 void getUnitVector(double x[3]) {x[0] = v[0]; x[1] = v[1]; x[2] = v[2];};
285
286 void setUnitVector(double x[3]) {v[0] = x[0]; v[1] = x[1]; v[2] = x[2];};
287
299 void move(const double &azimuth, const double &distance, Location &loc) const;
300
314 bool cross(const Location& x, Location& loc) const;
315
325 void rotate(Location& pole, double angle, Location& loc) const;
326
336 bool vectorTripleProduct(const Location& other, double vtp[]) const;
337
350 void move(const double vtp[], const double& a, Location& loc) const;
351
362 void move_north(const double x[], const double &distance, double z[]) const;
363
372 void move_north(const double &distance, Location &loc) const;
373
382 void rotate (const double x[], const double p[], const double &a,
383 double z[]) const;
384
385
392 double scalarTripleProduct(const Location& loc1, const Location& loc2) const;
393
400 double scalarTripleProduct(const double u[], const double w[]) const;
401
407 double scalarTripleProduct(const double u[], const double v[],
408 const double w[]) const;
409
416 bool vectorTripleProduct(const double u[], const double v[], double vtp[]) const;
417
425 bool vectorTripleProductNorthPole(const double u[], double w[]) const;
426
432 void move(const double v[], const double vtp[], const double& a,
433 double u[]) const;
434
442 double angle(const double u[], const double v[]) const;
443
450 double dot(const Location& other) const;
451
458 double dot(const double u[], const double v[]) const;
459
468 double cross(const double u[], const double v[], double w[]) const;
469
477 double crossNorth(const double u[], double w[]) const;
478
485 double normalize(double v[]) const;
486
491 double length(const double v[]) const;
492
502 static double EARTH_RADIUS;
503
504 static int getClassCount();
505
506 protected:
507
508 static int locationClassCount;
509
516 double v[3];
517
521 double radius;
522
523};
524
525// INLINE FUNCTIONS
526
527 inline void Location::setLocation(const double& latitude,
528 const double& longitude, const double& depth)
529{
530 // set x to the geocentric latitude in radians
531 double x = atan((1.0 - EARTH_E) * tan(latitude));
532 // z component of v is sin of geocentric latitude.
533 v[2] = sin(x);
534 // reuse x by setting it to cos of geocentric latitude
535 x = cos(x);
536 // compute x and y components of v
537 v[0] = x * cos(longitude);
538 v[1] = x * sin(longitude);
539 radius = getEarthRadius() - depth;
540}
541
542 inline double Location::getLat() const
543{
544 // Reference: Snyder, eq. 3-28, p. 17.
545 return atan(tan(asin(v[2]))/(1-EARTH_E));
546}
547
548 inline double Location::getLon() const
549{
550 return atan2(v[1],v[0]);
551}
552
553 inline double Location::getLatDegrees() const
554{
555 return RAD_TO_DEG * getLat();
556}
557
558 inline double Location::getLonDegrees() const
559{
560 return RAD_TO_DEG * getLon();
561}
562
563 inline double Location::getRadius() const
564{
565 return radius;
566}
567
568 inline void Location::setRadius(const double &r)
569{
570 radius = r;
571}
572
573 inline double Location::getDepth() const
574{
575 return getEarthRadius() - radius;
576}
577
578 inline void Location::setDepth(const double &depth)
579{
580 radius = getEarthRadius() - depth;
581}
582
583 inline double Location::distance(const Location& other) const
584{
585 return angle(v, other.v);
586}
587
588 inline double Location::distanceDegrees(const Location& other) const
589{
590 return angle(v, other.v) * RAD_TO_DEG;
591}
592
593 inline double Location::azimuth(const Location& other) const
594{
595 return azimuth(other, 0.);
596}
597
598 inline double Location::azimuth(const Location& other, const double& errorValue) const
599{
600 double azim=errorValue;
601 double c2[3];
602
603 // set c2 = the cross product of this x other.
604 if (cross(v, other.v, c2) > 0.)
605 // if the cross product has zero length then the two vectors are coincident
606 // (do nothing in that case; returns zero).
607 {
608 double c[3];
609 // set c = this x north pole
610 // if the cross product has zero length then this == north_pole
611 // or south pole and azimuth is indeterminant.
612 if (crossNorth(v, c) > 0.)
613 {
614 // set azimuth to the angle between (this x north pole) and
615 // (this x other).
616 azim = angle(c, c2);
617 // if the dot product of (this x other) . northPole < 0
618 if (c2[2] < 0.) azim = 2.0 * PI - azim;
619 }
620 }
621 return azim;
622}
623
624 inline double Location::azimuthDegrees(const Location& other) const
625{
626 return azimuth(other) * RAD_TO_DEG;
627}
628
629 inline double Location::azimuthDegrees(const Location& other, const double& errorValue) const
630{
631 double az = azimuth(other, errorValue);
632 if (az != errorValue)
633 az *= RAD_TO_DEG;
634 return az;
635}
636
637// find a location a specified distance and azimuth away from this.
638 inline void Location::move(const double &azimuth, const double &distance,
639 Location &loc) const
640{
641 //find a point distance radians due north of v
642 double x[3];
643 move_north(v, distance, x);
644 //x is now a point distance due north of v
645 //rotate x around v in direction of positive longitude
646 rotate(x, v, -azimuth, loc.v);
647 loc.radius=radius;
648}
649
650// given a vector x, find a vector z which is a specified distance north.
651 inline void Location::move_north(const double x[], const double &distance,
652 double z[]) const
653{
654 double c[3];
655 vectorTripleProductNorthPole(x, c);
656 move(x, c, distance, z);
657}
658
659// find a vector z which is a specified distance north of this.
660 inline void Location::move_north(const double &distance, Location &loc) const
661{
662 double vtp[3];
663 vectorTripleProductNorthPole(v, vtp);
664 move(v, vtp, distance, loc.v);
665}
666
667// given a precomputed normalized vector triple product, vtp, move this
668// a specified distance in the direction specified by vtp and return the
669// result in Location loc
670 inline void Location::move(const double vtp[], const double& a, Location& loc) const
671{
672 move(v, vtp, a, loc.v);
673 loc.radius = radius;
674}
675
676// given a vector w, and a precomputed normalized vector triple product, vtp,
677// move w a specified distance in the direction specified by vtp and return the
678// result in vector u.
679 inline void Location::move(const double w[], const double vtp[],
680 const double& a, double u[]) const
681{
682 double cosa = cos(a);
683 double sina = sin(a);
684 u[0] = cosa * w[0] + sina * vtp[0];
685 u[1] = cosa * w[1] + sina * vtp[1];
686 u[2] = cosa * w[2] + sina * vtp[2];
687}
688
689// find the normalized vector triple product using Location objects.
690// Compute vtp = (this x other) x this and return vtp.
691 inline bool Location::vectorTripleProduct(const Location& other, double vtp[]) const
692{
693 return vectorTripleProduct(v, other.v, vtp);
694}
695
696// find the normalized vector triple product using vectors.
697// Compute vtp = (u x v) x u and return vtp.
698 inline bool Location::vectorTripleProduct(const double u[], const double s[],
699 double w[]) const
700{
701 // set q = u cross s
702 double q0, q1, q2;
703 q0 = u[1] * s[2] - u[2] * s[1];
704 q1 = u[2] * s[0] - u[0] * s[2];
705 q2 = u[0] * s[1] - u[1] * s[0];
706 // set w = q cross u
707 w[0] = q1 * u[2] - q2 * u[1];
708 w[1] = q2 * u[0] - q0 * u[2];
709 w[2] = q0 * u[1] - q1 * u[0];
710 // normalize w to unit length. if length
711 // of w before normalization was zero, then
712 // w will = {0,0,0} and return false;
713 return normalize(w) != 0.;
714}
715
716// find the normalized vector triple product of vector u with north pole
717// Compute w = (u x n) x u and return w.
718 inline bool Location::vectorTripleProductNorthPole(const double u[],
719 double w[]) const
720{
721 w[0] = -u[0] * u[2];
722 w[1] = -u[1] * u[2];
723 w[2] = u[1] * u[1] + u[0] * u[0];
724 return normalize(w) != 0.;
725}
726
727 inline double Location::scalarTripleProduct(const Location& loc1, const Location& loc2) const
728{
729 return scalarTripleProduct(loc1.v, loc2.v, v);
730}
731
732 inline double Location::scalarTripleProduct(const double u[], const double p[],
733 const double w[]) const
734{
735 return u[0] * p[1] * w[2] + p[0] * w[1] * u[2] + w[0] * u[1] * p[2] -
736 w[0] * p[1] * u[2] - u[0] * w[1] * p[2] - p[0] * u[1] * w[2];
737}
738
739 inline double Location::scalarTripleProduct(const double v1[],
740 const double v2[]) const
741{
742 return scalarTripleProduct(v1, v2, v);
743}
744
745// rotate this Location around Location p by angle a and return result
746// in Location loc.
747 inline void Location::rotate(Location& p, double a, Location &loc) const
748{
749 rotate(v, p.v, a, loc.v);
750 loc.radius=radius;
751}
752
753// rotate vector x around vector p by angle a and return the result in vector z.
754 inline void Location::rotate (const double x[], const double p[],
755 const double &a, double z[]) const
756{
757 double d = x[0] * p[0] + x[1] * p[1] + x[2] * p[2]; // dot product
758 // if x and p are parallel, x needs no rotation.
759 if (d <= -1. || d >= 1.0)
760 {
761 z[0] = x[0]; z[1] = x[1]; z[2] = x[2];
762 }
763 else
764 {
765 double cosa = cos(a);
766 double sina = sin(a);
767 d *= (1-cosa);
768 z[0] = cosa * x[0] + d * p[0] + sina * (p[1] * x[2] - p[2] * x[1]);
769 z[1] = cosa * x[1] + d * p[1] + sina * (p[2] * x[0] - p[0] * x[2]);
770 z[2] = cosa * x[2] + d * p[2] + sina * (p[0] * x[1] - p[1] * x[0]);
771 double len = sqrt(z[0] * z[0] + z[1] * z[1] + z[2] * z[2]);
772 z[0] /= len;
773 z[1] /= len;
774 z[2] /= len;
775 }
776}
777
778// find angular separation of two unit vectors, in radians.
779 inline double Location::angle(const double u[], const double w[]) const
780{
781 return acos(max(min(dot(u, w), 1.0), -1.0));
782}
783
784 inline double Location::dot(const Location& other) const
785{
786 return v[0] * other.v[0] + v[1] * other.v[1] + v[2] * other.v[2];
787}
788
789// dot product
790 inline double Location::dot(const double u[], const double s[]) const
791{
792 return u[0] * s[0] + u[1] * s[1] + u[2] * s[2];
793}
794
795// cross product of this Location with other Location, result stored in
796// Location loc.
797// return true if length > zero.
798 inline bool Location::cross(const Location& x, Location &loc) const
799{
800 if (cross(v, x.v, loc.v) > 0.)
801 loc.radius = radius;
802 else
803 loc.radius = -1.;
804 return loc.radius > 0.;
805}
806
807// cross product of two unit vectors, normalized to unit length.
808// return the length of cross product before normalization.
809 inline double Location::cross(const double u[], const double s[],
810 double w[]) const
811{
812 w[0] = u[1] * s[2] - u[2] * s[1];
813 w[1] = u[2] * s[0] - u[0] * s[2];
814 w[2] = u[0] * s[1] - u[1] * s[0];
815 return normalize(w);
816}
817
818// cross product of vector u with north pole, normalized to unit length.
819// result returned in w. returns length prior to normalization.
820 inline double Location::crossNorth(const double u[], double w[]) const
821{
822 double len = u[0] * u[0] + u[1] * u[1];
823 if (len <= 0.)
824 {
825 len = 0.;
826 w[0] = 0.;
827 w[1] = 0.;
828 w[2] = 0.;
829 }
830 else
831 {
832 len = sqrt(len);
833 w[0] = u[1]/len;
834 w[1] = -u[0]/len;
835 w[2] = 0.;
836 }
837 return len;
838}
839
840// normalize vector to unit length. returns length prior to
841// normalization.
842 inline double Location::normalize(double u[]) const
843{
844 double len = length(u);
845 if (len == 0.)
846 {
847 u[0] = u[1] = u[2] = 0.0;
848 }
849 else
850 {
851 u[0] /= len;
852 u[1] /= len;
853 u[2] /= len;
854 }
855 return len;
856}
857
858// find the length of a vector.
859 inline double Location::length(const double u[]) const
860{
861 double l = dot(u, u);
862 if (l > 0.) return sqrt(l);
863 return 0.;
864}
865
866 inline string Location::toString() const
867{
868 ostringstream os;
869 os << std::fixed
870 << std::showpoint
871 << std::setprecision(4)
872 << std::setw(9) << getLatDegrees()
873 << " "
874 << std::setw(10) << getLonDegrees()
875 << " "
876 << std::setprecision(3)
877 << std::setw(10) << getDepth();
878 return os.str();
879}
880
881} // end slbm namespace
882
883#endif // Location_H
SLBMGlobals.h
SLBM_EXP_IMP
#define SLBM_EXP_IMP
Definition SLBMGlobals.h:180
slbm::Location
The Location Class manages a single point in/on the Earth, which is described by the GRS80 ellipsoid.
Definition Location.h:78
slbm::Location::getLon
double getLon() const
Retrieve the longitude of this Location. Value will be between -PI and PI radians.
Definition Location.h:548
slbm::Location::distanceDegrees
double distanceDegrees(const Location &other) const
Definition Location.h:588
slbm::Location::getGeocentricLatDegrees
double getGeocentricLatDegrees() const
Retrieve the geocentric latitude of this Location, degrees.
Definition Location.h:257
slbm::Location::setLocation
void setLocation(const double &lat, const double &lon, const double &depth)
Definition Location.h:527
slbm::Location::rotate
void rotate(Location &pole, double angle, Location &loc) const
Definition Location.h:747
slbm::Location::distance
double distance(const Location &other) const
Definition Location.h:583
slbm::Location::cross
bool cross(const Location &x, Location &loc) const
Definition Location.h:798
slbm::Location::getLonDegrees
double getLonDegrees() const
Retrieve the longitude of this Location. Value will be between -180 and 180 degrees.
Definition Location.h:558
slbm::Location::vectorTripleProductNorthPole
bool vectorTripleProductNorthPole(const double u[], double w[]) const
Compute the normalized vector triple product (u x northPole) x u and store result in w....
Definition Location.h:718
slbm::Location::azimuth
double azimuth(const Location &other) const
Find the azimuth from this Location to some other Location. Result will be between 0 and 2*PI radians...
Definition Location.h:593
slbm::Location::scalarTripleProduct
double scalarTripleProduct(const Location &loc1, const Location &loc2) const
Definition Location.h:727
slbm::Location::getUnitVector
const double * getUnitVector()
Definition Location.h:282
slbm::Location::crossNorth
double crossNorth(const double u[], double w[]) const
Definition Location.h:820
slbm::Location::Location
Location(const Location &location)
slbm::Location::Location
Location(const double v[], const double &radius)
slbm::Location::vectorTripleProduct
bool vectorTripleProduct(const Location &other, double vtp[]) const
Compute the vector triple product (this x other) x this, normalized to unit length....
Definition Location.h:691
slbm::Location::getClassCount
static int getClassCount()
slbm::Location::EARTH_RADIUS
static double EARTH_RADIUS
Definition Location.h:502
slbm::Location::getEarthRadius
double getEarthRadius() const
slbm::Location::move
void move(const double &azimuth, const double &distance, Location &loc) const
Retrieve a Location that is a specified distance away from this Location, in a specified direction.
Definition Location.h:638
slbm::Location::v
double v[3]
Definition Location.h:516
slbm::Location::azimuthDegrees
double azimuthDegrees(const Location &other) const
Find the azimuth from this Location to some other Location. Result will be between 0 and 360 degrees.
Definition Location.h:624
slbm::Location::~Location
virtual ~Location()
slbm::Location::setLocation
void setLocation(const double *u, const double &r)
Definition Location.h:133
slbm::Location::toString
string toString() const
Definition Location.h:866
slbm::Location::dot
double dot(const Location &other) const
Definition Location.h:784
slbm::Location::move_north
void move_north(const double x[], const double &distance, double z[]) const
Definition Location.h:651
slbm::Location::operator!=
bool operator!=(const Location &other)
Definition Location.h:113
slbm::Location::length
double length(const double v[]) const
Definition Location.h:859
slbm::Location::setDepth
void setDepth(const double &depth)
Definition Location.h:578
slbm::Location::locationClassCount
static int locationClassCount
Definition Location.h:508
slbm::Location::setUnitVector
void setUnitVector(double x[3])
Definition Location.h:286
slbm::Location::Location
Location()
slbm::Location::getDepth
double getDepth() const
Definition Location.h:573
slbm::Location::distanceKm
double distanceKm(Location &other) const
slbm::Location::getLatDegrees
double getLatDegrees() const
Retrieve the geographic latitude of this Location, degrees.
Definition Location.h:553
slbm::Location::operator==
bool operator==(const Location &other) const
slbm::Location::normalize
double normalize(double v[]) const
Definition Location.h:842
slbm::Location::Location
Location(const Location &loc1, const Location &loc2)
slbm::Location::getRadius
double getRadius() const
Definition Location.h:563
slbm::Location::operator=
Location & operator=(const Location &other)
slbm::Location::getGeocentricLat
double getGeocentricLat() const
Retrieve the geocentric latitude of this Location, radians.
Definition Location.h:250
slbm::Location::radius
double radius
Definition Location.h:521
slbm::Location::setRadius
void setRadius(const double &r)
Definition Location.h:568
slbm::Location::getLat
double getLat() const
Retrieve the geographic latitude of this Location, radians.
Definition Location.h:542
slbm::Location::Location
Location(const double &lat, const double &lon, const double &depth=0)
slbm::Location::getUnitVector
void getUnitVector(double x[3])
Definition Location.h:284
slbm::Location::angle
double angle(const double u[], const double v[]) const
Definition Location.h:779
slbm
Definition CrustalProfile.h:59