1// Source:  https://github.com/golang/geo
2/*
3Copyright 2014 Google Inc. All rights reserved.
4
5Licensed under the Apache License, Version 2.0 (the "License");
6you may not use this file except in compliance with the License.
7You may obtain a copy of the License at
8
9    http://www.apache.org/licenses/LICENSE-2.0
10
11Unless required by applicable law or agreed to in writing, software
12distributed under the License is distributed on an "AS IS" BASIS,
13WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14See the License for the specific language governing permissions and
15limitations under the License.
16*/
17
18package r3
19
20import (
21	"fmt"
22	"math"
23
24	"github.com/golang/geo/s1"
25)
26
27// Vector represents a point in ℝ³.
28type Vector struct {
29	X, Y, Z float64
30}
31
32// ApproxEqual reports whether v and ov are equal within a small epsilon.
33func (v Vector) ApproxEqual(ov Vector) bool {
34	const epsilon = 1e-16
35	return math.Abs(v.X-ov.X) < epsilon && math.Abs(v.Y-ov.Y) < epsilon && math.Abs(v.Z-ov.Z) < epsilon
36}
37
38func (v Vector) String() string { return fmt.Sprintf("(%0.24f, %0.24f, %0.24f)", v.X, v.Y, v.Z) }
39
40// Norm returns the vector's norm.
41func (v Vector) Norm() float64 { return math.Sqrt(v.Dot(v)) }
42
43// Norm2 returns the square of the norm.
44func (v Vector) Norm2() float64 { return v.Dot(v) }
45
46// Normalize returns a unit vector in the same direction as v.
47func (v Vector) Normalize() Vector {
48	if v == (Vector{0, 0, 0}) {
49		return v
50	}
51	return v.Mul(1 / v.Norm())
52}
53
54// IsUnit returns whether this vector is of approximately unit length.
55func (v Vector) IsUnit() bool {
56	const epsilon = 5e-14
57	return math.Abs(v.Norm2()-1) <= epsilon
58}
59
60// Abs returns the vector with nonnegative components.
61func (v Vector) Abs() Vector { return Vector{math.Abs(v.X), math.Abs(v.Y), math.Abs(v.Z)} }
62
63// Add returns the standard vector sum of v and ov.
64func (v Vector) Add(ov Vector) Vector { return Vector{v.X + ov.X, v.Y + ov.Y, v.Z + ov.Z} }
65
66// Sub returns the standard vector difference of v and ov.
67func (v Vector) Sub(ov Vector) Vector { return Vector{v.X - ov.X, v.Y - ov.Y, v.Z - ov.Z} }
68
69// Mul returns the standard scalar product of v and m.
70func (v Vector) Mul(m float64) Vector { return Vector{m * v.X, m * v.Y, m * v.Z} }
71
72// Dot returns the standard dot product of v and ov.
73func (v Vector) Dot(ov Vector) float64 { return v.X*ov.X + v.Y*ov.Y + v.Z*ov.Z }
74
75// Cross returns the standard cross product of v and ov.
76func (v Vector) Cross(ov Vector) Vector {
77	return Vector{
78		v.Y*ov.Z - v.Z*ov.Y,
79		v.Z*ov.X - v.X*ov.Z,
80		v.X*ov.Y - v.Y*ov.X,
81	}
82}
83
84// Distance returns the Euclidean distance between v and ov.
85func (v Vector) Distance(ov Vector) float64 { return v.Sub(ov).Norm() }
86
87// Angle returns the angle between v and ov.
88func (v Vector) Angle(ov Vector) s1.Angle {
89	return s1.Angle(math.Atan2(v.Cross(ov).Norm(), v.Dot(ov))) * s1.Radian
90}
91
92// Axis enumerates the 3 axes of ℝ³.
93type Axis int
94
95// The three axes of ℝ³.
96const (
97	XAxis Axis = iota
98	YAxis
99	ZAxis
100)
101
102// Ortho returns a unit vector that is orthogonal to v.
103// Ortho(-v) = -Ortho(v) for all v.
104func (v Vector) Ortho() Vector {
105	ov := Vector{0.012, 0.0053, 0.00457}
106	switch v.LargestComponent() {
107	case XAxis:
108		ov.Z = 1
109	case YAxis:
110		ov.X = 1
111	default:
112		ov.Y = 1
113	}
114	return v.Cross(ov).Normalize()
115}
116
117// LargestComponent returns the axis that represents the largest component in this vector.
118func (v Vector) LargestComponent() Axis {
119	t := v.Abs()
120
121	if t.X > t.Y {
122		if t.X > t.Z {
123			return XAxis
124		}
125		return ZAxis
126	}
127	if t.Y > t.Z {
128		return YAxis
129	}
130	return ZAxis
131}
132
133// SmallestComponent returns the axis that represents the smallest component in this vector.
134func (v Vector) SmallestComponent() Axis {
135	t := v.Abs()
136
137	if t.X < t.Y {
138		if t.X < t.Z {
139			return XAxis
140		}
141		return ZAxis
142	}
143	if t.Y < t.Z {
144		return YAxis
145	}
146	return ZAxis
147}
148
149// Cmp compares v and ov lexicographically and returns:
150//
151//   -1 if v <  ov
152//    0 if v == ov
153//   +1 if v >  ov
154//
155// This method is based on C++'s std::lexicographical_compare. Two entities
156// are compared element by element with the given operator. The first mismatch
157// defines which is less (or greater) than the other. If both have equivalent
158// values they are lexicographically equal.
159func (v Vector) Cmp(ov Vector) int {
160	if v.X < ov.X {
161		return -1
162	}
163	if v.X > ov.X {
164		return 1
165	}
166
167	// First elements were the same, try the next.
168	if v.Y < ov.Y {
169		return -1
170	}
171	if v.Y > ov.Y {
172		return 1
173	}
174
175	// Second elements were the same return the final compare.
176	if v.Z < ov.Z {
177		return -1
178	}
179	if v.Z > ov.Z {
180		return 1
181	}
182
183	// Both are equal
184	return 0
185}
186