diff --git a/GHA_triaxial/gha2_num.py b/GHA_triaxial/gha2_num.py index 8a9f313..c243f4a 100644 --- a/GHA_triaxial/gha2_num.py +++ b/GHA_triaxial/gha2_num.py @@ -17,124 +17,6 @@ def sph_azimuth(beta1, lam1, beta2, lam2): a += 2 * np.pi return a -def BETA_LAMBDA(ell, 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) - LAMBDA = (ell.ax**2 * np.sin(lamb)**2 + ell.ay**2 * np.cos(lamb)**2) / (ell.Ex**2 - ell.Ee**2 * np.cos(lamb)**2) - - # Erste Ableitungen von ΒETA und LAMBDA - BETA_ = (ell.ax**2 * ell.Ey**2 * np.sin(2*beta)) / (ell.Ex**2 - ell.Ey**2 * np.sin(beta)**2)**2 - LAMBDA_ = - (ell.b**2 * ell.Ee**2 * np.sin(2*lamb)) / (ell.Ex**2 - ell.Ee**2 * np.cos(lamb)**2)**2 - - # Zweite Ableitungen von ΒETA und LAMBDA - BETA__ = ((2 * ell.ax**2 * ell.Ey**4 * np.sin(2*beta)**2) / (ell.Ex**2 - ell.Ey**2 * np.sin(beta)**2)**3) + ((2 * ell.ax**2 * ell.Ey**2 * np.cos(2*beta)) / (ell.Ex**2 - ell.Ey**2 * np.sin(beta)**2)**2) - LAMBDA__ = (((2 * ell.b**2 * ell.Ee**4 * np.sin(2*lamb)**2) / (ell.Ex**2 - ell.Ee**2 * np.cos(lamb)**2)**3) - - ((2 * ell.b**2 * ell.Ee**2 * np.sin(2*lamb)) / (ell.Ex**2 - ell.Ee**2 * np.cos(lamb)**2)**2)) - - E = BETA * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2) - F = 0 - G = LAMBDA * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2) - - # Erste Ableitungen von E und G - E_beta = BETA_ * (ell.Ey**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(lamb)**2) - BETA * ell.Ey**2 * np.sin(2*beta) - E_lamb = BETA * ell.Ee**2 * np.sin(2*lamb) - - G_beta = - LAMBDA * ell.Ey**2 * np.sin(2*beta) - G_lamb = LAMBDA_ * (ell.Ey**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(lamb)**2) + LAMBDA * ell.Ee**2 * np.sin(2*lamb) - - # Zweite Ableitungen von E und G - E_beta_beta = BETA__ * (ell.Ey**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(lamb)**2) - 2 * BETA_ * ell.Ey**2 * np.sin(2*beta) - 2 * BETA * ell.Ey**2 * np.cos(2*beta) - E_beta_lamb = BETA_ * ell.Ee**2 * np.sin(2*lamb) - E_lamb_lamb = 2 * BETA * ell.Ee**2 * np.cos(2*lamb) - - G_beta_beta = - 2 * LAMBDA * ell.Ey**2 * np.cos(2*beta) - G_beta_lamb = - LAMBDA_ * ell.Ey**2 * np.sin(2*beta) - G_lamb_lamb = LAMBDA__ * (ell.Ey**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(lamb)**2) + 2 * LAMBDA_ * ell.Ee**2 * np.sin(2*lamb) + 2 * LAMBDA * ell.Ee**2 * np.cos(2*lamb) - - return (BETA, LAMBDA, E, G, - BETA_, LAMBDA_, BETA__, LAMBDA__, - E_beta, E_lamb, G_beta, G_lamb, - E_beta_beta, E_beta_lamb, E_lamb_lamb, - G_beta_beta, G_beta_lamb, G_lamb_lamb) - -def p_coef(beta, lamb): - - (BETA, LAMBDA, E, G, - BETA_, LAMBDA_, BETA__, LAMBDA__, - E_beta, E_lamb, G_beta, G_lamb, - E_beta_beta, E_beta_lamb, E_lamb_lamb, - G_beta_beta, G_beta_lamb, G_lamb_lamb) = BETA_LAMBDA(ell, beta, lamb) - - p_3 = - 0.5 * (E_lamb / G) - p_2 = (G_beta / G) - 0.5 * (E_beta / E) - p_1 = 0.5 * (G_lamb / G) - (E_lamb / E) - p_0 = 0.5 * (G_beta / E) - - p_33 = - 0.5 * ((E_beta_lamb * G - E_lamb * G_beta) / (G ** 2)) - p_22 = ((G * G_beta_beta - G_beta * G_beta) / (G ** 2)) - 0.5 * ((E * E_beta_beta - E_beta * E_beta) / (E ** 2)) - p_11 = 0.5 * ((G * G_beta_lamb - G_beta * G_lamb) / (G ** 2)) - ((E * E_beta_lamb - E_beta * E_lamb) / (E ** 2)) - p_00 = 0.5 * ((E * G_beta_beta - E_beta * G_beta) / (E ** 2)) - - return (BETA, LAMBDA, E, G, - p_3, p_2, p_1, p_0, - p_33, p_22, p_11, p_00) - -def buildODElamb(): - def ODE(lamb, v): - beta, beta_p, X3, X4 = v - - (BETA, LAMBDA, E, G, - p_3, p_2, p_1, p_0, - p_33, p_22, p_11, p_00) = p_coef(beta, lamb) - - dbeta = beta_p - dbeta_p = p_3 * beta_p ** 3 + p_2 * beta_p ** 2 + p_1 * beta_p + p_0 - dX3 = X4 - dX4 = (p_33 * beta_p ** 3 + p_22 * beta_p ** 2 + p_11 * beta_p + p_00) * X3 + \ - (3 * p_3 * beta_p ** 2 + 2 * p_2 * beta_p + p_1) * X4 - return np.array([dbeta, dbeta_p, dX3, dX4]) - - return ODE - -def q_coef(beta, lamb): - - (BETA, LAMBDA, E, G, - BETA_, LAMBDA_, BETA__, LAMBDA__, - E_beta, E_lamb, G_beta, G_lamb, - E_beta_beta, E_beta_lamb, E_lamb_lamb, - G_beta_beta, G_beta_lamb, G_lamb_lamb) = BETA_LAMBDA(ell, beta, lamb) - - q_3 = - 0.5 * (G_beta / E) - q_2 = (E_lamb / E) - 0.5 * (G_lamb / G) - q_1 = 0.5 * (E_beta / E) - (G_beta / G) - q_0 = 0.5 * (E_lamb / G) - - q_33 = - 0.5 * ((E * G_beta_lamb - E_lamb * G_lamb) / (E ** 2)) - q_22 = ((E * E_lamb_lamb - E_lamb * E_lamb) / (E ** 2)) - 0.5 * ((G * G_lamb_lamb - G_lamb * G_lamb) / (G ** 2)) - q_11 = 0.5 * ((E * E_beta_lamb - E_beta * E_lamb) / (E ** 2)) - ((G * G_beta_lamb - G_beta * G_lamb) / (G ** 2)) - q_00 = 0.5 * ((E_lamb_lamb * G - E_lamb * G_lamb) / (G ** 2)) - - return (BETA, LAMBDA, E, G, - q_3, q_2, q_1, q_0, - q_33, q_22, q_11, q_00) - -def buildODEbeta(): - def ODE(beta, v): - lamb, lamb_p, Y3, Y4 = v - - (BETA, LAMBDA, E, G, - q_3, q_2, q_1, q_0, - q_33, q_22, q_11, q_00) = q_coef(beta, lamb) - - dlamb = lamb_p - dlamb_p = q_3 * lamb_p ** 3 + q_2 * lamb_p ** 2 + q_1 * lamb_p + q_0 - dY3 = Y4 - dY4 = (q_33 * lamb_p ** 3 + q_22 * lamb_p ** 2 + q_11 * lamb_p + q_00) * Y3 + \ - (3 * q_3 * lamb_p ** 2 + 2 * q_2 * lamb_p + q_1) * Y4 - - return np.array([dlamb, dlamb_p, dY3, dY4]) - - return ODE # Panou 2013 def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float, lamb_2: float, @@ -155,6 +37,137 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float """ # h_x, h_y, h_e entsprechen E_x, E_y, E_e + 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) + LAMBDA = (ell.ax ** 2 * np.sin(lamb) ** 2 + ell.ay ** 2 * np.cos(lamb) ** 2) / ( + ell.Ex ** 2 - ell.Ee ** 2 * np.cos(lamb) ** 2) + + # Erste Ableitungen von ΒETA und LAMBDA + BETA_ = (ell.ax ** 2 * ell.Ey ** 2 * np.sin(2 * beta)) / (ell.Ex ** 2 - ell.Ey ** 2 * np.sin(beta) ** 2) ** 2 + LAMBDA_ = - (ell.b ** 2 * ell.Ee ** 2 * np.sin(2 * lamb)) / (ell.Ex ** 2 - ell.Ee ** 2 * np.cos(lamb) ** 2) ** 2 + + # Zweite Ableitungen von ΒETA und LAMBDA + BETA__ = ((2 * ell.ax ** 2 * ell.Ey ** 4 * np.sin(2 * beta) ** 2) / ( + ell.Ex ** 2 - ell.Ey ** 2 * np.sin(beta) ** 2) ** 3) + ( + (2 * ell.ax ** 2 * ell.Ey ** 2 * np.cos(2 * beta)) / ( + ell.Ex ** 2 - ell.Ey ** 2 * np.sin(beta) ** 2) ** 2) + LAMBDA__ = (((2 * ell.b ** 2 * ell.Ee ** 4 * np.sin(2 * lamb) ** 2) / ( + ell.Ex ** 2 - ell.Ee ** 2 * np.cos(lamb) ** 2) ** 3) - + ((2 * ell.b ** 2 * ell.Ee ** 2 * np.sin(2 * lamb)) / ( + ell.Ex ** 2 - ell.Ee ** 2 * np.cos(lamb) ** 2) ** 2)) + + E = BETA * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2) + F = 0 + G = LAMBDA * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2) + + # Erste Ableitungen von E und G + E_beta = BETA_ * ( + ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2) - BETA * ell.Ey ** 2 * np.sin( + 2 * beta) + E_lamb = BETA * ell.Ee ** 2 * np.sin(2 * lamb) + + G_beta = - LAMBDA * ell.Ey ** 2 * np.sin(2 * beta) + G_lamb = LAMBDA_ * ( + ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2) + LAMBDA * ell.Ee ** 2 * np.sin( + 2 * lamb) + + # Zweite Ableitungen von E und G + E_beta_beta = BETA__ * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin( + lamb) ** 2) - 2 * BETA_ * ell.Ey ** 2 * np.sin(2 * beta) - 2 * BETA * ell.Ey ** 2 * np.cos(2 * beta) + E_beta_lamb = BETA_ * ell.Ee ** 2 * np.sin(2 * lamb) + E_lamb_lamb = 2 * BETA * ell.Ee ** 2 * np.cos(2 * lamb) + + G_beta_beta = - 2 * LAMBDA * ell.Ey ** 2 * np.cos(2 * beta) + G_beta_lamb = - LAMBDA_ * ell.Ey ** 2 * np.sin(2 * beta) + G_lamb_lamb = LAMBDA__ * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin( + lamb) ** 2) + 2 * LAMBDA_ * ell.Ee ** 2 * np.sin(2 * lamb) + 2 * LAMBDA * ell.Ee ** 2 * np.cos(2 * lamb) + + return (BETA, LAMBDA, E, G, + BETA_, LAMBDA_, BETA__, LAMBDA__, + E_beta, E_lamb, G_beta, G_lamb, + E_beta_beta, E_beta_lamb, E_lamb_lamb, + G_beta_beta, G_beta_lamb, G_lamb_lamb) + + def p_coef(beta, lamb): + + (BETA, LAMBDA, E, G, + BETA_, LAMBDA_, BETA__, LAMBDA__, + E_beta, E_lamb, G_beta, G_lamb, + E_beta_beta, E_beta_lamb, E_lamb_lamb, + G_beta_beta, G_beta_lamb, G_lamb_lamb) = BETA_LAMBDA(beta, lamb) + + p_3 = - 0.5 * (E_lamb / G) + p_2 = (G_beta / G) - 0.5 * (E_beta / E) + p_1 = 0.5 * (G_lamb / G) - (E_lamb / E) + p_0 = 0.5 * (G_beta / E) + + p_33 = - 0.5 * ((E_beta_lamb * G - E_lamb * G_beta) / (G ** 2)) + p_22 = ((G * G_beta_beta - G_beta * G_beta) / (G ** 2)) - 0.5 * ((E * E_beta_beta - E_beta * E_beta) / (E ** 2)) + p_11 = 0.5 * ((G * G_beta_lamb - G_beta * G_lamb) / (G ** 2)) - ((E * E_beta_lamb - E_beta * E_lamb) / (E ** 2)) + p_00 = 0.5 * ((E * G_beta_beta - E_beta * G_beta) / (E ** 2)) + + return (BETA, LAMBDA, E, G, + p_3, p_2, p_1, p_0, + p_33, p_22, p_11, p_00) + + def buildODElamb(): + def ODE(lamb, v): + beta, beta_p, X3, X4 = v + + (BETA, LAMBDA, E, G, + p_3, p_2, p_1, p_0, + p_33, p_22, p_11, p_00) = p_coef(beta, lamb) + + dbeta = beta_p + dbeta_p = p_3 * beta_p ** 3 + p_2 * beta_p ** 2 + p_1 * beta_p + p_0 + dX3 = X4 + dX4 = (p_33 * beta_p ** 3 + p_22 * beta_p ** 2 + p_11 * beta_p + p_00) * X3 + \ + (3 * p_3 * beta_p ** 2 + 2 * p_2 * beta_p + p_1) * X4 + return np.array([dbeta, dbeta_p, dX3, dX4]) + + return ODE + + def q_coef(beta, lamb): + + (BETA, LAMBDA, E, G, + BETA_, LAMBDA_, BETA__, LAMBDA__, + E_beta, E_lamb, G_beta, G_lamb, + E_beta_beta, E_beta_lamb, E_lamb_lamb, + G_beta_beta, G_beta_lamb, G_lamb_lamb) = BETA_LAMBDA(beta, lamb) + + q_3 = - 0.5 * (G_beta / E) + q_2 = (E_lamb / E) - 0.5 * (G_lamb / G) + q_1 = 0.5 * (E_beta / E) - (G_beta / G) + q_0 = 0.5 * (E_lamb / G) + + q_33 = - 0.5 * ((E * G_beta_lamb - E_lamb * G_lamb) / (E ** 2)) + q_22 = ((E * E_lamb_lamb - E_lamb * E_lamb) / (E ** 2)) - 0.5 * ((G * G_lamb_lamb - G_lamb * G_lamb) / (G ** 2)) + q_11 = 0.5 * ((E * E_beta_lamb - E_beta * E_lamb) / (E ** 2)) - ((G * G_beta_lamb - G_beta * G_lamb) / (G ** 2)) + q_00 = 0.5 * ((E_lamb_lamb * G - E_lamb * G_lamb) / (G ** 2)) + + return (BETA, LAMBDA, E, G, + q_3, q_2, q_1, q_0, + q_33, q_22, q_11, q_00) + + def buildODEbeta(): + def ODE(beta, v): + lamb, lamb_p, Y3, Y4 = v + + (BETA, LAMBDA, E, G, + q_3, q_2, q_1, q_0, + q_33, q_22, q_11, q_00) = q_coef(beta, lamb) + + dlamb = lamb_p + dlamb_p = q_3 * lamb_p ** 3 + q_2 * lamb_p ** 2 + q_1 * lamb_p + q_0 + dY3 = Y4 + dY4 = (q_33 * lamb_p ** 3 + q_22 * lamb_p ** 2 + q_11 * lamb_p + q_00) * Y3 + \ + (3 * q_3 * lamb_p ** 2 + 2 * q_2 * lamb_p + q_1) * Y4 + + return np.array([dlamb, dlamb_p, dY3, dY4]) + + return ODE if lamb_1 != lamb_2: N = n @@ -164,7 +177,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float if abs(dlamb) < 1e-15: beta_0 = 0.0 else: - (_, _, E1, G1, *_) = BETA_LAMBDA(ell, beta_1, lamb_1) + (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1) beta_0 = np.sqrt(G1 / E1) * cot(alpha0_sph) ode_lamb = buildODElamb() @@ -196,7 +209,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float return False, None, None, None alpha0_sph = sph_azimuth(beta_1, lamb_1, beta_2, lamb_2) - (_, _, E1, G1, *_) = BETA_LAMBDA(ell, beta_1, lamb_1) + (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1) beta_p0_sph = np.sqrt(G1 / E1) * cot(alpha0_sph) guesses = [ @@ -220,7 +233,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float integrand = np.zeros(N + 1) for i in range(N + 1): - (_, _, Ei, Gi, *_) = BETA_LAMBDA(ell, beta_arr_c[i], lamb_arr_c[i]) + (_, _, 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 @@ -253,9 +266,9 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float beta_p_arr[i] = state[1] (_, _, E1, G1, - *_) = BETA_LAMBDA(ell, beta_arr[0], lamb_arr[0]) + *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) (_, _, E2, G2, - *_) = BETA_LAMBDA(ell, beta_arr[-1], lamb_arr[-1]) + *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1]) alpha_1 = arccot(np.sqrt(E1 / G1) * beta_p_arr[0]) alpha_2 = arccot(np.sqrt(E2 / G2) * beta_p_arr[-1]) @@ -263,7 +276,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float integrand = np.zeros(N + 1) for i in range(N + 1): (_, _, Ei, Gi, - *_) = BETA_LAMBDA(ell, beta_arr[i], lamb_arr[i]) + *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i]) integrand[i] = np.sqrt(Ei * beta_p_arr[i] ** 2 + Gi) h = abs(dlamb) / N @@ -341,9 +354,9 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float # Azimute (BETA1, LAMBDA1, E1, G1, - *_) = BETA_LAMBDA(ell, beta_arr[0], lamb_arr[0]) + *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) (BETA2, LAMBDA2, E2, G2, - *_) = BETA_LAMBDA(ell, beta_arr[-1], lamb_arr[-1]) + *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1]) alpha_1 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA1 / BETA1) * lambda_p_arr[0]) alpha_2 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA2 / BETA2) * lambda_p_arr[-1]) @@ -351,7 +364,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float integrand = np.zeros(N + 1) for i in range(N + 1): (_, _, Ei, Gi, - *_) = BETA_LAMBDA(ell, beta_arr[i], lamb_arr[i]) + *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i]) integrand[i] = np.sqrt(Ei + Gi * lambda_p_arr[i] ** 2) h = abs(dbeta) / N