LSST Applications g05325017a3+be04ab5bc8,g2079a07aa2+3e9fd84d81,g20cdd03214+600de29eb2,g2305ad1205+5e1b0c67c7,g2bbee38e9b+69e041a28d,g337abbeb29+69e041a28d,g3a166c0a6a+69e041a28d,g4322eb9e3a+9670871960,g4e26c3aaff+9125fe5495,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+608c162c81,g5c3423f6d4+600de29eb2,g607f77f49a+600de29eb2,g69c4c69b53+69e041a28d,g858d7b2824+600de29eb2,g99cad8db69+a7f8c65a2d,g9ddcbc5298+9a081db1e4,ga1e77700b3+2cbb763275,ga8d6e2f9c9+7a336fa28a,gae0086650b+585e252eca,gb0e22166c9+0e73c8378f,gb7c702c636+aec2f3b930,gbb8dafda3b+7d1409958d,gc120e1dc64+acd995d0d8,gc28159a63d+69e041a28d,gcdd4ae20e8+5ef4c697da,gcf0d15dbbd+5ef4c697da,gd1535ee943+1ef0950dc9,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+69e041a28d,gf905fbdb4b+471d9f2c03,gf907d663e3+01e3c6626b,gfbcc870c63+b310236976,w.2024.25
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.