from numpy import *
import scipy as sp
from ellipsoide import EllipsoidBiaxial
from typing import Tuple


def gha1(re: EllipsoidBiaxial, phi0: float, lamb0: float, alpha0:float, s: float) -> Tuple[float, float, float]:
    psi0 = re.phi2psi(phi0)
    clairant = arcsin(cos(psi0) * sin(alpha0))
    sigma0 = arcsin(sin(psi0) / cos(clairant))

    sqrt_sigma = lambda sigma: sqrt(1 + re.e_ ** 2 * cos(clairant) ** 2 * sin(sigma) ** 2)
    int_sqrt_sigma = lambda sigma: sp.integrate.quad(sqrt_sigma, sigma0, sigma)[0]

    f_sigma1_i = lambda sigma1_i: (int_sqrt_sigma(sigma1_i) - s / re.b)

    sigma1_0 = sigma0 + s / re.a
    sigma1 = sp.optimize.newton(f_sigma1_i, sigma1_0)
    psi1 = arcsin(cos(clairant) * sin(sigma1))
    phi1 = re.psi2phi(psi1)
    alpha1 = arcsin(sin(clairant) / cos(psi1))

    f_d_lambda = lambda sigma: sin(clairant) * sqrt_sigma(sigma) / (1 - cos(clairant) ** 2 * sin(sigma) ** 2)
    d_lambda = sqrt(1-re.e**2) * sp.integrate.quad(f_d_lambda, sigma0, sigma1)[0]
    lamb1 = lamb0 + d_lambda
    return phi1, lamb1, alpha1