LSST Applications g011c388f00+c533e7e796,g0265f82a02+71b8ded0d6,g16a3bce237+71b8ded0d6,g2079a07aa2+b9108c1c87,g2bbee38e9b+71b8ded0d6,g337abbeb29+71b8ded0d6,g3ddfee87b4+efeeb9cceb,g44050f54f7+816a004954,g4cf46543a9+ac9bf1619d,g50ff169b8f+8309cf5058,g52b1c1532d+43dac7135f,g858d7b2824+816a004954,g8a8a8dda67+43dac7135f,g99855d9996+36230435de,g9cb17f706c+56afb287bb,g9ddcbc5298+389b8f2b7e,ga1e77700b3+4bafba478f,ga8c6da7877+9a92598c84,gae46bcf261+71b8ded0d6,gb0e22166c9+713927f999,gb700894bec+0b6c79d923,gb8350603e9+a03320826c,gba4ed39666+e464e2e6f0,gbeb006f7da+8cd302297b,gc86a011abf+816a004954,gcf0d15dbbd+efeeb9cceb,gd0e876f1fb+efeeb9cceb,gd162630629+3eacc90e0c,gd44f2fa1a7+82b768a15d,gd8a3c24633+eaa4ad6639,gdaeeff99f8+6b435c3f92,ge2eec9bf53+fad368631c,ge79ae78c31+71b8ded0d6,gee10cc3b42+43dac7135f,gf1cff7945b+816a004954,w.2024.08
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
table::Key< int > a
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.