Berechnungen Biaxial

This commit is contained in:
2026-01-13 16:15:00 +01:00
parent 4b348533bb
commit 9dfb959eb0
4 changed files with 94 additions and 102 deletions

View File

@@ -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

View File

@@ -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

View File

@@ -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]