LSST Applications g1653933729+b86e4b8053,g2303c004a3+a99d6d42bd,g28da252d5a+3e0fd49a16,g2bbee38e9b+db5ba6d06c,g2bc492864f+db5ba6d06c,g2cdde0e794+b03b09fb0a,g3156d2b45e+b0ccf75d80,g347aa1857d+db5ba6d06c,g35bb328faa+b86e4b8053,g3a166c0a6a+db5ba6d06c,g3e281a1b8c+c9a36f98fc,g414038480c+70d5606f74,g42360bd172+a99d6d42bd,g461a3dce89+b86e4b8053,g5c9cb89fa5+b7dae1cdc6,g7f102674c4+d27c5e2769,g7f343241b9+516b07f5d7,g80478fca09+ee10125a54,g82479be7b0+11c35b2845,g858d7b2824+a99d6d42bd,g8cd86fa7b1+cb4c24add2,g9125e01d80+b86e4b8053,g979bb04a14+30ae7a28c8,g9ddcbc5298+b4cb6f9276,ga44b09e51c+390e0c76f2,ga5288a1d22+c32cd41bfc,gae0086650b+b86e4b8053,gb58c049af0+2fe46b4317,gc28159a63d+db5ba6d06c,gcd45df26be+a99d6d42bd,gcf0d15dbbd+be93eff7c7,gda6a2b7d83+be93eff7c7,gdaeeff99f8+b9455e0a82,ge5cf12406b+e20940066e,ge79ae78c31+db5ba6d06c,ge7c99ddb59+26aa1bc97f,gf048a9a2f4+f6d4943a12,gf0baf85859+e0282cf37d,w.2024.32
LSST Data Management Base Package
Loading...
Searching...
No Matches
utils.cc
Go to the documentation of this file.
1/*
2 * This file is part of sphgeom.
3 *
4 * Developed for the LSST Data Management System.
5 * This product includes software developed by the LSST Project
6 * (http://www.lsst.org).
7 * See the COPYRIGHT file at the top-level directory of this distribution
8 * for details of code ownership.
9 *
10 * This software is dual licensed under the GNU General Public License and also
11 * under a 3-clause BSD license. Recipients may choose which of these licenses
12 * to use; please see the files gpl-3.0.txt and/or bsd_license.txt,
13 * respectively. If you choose the GPL option then the following text applies
14 * (but note that there is still no warranty even if you opt for BSD instead):
15 *
16 * This program is free software: you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation, either version 3 of the License, or
19 * (at your option) any later version.
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program. If not, see <http://www.gnu.org/licenses/>.
28 */
29
32
33#include "lsst/sphgeom/utils.h"
34
35#include <cmath>
36
38
39
40namespace lsst {
41namespace sphgeom {
42
44 Vector3d const & a,
45 Vector3d const & b,
46 Vector3d const & n)
47{
48 Vector3d vxn = v.cross(n);
49 if (vxn.dot(a) > 0.0 && vxn.dot(b) < 0.0) {
50 // v is in the lune defined by the half great circle passing through
51 // n and a and the half great circle passing through n and b, so p
52 // is in the interior of the great circle segment from a to b. The
53 // angle θ between p and v satisfies ‖v‖ ‖n‖ sin θ = |v·n|,
54 // and ‖v‖ ‖n‖ cos θ = ‖v × n‖. The desired squared chord length is
55 // 4 sin²(θ/2).
56 double s = std::fabs(v.dot(n));
57 double c = vxn.getNorm();
58 double theta = (c == 0.0) ? 0.5 * PI : std::atan(s / c);
59 double d = std::sin(0.5 * theta);
60 return 4.0 * d * d;
61 }
62 return 4.0;
63}
64
66 Vector3d const & a,
67 Vector3d const & b,
68 Vector3d const & n)
69{
70 Vector3d vxn = v.cross(n);
71 if (vxn.dot(a) < 0.0 && vxn.dot(b) > 0.0) {
72 // v is in the lune defined by the half great circle passing through
73 // n and -a and the half great circle passing through n and -b, so p
74 // is in the interior of the great circle segment from a to b. The
75 // angle θ between p and v satisfies ‖v‖ ‖n‖ sin θ = |v·n|,
76 // and ‖v‖ ‖n‖ cos θ = -‖v × n‖. The desired squared chord length is
77 // 4 sin²(θ/2).
78 double s = std::fabs(v.dot(n));
79 double c = - vxn.getNorm();
80 double d = std::sin(0.5 * std::atan2(s, c));
81 return 4.0 * d * d;
82 }
83 return 0.0;
84}
85
87 UnitVector3d const & v1,
88 UnitVector3d const & v2)
89{
90 // For the details, see:
91 //
92 // The centroid and inertia tensor for a spherical triangle
93 // John E. Brock
94 // 1974, Naval Postgraduate School, Monterey Calif.
95 //
96 // https://openlibrary.org/books/OL25493734M/The_centroid_and_inertia_tensor_for_a_spherical_triangle
97
98 Vector3d x01 = v0.robustCross(v1); // twice the cross product of v0 and v1
99 Vector3d x12 = v1.robustCross(v2);
100 Vector3d x20 = v2.robustCross(v0);
101 double s01 = 0.5 * x01.normalize(); // sine of the angle between v0 and v1
102 double s12 = 0.5 * x12.normalize();
103 double s20 = 0.5 * x20.normalize();
104 double c01 = v0.dot(v1); // cosine of the angle between v0 and v1
105 double c12 = v1.dot(v2);
106 double c20 = v2.dot(v0);
107 double a0 = (s12 == 0.0 && c12 == 0.0) ? 0.0 : std::atan2(s12, c12);
108 double a1 = (s20 == 0.0 && c20 == 0.0) ? 0.0 : std::atan2(s20, c20);
109 double a2 = (s01 == 0.0 && c01 == 0.0) ? 0.0 : std::atan2(s01, c01);
110 return 0.5 * (x01 * a2 + x12 * a0 + x20 * a1);
111}
112
113}} // namespace lsst::sphgeom
table::Key< int > b
This file declares a class for representing unit vectors in ℝ³.
T atan2(T... args)
T atan(T... args)
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
double dot(Vector3d const &v) const
dot returns the inner product of this unit vector and v.
Vector3d robustCross(UnitVector3d const &v) const
a.robustCross(b) is (b + a).cross(b - a) - twice the cross product of a and b.
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition Vector3d.h:51
double dot(Vector3d const &v) const
dot returns the inner product of this vector and v.
Definition Vector3d.h:80
double getNorm() const
getNorm returns the L2 norm of this vector.
Definition Vector3d.h:88
double normalize()
normalize scales this vector to have unit norm and returns its norm prior to scaling.
Definition Vector3d.cc:48
Vector3d cross(Vector3d const &v) const
cross returns the cross product of this vector and v.
Definition Vector3d.h:108
T fabs(T... args)
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:65
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:86
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:43
constexpr double PI
Definition constants.h:43
T sin(T... args)
This file declares miscellaneous utility functions.