191 lines
6.2 KiB
Python
191 lines
6.2 KiB
Python
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
|
|
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_old(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, num):
|
|
phi, lamb, h = ell.cart2geod(point, "ligas3")
|
|
x, y, z = ell.geod2cart(phi, lamb, 0)
|
|
p, q = ell.p_q(x, y, z)
|
|
|
|
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 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
|
|
|
|
|
|
def checkLiouville(ell: ellipsoide.EllipsoidTriaxial, points):
|
|
constantValues = []
|
|
for point in points:
|
|
x = point[1]
|
|
dxds = point[2]
|
|
y = point[3]
|
|
dyds = point[4]
|
|
z = point[5]
|
|
dzds = point[6]
|
|
|
|
values = ell.p_q(x, y, z)
|
|
p = values["p"]
|
|
q = values["q"]
|
|
t1 = values["t1"]
|
|
t2 = values["t2"]
|
|
|
|
P = p[0]*dxds + p[1]*dyds + p[2]*dzds
|
|
Q = q[0]*dxds + q[1]*dyds + q[2]*dzds
|
|
alpha = arctan2(P, Q)
|
|
|
|
c = ell.ay**2 - (t1 * sin(alpha)**2 + t2 * cos(alpha)**2)
|
|
constantValues.append(c)
|
|
pass
|
|
|
|
|
|
def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM):
|
|
"""
|
|
Panou, Korakitits 2020, 5ff.
|
|
:param ell:
|
|
:param point:
|
|
:param alpha0:
|
|
:param s:
|
|
:param maxM:
|
|
:return:
|
|
"""
|
|
x, y, z = point
|
|
x_m = [x]
|
|
y_m = [y]
|
|
z_m = [z]
|
|
|
|
# erste Ableitungen (7-8)
|
|
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 = 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] * 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] +
|
|
1 / (1-ell.ee**2) ** 2 * y_m[p-i] * y_m[i] +
|
|
1 / (1-ell.ex**2) ** 2 * z_m[p-i] * z_m[i]) for i in range(0, p+1)])
|
|
|
|
# h Ableitungen (7)
|
|
h_ = lambda q: np.sum([comb(q, j) * (x_m[q-j+1] * x_m[j+1] +
|
|
1 / (1 - ell.ee ** 2) * y_m[q-j+1] * y_m[j+1] +
|
|
1 / (1 - ell.ex ** 2) * z_m[q-j+1] * z_m[j+1]) for j in range(0, q+1)])
|
|
|
|
# h/H Ableitungen (6)
|
|
hH_ = lambda t: 1/H_(0) * (h_(t) - fact(t) *
|
|
np.sum([H_(t+1-l) / (fact(t+1-l) * fact(l-1)) * hH_t[l-1] for l in range(1, t+1)]))
|
|
|
|
# xm, ym, zm Ableitungen (6)
|
|
x_ = lambda m: -np.sum([comb(m-2, k) * hH_t[m-2-k] * x_m[k] for k in range(0, m-2+1)])
|
|
y_ = lambda m: -1 / (1-ell.ee**2) * np.sum([comb(m-2, k) * hH_t[m-2-k] * y_m[k] for k in range(0, m-2+1)])
|
|
z_ = lambda m: -1 / (1-ell.ex**2) * np.sum([comb(m-2, k) * hH_t[m-2-k] * z_m[k] for k in range(0, m-2+1)])
|
|
|
|
hH_t = []
|
|
a_m = []
|
|
b_m = []
|
|
c_m = []
|
|
for m in range(0, maxM+1):
|
|
if m >= 2:
|
|
hH_t.append(hH_(m-2))
|
|
x_m.append(x_(m))
|
|
y_m.append(y_(m))
|
|
z_m.append(z_(m))
|
|
a_m.append(x_m[m] / fact(m))
|
|
b_m.append(y_m[m] / fact(m))
|
|
c_m.append(z_m[m] / fact(m))
|
|
|
|
# am, bm, cm (6)
|
|
x_s = 0
|
|
for a in reversed(a_m):
|
|
x_s = x_s * s + a
|
|
y_s = 0
|
|
for b in reversed(b_m):
|
|
y_s = y_s * s + b
|
|
z_s = 0
|
|
for c in reversed(c_m):
|
|
z_s = z_s * s + c
|
|
|
|
return x_s, y_s, z_s
|
|
pass
|
|
|
|
|
|
if __name__ == "__main__":
|
|
# ell = ellipsoide.EllipsoidTriaxial.init_name("Eitschberger1978")
|
|
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°
|
|
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)
|
|
|
|
# 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 |