LSSTApplications  19.0.0-11-g2ce9f25+2,20.0.0+1,20.0.0+11,20.0.0+2,20.0.0+3,20.0.0+4,20.0.0+5,20.0.0+8,20.0.0+9,20.0.0-1-g009f3de,20.0.0-1-g10df615+8,20.0.0-1-g253301a+4,20.0.0-1-g32a200e+8,20.0.0-1-g596936a+9,20.0.0-1-g8a53f90+1,20.0.0-1-gc96f8cb+10,20.0.0-1-gd1c87d7+1,20.0.0-15-g34741e2+3,20.0.0-2-g04cfba9+2,20.0.0-2-gec03fae+2,20.0.0-3-g082faa5+1,20.0.0-3-g2fa8bb8+6,20.0.0-3-gbdbfa727+1,20.0.0-4-gde602ef96+3,20.0.0-4-ge48a6ca+4,20.0.0-8-g1acaa36,20.0.0-8-g7eef53f7+5,w.2020.27
LSSTDataManagementBasePackage
utils.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2014-2015 AURA/LSST.
4  *
5  * This product includes software developed by the
6  * LSST Project (http://www.lsst.org/).
7  *
8  * This program is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the LSST License Statement and
19  * the GNU General Public License along with this program. If not,
20  * see <https://www.lsstcorp.org/LegalNotices/>.
21  */
22 
25 
26 #include "lsst/sphgeom/utils.h"
27 
28 #include <cmath>
29 
31 
32 
33 namespace lsst {
34 namespace sphgeom {
35 
37  Vector3d const & a,
38  Vector3d const & b,
39  Vector3d const & n)
40 {
41  Vector3d vxn = v.cross(n);
42  if (vxn.dot(a) > 0.0 && vxn.dot(b) < 0.0) {
43  // v is in the lune defined by the half great circle passing through
44  // n and a and the half great circle passing through n and b, so p
45  // is in the interior of the great circle segment from a to b. The
46  // angle θ between p and v satisfies ‖v‖ ‖n‖ sin θ = |v·n|,
47  // and ‖v‖ ‖n‖ cos θ = ‖v × n‖. The desired squared chord length is
48  // 4 sin²(θ/2).
49  double s = std::fabs(v.dot(n));
50  double c = vxn.getNorm();
51  double theta = (c == 0.0) ? 0.5 * PI : std::atan(s / c);
52  double d = std::sin(0.5 * theta);
53  return 4.0 * d * d;
54  }
55  return 4.0;
56 }
57 
59  Vector3d const & a,
60  Vector3d const & b,
61  Vector3d const & n)
62 {
63  Vector3d vxn = v.cross(n);
64  if (vxn.dot(a) < 0.0 && vxn.dot(b) > 0.0) {
65  // v is in the lune defined by the half great circle passing through
66  // n and -a and the half great circle passing through n and -b, so p
67  // is in the interior of the great circle segment from a to b. The
68  // angle θ between p and v satisfies ‖v‖ ‖n‖ sin θ = |v·n|,
69  // and ‖v‖ ‖n‖ cos θ = -‖v × n‖. The desired squared chord length is
70  // 4 sin²(θ/2).
71  double s = std::fabs(v.dot(n));
72  double c = - vxn.getNorm();
73  double d = std::sin(0.5 * std::atan2(s, c));
74  return 4.0 * d * d;
75  }
76  return 0.0;
77 }
78 
80  UnitVector3d const & v1,
81  UnitVector3d const & v2)
82 {
83  // For the details, see:
84  //
85  // The centroid and inertia tensor for a spherical triangle
86  // John E. Brock
87  // 1974, Naval Postgraduate School, Monterey Calif.
88  //
89  // https://openlibrary.org/books/OL25493734M/The_centroid_and_inertia_tensor_for_a_spherical_triangle
90 
91  Vector3d x01 = v0.robustCross(v1); // twice the cross product of v0 and v1
92  Vector3d x12 = v1.robustCross(v2);
93  Vector3d x20 = v2.robustCross(v0);
94  double s01 = 0.5 * x01.normalize(); // sine of the angle between v0 and v1
95  double s12 = 0.5 * x12.normalize();
96  double s20 = 0.5 * x20.normalize();
97  double c01 = v0.dot(v1); // cosine of the angle between v0 and v1
98  double c12 = v1.dot(v2);
99  double c20 = v2.dot(v0);
100  double a0 = (s12 == 0.0 && c12 == 0.0) ? 0.0 : std::atan2(s12, c12);
101  double a1 = (s20 == 0.0 && c20 == 0.0) ? 0.0 : std::atan2(s20, c20);
102  double a2 = (s01 == 0.0 && c01 == 0.0) ? 0.0 : std::atan2(s01, c01);
103  return 0.5 * (x01 * a2 + x12 * a0 + x20 * a1);
104 }
105 
106 }} // namespace lsst::sphgeom
lsst::utils::n
Definition: Demangle.cc:62
lsst::sphgeom::PI
constexpr double PI
Definition: constants.h:36
std::fabs
T fabs(T... args)
std::atan2
T atan2(T... args)
UnitVector3d.h
This file declares a class for representing unit vectors in ℝ³.
lsst::sphgeom::UnitVector3d::robustCross
Vector3d robustCross(UnitVector3d const &v) const
a.robustCross(b) is (b + a).cross(b - a) - twice the cross product of a and b.
Definition: UnitVector3d.h:161
lsst::sphgeom::UnitVector3d::dot
double dot(Vector3d const &v) const
dot returns the inner product of this unit vector and v.
Definition: UnitVector3d.h:152
lsst::sphgeom::Vector3d
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition: Vector3d.h:44
lsst::sphgeom::UnitVector3d
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
Definition: UnitVector3d.h:55
lsst::sphgeom::Vector3d::dot
double dot(Vector3d const &v) const
dot returns the inner product of this vector and v.
Definition: Vector3d.h:73
lsst::sphgeom::Vector3d::getNorm
double getNorm() const
getNorm returns the L2 norm of this vector.
Definition: Vector3d.h:81
lsst::sphgeom::Vector3d::normalize
double normalize()
normalize scales this vector to have unit norm and returns its norm prior to scaling.
Definition: Vector3d.cc:41
lsst::sphgeom::Vector3d::cross
Vector3d cross(Vector3d const &v) const
cross returns the cross product of this vector and v.
Definition: Vector3d.h:101
std::atan
T atan(T... args)
utils.h
This file declares miscellaneous utility functions.
lsst::sphgeom::getMaxSquaredChordLength
double getMaxSquaredChordLength(Vector3d const &v, Vector3d const &a, Vector3d const &b, Vector3d const &n)
Let p be the unit vector furthest from v that lies on the plane with normal n in the direction of the...
Definition: utils.cc:58
b
table::Key< int > b
Definition: TransmissionCurve.cc:467
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
std::sin
T sin(T... args)
a
table::Key< int > a
Definition: TransmissionCurve.cc:466
lsst::sphgeom::getMinSquaredChordLength
double getMinSquaredChordLength(Vector3d const &v, Vector3d const &a, Vector3d const &b, Vector3d const &n)
Let p be the unit vector closest to v that lies on the plane with normal n in the direction of the cr...
Definition: utils.cc:36
lsst::sphgeom::getWeightedCentroid
Vector3d getWeightedCentroid(UnitVector3d const &v0, UnitVector3d const &v1, UnitVector3d const &v2)
getWeightedCentroid returns the center of mass of the given spherical triangle (assuming a uniform ma...
Definition: utils.cc:79