@@ -15,166 +15,19 @@ which is nice.
1515*** DONE =legacy_write= read/write :unit:
1616** TODO New functionality :functionality:
1717
18- *** TODO Great Circle Arc calculation
18+ *** Incorporate Ellipticity in geometric calculations
1919
20- Use [[https://en.wikipedia.org/wiki/Haversine_formula][Haversine]] formula to start. Later can move on to Lambert's formula for long
21- lines ([[https://en.wikipedia.org/wiki/Geographical_distance#Lambert's_formula_for_long_lines][Wikipedia]])
22-
23- [[file:~/my_git/sac-format/src/sac_format.hpp::void calc_gcarc();]]
24-
25- ***** DONE n-vector for distance calculations
26-
27- [[https://en.wikipedia.org/wiki/N-vector][Wikipedia]], it has no singularities (unlike standard latitude/longitude).
28-
29- [[file:~/my_git/sac-format/src/sac_format.hpp::void calc_gcarc();]]
30-
31- ****** Turns out it is really quite simple to use
32-
33- [[https://www.ffi.no/en/research/n-vector/n-vector-explained][Norwegian Defense Research Establishment]]
34- [[https://github.com/FFI-no/n-vector][GitHub]]
35-
36- They provide an explanation (including using the ellipticity of the Earth) and
37- matlab codes. I can write up a C++ version.
38-
39- They also provide a reference paper: [[https://www.navlab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf][Gade (2010)]]
40-
41- ******* n-vector (normal)
42-
43- $\mathbf{n}_{AB}^{C}$ is the normal vector, from $A \rightarrow B$ in reference
44- frame $C$ ($C$ is often $E$ for Earth, $A$ is often $E$ for center of Earth).
45-
46- ******** Calculate n-vector from lat/lon
47-
48- $\mathbf{n}_{EA}^{E} = \left[ sin(\lambda), sin(\mu)\cdot cos(\lambda),
49- -cos(\mu)\cdot \cos(\lambda) \right]$
50-
51- where $\lambda \in \left[-\frac{\pi}{2},\frac{\pi}{2}\right]$
52- is latitude and $\mu \in \left(-\pi,\pi\right]$
53-
54- ******* p-vector (position)
55-
56- $\mathbf{p}_{EA}^{E}$ is the position vector, from the center of the Earth to
57- position $A$.
58-
59- ********* calculate p-vector from n-vector
60-
61- $\mathbf{p}_{EA}^{A} = \frac{b}{\sqrt{(n_{x}^{E})^{2} +
62- \frac{a^{2}}{b^{2}}(n_{y}^{E})^{2} +
63- \frac{a^{2}}{b^{2}}(n_{z}^{E})^{2}}}\left[n_{x}^{E},\frac{a^{2}}{b^{2}}n_{y}^{E},\frac{a^{2}}{b^{2}}n_{z}^{E}\right] +
64- h\mathbf{n}_{EA}^{A}$
65-
66- where $a,b$ are the semi-major and semi-minor axes of the reference ellipsoid
67- model being used.
68-
69- *NOTE* that usually for a reference ellipsoid we use $a,f$, where $a$ is the semi-major axis and $f$ is the inverse flattening and then to calculate $b$ use $b=a*(1-f)$.
70-
71- For *WGS84*: $a=6378137$ meters and $f=\frac{1}{298.257223563}$.
72-
73- ******* Distances between two points on surface of Earth
74-
75- $\Delta\mathbf{p}_{AB}^{E} = \mathbf{p}_{EB}^{E} - \mathbf{p}_{EA}^{E}$ is the distance vector
76-
77- Then, distance is then $\lvert \Delta\mathbf{p}_{AB}^{E} \rvert =
78- \sqrt{\Delta\mathbf{p}_{AB}^{E}\cdot\Delta\mathbf{p}_{AB}^{E}}$ and the
79- horizontal distance is the same, without the vertical component.
80-
81- ****** Simpler
82-
83- $\Delta\sigma = arctan\left(\frac{\lvert\mathbf{n}_{1}\times\mathbf{n}_{2}\rvert}{\mathbf{n}_{1}\cdot\mathbf{n}_{2}}\right)$ (assumes spherical Earth)
84-
85- ***** TODO Link paired headers
86- e.g., az and baz should be linked, if one is edited the other is updated
87- Location parameters (stla, stlo, evla, evlo) should be linked to dist and gcarc.
88-
89- ***** TODO Azimuth
90- function n_EA_E_and_EB_E2azimuth from
20+ Tried the n-vector stuff, didn't work. Moving on for now.
9121
92- [[https://github.com/pbrod/nvector/blob/v0.7.0/src/nvector/_core.py][GitHub]]
22+ *** DONE Great Circle Arc calculation
9323
94- ****** Summary
24+ Assumes spherical earth
9525
96- ******** Step 1
97-
98- Rotation matrix
99- $\mathbf{E} =
100- \begin{pmatrix}
101- 0 & 0 & 1\\
102- 0 & 1 & 0\\
103- -1 & 0 & 0
104- \end{pmatrix}$
105-
106- defines coordinate system where the x-direction points toward (latitude,
107- longitude) = (0, 0) and z-direction points to the north pole.
108-
109- $det(\mathbf{E}) = 1$
110-
111- ******** Step 2
112-
113- $\mathbf{n}_{EA}^{E}\rightarrow\mathbf{p}_{EA}^{E}$
114-
115- $\mathbf{n}_{EA}^{E}=\frac{1}{det(\mathbf{E})*\lvert\mathbf{n}_{EA}^{E}
116- \rvert}\mathbf{E}\cdot\mathbf{n}_{EA}^{E}$
117-
118- Note that both $\mathbf{E}$ and $\mathbf{n}_{EA}^{E}$ *are* normalized, I suspect
119- the renormalization is for the sake of numerical stability.
120-
121- $\mathbf{p}_{EL}^{E}=\frac{b}{\sqrt{(n_{x}^{E})^{2} +
122- \frac{a^{2}}{b^{2}}(n_{y}^{E})^{2} +
123- \frac{a^{2}}{b^{2}}(n_{z}^{E})^{2}}}\left[n_{x}^{E},\frac{a^{2}}{b^{2}}n_{y}^{E},\frac{a^{2}}{b^{2}}n_{z}^{E}\right]$
124-
125- $\mathbf{p}_{EA}^{E}=\mathbf{E}^{T}\cdot\left(\mathbf{p}_{EL}^{E} - h\mathbf{n}_{EA}^{E}\right)$ where $h$ is depth (below ellipsoid in meters).
126-
127- Calculate $\mathbf{p}_{EA}^{E}$ and $\mathbf{p}_{EB}^{E}$
128-
129- ******** Step 3
130-
131- $\mathbf{R}_{EN}$ from $\mathbf{n}_{EA}^{E}$ and $\mathbf{E}$:
132- $\mathbf{n}_{EA}^{E}=\frac{1}{det(\mathbf{E})*\lvert\mathbf{n}_{EA}^{E}\rvert}\mathbf{E}\cdot\mathbf{n}_{EA}^{E}$ (result is a matrix)
133- (renormalize for numerical stability)
134-
135- New coordinate frame $N$, where the z-axis is north (instead of E [east])
136-
137- $N_{z}^{E}=-\mathbf{n}_{EA}^{E}$ (points opposite direction of n-vector)
138-
139- $N_{y}^{E} = \begin{pmatrix}1\\0\\0\end{pmatrix}\times\mathbf{n}_{EA}^{E}$ (find
140- y-direction, (I think that this is the cross product with each column vector that forms $\mathbf{n}_{EA}^{E}$) equation 9) (points perpendicular
141- to plane formed by n-vector and Earth's spin axis)
142-
143- //
144- Normalize (I can't tell if this is a vector of a matrix at the moment,
145- everything /should/ be vectors, but their python code is very weird)
146- //
147-
148- $N_{y}^{E}=\frac{N_{y}^{E}}{\lvert N_{y}^{E}\rvert}$ normalize again
149-
150- This doesn't make sense to me:
151- $N_{y}^{E}[:, poles] = \begin{pmatrix}0\\1\\0\end{pmatrix}$ select y-axis direction
152-
153- Find x-axis direction of N
154- $N_{x}^{E} = N_{y}^{E}\times N_{z}^{E}$ just the right-hand rule (that is simple)
155-
156- Then some more shenanigans
157- $N_{xyz}^{E} = \left[N_{x}^{E}, N_{y}^{E}, N_{z}^{E}\right]$ where this is a full matrix (because each $N_{a}^{E}$ is a column vector)
158-
159- *Note*
160-
161-
162- $\mathbf{E}^{T} = \begin{pmatrix}0 & 0 & -1\\0 & 1 & 0\\1 & 0 & 0\end{pmatrix}$
163-
164- $\mathbf{R}_{EN}=\mathbf{E}^{T}\cdot N_{x}^{E}\cdot N_{y}^{E}\cdot N_{z}^{E}$
165-
166- ******** Step 4
167-
168- Then it seems you do a modified (?) version of
169- $\mathbf{R}_{EN}\cdot\mathbf{p}_{AB}^{E}$ to calculate $\mathbf{p}_{AB}^{N}$,
170- but they do weird python stuff with mdot, rollaxis (looks like they swap axis 0
171- and 1 [is axis a row or a column?]), and reshape. Their original version is
172- $\mathbf{p}_{AB}^{N}=\mathbf{R}_{EM}^{T}\cdot\mathbf{p}_{AB}^{E}$ (T means
173- transpose).
174-
175- ******** Step 5
26+ Use [[https://en.wikipedia.org/wiki/Haversine_formula][Haversine]] formula to start. Later can move on to Lambert's formula for long
27+ lines ([[https://en.wikipedia.org/wiki/Geographical_distance#Lambert's_formula_for_long_lines][Wikipedia]])
17628
177- Then the azimuth in radians is atan2$(\mathbf{p}_{AB}^{N}[1], \mathbf{p}_{AB}^{N}[0])$.
29+ *** DONE Azimuth calculation
30+ Assumes spherical earth
17831
17932*** TODO Read-inspect function
18033Like read, but it prints to std::cout for every header it readers (not data
0 commit comments