diff --git a/GHA_triaxial/numeric_examples_panou.py b/GHA_triaxial/numeric_examples_panou.py new file mode 100644 index 0000000..393d03a --- /dev/null +++ b/GHA_triaxial/numeric_examples_panou.py @@ -0,0 +1,111 @@ +import winkelumrechnungen as wu + +table1 = [ + (wu.deg2rad(0), wu.deg2rad(0), wu.deg2rad(0), wu.deg2rad(90), 1.00000000000, + wu.gms2rad([90,0,0.0000]), wu.gms2rad([90,0,0.0000]), 10018754.9569), + + (wu.deg2rad(1), wu.deg2rad(0), wu.deg2rad(-80), wu.deg2rad(5), 0.05883743460, + wu.gms2rad([179,7,12.2719]), wu.gms2rad([174,40,13.8487]), 8947130.7221), + + (wu.deg2rad(5), wu.deg2rad(0), wu.deg2rad(-60), wu.deg2rad(40), 0.34128138370, + wu.gms2rad([160,13,24.5001]), wu.gms2rad([137,26,47.0036]), 8004762.4330), + + (wu.deg2rad(30), wu.deg2rad(0), wu.deg2rad(-30), wu.deg2rad(175), 0.86632464962, + wu.gms2rad([91,7,30.9337]), wu.gms2rad([91,7,30.8672]), 19547128.7971), + + (wu.deg2rad(60), wu.deg2rad(0), wu.deg2rad(60), wu.deg2rad(175), 0.06207487624, + wu.gms2rad([2,52,26.2393]), wu.gms2rad([177,4,13.6373]), 6705715.1610), + + (wu.deg2rad(75), wu.deg2rad(0), wu.deg2rad(80), wu.deg2rad(120), 0.11708984898, + wu.gms2rad([23,20,34.7823]), wu.gms2rad([140,55,32.6385]), 2482501.2608), + + (wu.deg2rad(80), wu.deg2rad(0), wu.deg2rad(60), wu.deg2rad(90), 0.17478427424, + wu.gms2rad([72,26,50.4024]), wu.gms2rad([159,38,30.3547]), 3519745.1283) +] + +table2 = [ + (wu.deg2rad(0), wu.deg2rad(-90), wu.deg2rad(0), wu.deg2rad(89.5), 1.00000000000, + wu.gms2rad([90,0,0.0000]), wu.gms2rad([90,0,0.0000]), 19981849.8629), + + (wu.deg2rad(1), wu.deg2rad(-90), wu.deg2rad(1), wu.deg2rad(89.5), 0.18979826428, + wu.gms2rad([10,56,33.6952]), wu.gms2rad([169,3,26.4359]), 19776667.0342), + + (wu.deg2rad(5), wu.deg2rad(-90), wu.deg2rad(5), wu.deg2rad(89), 0.09398403161, + wu.gms2rad([5,24,48.3899]), wu.gms2rad([174,35,12.6880]), 18889165.0873), + + (wu.deg2rad(30), wu.deg2rad(-90), wu.deg2rad(30), wu.deg2rad(86), 0.06004022935, + wu.gms2rad([3,58,23.8038]), wu.gms2rad([176,2,7.2825]), 13331814.6078), + + (wu.deg2rad(60), wu.deg2rad(-90), wu.deg2rad(60), wu.deg2rad(78), 0.06076096484, + wu.gms2rad([6,56,46.4585]), wu.gms2rad([173,11,5.9592]), 6637321.6350), + + (wu.deg2rad(75), wu.deg2rad(-90), wu.deg2rad(75), wu.deg2rad(66), 0.05805851008, + wu.gms2rad([12,40,34.9009]), wu.gms2rad([168,20,26.7339]), 3267941.2812), + + (wu.deg2rad(80), wu.deg2rad(-90), wu.deg2rad(80), wu.deg2rad(55), 0.05817384452, + wu.gms2rad([18,35,40.7848]), wu.gms2rad([164,25,34.0017]), 2132316.9048) +] + +table3 = [ + (wu.deg2rad(0), wu.deg2rad(0.5), wu.deg2rad(80), wu.deg2rad(0.5), 0.05680316848, + wu.gms2rad([0,-0,16.0757]), wu.gms2rad([0,1,32.5762]), 8831874.3717), + + (wu.deg2rad(-1), wu.deg2rad(5), wu.deg2rad(75), wu.deg2rad(5), 0.05659149555, + wu.gms2rad([0,-1,47.2105]), wu.gms2rad([0,6,54.0958]), 8405370.4947), + + (wu.deg2rad(-5), wu.deg2rad(30), wu.deg2rad(60), wu.deg2rad(30), 0.04921108945, + wu.gms2rad([0,-4,22.3516]), wu.gms2rad([0,8,42.0756]), 7204083.8568), + + (wu.deg2rad(-30), wu.deg2rad(45), wu.deg2rad(30), wu.deg2rad(45), 0.04017812574, + wu.gms2rad([0,-3,41.2461]), wu.gms2rad([0,3,41.2461]), 6652788.1287), + + (wu.deg2rad(-60), wu.deg2rad(60), wu.deg2rad(5), wu.deg2rad(60), 0.02843082609, + wu.gms2rad([0,-8,40.4575]), wu.gms2rad([0,4,22.1675]), 7213412.4477), + + (wu.deg2rad(-75), wu.deg2rad(85), wu.deg2rad(1), wu.deg2rad(85), 0.00497802414, + wu.gms2rad([0,-6,44.6115]), wu.gms2rad([0,1,47.0474]), 8442938.5899), + + (wu.deg2rad(-80), wu.deg2rad(89.5), wu.deg2rad(0), wu.deg2rad(89.5), 0.00050178253, + wu.gms2rad([0,-1,27.9705]), wu.gms2rad([0,0,16.0490]), 8888783.7815) +] + +table4 = [ + (wu.deg2rad(0), wu.deg2rad(0), wu.deg2rad(0), wu.deg2rad(90), 1.00000000000, + wu.gms2rad([90,0,0.0000]), wu.gms2rad([90,0,0.0000]), 10018754.1714), + + (wu.deg2rad(1), wu.deg2rad(0), wu.deg2rad(0), wu.deg2rad(179.5), 0.30320665822, + wu.gms2rad([17,39,11.0942]), wu.gms2rad([162,20,58.9032]), 19884417.8083), + + (wu.deg2rad(5), wu.deg2rad(0), wu.deg2rad(-80), wu.deg2rad(170), 0.03104258442, + wu.gms2rad([178,12,51.5083]), wu.gms2rad([10,17,52.6423]), 11652530.7514), + + (wu.deg2rad(30), wu.deg2rad(0), wu.deg2rad(-75), wu.deg2rad(120), 0.24135347134, + wu.gms2rad([163,49,4.4615]), wu.gms2rad([68,49,50.9617]), 14057886.8752), + + (wu.deg2rad(60), wu.deg2rad(0), wu.deg2rad(-60), wu.deg2rad(40), 0.19408499032, + wu.gms2rad([157,9,33.5589]), wu.gms2rad([157,9,33.5589]), 13767414.8267), + + (wu.deg2rad(75), wu.deg2rad(0), wu.deg2rad(-30), wu.deg2rad(0.5), 0.00202789418, + wu.gms2rad([179,33,3.8613]), wu.gms2rad([179,51,57.0077]), 11661713.4496), + + (wu.deg2rad(80), wu.deg2rad(0), wu.deg2rad(-5), wu.deg2rad(120), 0.15201222384, + wu.gms2rad([61,5,33.9600]), wu.gms2rad([171,13,22.0148]), 11105138.2902), + + (wu.deg2rad(0), wu.deg2rad(0), wu.deg2rad(60), wu.deg2rad(0), 0.00000000000, + wu.gms2rad([0,0,0.0000]), wu.gms2rad([0,0,0.0000]), 6663348.2060) +] + +tables = [table1, table2, table3, table4] + +def get_example(table, example): + table -= 1 + example -= 1 + return tables[table][example] + +def get_tables(): + return tables + + +if __name__ == "__main__": + test = get_example(1, 4) + pass diff --git a/GHA_triaxial/panou.py b/GHA_triaxial/panou.py index 7d3b9c7..9f7c6ca 100644 --- a/GHA_triaxial/panou.py +++ b/GHA_triaxial/panou.py @@ -1,4 +1,5 @@ import numpy as np +from numpy import sin, cos, sqrt, arctan2 import ellipsoide import Numerische_Integration.num_int_runge_kutta as rk import winkelumrechnungen as wu @@ -6,33 +7,66 @@ import ausgaben as aus import GHA.rk as ghark from scipy.special import factorial as fact from math import comb +import GHA_triaxial.numeric_examples_panou as nep # Panou, Korakitits 2019 -def gha1_num(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, num): - phi, lamb, h = ell.cart2geod("ligas3", point) +def gha1_num_old(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, num): + phi, lamb, h = ell.cart2geod(point, "ligas3") x, y, z = ell.geod2cart(phi, lamb, 0) - values = ell.p_q(x, y, z) - H = values["H"] - p = values["p"] - q = values["q"] + p, q = ell.p_q(x, y, z) - 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) + dxds0 = p[0] * sin(alpha0) + q[0] * cos(alpha0) + dyds0 = p[1] * sin(alpha0) + q[1] * cos(alpha0) + dzds0 = p[2] * sin(alpha0) + q[2] * 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) + f1 = lambda x, dxds, y, dyds, z, dzds: dxds + f2 = lambda x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / ell.func_H(x, y, z) * x + f3 = lambda x, dxds, y, dyds, z, dzds: dyds + f4 = lambda x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / ell.func_H(x, y, z) * y/(1-ell.ee**2) + f5 = lambda x, dxds, y, dyds, z, dzds: dzds + f6 = lambda x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / ell.func_H(x, y, z) * z/(1-ell.ex**2) + + funktionswerte = rk.verfahren([f1, f2, f3, f4, f5, f6], [x, dxds0, y, dyds0, z, dzds0], s, num, fein=False) + P2 = funktionswerte[-1] + P2 = (P2[0], P2[2], P2[4]) + return P2 + +def buildODE(ell): + def ODE(v): + x, dxds, y, dyds, z, dzds = v + + H = ell.func_H(x, y, z) + h = dxds**2 + 1/(1-ell.ee**2)*dyds**2 + 1/(1-ell.ex**2)*dzds**2 + + ddx = -(h/H)*x + ddy = -(h/H)*y/(1-ell.ee**2) + ddz = -(h/H)*z/(1-ell.ex**2) + + return [dxds, ddx, dyds, ddy, dzds, ddz] + return ODE + +def gha1_num(ell, point, alpha0, s, num): + phi, lam, _ = ell.cart2geod(point, "ligas3") + x0, y0, z0 = ell.geod2cart(phi, lam, 0) + + p, q = ell.p_q(x0, y0, z0) + dxds0 = p[0] * sin(alpha0) + q[0] * cos(alpha0) + dyds0 = p[1] * sin(alpha0) + q[1] * cos(alpha0) + dzds0 = p[2] * sin(alpha0) + q[2] * cos(alpha0) + + v_init = [x0, dxds0, y0, dyds0, z0, dzds0] + + F = buildODE(ell) + + werte = rk.rk_chat(F, v_init, s, num) + x1, _, y1, _, z1, _ = werte[-1] + + return x1, y1, z1 - funktionswerte = rk.verfahren([f1, f2, f3, f4, f5, f6], [0, x, dxds0, y, dyds0, z, dzds0], s, num) - return funktionswerte def checkLiouville(ell: ellipsoide.EllipsoidTriaxial, points): constantValues = [] @@ -52,9 +86,9 @@ def checkLiouville(ell: ellipsoide.EllipsoidTriaxial, points): P = p[0]*dxds + p[1]*dyds + p[2]*dzds Q = q[0]*dxds + q[1]*dyds + q[2]*dzds - alpha = np.arctan(P/Q) + alpha = arctan2(P, Q) - c = ell.ay**2 - (t1 * np.sin(alpha)**2 + t2 * np.cos(alpha)**2) + c = ell.ay**2 - (t1 * sin(alpha)**2 + t2 * cos(alpha)**2) constantValues.append(c) pass @@ -63,9 +97,7 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM): """ Panou, Korakitits 2020, 5ff. :param ell: - :param x: - :param y: - :param z: + :param point: :param alpha0: :param s: :param maxM: @@ -77,21 +109,22 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM): z_m = [z] # erste Ableitungen (7-8) - sqrtH = np.sqrt(ell.p_q(x, y, z)["H"]) + H = x ** 2 + y ** 2 / (1 - ell.ee ** 2) ** 2 + z ** 2 / (1 - ell.ex ** 2) ** 2 + sqrtH = sqrt(H) n = np.array([x / sqrtH, y / ((1-ell.ee**2) * sqrtH), z / ((1-ell.ex**2) * sqrtH)]) u, v = ell.cart2para(np.array([x, y, z])) - G = np.sqrt(1 - ell.ex**2 * np.cos(u)**2 - ell.ee**2 * np.sin(u)**2 * np.sin(v)**2) - q = np.array([-1/G * np.sin(u) * np.cos(v), - -1/G * np.sqrt(1-ell.ee**2) * np.sin(u) * np.sin(v), - 1/G * np.sqrt(1-ell.ex**2) * np.cos(u)]) + G = sqrt(1 - ell.ex**2 * cos(u)**2 - ell.ee**2 * sin(u)**2 * sin(v)**2) + q = np.array([-1/G * sin(u) * cos(v), + -1/G * sqrt(1-ell.ee**2) * sin(u) * sin(v), + 1/G * sqrt(1-ell.ex**2) * cos(u)]) p = np.array([q[1]*n[2] - q[2]*n[1], q[2]*n[0] - q[0]*n[2], q[0]*n[1] - q[1]*n[0]]) - x_m.append(p[0] * np.sin(alpha0) + q[0] * np.cos(alpha0)) - y_m.append(p[1] * np.sin(alpha0) + q[1] * np.cos(alpha0)) - z_m.append(p[2] * np.sin(alpha0) + q[2] * np.cos(alpha0)) + x_m.append(p[0] * sin(alpha0) + q[0] * cos(alpha0)) + y_m.append(p[1] * sin(alpha0) + q[1] * cos(alpha0)) + z_m.append(p[2] * sin(alpha0) + q[2] * cos(alpha0)) # H Ableitungen (7) H_ = lambda p: np.sum([comb(p, i) * (x_m[p - i] * x_m[i] + @@ -143,32 +176,16 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM): if __name__ == "__main__": # ell = ellipsoide.EllipsoidTriaxial.init_name("Eitschberger1978") - ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980") - ellbi = ellipsoide.EllipsoidTriaxial.init_name("Bessel-biaxial") + ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980round") + # ellbi = ellipsoide.EllipsoidTriaxial.init_name("Bessel-biaxial") re = ellipsoide.EllipsoidBiaxial.init_name("Bessel") # Panou 2013, 7, Table 1, beta0=60° - beta1 = wu.deg2rad(60) - lamb1 = wu.deg2rad(0) - beta2 = wu.deg2rad(60) - lamb2 = wu.deg2rad(175) - P1 = ell.ell2cart(wu.deg2rad(60), wu.deg2rad(0)) - P2 = ell.ell2cart(wu.deg2rad(60), wu.deg2rad(175)) - para1 = ell.cart2para(P1) - para2 = ell.cart2para(P2) - cart1 = ell.para2cart(para1[0], para1[1]) - cart2 = ell.para2cart(para2[0], para2[1]) - ell11 = ell.cart2ell(P1) - ell21 = ell.cart2ell(P2) - ell1 = ell.cart2ell(cart1) - ell2 = ell.cart2ell(cart2) + beta0, lamb0, beta1, lamb1, c, alpha0, alpha1, s = nep.get_example(table=1, example=5) + P0 = ell.ell2cart(beta0, lamb0) + P1 = ell.ell2cart(beta1, lamb1) - c = 0.06207487624 - alpha0 = wu.gms2rad([2, 52, 26.2393]) - alpha1 = wu.gms2rad([177, 4, 13.6373]) - s = 6705715.1610 - pass - - P2_num = gha1_num(ell, P1, alpha0, s, 1000) - P2_ana = gha1_ana(ell, P1, alpha0, s, 70) + # P1_num = gha1_num(ell, P0, alpha0, s, 1000) + P1_num = gha1_num(ell, P0, alpha0, s, 10000) + P1_ana = gha1_ana(ell, P0, alpha0, s, 30) pass \ No newline at end of file diff --git a/GHA_triaxial/panou_2013_2GHA_num.py b/GHA_triaxial/panou_2013_2GHA_num.py index 680e335..ce2012c 100644 --- a/GHA_triaxial/panou_2013_2GHA_num.py +++ b/GHA_triaxial/panou_2013_2GHA_num.py @@ -319,13 +319,20 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep if __name__ == "__main__": ell = EllipsoidTriaxial.init_name("BursaSima1980round") - beta1 = np.deg2rad(75) - lamb1 = np.deg2rad(-90) - beta2 = np.deg2rad(75) - lamb2 = np.deg2rad(66) - a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2) - print(aus.gms("a1", a1, 4)) - print(aus.gms("a2", a2, 4)) + # beta1 = np.deg2rad(75) + # lamb1 = np.deg2rad(-90) + # beta2 = np.deg2rad(75) + # lamb2 = np.deg2rad(66) + # a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2) + # print(aus.gms("a1", a1, 4)) + # print(aus.gms("a2", a2, 4)) + # print(s) + cart1 = ell.para2cart(0, 0) + cart2 = ell.para2cart(0.4, 0.4) + beta1, lamb1 = ell.cart2ell(cart1) + beta2, lamb2 = ell.cart2ell(cart2) + + a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=2500) print(s) diff --git a/Numerische_Integration/num_int_runge_kutta.py b/Numerische_Integration/num_int_runge_kutta.py index d363a58..ca12ab8 100644 --- a/Numerische_Integration/num_int_runge_kutta.py +++ b/Numerische_Integration/num_int_runge_kutta.py @@ -1,4 +1,4 @@ -def verfahren(funktionen: list, startwerte: list, weite: float, schritte: int) -> list: +def verfahren(funktionen: list, startwerte: list, weite: float, schritte: int, fein: bool = True) -> list: """ Runge-Kutta-Verfahren für ein beliebiges DGLS :param funktionen: Liste mit allen Funktionen @@ -14,19 +14,21 @@ def verfahren(funktionen: list, startwerte: list, weite: float, schritte: int) - 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))] + if fein: + 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_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))] - 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_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) + werte.append(werte_korr) + else: + werte.append(werte_grob) return werte @@ -60,3 +62,23 @@ def zuschlaege(funktionen: list, startwerte: list, h: float) -> list: k_ = [(k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) / 6 for i in range(len(k1))] return k_ + +def rk_chat(F, v0: list, weite: float, schritte: int): + h = weite/schritte + v = v0 + werte = [v] + + for _ in range(schritte): + k1 = F(v) + k2 = F([v[i] + 0.5 * h * k1[i] for i in range(6)]) + k3 = F([v[i] + 0.5 * h * k2[i] for i in range(6)]) + k4 = F([v[i] + h * k3[i] for i in range(6)]) + + v = [ + v[i] + (h / 6) * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) + for i in range(6) + ] + + werte.append(v) + + return werte \ No newline at end of file diff --git a/dashborad.py b/dashboard.py similarity index 100% rename from dashborad.py rename to dashboard.py diff --git a/ellipsoide.py b/ellipsoide.py index 606e912..399f12c 100644 --- a/ellipsoide.py +++ b/ellipsoide.py @@ -1,4 +1,5 @@ import numpy as np +from numpy import sin, cos, arctan, arctan2, sqrt import winkelumrechnungen as wu import ausgaben as aus import jacobian_Ligas @@ -150,10 +151,15 @@ class EllipsoidTriaxial: b = 6356078.96290 return cls(ax, ay, b) elif name == "Fiction": - ax = 5500000 - ay = 4500000 + ax = 6000000 + ay = 5000000 b = 4000000 return cls(ax, ay, b) + elif name == "KarneyTest2024": + ax = np.sqrt(2) + ay = 1 + b = 1 / np.sqrt(2) + return cls(ax, ay, b) def point_on(self, point: np.ndarray) -> bool: """ @@ -170,7 +176,7 @@ class EllipsoidTriaxial: def ellu2cart(self, beta: float, lamb: float, u: float) -> np.ndarray: """ Panou 2014 12ff. - Ellipsoidische Breite+Länge sind nicht gleich der geodätischen + Elliptische Breite+Länge sind nicht gleich der geodätischen Verhältnisse des Ellipsoids bekannt, Größe verändern bis Punkt erreicht, dann ist u die Größe entlang der z-Achse :param beta: ellipsoidische Breite [rad] @@ -196,38 +202,47 @@ class EllipsoidTriaxial: return np.array([x, y, z]) - def ell2cart(self, beta: float, lamb: float) -> np.ndarray: + def ell2cart(self, beta: float | np.ndarray, lamb: float | np.ndarray) -> np.ndarray: """ Panou, Korakitis 2019 2 - :param beta: ellipsoidische Breite [rad] - :param lamb: ellipsoidische Länge [rad] + :param beta: elliptische Breite [rad] + :param lamb: elliptische Länge [rad] :return: Punkt in kartesischen Koordinaten """ - if beta == -np.pi/2: - return np.array([0, 0, -self.b]) - elif beta == np.pi/2: - return np.array([0, 0, self.b]) - elif beta == 0 and lamb == -np.pi/2: - return np.array([0, -self.ay, 0]) - elif beta == 0 and lamb == np.pi/2: - return np.array([0, self.ay, 0]) - elif beta == 0 and lamb == 0: - return np.array([self.ax, 0, 0]) - elif beta == 0 and lamb == np.pi: - return np.array([-self.ax, 0, 0]) - else: - B = self.Ex**2 * np.cos(beta)**2 + self.Ee**2 * np.sin(beta)**2 - L = self.Ex**2 - self.Ee**2 * np.cos(lamb)**2 - x = self.ax / self.Ex * np.sqrt(B) * np.cos(lamb) - y = self.ay * np.cos(beta) * np.sin(lamb) - z = self.b / self.Ex * np.sin(beta) * np.sqrt(L) - return np.array([x, y, z]) + beta = np.asarray(beta, dtype=float) + lamb = np.asarray(lamb, dtype=float) + + beta, lamb = np.broadcast_arrays(beta, lamb) + + B = self.Ex ** 2 * np.cos(beta) ** 2 + self.Ee ** 2 * np.sin(beta) ** 2 + L = self.Ex ** 2 - self.Ee ** 2 * np.cos(lamb) ** 2 + + x = self.ax / self.Ex * np.sqrt(B) * np.cos(lamb) + y = self.ay * np.cos(beta) * np.sin(lamb) + z = self.b / self.Ex * np.sin(beta) * np.sqrt(L) + + xyz = np.stack((x, y, z), axis=-1) + + # Pole + mask_south = beta == -np.pi / 2 + mask_north = beta == np.pi / 2 + xyz[mask_south] = np.array([0, 0, -self.b]) + xyz[mask_north] = np.array([0, 0, self.b]) + + # Äquator + mask_eq = beta == 0 + xyz[mask_eq & (lamb == -np.pi / 2)] = np.array([0, -self.ay, 0]) + xyz[mask_eq & (lamb == np.pi / 2)] = np.array([0, self.ay, 0]) + xyz[mask_eq & (lamb == 0)] = np.array([self.ax, 0, 0]) + xyz[mask_eq & (lamb == np.pi)] = np.array([-self.ax, 0, 0]) + + return xyz def cart2ellu(self, point: np.ndarray) -> tuple[float, float, float]: """ Panou 2014 15ff. :param point: Punkt in kartesischen Koordinaten - :return: ellipsoidische Breite, ellipsoidische Länge, Größe entlang der z-Achse + :return: elliptische Breite, elliptische Länge, Größe entlang der z-Achse """ x, y, z = point c2 = self.ax**2 + self.ay**2 + self.b**2 - x**2 - y**2 - z**2 @@ -246,7 +261,10 @@ class EllipsoidTriaxial: # print(s1, s2, s3) beta = np.arctan(np.sqrt((-self.b**2 - s2) / (self.ay**2 + s2))) - lamb = np.arctan(np.sqrt((-self.ay**2 - s3) / (self.ax**2 + s3))) + if abs((-self.ay**2 - s3) / (self.ax**2 + s3)) > 1e-7: + lamb = np.arctan(np.sqrt((-self.ay**2 - s3) / (self.ax**2 + s3))) + else: + lamb = 0 u = np.sqrt(self.b**2 + s1) return beta, lamb, u @@ -255,7 +273,7 @@ class EllipsoidTriaxial: """ Panou, Korakitis 2019 2f. :param point: Punkt in kartesischen Koordinaten - :return: ellipsoidische Breite, ellipsoidische Länge + :return: elliptische Breite, elliptische Länge """ x, y, z = point @@ -312,7 +330,7 @@ class EllipsoidTriaxial: return beta, lamb - def cart2geod(self, mode: str, point: np.ndarray, maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]: + def cart2geod(self, point: np.ndarray, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]: """ Ligas 2012 :param mode: ligas1, ligas2, oder ligas3 @@ -328,39 +346,21 @@ class EllipsoidTriaxial: if abs(xG) < eps and abs(yG) < eps: # Punkt in der z-Achse phi = np.pi / 2 if zG > 0 else -np.pi / 2 lamb = 0.0 - h = abs(zG) - ell.b + h = abs(zG) - self.b return phi, lamb, h elif abs(xG) < eps and abs(zG) < eps: # Punkt in der y-Achse phi = 0.0 lamb = np.pi / 2 if yG > 0 else -np.pi / 2 - h = abs(yG) - ell.ay + h = abs(yG) - self.ay return phi, lamb, h elif abs(yG) < eps and abs(zG) < eps: # Punkt in der x-Achse phi = 0.0 lamb = 0.0 if xG > 0 else np.pi - h = abs(xG) - ell.ax + h = abs(xG) - self.ax return phi, lamb, h - # elif abs(zG) < eps: # Punkt in der xy-Ebene - # phi = 0 - # lamb = np.arctan2(yG / ell.ay**2, xG / ell.ax**2) - # rG = np.sqrt(xG ** 2 + yG ** 2) - # pE = np.array([self.ax * xG / rG, self.ax * yG / rG, self.ax * zG / rG], dtype=np.float64) - # rE = np.sqrt(pE[0] ** 2 + pE[1] ** 2) - # h = rG - rE - # return phi, lamb, h - # - # elif abs(yG) < eps: # Punkt in der xz-Ebene - # phi = np.arctan2(zG / ell.b**2, xG / ell.ax**2) - # lamb = 0 if xG > 0 else np.pi - # rG = np.sqrt(xG ** 2 + zG ** 2) - # pE = np.array([self.ax * xG / rG, self.ax * yG / rG, self.ax * zG / rG], dtype=np.float64) - # rE = np.sqrt(pE[0] ** 2 + pE[2] ** 2) - # h = rG - rE - # return phi, lamb, h - rG = np.sqrt(xG ** 2 + yG ** 2 + zG ** 2) pE = np.array([self.ax * xG / rG, self.ax * yG / rG, self.ax * zG / rG], dtype=np.float64) @@ -399,7 +399,7 @@ class EllipsoidTriaxial: return phi, lamb, h - def geod2cart(self, phi: float, lamb: float, h: float) -> np.ndarray: + def geod2cart(self, phi: float | np.ndarray, lamb: float | np.ndarray, h: float) -> np.ndarray: """ Ligas 2012, 250 :param phi: geodätische Breite [rad] @@ -419,7 +419,7 @@ class EllipsoidTriaxial: :param point: Punkt in kartesischen Koordinaten, der gelotet werden soll :return: Lotpunkt in kartesischen Koordinaten, geodätische Koordinaten des Punktes """ - phi, lamb, h = self.cart2geod("ligas3", point) + phi, lamb, h = self.cart2geod(point, "ligas3") x, y, z = self. geod2cart(phi, lamb, 0) return np.array([x, y, z]), phi, lamb, h @@ -430,11 +430,11 @@ class EllipsoidTriaxial: :param h: Höhe über dem Ellipsoid :return: hochgeloteter Punkt """ - phi, lamb, _ = self.cart2geod("ligas3", point) + phi, lamb, _ = self.cart2geod(point, "ligas3") pointH = self. geod2cart(phi, lamb, h) return pointH - def para2cart(self, u: float, v: float) -> np.ndarray: + def para2cart(self, u: float | np.ndarray, v: float | np.ndarray) -> np.ndarray: """ Panou, Korakitits 2020, 4 :param u: Parameter u @@ -444,6 +444,7 @@ class EllipsoidTriaxial: x = self.ax * np.cos(u) * np.cos(v) y = self.ay * np.cos(u) * np.sin(v) z = self.b * np.sin(u) + z = np.broadcast_to(z, np.shape(x)) return np.array([x, y, z]) def cart2para(self, point: np.ndarray) -> tuple[float, float]: @@ -471,7 +472,37 @@ class EllipsoidTriaxial: return u, v - def p_q(self, x, y, z) -> dict: + def ell2para(self, beta, lamb) -> tuple[float, float]: + cart = self.ell2cart(beta, lamb) + return self.cart2para(cart) + + def para2ell(self, u, v) -> tuple[float, float]: + cart = self.para2cart(u, v) + return self.cart2ell(cart) + + def para2geod(self, u: float, v: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]: + cart = self.para2cart(u, v) + return self.cart2geod(cart, mode, maxIter, maxLoa) + + def geod2para(self, phi, lamb, h) -> tuple[float, float]: + cart = self.geod2cart(phi, lamb, h) + return self.cart2para(cart) + + def ell2geod(self, beta, lamb, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]: + cart = self.ell2cart(beta, lamb) + return self.cart2geod(cart, mode, maxIter, maxLoa) + + def func_H(self, x, y, z): + return x ** 2 + y ** 2 / (1 - self.ee ** 2) ** 2 + z ** 2 / (1 - self.ex ** 2) ** 2 + + def func_n(self, x, y, z, H=None): + if H is None: + H = self.func_H(x, y, z) + return np.array([x / sqrt(H), + y / ((1 - self.ee ** 2) * sqrt(H)), + z / ((1 - self.ex ** 2) * sqrt(H))]) + + def p_q(self, x, y, z) -> tuple[np.ndarray, np.ndarray]: """ Berechnung sämtlicher Größen :param x: x @@ -479,11 +510,9 @@ class EllipsoidTriaxial: :param z: z :return: Dictionary sämtlicher Größen """ - H = x ** 2 + y ** 2 / (1 - self.ee ** 2) ** 2 + z ** 2 / (1 - self.ex ** 2) ** 2 + n = self.func_n(x, y, z) - n = np.array([x / np.sqrt(H), y / ((1 - self.ee ** 2) * np.sqrt(H)), z / ((1 - self.ex ** 2) * np.sqrt(H))]) - - beta, lamb, u = self.cart2ellu(np.array([x, y, z])) + beta, lamb = self.cart2ell(np.array([x, y, z])) B = self.Ex ** 2 * np.cos(beta) ** 2 + self.Ee ** 2 * np.sin(beta) ** 2 L = self.Ex ** 2 - self.Ee ** 2 * np.cos(lamb) ** 2 @@ -507,11 +536,9 @@ class EllipsoidTriaxial: p = np.array([p1, p2, p3]) q = np.array([n[1] * p[2] - n[2] * p[1], n[2] * p[0] - n[0] * p[2], - n[1] * p[1] - n[1] * p[0]]) + n[0] * p[1] - n[1] * p[0]]) - return {"H": H, "n": n, "beta": beta, "lamb": lamb, "u": u, "B": B, "L": L, "c1": c1, "c0": c0, "t1": t1, - "t2": t2, - "F": F, "p": p, "q": q} + return p, q if __name__ == "__main__": @@ -531,19 +558,20 @@ if __name__ == "__main__": cart_para = ell.para2cart(para[0], para[1]) diff_para = np.sum(np.abs(point-cart_para)) - geod = ell.cart2geod("ligas1", point) - cart_geod = ell.geod2cart(geod[0], geod[1], geod[2]) - diff_geod1 = np.sum(np.abs(point-cart_geod)) + # geod = ell.cart2geod(point, "ligas1") + # cart_geod = ell.geod2cart(geod[0], geod[1], geod[2]) + # diff_geod1 = np.sum(np.abs(point-cart_geod)) + # + # geod = ell.cart2geod(point, "ligas2") + # cart_geod = ell.geod2cart(geod[0], geod[1], geod[2]) + # diff_geod2 = np.sum(np.abs(point-cart_geod)) - geod = ell.cart2geod("ligas2", point) - cart_geod = ell.geod2cart(geod[0], geod[1], geod[2]) - diff_geod2 = np.sum(np.abs(point-cart_geod)) - - geod = ell.cart2geod("ligas3", point) + geod = ell.cart2geod(point, "ligas3") cart_geod = ell.geod2cart(geod[0], geod[1], geod[2]) diff_geod3 = np.sum(np.abs(point-cart_geod)) - diff_list.append([beta_deg, lamb_deg, diff_ell, diff_para, diff_geod1, diff_geod2, diff_geod3]) + diff_list.append([beta_deg, lamb_deg, diff_ell, diff_para, diff_geod3]) + diff_list.append([diff_ell]) diff_list = np.array(diff_list) pass \ No newline at end of file diff --git a/show_constant_lines.py b/show_constant_lines.py new file mode 100644 index 0000000..d20c6e6 --- /dev/null +++ b/show_constant_lines.py @@ -0,0 +1,30 @@ +import numpy as np +import plotly.graph_objects as go +from ellipsoide import EllipsoidTriaxial +import winkelumrechnungen as wu +from dashboard import ellipsoid_figure + +u = np.linspace(0, 2*np.pi, 51) +v = np.linspace(0, np.pi, 51) +ell = EllipsoidTriaxial.init_name("BursaSima1980round") +points = [] +lines = [] +for u_i, u_value in enumerate(u): + for v_i, v_value in enumerate(v): + cart = ell.ell2cart(u_value, v_value) + if u_i != 0 and v_i != 0: + lines.append((points[-1], cart, "red")) + points.append(cart) +points = [] +for v_i, v_value in enumerate(v): + for u_i, u_value in enumerate(u): + cart = ell.ell2cart(u_value, v_value) + if u_i != 0 and v_i != 0: + lines.append((points[-1], cart, "blue")) + points.append(cart) +ax = ell.ax +ay = ell.ay +b = ell.b + +figu = ellipsoid_figure(ax, ay, b, lines=lines) +figu.show() \ No newline at end of file diff --git a/test_algorithms.py b/test_algorithms.py new file mode 100644 index 0000000..f48298b --- /dev/null +++ b/test_algorithms.py @@ -0,0 +1,84 @@ +import GHA_triaxial.numeric_examples_panou as nep +import ellipsoide +from GHA_triaxial.panou_2013_2GHA_num import gha2_num +from GHA_triaxial.panou import gha1_ana, gha1_num +import numpy as np +import time + +def test(): + ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980round") + + tables = nep.get_tables() + + diffs_gha1_num = [] + diffs_gha1_ana = [] + diffs_gha2_num = [] + times_gha1_num = [] + times_gha1_ana = [] + times_gha2_num = [] + + for table in tables: + diffs_gha1_num.append([]) + diffs_gha1_ana.append([]) + diffs_gha2_num.append([]) + times_gha1_num.append([]) + times_gha1_ana.append([]) + times_gha2_num.append([]) + + for example in table: + beta0, lamb0, beta1, lamb1, c, alpha0, alpha1, s = example + P0 = ell.ell2cart(beta0, lamb0) + P1 = ell.ell2cart(beta1, lamb1) + + start = time.perf_counter() + try: + P1_num = gha1_num(ell, P0, alpha0, s, 10000) + end = time.perf_counter() + diff_P1_num = np.linalg.norm(P1 - P1_num) + except: + end = time.perf_counter() + diff_P1_num = None + time_gha1_num = end - start + + start = time.perf_counter() + try: + P1_ana = gha1_ana(ell, P0, alpha0, s, 50) + end = time.perf_counter() + diff_P1_ana = np.linalg.norm(P1 - P1_ana) + except: + end = time.perf_counter() + diff_P1_ana = None + time_gha1_ana = end - start + + start = time.perf_counter() + try: + alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=1000) + end = time.perf_counter() + diff_s_num = abs(s - s_num) + except: + end = time.perf_counter() + diff_s_num = None + time_gha2_num = None + time_gha2_num = end - start + + diffs_gha1_num[-1].append(diff_P1_num) + diffs_gha1_ana[-1].append(diff_P1_ana) + diffs_gha2_num[-1].append(diff_s_num) + times_gha1_num[-1].append(time_gha1_num) + times_gha1_ana[-1].append(time_gha1_ana) + times_gha2_num[-1].append(time_gha2_num) + print(diffs_gha1_num, diffs_gha1_ana, diffs_gha2_num) + print(times_gha1_num, times_gha1_ana, times_gha2_num) + +def display(): + diffs = [[{'gha1_num': np.float64(3.410763124264611e-05), 'gha1_ana': np.float64(3.393273802112796e-05), 'gha2_num': np.float64(3.3931806683540344e-05)}, {'gha1_num': np.float64(0.0008736425000530604), 'gha1_ana': np.float64(0.0008736458415010259), 'gha2_num': None}, {'gha1_num': np.float64(0.0007739730058338136), 'gha1_ana': np.float64(0.0007739621469802854), 'gha2_num': np.float64(1.5832483768463135e-07)}, {'gha1_num': np.float64(0.00010554956741100295), 'gha1_ana': np.float64(8.814246009944831), 'gha2_num': np.float64(4.864111542701721e-05)}, {'gha1_num': np.float64(0.0002135908394614854), 'gha1_ana': np.float64(0.0002138610897967267), 'gha2_num': np.float64(5.0179407158866525)}, {'gha1_num': np.float64(0.00032727226891456654), 'gha1_ana': np.float64(0.00032734569198545905), 'gha2_num': np.float64(9.735533967614174e-05)}, {'gha1_num': np.float64(0.0005195973303787956), 'gha1_ana': np.float64(0.0005197766935509641), 'gha2_num': None}], [{'gha1_num': np.float64(1.780250537652368e-05), 'gha1_ana': np.float64(1.996805145339501e-05), 'gha2_num': np.float64(1.8164515495300293e-05)}, {'gha1_num': np.float64(4.8607540473363564e-05), 'gha1_ana': np.float64(2205539.954949392), 'gha2_num': None}, {'gha1_num': np.float64(0.00017376854985685854), 'gha1_ana': np.float64(328124.1513636429), 'gha2_num': np.float64(0.17443156614899635)}, {'gha1_num': np.float64(5.83429352558999e-05), 'gha1_ana': np.float64(0.01891628037258558), 'gha2_num': np.float64(1.4207654744386673)}, {'gha1_num': np.float64(0.0006421087024666934), 'gha1_ana': np.float64(0.0006420400127297228), 'gha2_num': np.float64(0.12751091085374355)}, {'gha1_num': np.float64(0.0004456207867164434), 'gha1_ana': np.float64(0.0004455649707698245), 'gha2_num': np.float64(0.00922046648338437)}, {'gha1_num': np.float64(0.0002340879908275419), 'gha1_ana': np.float64(0.00023422217242111216), 'gha2_num': np.float64(0.001307751052081585)}], [{'gha1_num': np.float64(976.6580096633622), 'gha1_ana': np.float64(976.6580096562798), 'gha2_num': np.float64(6.96033239364624e-05)}, {'gha1_num': np.float64(2825.2936643258527), 'gha1_ana': np.float64(2794.954866417055), 'gha2_num': np.float64(1.3615936040878296e-05)}, {'gha1_num': np.float64(1248.8942058074501), 'gha1_ana': np.float64(538.5550561841195), 'gha2_num': np.float64(3.722589462995529e-05)}, {'gha1_num': np.float64(2201.1793359793814), 'gha1_ana': np.float64(3735.376499414938), 'gha2_num': np.float64(1.4525838196277618e-05)}, {'gha1_num': np.float64(2262.134819997246), 'gha1_ana': np.float64(25549.567793410763), 'gha2_num': np.float64(9.328126907348633e-06)}, {'gha1_num': np.float64(2673.219788119847), 'gha1_ana': np.float64(21760.866677295206), 'gha2_num': np.float64(8.635222911834717e-06)}, {'gha1_num': np.float64(1708.758419275875), 'gha1_ana': np.float64(3792.1128807063437), 'gha2_num': np.float64(2.4085864424705505e-05)}], [{'gha1_num': np.float64(0.7854659044152204), 'gha1_ana': np.float64(0.785466068424286), 'gha2_num': np.float64(0.785466069355607)}, {'gha1_num': np.float64(237.79878717216718), 'gha1_ana': np.float64(1905080.064324282), 'gha2_num': None}, {'gha1_num': np.float64(55204.601699830164), 'gha1_ana': np.float64(55204.60175211949), 'gha2_num': None}, {'gha1_num': np.float64(12766.348063015519), 'gha1_ana': np.float64(12766.376619517901), 'gha2_num': np.float64(12582.786206113175)}, {'gha1_num': np.float64(29703.049988324146), 'gha1_ana': np.float64(29703.056427749252), 'gha2_num': np.float64(28933.668131249025)}, {'gha1_num': np.float64(43912.03007182513), 'gha1_ana': np.float64(43912.03007528712), 'gha2_num': None}, {'gha1_num': np.float64(28522.29828970693), 'gha1_ana': np.float64(28522.29830145182), 'gha2_num': None}, {'gha1_num': np.float64(17769.115549537233), 'gha1_ana': np.float64(17769.115549483362), 'gha2_num': np.float64(17769.121286311187)}]] + arr = [] + for table in diffs: + for example in table: + arr.append([example['gha1_num'], example['gha1_ana'], example['gha2_num']]) + arr = np.array(arr) + pass + +if __name__ == "__main__": + # test() + display() diff --git a/winkelumrechnungen.py b/winkelumrechnungen.py index 0b0b7f8..8985c21 100644 --- a/winkelumrechnungen.py +++ b/winkelumrechnungen.py @@ -1,4 +1,5 @@ from numpy import * +import numpy as np def deg2gms(deg: float) -> list: @@ -10,13 +11,13 @@ def deg2gms(deg: float) -> list: :rtype: list """ gra = deg // 1 - min = gra % 1 + minu = gra % 1 gra = gra // 1 - min *= 60 - sek = min % 1 - min = min // 1 + minu *= 60 + sek = minu % 1 + minu = minu // 1 sek *= 60 - return [gra, min, sek] + return [gra, minu, sek] def deg2gra(deg: float) -> float: @@ -30,13 +31,13 @@ def deg2gra(deg: float) -> float: return deg * 10/9 -def deg2rad(deg: float) -> float: +def deg2rad(deg: float | np.ndarray) -> float | np.ndarray: """ Umrechnung von Grad in Radiant :param deg: Winkel in Grad - :type deg: float + :type deg: float or np.ndarray :return: Winkel in Radiant - :rtype: float + :rtype: float or np.ndarray """ return deg * pi / 180 @@ -51,13 +52,13 @@ def gra2gms(gra: float) -> list: """ deg = gra2deg(gra) gra = deg // 1 - min = gra % 1 + minu = gra % 1 gra = gra // 1 - min *= 60 - sek = min % 1 - min = min // 1 + minu *= 60 + sek = minu % 1 + minu = minu // 1 sek *= 60 - return [gra, min, sek] + return [gra, minu, sek] def gra2rad(gra: float) -> float: @@ -113,13 +114,13 @@ def rad2gms(rad: float) -> list: :rtype: list """ deg = rad2deg(rad) - min = deg % 1 + minu = deg % 1 gra = deg // 1 - min *= 60 - sek = min % 1 - min = min // 1 + minu *= 60 + sek = minu % 1 + minu = minu // 1 sek *= 60 - return [gra, min, sek] + return [gra, minu, sek] def gms2rad(gms: list) -> float: