GHA triaxial nach Panou Korakitis 2019 läuft

bei ax=ay Fehler
irgendwas bei gha1 biaxial noch falsch, evtl Umrechnung

GHA/rk.py

@@ -0,0 +1,17 @@
+import Numerische_Integration.num_int_runge_kutta as rk
+from numpy import sin, cos, tan
+import winkelumrechnungen as wu
+from ellipsoide import EllipsoidBiaxial
+
+def gha1(re, x0, y0, z0, A0, s, num):
+    phi0, lamb0, h0 = re.cart2ell(0.001, wu.gms2rad([0, 0, 0.001]), x0, y0, z0)
+
+    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, phi0, lamb0, A0],
+                                  s, num)
+    coords = re.ell2cart(funktionswerte[-1][1], funktionswerte[-1][2], h0)
+    return coords

GHA_triaxial/panou.py

@@ -0,0 +1,62 @@
+import numpy as np
+import ellipsoide
+import Numerische_Integration.num_int_runge_kutta as rk
+import winkelumrechnungen as wu
+import ausgaben as aus
+import GHA.rk as ghark
+
+# Panou, Korakitits 2019
+
+def gha1(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, num):
+    H = x**2 + y**2 / (1-ell.ee**2)**2 + z**2/(1-ell.ex**2)**2
+
+    n = np.array([x/np.sqrt(H), y/((1-ell.ee**2)*np.sqrt(H)), z/((1-ell.ex**2)*np.sqrt(H))])
+
+    beta, lamb, u = ell.cart2ell(x, y, z)
+    B = ell.Ex**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(beta)**2
+    L = ell.Ex**2 - ell.Ee**2 * np.cos(lamb)**2
+
+    c1 = x**2 + y**2 + z**2 - (ell.ax**2 + ell.ay**2 + ell.b**2)
+    c0 = (ell.ax**2*ell.ay**2 + ell.ax**2*ell.b**2+ell.ay**2*ell.b**2 -
+          (ell.ay**2+ell.b**2)*x**2 - (ell.ax**2+ell.b**2)*y**2 - (ell.ax**2+ell.ay**2)*z**2)
+    t2 = (-c1 + np.sqrt(c1**2 - 4*c0)) / 2
+    t1 = c0 / t2
+
+    F = ell.Ey**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(lamb)**2
+    p1 = -np.sqrt(L/(F*t2)) * ell.ax/ell.Ex * np.sqrt(B) * np.sin(lamb)
+    p2 = np.sqrt(L/(F*t2)) * ell.ay * np.cos(beta) * np.cos(lamb)
+    p3 = 1 / np.sqrt(F*t2) * (ell.b*ell.Ee**2)/(2*ell.Ex) * np.sin(beta) * np.sin(2*lamb)
+    p = np.array([p1, p2, p3])
+    q = np.cross(n, p)
+
+    dxds0 = p[0] * np.sin(alpha0) + q[0] * np.cos(alpha0)
+    dyds0 = p[1] * np.sin(alpha0) + q[1] * np.cos(alpha0)
+    dzds0 = p[2] * np.sin(alpha0) + q[2] * np.cos(alpha0)
+
+    h = lambda dxds, dyds, dzds: dxds**2 + 1/(1-ell.ee**2)*dyds**2 + 1/(1-ell.ex**2)*dzds**2
+
+    f1 = lambda s, x, dxds, y, dyds, z, dzds: dxds
+    f2 = lambda s, x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / H * x
+    f3 = lambda s, x, dxds, y, dyds, z, dzds: dyds
+    f4 = lambda s, x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / H * y/(1-ell.ee**2)
+    f5 = lambda s, x, dxds, y, dyds, z, dzds: dzds
+    f6 = lambda s, x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / H * z/(1-ell.ex**2)
+
+    funktionswerte = rk.verfahren([f1, f2, f3, f4, f5, f6], [0, x, dxds0, y, dyds0, z, dzds0], s, num)
+    return funktionswerte
+
+
+if __name__ == "__main__":
+    ell = ellipsoide.EllipsoidTriaxial.init_name("Eitschberger1978")
+    ellbi = ellipsoide.EllipsoidTriaxial.init_name("Bessel-biaxial")
+    re = ellipsoide.EllipsoidBiaxial.init_name("Bessel")
+    x0 = 5672455.1954766
+    y0 = 2698193.7242382686
+    z0 = 1103177.6450055107
+    alpha0 = wu.gms2rad([20, 0, 0])
+    s = 10000
+    num = 10000
+    werteTri = gha1(ellbi, x0, y0, z0, alpha0, s, num)
+    print(aus.xyz(werteTri[-1][1], werteTri[-1][3], werteTri[-1][5], 8))
+    werteBi = ghark.gha1(re, x0, y0, z0, alpha0, s, num)
+    print(aus.xyz(werteBi[0], werteBi[1], werteBi[2], 8))

ellipsoide.py

@@ -54,9 +54,9 @@ class EllipsoidBiaxial:
 
     phi2p = lambda self, phi: self.N(phi) * np.cos(phi)
 
-    def ellipsoidische_Koords (self, Eh, Ephi, x, y, z):
+    def cart2ell(self, Eh, Ephi, x, y, z):
         p = np.sqrt(x**2+y**2)
-        print(f"p = {round(p, 5)} m")
+        # print(f"p = {round(p, 5)} m")
 
         lamb = np.arctan(y/x)
 
@@ -80,9 +80,18 @@ class EllipsoidBiaxial:
                 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")
+            # 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):
+        W = np.sqrt(1 - self.e**2 * np.sin(phi)**2)
+        N = self.a / W
+        x = (N+h) * np.cos(phi) * np.cos(lamb)
+        y = (N+h) * np.cos(phi) * np.sin(lamb)
+        z = (N * (1-self.e**2) + h) * np.sin(lamb)
+        return x, y, z
+
 class EllipsoidTriaxial:
     def __init__(self, ax: float, ay: float, b: float):
         self.ax = ax
@@ -94,6 +103,9 @@ class EllipsoidTriaxial:
         self.ex_ = np.sqrt((self.ax**2 - self.b**2) / self.b**2)
         self.ey_ = np.sqrt((self.ay**2 - self.b**2) / self.b**2)
         self.ee_ = np.sqrt((self.ax**2 - self.ay**2) / self.ay**2)
+        self.Ex = np.sqrt(self.ax**2 - self.b**2)
+        self.Ey = np.sqrt(self.ay**2 - self.b**2)
+        self.Ee = np.sqrt(self.ax**2 - self.ay**2)
 
     @classmethod
     def init_name(cls, name: str):
@@ -122,6 +134,11 @@ class EllipsoidTriaxial:
             ay = 6378105
             b = 6356754
             return cls(ax, ay, b)
+        elif name == "Bessel-biaxial":
+            ax = 6377397.155085
+            ay = 6377397.15508
+            b = 6356078.96290
+            return cls(ax, ay, b)
 
     def ell2cart(self, beta, lamb, u):
         """
@@ -233,9 +250,9 @@ if __name__ == "__main__":
     stellen = 20
     geod1 = ellips.cart2geod("ligas1", 5712200, 2663400, 1106000)
     print(aus.gms("phi", geod1[0], stellen), aus.gms("lambda", geod1[1], stellen), "h =", geod1[2])
-    geod2 = ellips.cart2geod("ligas2", 5712216.95426783, 2663487.024865021, 1106098.8415910944)
+    geod2 = ellips.cart2geod("ligas2", 5712200, 2663400, 1106000)
     print(aus.gms("phi", geod2[0], stellen), aus.gms("lambda", geod2[1], stellen), "h =", geod2[2])
-    geod3 = ellips.cart2geod("ligas3", 5712216.95426783, 2663487.024865021, 1106098.8415910944)
+    geod3 = ellips.cart2geod("ligas3", 5712200, 2663400, 1106000)
     print(aus.gms("phi", geod3[0], stellen), aus.gms("lambda", geod3[1], stellen), "h =", geod3[2])
     cart1 = ellips.geod2cart(geod1[0], geod1[1], geod1[2])
     print(aus.xyz(cart1[0], cart1[1], cart1[2], 10))
@@ -243,4 +260,7 @@ if __name__ == "__main__":
     print(aus.xyz(cart2[0], cart2[1], cart2[2], 10))
     cart3 = ellips.geod2cart(geod3[0], geod3[1], geod3[2])
     print(aus.xyz(cart3[0], cart3[1], cart3[2], 10))
+
+    test_cart = ellips.geod2cart(0.175, 0.444, 100)
+    print(aus.xyz(test_cart[0], test_cart[1], test_cart[2], 10))
     pass

hausarbeit.py

@@ -74,5 +74,5 @@ 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)
+phi_p3, lambda_p3, h_p3 = re.cart2ell(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")