Benutzer:Dirk Huenniger/geo2

Erstmal müssen wir von geographischen Koordinaten in Kartesische Kooridinaten.

a) Kirche Nieukerk  N 51 27.359 E 6 22.345  H= 95 ü.NN
b) Kirche Aldekerk  N 51 26.377 E 6 24.947  H= 133 ü.NN
c) Baggersee Meerendonker Straße  N 51 26.080 E 6 20.488  26m ü.NN

Wikipedia sagt (w:Referenzellipsoid):

Umrechnung in geozentrische Kartesische Koordinaten

In einem geozentrischen rechtwinkligen Bezugssystem, dessen Ursprung im Mittelpunkt des Rotationsellipsoids liegt und in Richtung der Rotationsachse ($Z$) sowie des Nullmeridians ($X$) ausgerichtet ist, gilt dann

$X=(N_{\varphi }+h)\cos \varphi \cdot \cos \lambda ;\,\!$

$Y=(N_{\varphi }+h)\cos \varphi \cdot \sin \lambda ;\,\!$

$Z=(N_{\varphi }(1-\varepsilon ^{2})+h)\sin \varphi ;\,\!$

mit

$a$ - große Halbachse (Parameter des Referenzellipsoids)

$b$ - kleine Halbachse (Parameter des Referenzellipsoids)

$\varepsilon ={\frac {\sqrt {a^{2}-b^{2}}}{a}}$ - numerische Exzentrizität

$N_{\varphi }={\frac {a}{\sqrt {1-\varepsilon ^{2}\sin ^{2}\varphi }}}\,\!$ - Krümmungsradius des Ersten Vertikals, d.h. der Abstand des Lotfußpunktes vom Schnittpunkt des verlängerten Lots mit der Z-Achse.

Für GPS gilt

a=6378137.0

b= 6356752.3142

Er mal brauchen wir die Winkel im Bogenmaß

$\alpha =\left(\mathrm {Grad} +{\frac {\mathrm {Minuten} }{60}}\right){\frac {\pi }{180}}$

Das können wir nun in die Wikipedia Formeln einsetzen und erhalten 3 Vektoren in Kartesische Koordinaten $\mathbf {a} ,\mathbf {b} ,\mathbf {c}$ im $\mathbb {R} ^{3}$. Durch den Fettdruck wird angedeute das jeweils ein Vektor mit 3 Komponenten durch einen fetten Buchstaben bezeichnet wird. Der Mittelpunkt zwischen a und b ergibt sich nach.

$\mathbf {d} ={\frac {\mathbf {a} +\mathbf {b} }{2}}$

Nun suchen wir die Ebene die von den Vektoren $\mathbf {a} ,\mathbf {b} ,\mathbf {c}$ aufgespannt wird. Wir berechnen nicht diese Ebene sonderen eine hierzu paralle Ebene durch den Koordianten Usprung sie wird von den beiden Vektoren $\mathbf {ab} =\mathbf {a} -\mathbf {b}$ und $\mathbf {ac} =\mathbf {a} -\mathbf {c}$ aufgespannt. $\mathbf {a}$ heirbei der name eines Vektor und zwar des Differenzvektors aus den Vektoren $\mathbf {a}$ und $\mathbf {b}$. Insbesondere bezeichnet $\mathbf {a}$ nicht das Skalarprodukt von $\mathbf {a}$ und $\mathbf {b}$. Nun bestimmen wir einen vektor der auf dieser Ebene senkrecht steht. Er ergibt sich mit Hilfe des Kreuzproduktes nach:

$\mathbf {n} =\mathbf {ab} \times \mathbf {ac}$

Wir normieren diesen nun das. Heist wir dividieren ihn durch seine Länge. Damit er die Länge eins bekommt und bezeichnen das Ergebnis mit $\mathbf {un}$

$\mathbf {un} ={\frac {\mathbf {n} }{\sqrt {\mathbf {n} \cdot \mathbf {n} }}}$

Der Vektor $\mathbf {d}$ zeigt nun wenn man ihm in Ursprung anheftet und seinen endpunkt betrachtet genau in diesem Endpunkt exakt entgegen der Schwerkraft, sofern man die Erde hier als Kugel annimmt. Daher ergibt sich ein normierter Vektor in Richtung des Schwerkraft durch

$\mathbf {dn} =-{\frac {\mathbf {d} }{\sqrt {\mathbf {d} \cdot \mathbf {d} }}}$

Nun ziehen wir von $mathbf{dn}$ seinen anteil in richtung von $\mathbf {un}$ ab. Das ergibt die Richtung in die eine auf der Ebene ruhende Kugel aufgrund des Schwerkraft zu rollen beginnen würde.

$\mathbf {mden} =(\mathbf {dn} \cdot \mathbf {un} )\mathbf {dn}$

Diesen bringen wir noch auf die Länge von 1983 Metern.

$\mathbf {g} =-1983{\frac {\mathbf {mden} }{\sqrt {\mathbf {mden} \cdot \mathbf {mden} }}}$

und addieren Ihn zum Punkt d.

$\mathbf {p} =\mathbf {g} +\mathbf {d}$

Nun beginnt die Rückrechnung in geographsiche Korrdinaten. Die geht nicht exakt, weil die Formel nicht umkehrbar sind, bzw. noch niemand eine Umkehrung gefunden hat. Der Längengrad geht noch exakt.

$\lambda =\mathrm {arctan} \left({\frac {y}{x}}\right)$

Für deb Bereiten grad ergibt sich mit sofern man die höhe erst einmal 0 setzt, was ja am Niederrhein vertretbar ist.

$\phi =\mathrm {arctan} \left({\frac {z}{\sqrt {x^{2}+y^{2}}}}{\frac {1}{1-\epsilon ^{2}}}\right)$

Nun können wir nach der Wikipedia Formel eine Schätzung für h berechnen.

$h={\frac {\sqrt {x^{2}+y^{2}}}{\mathrm {cos} \phi }}-N_{\phi }$

Die Formel

$h\approx {\sqrt {X^{2}+Y^{2}}}\cdot \cos \varphi +Z\sin \varphi -a^{2}/N_{\varphi }$

ist aus Wikipedia übernommen.

Die Formel

$\phi =\mathrm {arctan} \left({\frac {z}{\sqrt {x^{2}+y^{2}}}}{\frac {h+N_{\varphi }}{N_{\varphi }+(1-\epsilon ^{2})}}\right)$ ergibt sich durch umstellen der oben genannten Formel aus Wikipedia

Zur erhöhnung der Genauigkeit wird dieses Verfahren noch zwei mal iteriert. Das heist mit den Schätzung von $h$ und $phi$ werden jeweils neue bessere Schätzungen von $h$ und $phi$ berechnet