import numpy as np


class Lambert:
    def __init__(
        self,
        X0=1.7e6,
        Y0=2.2e6,
        lambda0=3,
        phi0=43,
        phi1=42.25,
        phi2=43.75,
        a=6378137,
        e=0.081819190842622,
    ):
        self._X0 = X0
        self._Y0 = Y0

        self._lambda0 = np.radians(lambda0)
        self._phi0 = np.radians(phi0)
        self._phi1 = np.radians(phi1)
        self._phi2 = np.radians(phi2)

        self._a = a
        self._e = e

        self._set_n()
        self._set_rho0()

    def cartesian(self, lam, phi):
        lam = np.radians(lam)
        phi = np.radians(phi)

        theta = self._n * (lam - self._lambda0)
        rho = self._rho(phi)
        rho0 = self._rho(self._phi0)

        X = self._X0 + rho * np.sin(theta)
        Y = self._Y0 + rho0 - rho * np.cos(theta)

        return X, Y

    def _set_n(self):
        self._n = (
            np.log(np.cos(self._phi2) / np.cos(self._phi1))
            + (
                np.log(
                    (1 - (self._e * np.sin(self._phi1) ** 2))
                    / (1 - (self._e * np.sin(self._phi2) ** 2))
                )
                / 2
            )
        ) / (
            np.log(
                (
                    np.tan(self._phi1 / 2 + np.pi / 4)
                    * (
                        (1 - self._e * np.sin(self._phi2))
                        * (1 + self._e * np.sin(self._phi1))
                    )
                    ** (self._e / 2)
                )
                / (
                    np.tan(self._phi2 / 2 + np.pi / 4)
                    * (
                        (1 + self._e * np.sin(self._phi2))
                        * (1 - self._e * np.sin(self._phi1))
                    )
                    ** (self._e / 2)
                )
            )
        )

    def _set_rho0(self):
        self._rho0 = (
            self._a
            * np.cos(self._phi1)
            / (self._n * np.sqrt(1 - self._e**2 * np.sin(self._phi1) ** 2))
            * (
                np.tan(self._phi1 / 2 + np.pi / 4)
                * (
                    (1 - self._e * np.sin(self._phi1))
                    / (1 + self._e * np.sin(self._phi1))
                )
                ** (self._e / 2)
            )
            ** self._n
        )

    def _rho(self, phi):
        return (
            self._rho0
            * (
                1
                / np.tan(phi / 2 + np.pi / 4)
                * ((1 + self._e * np.sin(phi)) / (1 - self._e * np.sin(phi)))
                ** (self._e / 2)
            )
            ** self._n
        )