diff --git a/GHA_biaxial/bessel.py b/GHA_biaxial/bessel.py index 931999e..898124f 100644 --- a/GHA_biaxial/bessel.py +++ b/GHA_biaxial/bessel.py @@ -1,25 +1,26 @@ from numpy import * import scipy as sp +from ellipsoide import EllipsoidBiaxial +from typing import Tuple -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)) +def gha1(re: EllipsoidBiaxial, phi0: float, lamb0: float, alpha0:float, s: float) -> Tuple[float, float, float]: + psi0 = re.phi2psi(phi0) + clairant = arcsin(cos(psi0) * sin(alpha0)) + sigma0 = arcsin(sin(psi0) / cos(clairant)) - 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] + sqrt_sigma = lambda sigma: sqrt(1 + re.e_ ** 2 * cos(clairant) ** 2 * sin(sigma) ** 2) + int_sqrt_sigma = lambda sigma: sp.integrate.quad(sqrt_sigma, sigma0, sigma)[0] - f_sigma_p2i = lambda sigma_p2i: (int_sqrt_sigma(sigma_p2i) - s / re.b) + f_sigma1_i = lambda sigma1_i: (int_sqrt_sigma(sigma1_i) - 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)) + sigma1_0 = sigma0 + s / re.a + sigma1 = sp.optimize.newton(f_sigma1_i, sigma1_0) + psi1 = arcsin(cos(clairant) * sin(sigma1)) + phi1 = re.psi2phi(psi1) + alpha1 = arcsin(sin(clairant) / cos(psi1)) - 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 + f_d_lambda = lambda sigma: sin(clairant) * sqrt_sigma(sigma) / (1 - cos(clairant) ** 2 * sin(sigma) ** 2) + d_lambda = sqrt(1-re.e**2) * sp.integrate.quad(f_d_lambda, sigma0, sigma1)[0] + lamb1 = lamb0 + d_lambda + return phi1, lamb1, alpha1 diff --git a/GHA_biaxial/gauss.py b/GHA_biaxial/gauss.py index e89c264..834fb4d 100644 --- a/GHA_biaxial/gauss.py +++ b/GHA_biaxial/gauss.py @@ -1,89 +1,84 @@ from numpy import sin, cos, pi, sqrt, tan, arcsin, arccos, arctan import ausgaben as aus +from ellipsoide import EllipsoidBiaxial +from typing import Tuple -def gha1(re, phi_p1, lambda_p1, A_p1, s, eps): +def gha1(re: EllipsoidBiaxial, phi0: float, lamb0: float, alpha0: float, s: float, eps: float = 1e-12) -> Tuple[float, float, float]: """ 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 phi0: Breite Punkt 0 + :param lamb0: Länge Punkt 0 + :param alpha0: Azimut der geodätischen Linie in Punkt 1 + :param s: Strecke zu Punkt 1 :param eps: Abbruchkriterium für Winkelgrößen - :return: Breite, Länge, Azimut von Punkt 2 + :return: Breite, Länge, Azimut von Punkt 21 """ - 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 + F1 = lambda alpha, phi, s: 1 + (2 + 3 * t(phi) ** 2 + 2 * eta(phi) ** 2) / (24 * re.N(phi) ** 2) * sin(alpha) ** 2 * s ** 2 \ + + (t(phi)**2 - 1) * eta(phi) ** 2 / (8 * re.N(phi) ** 2) * cos(alpha) ** 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 + F2 = lambda alpha, phi, s: 1 + t(phi) ** 2 / (24 * re.N(phi) ** 2) * sin(alpha) ** 2 * s ** 2 \ + - (1 + eta(phi)**2 - 9 * eta(phi)**2 * t(phi)**2) / (24 * re.N(phi) ** 2) * cos(alpha) ** 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 + F3 = lambda alpha, phi, s: 1 + (1 + eta(phi) ** 2) / (12 * re.N(phi) ** 2) * sin(alpha) ** 2 * s ** 2 \ + + (3 + 8 * eta(phi)**2) / (24 * re.N(phi) ** 2) * cos(alpha) ** 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) + phi1_i = lambda alpha, phi: phi0 + cos(alpha) / re.M(phi) * s * F1(alpha, phi, s) + lamb1_i = lambda alpha, phi: lamb0 + sin(alpha) / (re.N(phi) * cos(phi)) * s * F2(alpha, phi, s) + alpha1_i = lambda alpha, phi: alpha0 + sin(alpha) * tan(phi) / re.N(phi) * s * F3(alpha, phi, s) - phi_p0_i = lambda phi2: (phi_p1 + phi2) / 2 - A_p1_i = lambda A2: (A_p1 + A2) / 2 + phi_m_i = lambda phi1: (phi0 + phi1) / 2 + alpha_m_i = lambda alpha1: (alpha0 + alpha1) / 2 - phi_p0 = [] - A_p0 = [] - phi_p2 = [] - lambda_p2 = [] - A_p2 = [] + phi_m = [] + alpha_m = [] + phi1 = [] + lamb1 = [] + alpha1 = [] - # 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) + # 1. Näherung für P1 + phi1.append(phi0 + cos(alpha0) / re.M(phi0) * s) + lamb1.append(lamb0 + sin(alpha0) / (re.N(phi0) * cos(phi0)) * s) + alpha1.append(alpha0 + sin(alpha0) * tan(phi0) / re.N(phi0) * 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])) + # Berechnug P_m durch Mittelbildung + phi_m.append(phi_m_i(phi1[-1])) + alpha_m.append(alpha_m_i(alpha1[-1])) + # Berechnung P1 + phi1.append(phi1_i(alpha_m[-1], phi_m[-1])) + lamb1.append(lamb1_i(alpha_m[-1], phi_m[-1])) + alpha1.append(alpha1_i(alpha_m[-1], phi_m[-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: + if abs(phi1[-2] - phi1[-1]) < eps and \ + abs(lamb1[-2] - lamb1[-1]) < eps and \ + abs(alpha1[-2] - alpha1[-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 + return phi1[-1], lamb1[-1], alpha1[-1] -def gha2(re, phi_p1, lambda_p1, phi_p2, lambda_p2): +def gha2(re: EllipsoidBiaxial, phi0: float, lamb0: float, phi1: float, lamb1: float): """ 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 + :param phi0: Breite Punkt 1 + :param lamb0: Länge Punkt 1 + :param phi1: Breite Punkt 2 + :param lamb1: 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 + phi_0 = (phi0 + phi1) / 2 + d_phi = phi1 - phi0 + d_lambda = lamb1 - lamb0 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 @@ -100,8 +95,7 @@ def gha2(re, phi_p1, lambda_p1, phi_p2, lambda_p2): 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 + alpha0 = A_0 - d_A / 2 + alpha1 = A_0 + d_A / 2 - return s, A_p1, A_p2_p1 + return alpha0, alpha1, s diff --git a/GHA_biaxial/rk.py b/GHA_biaxial/rk.py index 85101c8..aca6631 100644 --- a/GHA_biaxial/rk.py +++ b/GHA_biaxial/rk.py @@ -3,19 +3,18 @@ from numpy import sin, cos, tan import winkelumrechnungen as wu from ellipsoide import EllipsoidBiaxial import numpy as np +from numpy.typing import NDArray -def gha1(re: EllipsoidBiaxial, x0, y0, z0, A0, s, num): - phi0, lamb0, h0 = re.cart2ell(0.001, wu.gms2rad([0, 0, 0.001]), x0, y0, z0) - +def gha1(re: EllipsoidBiaxial, phi0: float, lamb0: float, alpha0: float, s: float, num: int) -> Tuple[float, float, float]: def buildODE(): def ODE(s, v): phi, lam, A = v - dphi = cos(A) * re.V(phi) ** 3 / re.c - dlam = sin(A) * re.V(phi) / (cos(phi) * re.c) - dA = tan(phi) * sin(A) * re.V(phi) / re.c + V = re.V(phi) + dphi = cos(A) * V ** 3 / re.c + dlam = sin(A) * V / (cos(phi) * re.c) + dA = tan(phi) * sin(A) * V / re.c return np.array([dphi, dlam, dA]) return ODE - _, funktionswerte = rk.rk4(buildODE(), 0, np.array([phi0, lamb0, A0]), s, num) - coords = re.ell2cart(funktionswerte[-1][0], funktionswerte[-1][1], h0) - return coords + _, funktionswerte = rk.rk4(buildODE(), 0, np.array([phi0, lamb0, alpha0]), s, num) + return funktionswerte[-1][0], funktionswerte[-1][1], funktionswerte[-1][2] diff --git a/ellipsoide.py b/ellipsoide.py index 615bd42..96b59e7 100644 --- a/ellipsoide.py +++ b/ellipsoide.py @@ -46,24 +46,25 @@ class EllipsoidBiaxial: M = lambda self, phi: self.c / self.V(phi) ** 3 N = lambda self, phi: self.c / self.V(phi) - beta2psi = lambda self, beta: arctan(self.a / self.b * np.tan(beta)) - beta2phi = lambda self, beta: arctan(self.a ** 2 / self.b ** 2 * np.tan(beta)) + beta2psi = lambda self, beta: np.arctan2(self.a * np.sin(beta), self.b * np.cos(beta)) + beta2phi = lambda self, beta: np.arctan2(self.a ** 2 * np.sin(beta), self.b ** 2 * np.cos(beta)) - psi2beta = lambda self, psi: arctan(self.b / self.a * np.tan(psi)) - psi2phi = lambda self, psi: arctan(self.a / self.b * np.tan(psi)) + psi2beta = lambda self, psi: np.arctan2(self.b * np.sin(psi), self.a * np.cos(psi)) + psi2phi = lambda self, psi: np.arctan2(self.a * np.sin(psi), self.b * np.cos(psi)) - phi2beta = lambda self, phi: arctan(self.b ** 2 / self.a ** 2 * np.tan(phi)) - phi2psi = lambda self, phi: arctan(self.b / self.a * np.tan(phi)) + phi2beta = lambda self, phi: np.arctan2(self.b**2 * np.sin(phi), self.a**2 * np.cos(phi)) + phi2psi = lambda self, phi: np.arctan2(self.b * np.sin(phi), self.a * np.cos(phi)) phi2p = lambda self, phi: self.N(phi) * cos(phi) - def cart2ell(self, Eh, Ephi, x, y, z): + def bi_cart2ell(self, point: NDArrayself, Eh: float = 0.001, Ephi: float = wu.gms2rad([0, 0, 0.001])) -> Tuple[float, float, float]: + x, y, z = point + + lamb = arctan2(y, x) + p = sqrt(x**2+y**2) - # print(f"p = {round(p, 5)} m") - lamb = arctan(y/x) - - phi_null = arctan(z/p*(1-self.e**2)**-1) + phi_null = arctan2(z, p*(1 - self.e**2)) hi = [0] phii = [phi_null] @@ -71,9 +72,9 @@ class EllipsoidBiaxial: i = 0 while True: - N = self.a*(1-self.e**2*sin(phii[i])**2)**(-1/2) - h = p/cos(phii[i])-N - phi = arctan(z/p*(1-(self.e**2*N)/(N+h))**(-1)) + N = self.a / sqrt(1 - self.e**2 * sin(phii[i])**2) + h = p / cos(phii[i]) - N + phi = arctan2(z, p * (1-(self.e**2*N) / (N+h))) hi.append(h) phii.append(phi) dh = abs(hi[i]-h) @@ -82,18 +83,15 @@ class EllipsoidBiaxial: 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") - pass return phi, lamb, h - def ell2cart(self, phi, lamb, h): + def bi_ell2cart(self, phi: float, lamb: float, h: float) -> NDArray: W = sqrt(1 - self.e**2 * sin(phi)**2) N = self.a / W x = (N+h) * cos(phi) * cos(lamb) y = (N+h) * cos(phi) * sin(lamb) - z = (N * (1-self.e**2) + h) * sin(lamb) - return x, y, z + z = (N * (1-self.e**2) + h) * sin(phi) + return np.array([x, y, z]) class EllipsoidTriaxial: def __init__(self, ax: float, ay: float, b: float):