diff --git a/GHA_triaxial/ES_gha2.py b/GHA_triaxial/ES_gha2.py index 71d8844..0c851e3 100644 --- a/GHA_triaxial/ES_gha2.py +++ b/GHA_triaxial/ES_gha2.py @@ -14,7 +14,7 @@ P_next: NDArray = None P_end: NDArray = None stepLen: float = None -def Bogenlaenge(P1, P2): +def Bogenlaenge(P1: NDArray, P2: NDArray) -> float: """ Berechnung der mittleren Bogenlänge zwischen zwei kartesischen Punkten :param P1: kartesische Koordinate Punkt 1 @@ -23,13 +23,13 @@ def Bogenlaenge(P1, P2): """ R1 = np.linalg.norm(P1) R2 = np.linalg.norm(P2) - R = 0.5*(R1 + R2) + R = 0.5 * (R1 + R2) theta = arccos(P1 @ P2 / (R1 * R2)) - s = R * theta + s = float(R * theta) return s -def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, stepLenTarget: float = None, sigmaStep: float = 1e-7, stopeval: int = 1000, maxSteps: int = 10000, all_points: bool = False): +def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, stepLenTarget: float = None, sigmaStep: float = 1e-5, stopeval: int = 1000, maxSteps: int = 10000, all_points: bool = False): """ Berechnen der 2. GHA mithilfe der CMA-ES. Die CMA-ES optimiert sukzessive einzelne Punkte, die einen definierten Abstand (stepLenTarget) zum vorherigen und den kürzesten @@ -49,7 +49,7 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, stepLenTarget: flo ell_ES = ell P_start = P0 P_end = Pk - if stepLenTarget == None: + if stepLenTarget is None: R0 = (ell.ax + ell.ay + ell.b) / 3 stepLenTarget = R0 * 1 / 600 stepLen = stepLenTarget @@ -92,14 +92,13 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, stepLenTarget: flo # sin(v); # sin(u)] auf # Einheitskugel - # % arg_u = max(-1, min(1, q(3))); + #arg_u = max(-1, min(1, q(3))); # # %Quadrantenabfrage - # % u0 = mod(asin(arg_u) + pi / 2, pi) - pi / 2; - # % v0 = atan2(q(2), q(1)); - # % xmean_init = [u0; - # v0]; - xmean_init = ell.cartonell(P_prev + stepLen * (P_end - P_prev) / np.linalg.norm(P_end - P_prev)) + #u0 = mod(asin(arg_u) + pi / 2, pi) - pi / 2; + #v0 = atan2(q(2), q(1)); + #xmean_init = [u0;v0]; + xmean_init = ell.point_onto_ellipsoid(P_prev + stepLen * (P_end - P_prev) / np.linalg.norm(P_end - P_prev)) # [~, ~, aux] = geoLength(xmean_init); # print('Startguess: d_step=%.3f (soll %.3f), d_to_target=%.3f\n', aux(1), stepLen, aux(2)); @@ -107,7 +106,7 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, stepLenTarget: flo print(f'[Punkt {i}] Optimiere nächsten Punkt: Restdistanz = {round(d_remain, 3)} m') xmean_init = np.array(ell_ES.cart2para(xmean_init)) - u, v = escma(geoLength,2, xmean_init, sigmaStep, -np.inf, stopeval) + u, v = escma(geoLength, N=2, xmean=xmean_init, sigma=sigmaStep, stopfitness=-np.inf, stopeval=stopeval) P_next = ell.para2cart(u, v) @@ -128,26 +127,30 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, stepLenTarget: flo totalLen += d_step P_prev = P_next - print('Maximale Schrittanzahl erreicht.') P_all.append(P_end) totalLen += Bogenlaenge(P_prev, P_end) - p0i = ell.cartonell(P0 + 10 * (P_all[1] - P0) / np.linalg.norm(P_all[1] - P0)) + p0i = ell.point_onto_ellipsoid(P0 + 10 * (P_all[1] - P0) / np.linalg.norm(P_all[1] - P0)) sigma0 = (p0i - P0) / np.linalg.norm(p0i - P0) alpha0 = sigma2alpha(ell_ES, sigma0, P0) - p1i = ell.cartonell(Pk - 10 * (Pk - P_all[-2]) / np.linalg.norm(Pk - P_all[-2])) + p1i = ell.point_onto_ellipsoid(Pk - 10 * (Pk - P_all[-2]) / np.linalg.norm(Pk - P_all[-2])) sigma1 = (Pk - p1i) / np.linalg.norm(Pk - p1i) alpha1 = sigma2alpha(ell_ES, sigma1, Pk) if all_points: return alpha0, alpha1, totalLen, np.array(P_all) else: - return alpha0, alpha1, totalLen + return alpha0, alpha1, totalLen -def geoLength(P_candidate): +def geoLength(P_candidate: Tuple) -> float: + """ + Berechung der Fitness eines Kandidaten anhand der Strecken + :param P_candidate: Kandidat in parametrischen Koordinaten + :return: Fitness-Wert + """ # P_candidate = [u;v] des naechsten Punktes. # Ziel: Distanz zum Ziel minimieren, aber Schrittlaenge ~ stepLenTarget erzwingen. u, v = P_candidate @@ -165,7 +168,7 @@ def geoLength(P_candidate): pen_step = ((d_step - stepLen) / stepLen)**2 # falls Punkt "weg" vom Ziel geht, extra bestrafen - pen_away = max(0, (d_to_target - d_prev_to_target) / stepLen)**2 + pen_away = max(0.0, (d_to_target - d_prev_to_target) / stepLen)**2 # Gewichtungen alpha = 1e2 @@ -175,13 +178,20 @@ def geoLength(P_candidate): f = d_to_target * (1 + alpha * pen_step + gamma * pen_away) # Für Debug / Extraktion - aux = [d_step, d_to_target] + # aux = [d_step, d_to_target] return f # , P_candidate, aux -def show_points(points: NDArray, pointsNum:NDArray, p0: NDArray, p1: NDArray): +def show_points(points: NDArray, pointsES: NDArray, p0: NDArray, p1: NDArray): + """ + Anzeigen der Punkte + :param points: wahre Punkte der Linie + :param pointsES: Punkte der Linie aus ES + :param p0: wahrer Startpunkt + :param p1: wahrer Endpunkt + """ fig = go.Figure() - fig.add_scatter3d(x=pointsNum[:, 0], y=pointsNum[:, 1], z=pointsNum[:, 2], + fig.add_scatter3d(x=pointsES[:, 0], y=pointsES[:, 1], z=pointsES[:, 2], mode='lines', line=dict(color="green", width=3), name="Numerisch") fig.add_scatter3d(x=points[:, 0], y=points[:, 1], z=points[:, 2], mode='lines', line=dict(color="red", width=3), name="ES") @@ -207,11 +217,12 @@ if __name__ == '__main__': beta1, lamb1 = (0.7, 0.3) P1 = ell.ell2cart(beta1, lamb1) - alpha0, alpha1, s, betas, lambs = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=5000, all_points=True) + alpha0, alpha1, s_num, betas, lambs = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=10000, all_points=True) points_num = [] for beta, lamb in zip(betas, lambs): points_num.append(ell.ell2cart(beta, lamb)) points_num = np.array(points_num) - alpha0, alpha1, s, points = gha2_ES(ell, P0, P1) + alpha0, alpha1, s, points = gha2_ES(ell, P0, P1, all_points=True, sigmaStep=1e-5) + print(s - s_num) show_points(points, points_num, P0, P1) diff --git a/GHA_triaxial/approx_gha1.py b/GHA_triaxial/approx_gha1.py index f87180d..379e3ae 100644 --- a/GHA_triaxial/approx_gha1.py +++ b/GHA_triaxial/approx_gha1.py @@ -4,7 +4,17 @@ from panou import louville_constant, func_sigma_ell, gha1_ana import plotly.graph_objects as go import winkelumrechnungen as wu -def gha1(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: float, all_points: bool = False) -> Tuple[NDArray, float] | Tuple[NDArray, float, NDArray]: +def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: float, all_points: bool = False) -> Tuple[NDArray, float] | Tuple[NDArray, float, NDArray]: + """ + Berechung einer Näherungslösung der ersten Hauptaufgabe + :param ell: Ellipsoid + :param p0: Anfangspunkt + :param alpha0: Azimut im Anfangspunkt + :param s: Strecke bis zum Endpunkt + :param ds: Länge einzelner Streckenelemente + :param all_points: Ausgabe aller Punkte als Array? + :return: Endpunkt, Azimut im Endpunkt, optional alle Punkte + """ l0 = louville_constant(ell, p0, alpha0) points = [p0] alphas = [alpha0] @@ -17,7 +27,7 @@ def gha1(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: fl alpha1 = alphas[-1] sigma = func_sigma_ell(ell, p1, alpha1) p2 = p1 + ds_step * sigma - p2 = ell.cartonell(p2) + p2 = ell.point_onto_ellipsoid(p2) ds_step = np.linalg.norm(p2 - p1) points.append(p2) @@ -33,7 +43,13 @@ def gha1(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: fl else: return points[-1], alphas[-1] -def show_points(points, p0, p1): +def show_points(points: NDArray, p0: NDArray, p1: NDArray): + """ + Anzeigen der Punkte + :param points: Array aller approximierten Punkte + :param p0: Startpunkt + :param p1: wahrer Endpunkt + """ fig = go.Figure() fig.add_scatter3d(x=points[:, 0], y=points[:, 1], z=points[:, 2], @@ -59,6 +75,6 @@ if __name__ == '__main__': alpha0 = wu.deg2rad(90) s = 1000000 P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=60, maxPartCircum=32) - P1_app, alpha1_app, points = gha1(ell, P0, alpha0, s, ds=5000, all_points=True) + P1_app, alpha1_app, points = gha1_approx(ell, P0, alpha0, s, ds=5000, all_points=True) show_points(points, P0, P1_ana) print(np.linalg.norm(P1_app - P1_ana)) diff --git a/GHA_triaxial/approx_gha2.py b/GHA_triaxial/approx_gha2.py index b47ac5d..70879a1 100644 --- a/GHA_triaxial/approx_gha2.py +++ b/GHA_triaxial/approx_gha2.py @@ -27,7 +27,7 @@ def gha2(ell: EllipsoidTriaxial, p0: NDArray, p1: NDArray, ds: float, all_points new_points.append(points[i]) pi = points[i] + 1/2 * (points[i+1] - points[i]) - pi = ell.cartonell(pi) + pi = ell.point_onto_ellipsoid(pi) new_points.append(pi) @@ -39,11 +39,11 @@ def gha2(ell: EllipsoidTriaxial, p0: NDArray, p1: NDArray, ds: float, all_points if np.average(elements) < ds: break - p0i = ell.cartonell(p0 + ds/100 * (points[1] - p0) / np.linalg.norm(points[1] - p0)) + p0i = ell.point_onto_ellipsoid(p0 + ds / 100 * (points[1] - p0) / np.linalg.norm(points[1] - p0)) sigma0 = (p0i - p0) / np.linalg.norm(p0i - p0) alpha0 = sigma2alpha(ell, sigma0, p0) - p1i = ell.cartonell(p1 - ds/100 * (p1 - points[-2]) / np.linalg.norm(p1 - points[-2])) + p1i = ell.point_onto_ellipsoid(p1 - ds / 100 * (p1 - points[-2]) / np.linalg.norm(p1 - points[-2])) sigma1 = (p1 - p1i) / np.linalg.norm(p1 - p1i) alpha1 = sigma2alpha(ell, sigma1, p1) @@ -55,6 +55,13 @@ def gha2(ell: EllipsoidTriaxial, p0: NDArray, p1: NDArray, ds: float, all_points return alpha0, alpha1, s def show_points(points: NDArray, points_app: NDArray, p0: NDArray, p1: NDArray): + """ + Anzeigen der Punkte + :param points: wahre Punkte der Linie + :param points_app: approximierte Punkte der Linie + :param p0: wahrer Startpunkt + :param p1: wahrer Endpunkt + """ fig = go.Figure() fig.add_scatter3d(x=points[:, 0], y=points[:, 1], z=points[:, 2], diff --git a/GHA_triaxial/numeric_examples_karney.py b/GHA_triaxial/numeric_examples_karney.py index 5303b10..173683d 100644 --- a/GHA_triaxial/numeric_examples_karney.py +++ b/GHA_triaxial/numeric_examples_karney.py @@ -1,7 +1,13 @@ import random import winkelumrechnungen as wu +from typing import List, Tuple -def line2example(line): +def line2example(line: str) -> List: + """ + Line-String in Liste umwandeln + :param line: Line-String + :return: Liste mit Zahlenwerten + """ split = line.split() example = [float(value) for value in split[:7]] for i, value in enumerate(example): @@ -10,13 +16,16 @@ def line2example(line): # example[i] = value return example -def get_random_examples(num): +def get_random_examples(num: int, seed: int = None) -> List: """ + Rückgabe zufälliger Beispiele beta0, lamb0, alpha0, beta1, lamb1, alpha1, s12 - :param num: - :return: + :param num: Anzahl zufälliger Beispiele + :param seed: Random-Seed + :return: Liste mit Beispielen """ - # random.seed(42) + if seed is not None: + random.seed(seed) with open("Karney_2024_Testset.txt") as datei: lines = datei.readlines() examples = [] @@ -25,11 +34,12 @@ def get_random_examples(num): examples.append(example) return examples -def get_examples(l_i): +def get_examples(l_i: List) -> List: """ + Rückgabe ausgewählter Beispiele beta0, lamb0, alpha0, beta1, lamb1, alpha1, s12 - :param num: - :return: + :param l_i: Liste von Indizes + :return: Liste mit Beispielen """ with open("Karney_2024_Testset.txt") as datei: lines = datei.readlines() @@ -39,5 +49,6 @@ def get_examples(l_i): examples.append(example) return examples + if __name__ == "__main__": get_random_examples(10) \ No newline at end of file diff --git a/GHA_triaxial/numeric_examples_panou.py b/GHA_triaxial/numeric_examples_panou.py index 5f025cc..3040aa8 100644 --- a/GHA_triaxial/numeric_examples_panou.py +++ b/GHA_triaxial/numeric_examples_panou.py @@ -1,113 +1,134 @@ +from typing import Tuple, List + import winkelumrechnungen as wu import random 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(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(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(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(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(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(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) + (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(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(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(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(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(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(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) + (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(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(-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(-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(-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(-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(-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) + (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(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(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(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(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(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(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(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) + (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): +def get_example(table: int, example: int) -> Tuple: + """ + Rückgabe eines Beispiels + :param table: Tabellen-Nummer + :param example: Beispiel-Nummer + :return: Bespiel + """ table -= 1 example -= 1 + tables = get_tables() return tables[table][example] -def get_tables(): +def get_tables() -> List: + """ + Rückgabe aller Tabellen + :return: Alle Tabellen + """ return tables -def get_random_examples(num): - random.seed(42) +def get_random_examples(num: int, seed: int = None) -> List: + """ + Rückgabe zufäliger Beispiele + :param num: Anzahl Beispiele + :param seed: Random-Seed + :return: + """ + if seed is not None: + random.seed(seed) + examples = [] for i in range(num): table = random.randint(1, 4) @@ -120,7 +141,6 @@ def get_random_examples(num): return examples - if __name__ == "__main__": # test = get_example(1, 4) examples = get_random_examples(5) diff --git a/GHA_triaxial/panou_2013_2GHA_num.py b/GHA_triaxial/panou_2013_2GHA_num.py index 694b300..deb9212 100644 --- a/GHA_triaxial/panou_2013_2GHA_num.py +++ b/GHA_triaxial/panou_2013_2GHA_num.py @@ -10,7 +10,7 @@ from numpy.typing import NDArray # Panou 2013 def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float, lamb_2: float, n: int = 16000, epsilon: float = 10**-12, iter_max: int = 30, all_points: bool = False - ) -> Tuple[float, float, float]| Tuple[float, float, float, NDArray, NDArray]: + ) -> Tuple[float, float, float] | Tuple[float, float, float, NDArray, NDArray]: """ :param ell: triaxiales Ellipsoid @@ -252,7 +252,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float else: return alpha_1, alpha_2, s - else: # lamb_1 == lamb_2 + else: # lamb_1 == lamb_2 N = n dbeta = beta_2 - beta_1 diff --git a/Hansen_ES_CMA.py b/Hansen_ES_CMA.py index 290cc76..d90b602 100644 --- a/Hansen_ES_CMA.py +++ b/Hansen_ES_CMA.py @@ -11,7 +11,7 @@ def felli(x): def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-10, stopeval=None, func_args=(), func_kwargs=None, seed=None, - bestEver = np.inf, noImproveGen = 0, absTolImprove = 1e-10, maxNoImproveGen = 100, sigmaImprove = 1e-12): + bestEver = np.inf, noImproveGen = 0, absTolImprove = 1e-12, maxNoImproveGen = 100, sigmaImprove = 1e-12): if func_kwargs is None: func_kwargs = {} diff --git a/ellipsoide.py b/ellipsoide.py index 76bbc24..615bd42 100644 --- a/ellipsoide.py +++ b/ellipsoide.py @@ -154,8 +154,8 @@ class EllipsoidTriaxial: return cls(ax, ay, b) elif name == "Fiction": ax = 6000000 - ay = 5000000 - b = 4000000 + ay = 4000000 + b = 2000000 return cls(ax, ay, b) elif name == "KarneyTest2024": ax = sqrt(2) @@ -294,9 +294,9 @@ class EllipsoidTriaxial: def ell2cart_bektas(self, beta: float | NDArray, omega: float | NDArray) -> NDArray: """ Bektas 2015 - :param beta: - :param omega: - :return: + :param beta: elliptische Breite [rad] + :param omega: elliptische Länge [rad] + :return: Punkt in kartesischen Koordinaten """ x = self.ax * cos(omega) * sqrt((self.ax**2 - self.ay**2 * sin(beta)**2 - self.b**2 * cos(beta)**2) / (self.ax**2 - self.b**2)) y = self.ay * cos(beta) * sin(omega) @@ -307,9 +307,9 @@ class EllipsoidTriaxial: def ell2cart_karney(self, beta: float | NDArray, lamb: float | NDArray) -> NDArray: """ Karney 2025 Geographic Lib - :param beta: - :param lamb: - :return: + :param beta: elliptische Breite [rad] + :param lamb: elliptische Länge [rad] + :return: Punkt in kartesischen Koordinaten """ k = sqrt(self.ay**2 - self.b**2) / sqrt(self.ax**2 - self.b**2) k_ = sqrt(self.ax**2 - self.ay**2) / sqrt(self.ax**2 - self.b**2) @@ -321,10 +321,10 @@ class EllipsoidTriaxial: def cart2ell(self, point: NDArray, eps: float = 1e-12, maxI: int = 100) -> Tuple[float, float]: """ Panou, Korakitis 2019 3f. (num) - :param point: - :param eps: - :param maxI: - :return: + :param point: Punkt in kartesischen Koordinaten + :param eps: zu erreichende Genauigkeit + :param maxI: maximale Anzahl Iterationen + :return: elliptische Breite und Länge [rad] """ x, y, z = point beta, lamb = self.cart2ell_panou(point) @@ -427,10 +427,10 @@ class EllipsoidTriaxial: def cart2ell_bektas(self, point: NDArray, eps: float = 1e-12, maxI: int = 100) -> Tuple[float, float]: """ Bektas 2015 - :param point: - :param eps: - :param maxI: - :return: + :param point: Punkt in kartesischen Koordinaten + :param eps: zu erreichende Genauigkeit + :param maxI: maximale Anzahl Iterationen + :return: elliptische Breite und Länge [rad] """ x, y, z = point phi, lamb = self.cart2para(point) @@ -575,26 +575,70 @@ class EllipsoidTriaxial: return u, v def ell2para(self, beta: float, lamb: float) -> Tuple[float, float]: + """ + Umrechung von elliptischen in parametrische Koordinaten (über kartesische Koordinaten) + :param beta: elliptische Breite + :param lamb: elliptische Länge + :return: parametrische Koordinaten + """ cart = self.ell2cart(beta, lamb) return self.cart2para(cart) def para2ell(self, u: float, v: float) -> Tuple[float, float]: + """ + Umrechung von parametrischen in elliptische Koordinaten (über kartesische Koordinaten) + :param u: u + :param v: v + :return: elliptische Koordinaten + """ 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]: + """ + Umrechung von parametrischen in geodätische Koordinaten (über kartesische Koordinaten) + :param u: u + :param v: v + :param mode: ligas1, ligas2, oder ligas3 + :param maxIter: maximale Anzahl Iterationen + :param maxLoa: Level of Accuracy, das erreicht werden soll + :return: geodätische Koordinaten + """ cart = self.para2cart(u, v) return self.cart2geod(cart, mode, maxIter, maxLoa) def geod2para(self, phi: float, lamb: float, h: float) -> Tuple[float, float]: + """ + Umrechung von geodätischen in parametrische Koordinaten (über kartesische Koordinaten) + :param phi: u + :param lamb: v + :param h: geodätische Höhe + :return: parametrische Koordinaten + """ cart = self.geod2cart(phi, lamb, h) return self.cart2para(cart) def ell2geod(self, beta: float, lamb: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> Tuple[float, float, float]: + """ + Umrechung von elliptischen in geodätische Koordinaten (über kartesische Koordinaten) + :param beta: elliptische Breite + :param lamb: eliptische Länge + :param mode: ligas1, ligas2, oder ligas3 + :param maxIter: maximale Anzahl Iterationen + :param maxLoa: Level of Accuracy, das erreicht werden soll + :return: geodätische Koordinaten + """ cart = self.ell2cart(beta, lamb) return self.cart2geod(cart, mode, maxIter, maxLoa) def geod2ell(self, phi: float, lamb: float, h: float) -> Tuple[float, float]: + """ + Umrechung von geodätischen in elliptische Koordinaten (über kartesische Koordinaten) + :param phi: u + :param lamb: v + :param h: geodätische Höhe + :return: elliptische Koordinaten + """ cart = self.geod2cart(phi, lamb, h) return self.cart2ell(cart) @@ -610,9 +654,9 @@ class EllipsoidTriaxial: else: return False - def cartonell(self, point: NDArray) -> NDArray: + def point_onto_ellipsoid(self, point: NDArray) -> NDArray: """ - Berechnung des Lotpunktes auf einem Ellipsoiden + Berechnung des Lotpunktes entlang der Normalkrümmung auf einem Ellipsoiden :param point: Punkt in kartesischen Koordinaten, der gelotet werden soll :return: Lotpunkt in kartesischen Koordinaten """ @@ -620,7 +664,7 @@ class EllipsoidTriaxial: p = self. geod2cart(phi, lamb, 0) return p - def cartellh(self, point: NDArray, h: float) -> NDArray: + def cart_ellh(self, point: NDArray, h: float) -> NDArray: """ Punkt auf Ellipsoid hoch loten :param point: Punkt auf dem Ellipsoid diff --git a/utils.py b/utils.py index c8b6765..78403d3 100644 --- a/utils.py +++ b/utils.py @@ -13,7 +13,6 @@ def sigma2alpha(ell: EllipsoidTriaxial, sigma: NDArray, point: NDArray) -> float :param point: Punkt :return: Richtungswinkel """ - "" p, q = pq_ell(ell, point) P = float(p @ sigma) Q = float(q @ sigma) diff --git a/winkelumrechnungen.py b/winkelumrechnungen.py index 8985c21..eb877c2 100644 --- a/winkelumrechnungen.py +++ b/winkelumrechnungen.py @@ -10,13 +10,10 @@ def deg2gms(deg: float) -> list: :return: Winkel in Grad-Minuten-Sekunden :rtype: list """ - gra = deg // 1 - minu = gra % 1 - gra = gra // 1 - minu *= 60 - sek = minu % 1 - minu = minu // 1 - sek *= 60 + gra = int(deg) + minu_f = (deg - gra) * 60 + minu = int(minu_f) + sek = (minu_f - minu) * 60 return [gra, minu, sek] @@ -51,13 +48,10 @@ def gra2gms(gra: float) -> list: :rtype: list """ deg = gra2deg(gra) - gra = deg // 1 - minu = gra % 1 - gra = gra // 1 - minu *= 60 - sek = minu % 1 - minu = minu // 1 - sek *= 60 + gra = int(deg) + minu_f = (deg - gra) * 60 + minu = int(minu_f) + sek = (minu_f - minu) * 60 return [gra, minu, sek] @@ -114,12 +108,10 @@ def rad2gms(rad: float) -> list: :rtype: list """ deg = rad2deg(rad) - minu = deg % 1 - gra = deg // 1 - minu *= 60 - sek = minu % 1 - minu = minu // 1 - sek *= 60 + gra = int(deg) + minu_f = (deg - gra) * 60 + minu = int(minu_f) + sek = (minu_f - minu) * 60 return [gra, minu, sek]