From c89b932d551659f4de52eb637e0cf19a5896b1bb Mon Sep 17 00:00:00 2001 From: Hendrik Date: Wed, 15 Oct 2025 11:32:56 +0200 Subject: [PATCH] Meine alten Sachen aus RALV --- GHA/__init__.py | 0 GHA/bessel.py | 25 +++ GHA/gauss.py | 107 ++++++++++++ Numerische_Integration/__init__.py | 0 Numerische_Integration/num_int_runge_kutta.py | 62 +++++++ .../num_int_trapezformel.py | 46 +++++ Numerische_Integration/rk_DGLS_1Ordnung.py | 7 + Numerische_Integration/rk_DGL_1Ordnung.py | 6 + Numerische_Integration/rk_DGL_2Ordnung.py | 10 ++ Numerische_Integration/rk_DGL_2Ordnung_2.py | 7 + Richtungsreduktion.py | 14 ++ ausgaben.py | 27 +++ ellipsoide.py | 83 +++++++++ hausarbeit.py | 78 +++++++++ phi_ellipse.py | 55 ++++++ s_ellipse.py | 78 +++++++++ winkelumrechnungen.py | 158 ++++++++++++++++++ 17 files changed, 763 insertions(+) create mode 100644 GHA/__init__.py create mode 100644 GHA/bessel.py create mode 100644 GHA/gauss.py create mode 100644 Numerische_Integration/__init__.py create mode 100644 Numerische_Integration/num_int_runge_kutta.py create mode 100644 Numerische_Integration/num_int_trapezformel.py create mode 100644 Numerische_Integration/rk_DGLS_1Ordnung.py create mode 100644 Numerische_Integration/rk_DGL_1Ordnung.py create mode 100644 Numerische_Integration/rk_DGL_2Ordnung.py create mode 100644 Numerische_Integration/rk_DGL_2Ordnung_2.py create mode 100644 Richtungsreduktion.py create mode 100644 ausgaben.py create mode 100644 ellipsoide.py create mode 100644 hausarbeit.py create mode 100644 phi_ellipse.py create mode 100644 s_ellipse.py create mode 100644 winkelumrechnungen.py diff --git a/GHA/__init__.py b/GHA/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/GHA/bessel.py b/GHA/bessel.py new file mode 100644 index 0000000..931999e --- /dev/null +++ b/GHA/bessel.py @@ -0,0 +1,25 @@ +from numpy import * +import scipy as sp + + +def gha1(re, phi_p1, lambda_p1, A_p1, s): + psi_p1 = re.phi2psi(phi_p1) + A_0 = arcsin(cos(psi_p1) * sin(A_p1)) + temp = sin(psi_p1) / cos(A_0) + sigma_p1 = arcsin(sin(psi_p1) / cos(A_0)) + + sqrt_sigma = lambda sigma: sqrt(1 + re.e_ ** 2 * cos(A_0) ** 2 * sin(sigma) ** 2) + int_sqrt_sigma = lambda sigma: sp.integrate.quad(sqrt_sigma, sigma_p1, sigma)[0] + + f_sigma_p2i = lambda sigma_p2i: (int_sqrt_sigma(sigma_p2i) - s / re.b) + + sigma_p2_0 = sigma_p1 + s / re.a + sigma_p2 = sp.optimize.newton(f_sigma_p2i, sigma_p2_0) + psi_p2 = arcsin(cos(A_0) * sin(sigma_p2)) + phi_p2 = re.psi2phi(psi_p2) + A_p2 = arcsin(sin(A_0) / cos(psi_p2)) + + f_d_lambda = lambda sigma: sin(A_0) * sqrt_sigma(sigma) / (1 - cos(A_0)**2 * sin(sigma)**2) + d_lambda = sqrt(1-re.e**2) * sp.integrate.quad(f_d_lambda, sigma_p1, sigma_p2)[0] + lambda_p2 = lambda_p1 + d_lambda + return phi_p2, lambda_p2, A_p2 diff --git a/GHA/gauss.py b/GHA/gauss.py new file mode 100644 index 0000000..e89c264 --- /dev/null +++ b/GHA/gauss.py @@ -0,0 +1,107 @@ +from numpy import sin, cos, pi, sqrt, tan, arcsin, arccos, arctan +import ausgaben as aus + + +def gha1(re, phi_p1, lambda_p1, A_p1, s, eps): + """ + Berechnung der 1. Geodätische Hauptaufgabe nach Gauß´schen Mittelbreitenformeln + :param re: Klasse Ellipsoid + :param phi_p1: Breite Punkt 1 + :param lambda_p1: Länge Punkt 1 + :param A_p1: Azimut der geodätischen Linie in Punkt 1 + :param s: Strecke zu Punkt 2 + :param eps: Abbruchkriterium für Winkelgrößen + :return: Breite, Länge, Azimut von Punkt 2 + """ + + t = lambda phi: tan(phi) + eta = lambda phi: sqrt(re.e_ ** 2 * cos(phi) ** 2) + + F1 = lambda A, phi, s: 1 + (2 + 3 * t(phi)**2 + 2 * eta(phi)**2) / (24 * re.N(phi) ** 2) * sin(A) ** 2 * s ** 2 \ + + (t(phi)**2 - 1) * eta(phi) ** 2 / (8 * re.N(phi) ** 2) * cos(A) ** 2 * s ** 2 + + F2 = lambda A, phi, s: 1 + t(phi) ** 2 / (24 * re.N(phi) ** 2) * sin(A) ** 2 * s ** 2 \ + - (1 + eta(phi)**2 - 9 * eta(phi)**2 * t(phi)**2) / (24 * re.N(phi) ** 2) * cos(A) ** 2 * s ** 2 + + F3 = lambda A, phi, s: 1 + (1 + eta(phi)**2) / (12 * re.N(phi) ** 2) * sin(A) ** 2 * s ** 2 \ + + (3 + 8 * eta(phi)**2) / (24 * re.N(phi) ** 2) * cos(A) ** 2 * s ** 2 + + phi_p2_i = lambda A, phi: phi_p1 + cos(A) / re.M(phi) * s * F1(A, phi, s) + lambda_p2_i = lambda A, phi: lambda_p1 + sin(A) / (re.N(phi) * cos(phi)) * s * F2(A, phi, s) + A_p2_i = lambda A, phi: A_p1 + sin(A) * tan(phi) / re.N(phi) * s * F3(A, phi, s) + + phi_p0_i = lambda phi2: (phi_p1 + phi2) / 2 + A_p1_i = lambda A2: (A_p1 + A2) / 2 + + phi_p0 = [] + A_p0 = [] + phi_p2 = [] + lambda_p2 = [] + A_p2 = [] + + # 1. Näherung für P2 + phi_p2.append(phi_p1 + cos(A_p1) / re.M(phi_p1) * s) + lambda_p2.append(lambda_p1 + sin(A_p1) / (re.N(phi_p1) * cos(phi_p1)) * s) + A_p2.append(A_p1 + sin(A_p1) * tan(phi_p1) / re.N(phi_p1) * s) + + while True: + # Berechnug P0 durch Mittelbildung + phi_p0.append(phi_p0_i(phi_p2[-1])) + A_p0.append(A_p1_i(A_p2[-1])) + # Berechnung P2 + phi_p2.append(phi_p2_i(A_p0[-1], phi_p0[-1])) + lambda_p2.append(lambda_p2_i(A_p0[-1], phi_p0[-1])) + A_p2.append(A_p2_i(A_p0[-1], phi_p0[-1])) + + # Abbruchkriterium + if abs(phi_p2[-2] - phi_p2[-1]) < eps and \ + abs(lambda_p2[-2] - lambda_p2[-1]) < eps and \ + abs(A_p2[-2] - A_p2[-1]) < eps: + break + + nks = 5 + for i in range(len(phi_p2)): + print(f"P2[{i}]: {aus.gms('phi', phi_p2[i], nks)}\t{aus.gms('lambda', lambda_p2[i], nks)}\t{aus.gms('A', A_p2[i], nks)}") + if i != len(phi_p2)-1: + print(f"P0[{i}]: {aus.gms('phi', phi_p0[i], nks)}\t\t\t\t\t\t\t\t{aus.gms('A', A_p0[i], nks)}") + return phi_p2, lambda_p2, A_p2 + + +def gha2(re, phi_p1, lambda_p1, phi_p2, lambda_p2): + """ + Berechnung der 2. Geodätische Hauptaufgabe nach Gauß´schen Mittelbreitenformeln + :param re: Klasse Ellipsoid + :param phi_p1: Breite Punkt 1 + :param lambda_p1: Länge Punkt 1 + :param phi_p2: Breite Punkt 2 + :param lambda_p2: Länge Punkt 2 + :return: Länge der geodätischen Linie, Azimut von P1 nach P2, Azimut von P2 nach P1 + """ + + t = lambda phi: tan(phi) + eta = lambda phi: sqrt(re.e_ ** 2 * cos(phi) ** 2) + + phi_0 = (phi_p1 + phi_p2) / 2 + d_phi = phi_p2 - phi_p1 + d_lambda = lambda_p2 - lambda_p1 + + f_A = lambda phi: (2 + 3*t(phi)**2 + 2*eta(phi)**2) / 24 + f_B = lambda phi: ((t(phi)**2 - 1) * eta(phi)**2) / 8 + # f_C = lambda phi: (t(phi)**2) / 24 + f_D = lambda phi: (1 + eta(phi)**2 - 9 * eta(phi)**2 * t(phi)**2) / 24 + + F1 = lambda phi: d_phi * re.M(phi) * (1 - f_A(phi) * d_lambda ** 2 * cos(phi) ** 2 - + f_B(phi) * d_phi ** 2 / re.V(phi) ** 4) + F2 = lambda phi: d_lambda * re.N(phi) * cos(phi) * (1 - 1 / 24 * d_lambda ** 2 * sin(phi) ** 2 + + f_D(phi) * d_phi ** 2 / re.V(phi) ** 4) + + s = sqrt(F1(phi_0) ** 2 + F2(phi_0) ** 2) + A_0 = arctan(F2(phi_0) / F1(phi_0)) + d_A = d_lambda * sin(phi_0) * (1 + (1 + eta(phi_0) ** 2) / 12 * s ** 2 * sin(A_0) ** 2 / re.N(phi_0) ** 2 + + (3 + 8 * eta(phi_0) ** 2) / 24 * s ** 2 * cos(A_0) ** 2 / re.N(phi_0) ** 2) + + A_p1 = A_0 - d_A / 2 + A_p2 = A_0 + d_A / 2 + A_p2_p1 = A_p2 + pi + + return s, A_p1, A_p2_p1 diff --git a/Numerische_Integration/__init__.py b/Numerische_Integration/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/Numerische_Integration/num_int_runge_kutta.py b/Numerische_Integration/num_int_runge_kutta.py new file mode 100644 index 0000000..d363a58 --- /dev/null +++ b/Numerische_Integration/num_int_runge_kutta.py @@ -0,0 +1,62 @@ +def verfahren(funktionen: list, startwerte: list, weite: float, schritte: int) -> list: + """ + Runge-Kutta-Verfahren für ein beliebiges DGLS + :param funktionen: Liste mit allen Funktionen + :param startwerte: Liste mit allen Startwerten der Variablen + :param weite: gesamte Weite über die integriert werden soll + :param schritte: Anzahl der Schritte über die gesamte Weite + :return: Liste mit Listen für alle Wertepaare + """ + h = weite / schritte + werte = [startwerte] + for i in range(schritte): + + zuschlaege_grob = zuschlaege(funktionen, werte[-1], h) + werte_grob = [werte[-1][j] if j == 0 else werte[-1][j] + zuschlaege_grob[j - 1] + for j in range(len(startwerte))] + + zuschlaege_fein_1 = zuschlaege(funktionen, werte[-1], h / 2) + werte_fein_1 = [werte[-1][j] + h/2 if j == 0 else werte[-1][j]+zuschlaege_fein_1[j-1] + for j in range(len(startwerte))] + + zuschlaege_fein_2 = zuschlaege(funktionen, werte_fein_1, h / 2) + werte_fein_2 = [werte_fein_1[j] + h/2 if j == 0 else werte_fein_1[j]+zuschlaege_fein_2[j-1] + for j in range(len(startwerte))] + + werte_korr = [werte_fein_2[j] if j == 0 else werte_fein_2[j] + 1/15 * (werte_fein_2[j] - werte_grob[j]) + for j in range(len(startwerte))] + + werte.append(werte_korr) + return werte + + +def zuschlaege(funktionen: list, startwerte: list, h: float) -> list: + """ + Berechnung der Zuschläge eines einzelnen Schritts + :param funktionen: Liste mit allen Funktionen + :param startwerte: Liste mit allen Startwerten der Variablen + :param h: Schrittweite + :return: Liste mit Zuschlägen für die einzelnen Variablen + """ + werte = [wert for wert in startwerte] + + k1 = [h * funktion(*werte) for funktion in funktionen] + + werte = [startwerte[i] + (h / 2 if i == 0 else k1[i - 1] / 2) + for i in range(len(startwerte))] + + k2 = [h * funktion(*werte) for funktion in funktionen] + + werte = [startwerte[i] + (h / 2 if i == 0 else k2[i - 1] / 2) + for i in range(len(startwerte))] + + k3 = [h * funktion(*werte) for funktion in funktionen] + + werte = [startwerte[i] + (h if i == 0 else k3[i - 1]) + for i in range(len(startwerte))] + + k4 = [h * funktion(*werte) for funktion in funktionen] + + k_ = [(k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) / 6 for i in range(len(k1))] + + return k_ diff --git a/Numerische_Integration/num_int_trapezformel.py b/Numerische_Integration/num_int_trapezformel.py new file mode 100644 index 0000000..4502d22 --- /dev/null +++ b/Numerische_Integration/num_int_trapezformel.py @@ -0,0 +1,46 @@ +def f(xi, yi): + ys = xi + yi + return ys + + +# Gegeben: +x0 = 0 +y0 = 0 +h = 0.2 +xmax = 0.4 +Ey = 0.0001 + +x = [x0] +y = [[y0]] + +n = 0 +while x[-1] < xmax: + x.append(x[-1]+h) + fn = f(x[n], y[n][-1]) + if n == 0: + y_neu = y[n][-1] + h * fn + else: + y_neu = y[n-1][-1] + 2*h * fn + y.append([y_neu]) + dy = 1 + while dy > Ey: + y_neu = y[n][-1] + h/2 * (fn + f(x[n+1], y[n+1][-1])) + y[-1].append(y_neu) + dy = abs(y[-1][-2]-y[-1][-1]) + n += 1 + +print(x) +print(y) + +werte = [] +for i in range(len(x)): + werte.append((x[i], y[i][-1])) + +for paar in werte: + print(f"({round(paar[0], 5)}, {round(paar[1], 5)})") + +integral = 0 +for i in range(len(werte)-1): + integral += (werte[i+1][1]+werte[i][1])/2 * (werte[i+1][0]-werte[i][0]) + +print(f"Integral = {round(integral, 5)}") diff --git a/Numerische_Integration/rk_DGLS_1Ordnung.py b/Numerische_Integration/rk_DGLS_1Ordnung.py new file mode 100644 index 0000000..e9b8643 --- /dev/null +++ b/Numerische_Integration/rk_DGLS_1Ordnung.py @@ -0,0 +1,7 @@ +import num_int_runge_kutta as rk + +f = lambda ti, xi, yi: yi - ti +g = lambda ti, xi, yi: xi + yi + +funktionswerte = rk.verfahren([f, g], [0, 0, 1], 0.6, 3) +print(funktionswerte) diff --git a/Numerische_Integration/rk_DGL_1Ordnung.py b/Numerische_Integration/rk_DGL_1Ordnung.py new file mode 100644 index 0000000..253cb85 --- /dev/null +++ b/Numerische_Integration/rk_DGL_1Ordnung.py @@ -0,0 +1,6 @@ +import num_int_runge_kutta as rk + +f = lambda xi, yi: xi + yi + +funktionswerte = rk.verfahren([f], [0, 0], 1, 5) +print(funktionswerte) diff --git a/Numerische_Integration/rk_DGL_2Ordnung.py b/Numerische_Integration/rk_DGL_2Ordnung.py new file mode 100644 index 0000000..37a86e4 --- /dev/null +++ b/Numerische_Integration/rk_DGL_2Ordnung.py @@ -0,0 +1,10 @@ +import numpy as np +import num_int_runge_kutta as rk + +f = lambda ti, ui, phii: -4 * np.sin(phii) +g = lambda ti, ui, phii: ui + +funktionswerte = rk.verfahren([f, g], [0, 1, 0], 0.6, 3) + +for wert in funktionswerte: + print(f"t = {round(wert[0],1)}s -> phi = {round(wert[2],5)}, phip = {round(wert[1],5)}, v = {round(2.45 * wert[1],5)}") diff --git a/Numerische_Integration/rk_DGL_2Ordnung_2.py b/Numerische_Integration/rk_DGL_2Ordnung_2.py new file mode 100644 index 0000000..4d8f4ab --- /dev/null +++ b/Numerische_Integration/rk_DGL_2Ordnung_2.py @@ -0,0 +1,7 @@ +import num_int_runge_kutta as rk + +f = lambda xi, yi, ui: ui +g = lambda xi, yi, ui: 4 * ui - 4 * yi + +funktionswerte = rk.verfahren([f, g], [0, 0, 1], 0.6, 2) +print(funktionswerte) diff --git a/Richtungsreduktion.py b/Richtungsreduktion.py new file mode 100644 index 0000000..6b7f031 --- /dev/null +++ b/Richtungsreduktion.py @@ -0,0 +1,14 @@ +from numpy import sqrt, cos, sin +import ellipsoide +import ausgaben as aus +import winkelumrechnungen as wu + + +re = ellipsoide.Ellipsoid.init_name("Bessel") +eta = lambda phi: sqrt(re.e_ ** 2 * cos(phi) ** 2) + +A = wu.gms2rad([327,0,0]) +phi = wu.gms2rad([35,0,0]) +h = 3500 +dA = eta(phi)**2 * h * sin(A)*cos(A) / re.N(phi) +print(aus.gms("dA", dA, 10)) \ No newline at end of file diff --git a/ausgaben.py b/ausgaben.py new file mode 100644 index 0000000..8c3e55f --- /dev/null +++ b/ausgaben.py @@ -0,0 +1,27 @@ +import winkelumrechnungen as wu + + +def xyz(x: float, y: float, z: float, stellen: int) -> str: + """ + Erzeugen eines mehrzeiligen Strings zur Ausgabe von Koordinaten. + :param x: x-Koordinate + :param y: y-Koordinate + :param z: z-Koordinate + :param stellen: Anzahl Nachkommastellen + :return: String zur Ausgabe der Koordinaten + """ + return f"""x = {(round(x,stellen))} m +y = {(round(y,stellen))} m +z = {(round(z,stellen))} m""" + + +def gms(name: str, rad: float, stellen: int) -> str: + """ + Erzeugen eines Strings zur Ausgabe eines Winkel in Grad-Minuten-Sekunden. + :param name: Bezeichnung des Winkels + :param rad: Winkel in Radiant + :param stellen: Anzahl Nachkommastellen + :return: String zur Ausgabe des Winkels + """ + gms = wu.rad2gms(rad) + return f"{name} = {int(gms[0])}° {int(gms[1])}' {round(gms[2],stellen):.{stellen}f}''" diff --git a/ellipsoide.py b/ellipsoide.py new file mode 100644 index 0000000..2de6f6b --- /dev/null +++ b/ellipsoide.py @@ -0,0 +1,83 @@ +import numpy as np +import winkelumrechnungen as wu +import ausgaben as aus + + +class Ellipsoid: + def __init__(self, a: float, b: float): + self.a = a + self.b = b + self.c = a ** 2 / b + self.e = np.sqrt(a ** 2 - b ** 2) / a + self.e_ = np.sqrt(a ** 2 - b ** 2) / b + + @classmethod + def init_name(cls, name: str): + if name == "Bessel": + a = 6377397.15508 + b = 6356078.96290 + return cls(a, b) + elif name == "Hayford": + a = 6378388 + f = 1/297 + b = a - a * f + return cls(a, b) + elif name == "Krassowski": + a = 6378245 + f = 298.3 + b = a - a * f + return cls(a, b) + elif name == "WGS84": + a = 6378137 + f = 298.257223563 + b = a - a * f + return cls(a, b) + + @classmethod + def init_af(cls, a: float, f: float): + b = a - a * f + return cls(a, b) + + V = lambda self, phi: np.sqrt(1 + self.e_ ** 2 * np.cos(phi) ** 2) + M = lambda self, phi: self.c / self.V(phi) ** 3 + N = lambda self, phi: self.c / self.V(phi) + + beta2psi = lambda self, beta: np.arctan(self.a / self.b * np.tan(beta)) + beta2phi = lambda self, beta: np.arctan(self.a ** 2 / self.b ** 2 * np.tan(beta)) + + psi2beta = lambda self, psi: np.arctan(self.b / self.a * np.tan(psi)) + psi2phi = lambda self, psi: np.arctan(self.a / self.b * np.tan(psi)) + + phi2beta = lambda self, phi: np.arctan(self.b ** 2 / self.a ** 2 * np.tan(phi)) + phi2psi = lambda self, phi: np.arctan(self.b / self.a * np.tan(phi)) + + phi2p = lambda self, phi: self.N(phi) * np.cos(phi) + + def ellipsoidische_Koords (self, Eh, Ephi, x, y, z): + p = np.sqrt(x**2+y**2) + print(f"p = {round(p, 5)} m") + + lamb = np.arctan(y/x) + + phi_null = np.arctan(z/p*(1-self.e**2)**-1) + + hi = [0] + phii = [phi_null] + + i = 0 + + while True: + N = self.a*(1-self.e**2*np.sin(phii[i])**2)**(-1/2) + h = p/np.cos(phii[i])-N + phi = np.arctan(z/p*(1-(self.e**2*N)/(N+h))**(-1)) + hi.append(h) + phii.append(phi) + dh = abs(hi[i]-h) + dphi = abs(phii[i]-phi) + i = i+1 + if dh < Eh: + if dphi < Ephi: + break + for i in range(len(phii)): + print(f"P3[{i}]: {aus.gms('phi', phii[i], 5)}\th = {round(hi[i], 5)} m") + return phi, lamb, h diff --git a/hausarbeit.py b/hausarbeit.py new file mode 100644 index 0000000..95089b2 --- /dev/null +++ b/hausarbeit.py @@ -0,0 +1,78 @@ +from numpy import cos, sin, tan +import winkelumrechnungen as wu +import s_ellipse as s_ell +import ausgaben as aus +from ellipsoide import Ellipsoid +import Numerische_Integration.num_int_runge_kutta as rk +import GHA.gauss as gauss +import GHA.bessel as bessel + +matrikelnummer = "6044051" +print(f"Matrikelnummer: {matrikelnummer}") + +m4 = matrikelnummer[-4] +m3 = matrikelnummer[-3] +m2 = matrikelnummer[-2] +m1 = matrikelnummer[-1] +print(f"m1={m1}\tm2={m2}\tm3={m3}\tm4={m4}") + +re = Ellipsoid.init_name("Bessel") +nks = 3 + +print(f"\na = {re.a} m\nb = {re.b} m\nc = {re.c} m\ne' = {re.e_}") + +print("\n\nAufgabe 1 (via Reihenentwicklung)") +phi1 = re.beta2phi(wu.gms2rad([int("1"+m1), int("1"+m2), int("1"+m3)])) +print(f"{aus.gms('phi_P1', phi1, 5)}") +s = s_ell.reihenentwicklung(re.e_, re.c, phi1) +print(f'\ns_P1 = {round(s,nks)} m') +# s_poly = s_ell.polyapp_tscheby_bessel(phi1) +# print(f'Meridianbogenlänge: {round(s_poly,nks)} m (Polynomapproximation)') + +print("\n\nAufgabe 2 (via Gauß´schen Mittelbreitenformeln)") + +lambda1 = wu.gms2rad([int("1"+m3), int("1"+m2), int("1"+m4)]) +A12 = wu.gms2rad([90, 0, 0]) +s12 = float("1"+m4+m3+m2+"."+m1+"0") + +# print("\nvia Runge-Kutta-Verfahren") +# re = Ellipsoid.init_name("Bessel") +# f_phi = lambda s, phi, lam, A: cos(A) * re.V(phi) ** 3 / re.c +# f_lam = lambda s, phi, lam, A: sin(A) * re.V(phi) / (cos(phi) * re.c) +# f_A = lambda s, phi, lam, A: tan(phi) * sin(A) * re.V(phi) / re.c +# +# funktionswerte = rk.verfahren([f_phi, f_lam, f_A], +# [0, phi1, lambda1, A12], +# s12, 1) +# +# for s, phi, lam, A in funktionswerte: +# print(f"{s} m, {aus.gms('phi', phi, nks)}, {aus.gms('lambda', lam, nks)}, {aus.gms('A', A, nks)}") + +# print("via Gauß´schen Mittelbreitenformeln") +phi_p2, lambda_p2, A_p2 = gauss.gha1(Ellipsoid.init_name("Bessel"), + phi_p1=phi1, + lambda_p1=lambda1, + A_p1=A12, + s=s12, eps=wu.gms2rad([1*10**-8, 0, 00])) + +print(f"\nP2: {aus.gms('phi', phi_p2[-1], nks)}\t\t{aus.gms('lambda', lambda_p2[-1], nks)}\t{aus.gms('A', A_p2[-1], nks)}") + +# print("\nvia Verfahren nach Bessel") +# phi_p2, lambda_p2, A_p2 = bessel.gha1(Ellipsoid.init_name("Bessel"), +# phi_p1=phi1, +# lambda_p1=lambda1, +# A_p1=A12, +# s=s12) +# print(f"P2: {aus.gms('phi', phi_p2, nks)}, {aus.gms('lambda', lambda_p2, nks)}, {aus.gms('A', A_p2, nks)}") + +print("\n\nAufgabe 3") +p = re.phi2p(phi_p2[-1]) +print(f"p = {round(p,2)} m") + + +print("\n\nAufgabe 4") +x = float("4308"+m3+"94.556") +y = float("1214"+m2+"88.242") +z = float("4529"+m4+"03.878") +phi_p3, lambda_p3, h_p3 = re.ellipsoidische_Koords(0.001, wu.gms2rad([0, 0, 0.001]), x, y, z) +print(f"\nP3: {aus.gms('phi', phi_p3, nks)}, {aus.gms('lambda', lambda_p3, 5)}, h = {round(h_p3,nks)} m") diff --git a/phi_ellipse.py b/phi_ellipse.py new file mode 100644 index 0000000..29b1a37 --- /dev/null +++ b/phi_ellipse.py @@ -0,0 +1,55 @@ +import winkelumrechnungen as wu + + +def polyapp_tscheby_hayford(s: float) -> float: + """ + Berechnung der ellipsoidisch geodätischen Breite. + Polynomapproximation mittels Tschebyscheff-Polynomen. + Auf dem Hayford-Ellipsoid. + :param s: Strecke auf einer Ellipse vom Äquator aus + :type s: float + :return: ellipsoidisch geodätische Breite + :rtype: float + """ + c0 = 1 + c1 = -0.00837809325 + c2 = 0.00428127367 + c3 = -0.00114523986 + c4 = 0.00023219707 + c5 = -0.00004421222 + c6 = 0.00000570244 + alpha = wu.gms2rad([0, 0, 325643.97199]) + s90 = 10002288.2990 + + xi = s/s90 + + phi = alpha * xi * (c0*xi**(2*0) + c1*xi**(2*1) + c2*xi**(2*2) + c3*xi**(2*3) + + c4*xi**(2*4) + c5*xi**(2*5) + c6*xi**(2*6)) + return phi + + +def polyapp_tscheby_bessel(s: float) -> float: + """ + Berechnung der ellipsoidisch geodätischen Breite. + Polynomapproximation mittels Tschebyscheff-Polynomen. + Auf dem Bessel-Ellipsoid. + :param s: Strecke auf einer Ellipse vom Äquator aus + :type s: float + :return: ellipsoidisch geodätische Breite + :rtype: float + """ + c0 = 1 + c1 = -0.00831729565 + c2 = 0.00424914906 + c3 = -0.00113566119 + c4 = 0.00022976983 + c5 = -0.00004363980 + c6 = 0.00000562025 + alpha = wu.gms2rad([0, 0, 325632.08677]) + s90 = 10000855.7644 + + xi = s/s90 + + phi = alpha * xi * (c0*xi**(2*0) + c1*xi**(2*1) + c2*xi**(2*2) + c3*xi**(2*3) + + c4*xi**(2*4) + c5*xi**(2*5) + c6*xi**(2*6)) + return phi diff --git a/s_ellipse.py b/s_ellipse.py new file mode 100644 index 0000000..b9ab052 --- /dev/null +++ b/s_ellipse.py @@ -0,0 +1,78 @@ +import numpy as np + + +def reihenentwicklung(es: float, c: float, phi: float) -> float: + """ + Berechnung der Strecke auf einer Ellipse. + Reihenentwicklung. + :param es: zweite numerische Exzentrizität + :type es: float + :param c: Polkrümmungshalbmesser + :type c: float + :param phi: ellipsoidisch geodästische Breite in Radiant + :type phi: float + :return: Strecke auf einer Ellipse vom Äquator aus + :rtype: float + """ + Ass = 1 - 3/4*es**2 + 45/64*es**4 - 175/256*es**6 + 11025/16384*es**8 + Bss = - 3/4*es**2 + 15/16*es**4 - 525/512*es**6 + 2205/2048*es**8 + Css = 15/64*es**4 - 105/256*es**6 + 2205/4096*es**8 + Dss = - 35/512*es**6 + 315/2048*es**8 + print(f"A'' = {round(Ass, 10):.10f}\nB'' = {round(Bss, 10):.10f}\nC'' = {round(Css, 10):.10f}\nD'' = {round(Dss, 10):.10f}") + + s = c * (Ass*phi + 1/2*Bss * np.sin(2*phi) + 1/4*Css * np.sin(4*phi) + 1/6*Dss * np.sin(6*phi)) + return s + + +def polyapp_tscheby_hayford(phi: float) -> float: + """ + Berechnung der Strecke auf einer Ellipse. + Polynomapproximation mittels Tschebyscheff-Polynomen. + Auf dem Hayford-Ellipsoid. + :param phi: ellipsoidisch geodästische Breite in Radiant + :type phi: float + :return: Strecke auf einer Ellipse vom Äquator aus + :rtype: float + """ + c1s = 0.00829376218 + c2s = -0.00398963425 + c3s = 0.00084200710 + c4s = -0.0000648906 + c5s = -0.00001075680 + c6s = 0.00000396474 + c7s = -0.00000046347 + alpha = 9951793.0123 + + xi = 2 / np.pi * phi + xis = xi ^ 2 + + s = alpha * xi * (1 + c1s * xis + c2s * xis**2 + c3s * xis**3 + + c4s * xis**4 + c5s * xis**5 + c6s * xis**6 + c7s * xis**7) + return s + + +def polyapp_tscheby_bessel(phi: float) -> float: + """ + Berechnung der Strecke auf einer Ellipse. + Polynomapproximation mittels Tschebyscheff-Polynomen. + Auf dem Bessel-Ellipsoid. + :param phi: ellipsoidisch geodästische Breite in Radiant + :type phi: float + :return: Strecke auf einer Ellipse vom Äquator aus + :rtype: float + """ + c1s = 0.00823417717 + c2s = -0.00396170744 + c3s = 0.00083680249 + c4s = -0.00006488462 + c5s = -0.00001053242 + c6s = 0.00000390854 + c7s = -0.00000045768 + alpha = 9950730.8876 + + xi = 2 / np.pi * phi + xis = xi ** 2 + + s = alpha * xi * (1 + c1s * xis + c2s * xis**2 + c3s * xis**3 + + c4s * xis**4 + c5s * xis**5 + c6s * xis**6 + c7s * xis**7) + return s diff --git a/winkelumrechnungen.py b/winkelumrechnungen.py new file mode 100644 index 0000000..0b0b7f8 --- /dev/null +++ b/winkelumrechnungen.py @@ -0,0 +1,158 @@ +from numpy import * + + +def deg2gms(deg: float) -> list: + """ + Umrechnung von Grad in Grad-Minuten-Sekunden + :param deg: Winkel in Grad + :type deg: float + :return: Winkel in Grad-Minuten-Sekunden + :rtype: list + """ + gra = deg // 1 + min = gra % 1 + gra = gra // 1 + min *= 60 + sek = min % 1 + min = min // 1 + sek *= 60 + return [gra, min, sek] + + +def deg2gra(deg: float) -> float: + """ + Umrechnung von Grad in Gon + :param deg: Winkel in Grad + :type deg: float + :return: Winkel in Gon + :rtype: float + """ + return deg * 10/9 + + +def deg2rad(deg: float) -> float: + """ + Umrechnung von Grad in Radiant + :param deg: Winkel in Grad + :type deg: float + :return: Winkel in Radiant + :rtype: float + """ + return deg * pi / 180 + + +def gra2gms(gra: float) -> list: + """ + Umrechnung von Gon in Grad-Minuten-Sekunden + :param gra: Winkel in Gon + :type gra: float + :return: Winkel in Grad-Minuten-Sekunden + :rtype: list + """ + deg = gra2deg(gra) + gra = deg // 1 + min = gra % 1 + gra = gra // 1 + min *= 60 + sek = min % 1 + min = min // 1 + sek *= 60 + return [gra, min, sek] + + +def gra2rad(gra: float) -> float: + """ + Umrechnung von Gon in Radiant + :param gra: Winkel in Gon + :type gra: float + :return: Winkel in Radiant + :rtype: float + """ + return gra * pi / 200 + + +def gra2deg(gra: float) -> float: + """ + Umrechnung von Gon in Grad + :param gra: Winkel in Gon + :type gra: float + :return: Winkel in Grad + :rtype: float + """ + return gra * 9/10 + + +def rad2deg(rad: float) -> float: + """ + Umrechnung von Radiant in Grad + :param rad: Winkel in Radiant + :type rad: float + :return: Winkel in Grad + :rtype: float + """ + return rad * 180 / pi + + +def rad2gra(rad: float) -> float: + """ + Umrechnung von Radiant in Gon + :param rad: Winkel in Radiant + :type rad: float + :return: Winkel in Gon + :rtype: float + """ + return rad * 200 / pi + + +def rad2gms(rad: float) -> list: + """ + Umrechnung von Radiant in Grad-Minuten-Sekunden + :param rad: Winkel in Radiant + :type rad: float + :return: Winkel in Grad-Minuten-Sekunden + :rtype: list + """ + deg = rad2deg(rad) + min = deg % 1 + gra = deg // 1 + min *= 60 + sek = min % 1 + min = min // 1 + sek *= 60 + return [gra, min, sek] + + +def gms2rad(gms: list) -> float: + """ + Umrechnung von Grad-Minuten-Sekunden in Radiant + :param gms: Winkel in Grad-Minuten-Sekunden + :type gms: list + :return: Winkel in Radiant + :rtype: float + """ + deg = gms[0] + gms[1] / 60 + gms[2] / 3600 + return deg2rad(deg) + + +def gms2deg(gms: list) -> float: + """ + Umrechnung von Grad-Minuten-Sekunden in Grad + :param gms: Winkel in Grad-Minuten-Sekunden + :type gms: list + :return: Winkel in Grad + :rtype: float + """ + deg = gms[0] + gms[1] / 60 + gms[2] / 3600 + return deg + + +def gms2gra(gms: list) -> float: + """ + Umrechnung von Grad-Minuten-Sekunden in Gon + :param gms: Winkel in Grad-Minuten-Sekunden + :type gms: list + :return: Winkel in Gon + :rtype: float + """ + deg = gms[0] + gms[1] / 60 + gms[2] / 3600 + return deg2gra(deg)