diff --git a/GHA_triaxial/panou_2013_2GHA_num.py b/GHA_triaxial/panou_2013_2GHA_num.py index deb9212..4d4d229 100644 --- a/GHA_triaxial/panou_2013_2GHA_num.py +++ b/GHA_triaxial/panou_2013_2GHA_num.py @@ -30,6 +30,22 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float def arccot(x): return np.arctan2(1.0, x) + def cot(a): + return np.cos(a) / np.sin(a) + + def wrap_to_pi(x): + return (x + np.pi) % (2 * np.pi) - np.pi + + def sph_azimuth(beta1, lam1, beta2, lam2): + # sphärischer Anfangsazimut (von Norden/meridian, im Bogenmaß) + dlam = wrap_to_pi(lam2 - lam1) + y = np.sin(dlam) * np.cos(beta2) + x = np.cos(beta1) * np.sin(beta2) - np.sin(beta1) * np.cos(beta2) * np.cos(dlam) + a = np.arctan2(y, x) # (-pi, pi] + if a < 0: + a += 2 * np.pi + return a + def BETA_LAMBDA(beta, lamb): BETA = (ell.ay**2 * np.sin(beta)**2 + ell.b**2 * np.cos(beta)**2) / (ell.Ex**2 - ell.Ey**2 * np.sin(beta)**2) @@ -158,11 +174,13 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float N = n dlamb = lamb_2 - lamb_1 + alpha0_sph = sph_azimuth(beta_1, lamb_1, beta_2, lamb_2) if abs(dlamb) < 1e-15: beta_0 = 0.0 else: - beta_0 = (beta_2 - beta_1) / (lamb_2 - lamb_1) + (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1) + beta_0 = np.sqrt(G1 / E1) * cot(alpha0_sph) converged = False iterations = 0 @@ -170,40 +188,76 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float # funcs = functions() ode_lamb = buildODElamb() - for i in range(iter_max): - iterations = i + 1 + def solve_newton(beta_p0_init: float): + beta_p0 = float(beta_p0_init) - # startwerte = [lamb_1, beta_1, beta_0, 0.0, 1.0] - startwerte = np.array([beta_1, beta_0, 0.0, 1.0]) + for _ in range(iter_max): + startwerte = np.array([beta_1, beta_p0, 0.0, 1.0], dtype=float) + lamb_list, states = rk.rk4(ode_lamb, lamb_1, startwerte, dlamb, N, False) - # werte = rk.verfahren(funcs, startwerte, dlamb, N) - lamb_list, werte = rk.rk4(ode_lamb, lamb_1, startwerte, dlamb, N, False) - # lamb_end, beta_end, beta_p_end, X3_end, X4_end = werte[-1] - lamb_end = lamb_list[-1] - beta_end, beta_p_end, X3_end, X4_end = werte[-1] + beta_end, beta_p_end, X3_end, X4_end = states[-1] + delta = beta_end - beta_2 - d_beta_end_d_beta0 = X3_end - delta = beta_end - beta_2 + if abs(delta) < epsilon: + return True, beta_p0, lamb_list, states - if abs(delta) < epsilon: - converged = True - break + d_beta_end_d_beta0 = X3_end + if abs(d_beta_end_d_beta0) < 1e-20: + return False, None, None, None - if abs(d_beta_end_d_beta0) < 1e-20: - raise RuntimeError("Abbruch.") + step = delta / d_beta_end_d_beta0 + max_step = 0.5 + if abs(step) > max_step: + step = np.sign(step) * max_step - max_step = 0.5 - step = delta / d_beta_end_d_beta0 - if abs(step) > max_step: - step = np.sign(step) * max_step - beta_0 = beta_0 - step + beta_p0 = beta_p0 - step - if not converged: - raise RuntimeError("konvergiert nicht.") + return False, None, None, None - # Z - # werte = rk.verfahren(funcs, [lamb_1, beta_1, beta_0, 0.0, 1.0], dlamb, N, False) - lamb_list, werte = rk.rk4(ode_lamb, lamb_1, np.array([beta_1, beta_0, 0.0, 1.0]), dlamb, N, False) + alpha0_sph = sph_azimuth(beta_1, lamb_1, beta_2, lamb_2) + (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1) + beta_p0_sph = np.sqrt(G1 / E1) * cot(alpha0_sph) + + guesses = [ + beta_p0_sph, + 0.5 * beta_p0_sph, + 2.0 * beta_p0_sph, + -beta_p0_sph, + -0.5 * beta_p0_sph, + ] + + best = None + + for g in guesses: + ok, beta_p0_sol, lamb_list_cand, states_cand = solve_newton(g) + if not ok: + continue + + beta_arr_c = np.array([st[0] for st in states_cand], dtype=float) + beta_p_arr_c = np.array([st[1] for st in states_cand], dtype=float) + lamb_arr_c = np.array(lamb_list_cand, dtype=float) + + integrand = np.zeros(N + 1) + for i in range(N + 1): + (_, _, Ei, Gi, *_) = BETA_LAMBDA(beta_arr_c[i], lamb_arr_c[i]) + integrand[i] = np.sqrt(Ei * beta_p_arr_c[i] ** 2 + Gi) + + h = abs(dlamb) / N + if N % 2 == 0: + S = integrand[0] + integrand[-1] \ + + 4.0 * np.sum(integrand[1:-1:2]) \ + + 2.0 * np.sum(integrand[2:-1:2]) + s_cand = h / 3.0 * S + else: + s_cand = np.trapz(integrand, dx=h) + + if (best is None) or (s_cand < best[0]): + best = (s_cand, beta_p0_sol, lamb_list_cand, states_cand) + + if best is None: + raise RuntimeError("Keine Multi-Start-Variante konvergiert.") + + s_best, beta_0, lamb_list, werte = best beta_arr = np.zeros(N + 1) # lamb_arr = np.zeros(N + 1) diff --git a/dashboard.py b/dashboard.py index 7caa4b5..6eca749 100644 --- a/dashboard.py +++ b/dashboard.py @@ -13,7 +13,7 @@ from GHA_triaxial.approx_gha1 import gha1_approx from GHA_triaxial.panou_2013_2GHA_num import gha2_num from GHA_triaxial.ES_gha2 import gha2_ES -from GHA_triaxial.approx_gha2 import gha2 +from GHA_triaxial.approx_gha2 import gha2_approx app = Dash(__name__, suppress_callback_exceptions=True) @@ -194,7 +194,9 @@ app.layout = html.Div( {"label": "Bursa1972", "value": "Bursa1972"}, {"label": "Bursa1970", "value": "Bursa1970"}, {"label": "BesselBiaxial", "value": "BesselBiaxial"}, + {"label": "KarneyTest2024", "value": "KarneyTest2024"}, {"label": "Fiction", "value": "Fiction"}, + ], value="", style={"width": "300px", "marginBottom": "16px"}, @@ -241,12 +243,25 @@ app.layout = html.Div( "minWidth": "520px", "position": "sticky", "top": "0", - "marginTop": "-150px", + "marginTop": "-50px", }, children=[ + html.Label("Koordinatenart wählen:"), + dcc.Dropdown( + id="dropdown-coors-type", + options=[ + {"label": "Ellipsoidisch", "value": "ell"}, + {"label": "Parametrisch", "value": "para"}, + {"label": "Geodätisch", "value": "geod"}, + + ], + value="ell", + clearable=False, + style={"width": "300px"}, + ), dcc.Graph( id="ellipsoid-plot", - style={"height": "90vh", "width": "100%"}, + style={"height": "85vh", "width": "100%"}, config={"responsive": True} ) ], @@ -297,7 +312,7 @@ def render_content(tab): options=[ {"label": "Analytisch", "value": "analytisch"}, {"label": "Numerisch", "value": "numerisch"}, - {"label": "Stochastisch (ES)", "value": "stochastisch"}, + #{"label": "Stochastisch (ES)", "value": "stochastisch"}, {"label": "Approximiert", "value": "approx"}, ], value=[], @@ -602,7 +617,7 @@ def compute_gha2_num(n2, beta0, lamb0, beta1, lamb1, ax, ay, b, method2): P0 = ell.ell2cart(beta0_rad, lamb0_rad) P1 = ell.ell2cart(beta1_rad, lamb1_rad) - a0_num, a1_num, s_num, beta_arr, lamb_arr = gha2_num(ell, beta0_rad, lamb0_rad, beta1_rad, lamb1_rad, all_points=True) + a0_num, a1_num, s_num, beta_arr, lamb_arr = gha2_num(ell, beta0_rad, lamb0_rad, beta1_rad, lamb1_rad, all_points=True, n=3000) polyline = [] for b_rad, l_rad in zip(beta_arr, lamb_arr): @@ -652,7 +667,7 @@ def compute_gha2_stoch(n2, beta0, lamb0, beta1, lamb1, ax, ay, b, method2): P0 = ell.ell2cart(beta0_rad, lamb0_rad) P1 = ell.ell2cart(beta1_rad, lamb1_rad) - a0_stoch, a1_stoch, s_stoch, points = gha2_ES(ell, P0, P1, all_points=True, sigmaStep=1e-5) + a0_stoch, a1_stoch, s_stoch, points = gha2_ES(ell, P0, P1, all_points=True) out = html.Div([ html.Strong("Stochastisch (ES): "), @@ -691,7 +706,7 @@ def compute_gha2_approx(n2, beta0, lamb0, beta1, lamb1, ax, ay, b, method2): P0 = ell.ell2cart(beta0_rad, lamb0_rad) P1 = ell.ell2cart(beta1_rad, lamb1_rad) - a0_app, a1_app, s_app, points = gha2(ell, P0, P1, ds=1e-4, all_points=True) + a0_app, a1_app, s_app, points = gha2_approx(ell, P0, P1, ds=1e-4, all_points=True) out = html.Div([ html.Strong("Approximiert: "), @@ -706,6 +721,7 @@ def compute_gha2_approx(n2, beta0, lamb0, beta1, lamb1, ax, ay, b, method2): Input("input-ax", "value"), Input("input-ay", "value"), Input("input-b", "value"), + Input("dropdown-coors-type", "value"), Input("store-gha1-ana", "data"), Input("store-gha1-num", "data"), Input("store-gha1-stoch", "data"), @@ -714,13 +730,13 @@ def compute_gha2_approx(n2, beta0, lamb0, beta1, lamb1, ax, ay, b, method2): Input("store-gha2-stoch", "data"), Input("store-gha2-approx", "data"), ) -def render_all(ax, ay, b, store_gha1_ana, store_gha1_num, store_gha1_stoch, store_gha1_approx, store_gha2_num, store_gha2_stoch, store_gha2_approx): +def render_all(ax, ay, b, coords_type, store_gha1_ana, store_gha1_num, store_gha1_stoch, store_gha1_approx, store_gha2_num, store_gha2_stoch, store_gha2_approx): if None in (ax, ay, b): return go.Figure() ell = EllipsoidTriaxial(ax, ay, b) fig = ellipsoid_figure(ell, title="") - fig = figure_constant_lines(fig, ell, "ell") + fig = figure_constant_lines(fig, ell, coords_type) def add_from_store(fig, store): if not store: