LSST Applications g0d97872fb5+76bdf3f114,g1653933729+34a971ddd9,g28da252d5a+072f89fe25,g2bbee38e9b+a99b0ab4cd,g2bc492864f+a99b0ab4cd,g2ca4be77d2+cae847d7e5,g2cdde0e794+704103fe75,g3156d2b45e+6e87dc994a,g347aa1857d+a99b0ab4cd,g35bb328faa+34a971ddd9,g3a166c0a6a+a99b0ab4cd,g3e281a1b8c+8ec26ec694,g4005a62e65+ba0306790b,g414038480c+9ed5ed841a,g569e0e2b34+cb4faa46ad,g5a97de2502+520531a62c,g717e5f8c0f+29153700a5,g7ede599f99+924b0fe0eb,g80478fca09+17051a22cc,g82479be7b0+061a565528,g858d7b2824+29153700a5,g8b782ad322+29153700a5,g8cd86fa7b1+05420e7f7d,g9125e01d80+34a971ddd9,ga5288a1d22+e7f674aaf3,gae0086650b+34a971ddd9,gae74b0b5c6+45ef5cdc51,gb58c049af0+ace264a4f2,gc28159a63d+a99b0ab4cd,gcf0d15dbbd+8051a81198,gda6a2b7d83+8051a81198,gdaeeff99f8+7774323b41,gdf4d240d4a+34a971ddd9,ge2409df99d+de9c54f8ec,ge33fd446bb+29153700a5,ge79ae78c31+a99b0ab4cd,gf0baf85859+890af219f9,gf5289d68f6+234e029470,w.2024.36
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:49
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.