```"""geodetic transforms to auxilary coordinate systems involving latitude"""
import typing

from .ellipsoid import Ellipsoid
from .utils import sanitize
from .rcurve import rcurve_transverse

try:
from numpy import radians, degrees, tan, sin, exp, pi, sqrt, inf, vectorize
from numpy import arctan as atan, arcsinh as asinh, arctanh as atanh  # noqa: A001

use_numpy = True
except ImportError:
from math import atan, radians, degrees, tan, sin, asinh, atanh, exp, pi, sqrt, inf

use_numpy = False

__all__ = [
"geodetic2isometric",
"isometric2geodetic",
"geodetic2rectifying",
"rectifying2geodetic",
"geodetic2conformal",
"conformal2geodetic",
"geodetic2parametric",
"parametric2geodetic",
"geodetic2geocentric",
"geocentric2geodetic",
"geodetic2authalic",
"authalic2geodetic",
"geod2geoc",
"geoc2geod",
]

if typing.TYPE_CHECKING:
from numpy import ndarray

def geoc2geod(geocentric_lat: "ndarray", geocentric_distance: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
"""
convert geocentric latitude to geodetic latitude, consider mean sea level altitude

like Matlab geoc2geod()

Parameters
----------
geocentric_lat : "ndarray"
geocentric latitude
geocentric_distance: "ndarray"
distance from planet center, meters (NOT altitude above ground!)
ell : Ellipsoid, optional
reference ellipsoid (default WGS84)
deg : bool, optional

Returns
-------
geodetic_lat : "ndarray"
geodetic latiude

References
----------
Long, S.A.T. "General-Altitude Transformation between Geocentric
and Geodetic Coordinates. Celestial Mechanics (12), 2, p. 225-230 (1975)
doi: 10.1007/BF01230214"
"""
geocentric_lat, ell = sanitize(geocentric_lat, ell, deg)

r = geocentric_distance / ell.semimajor_axis

geodetic_lat = (
geocentric_lat
+ (sin(2 * geocentric_lat) / r) * ell.flattening
+ ((1 / r ** 2 + 1 / (4 * r)) * sin(4 * geocentric_lat)) * ell.flattening ** 2
)

return degrees(geodetic_lat) if deg else geodetic_lat

def geodetic2geocentric(geodetic_lat: "ndarray", alt_m: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
"""
convert geodetic latitude to geocentric latitude on spheroid surface

like Matlab geocentricLatitude() with alt_m = 0
like Matlab geod2geoc()

Parameters
----------
geodetic_lat : "ndarray"
geodetic latitude
alt_m: "ndarray"
altitude above ellipsoid
ell : Ellipsoid, optional
reference ellipsoid (default WGS84)
deg : bool, optional

Returns
-------
geocentric_lat : "ndarray"
geocentric latiude

Notes
-----
Equations from J. P. Snyder, "Map Projections - A Working Manual",
US Geological Survey Professional Paper 1395, US Government Printing
Office, Washington, DC, 1987, pp. 13-18.
"""
geodetic_lat, ell = sanitize(geodetic_lat, ell, deg)
r = rcurve_transverse(geodetic_lat, ell, deg=False)
geocentric_lat = atan((1 - ell.eccentricity ** 2 * (r / (r + alt_m))) * tan(geodetic_lat))

return degrees(geocentric_lat) if deg else geocentric_lat

geod2geoc = geodetic2geocentric

def geocentric2geodetic(geocentric_lat: "ndarray", alt_m: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
"""
converts from geocentric latitude to geodetic latitude

like Matlab geodeticLatitudeFromGeocentric() when alt_m = 0
like Matlab geod2geoc() but with sea level altitude rather than planet center distance

Parameters
----------
geocentric_lat : "ndarray"
geocentric latitude
alt_m: "ndarray"
altitude above ellipsoid
ell : Ellipsoid, optional
reference ellipsoid (default WGS84)
deg : bool, optional

Returns
-------
geodetic_lat : "ndarray"
geodetic latiude

Notes
-----
Equations from J. P. Snyder, "Map Projections - A Working Manual",
US Geological Survey Professional Paper 1395, US Government Printing
Office, Washington, DC, 1987, pp. 13-18.
"""
geocentric_lat, ell = sanitize(geocentric_lat, ell, deg)
r = rcurve_transverse(geocentric_lat, ell, deg=False)
geodetic_lat = atan(tan(geocentric_lat) / (1 - ell.eccentricity ** 2 * (r / (r + alt_m))))

return degrees(geodetic_lat) if deg else geodetic_lat

def geodetic2isometric_point(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
geodetic_lat, ell = sanitize(geodetic_lat, ell, deg)

e = ell.eccentricity

if abs(geodetic_lat - pi / 2) <= 1e-9:
isometric_lat = inf
elif abs(-geodetic_lat - pi / 2) <= 1e-9:
isometric_lat = -inf
else:
isometric_lat = asinh(tan(geodetic_lat)) - e * atanh(e * sin(geodetic_lat))
# same results
# a1 = e * sin(geodetic_lat)
# y = (1 - a1) / (1 + a1)
# a2 = pi / 4 + geodetic_lat / 2
# isometric_lat = log(tan(a2) * (y ** (e / 2)))
# isometric_lat = log(tan(a2)) + e/2 * log((1-e*sin(geodetic_lat)) / (1+e*sin(geodetic_lat)))

return degrees(isometric_lat) if deg else isometric_lat

def geodetic2isometric(geodetic_lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
"""
computes isometric latitude on an ellipsoid

like Matlab map.geodesy.IsometricLatitudeConverter.forward()

Parameters
----------
lat : "ndarray"
geodetic latitude
ell : Ellipsoid, optional
reference ellipsoid (default WGS84)
deg : bool, optional

Returns
-------
isolat : "ndarray"
isometric latiude

Notes
-----
Isometric latitude is an auxiliary latitude proportional to the spacing
of parallels of latitude on an ellipsoidal mercator projection.
Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes,
School of Mathematical and Geospatial Sciences, RMIT University,
January 2010
"""
if use_numpy:
fun = vectorize(geodetic2isometric_point)
return fun(geodetic_lat, ell, deg)[()]
else:
return geodetic2isometric_point(geodetic_lat, ell, deg)

def isometric2geodetic(isometric_lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
"""
converts from isometric latitude to geodetic latitude

like Matlab map.geodesy.IsometricLatitudeConverter.inverse()

Parameters
----------
isometric_lat : "ndarray"
isometric latitude
ell : Ellipsoid, optional
reference ellipsoid (default WGS84)
deg : bool, optional

Returns
-------
geodetic_lat : "ndarray"
geodetic latiude

Notes
-----
Equations from J. P. Snyder, "Map Projections - A Working Manual",
US Geological Survey Professional Paper 1395, US Government Printing
Office, Washington, DC, 1987, pp. 13-18.
"""
# NOT sanitize for isometric2geo
if deg:

conformal_lat = 2 * atan(exp(isometric_lat)) - (pi / 2)
geodetic_lat = conformal2geodetic(conformal_lat, ell, deg=False)

return degrees(geodetic_lat) if deg else geodetic_lat

def conformal2geodetic(conformal_lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
"""
converts from conformal latitude to geodetic latitude

like Matlab map.geodesy.ConformalLatitudeConverter.inverse()

Parameters
----------
conformal_lat : "ndarray"
conformal latitude
ell : Ellipsoid, optional
reference ellipsoid (default WGS84)
deg : bool, optional

Returns
-------
geodetic_lat : "ndarray"
geodetic latiude

Notes
-----
Equations from J. P. Snyder, "Map Projections - A Working Manual",
US Geological Survey Professional Paper 1395, US Government Printing
Office, Washington, DC, 1987, pp. 13-18.
"""
conformal_lat, ell = sanitize(conformal_lat, ell, deg)

e = ell.eccentricity
f1 = e ** 2 / 2 + 5 * e ** 4 / 24 + e ** 6 / 12 + 13 * e ** 8 / 360
f2 = 7 * e ** 4 / 48 + 29 * e ** 6 / 240 + 811 * e ** 8 / 11520
f3 = 7 * e ** 6 / 120 + 81 * e ** 8 / 1120
f4 = 4279 * e ** 8 / 161280

geodetic_lat = (
conformal_lat
+ f1 * sin(2 * conformal_lat)
+ f2 * sin(4 * conformal_lat)
+ f3 * sin(6 * conformal_lat)
+ f4 * sin(8 * conformal_lat)
)

return degrees(geodetic_lat) if deg else geodetic_lat

def geodetic2conformal(geodetic_lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
"""
converts from geodetic latitude to conformal latitude

like Matlab map.geodesy.ConformalLatitudeConverter.forward()

Parameters
----------
geodetic_lat : "ndarray"
geodetic latitude
ell : Ellipsoid, optional
reference ellipsoid (default WGS84)
deg : bool, optional

Returns
-------
conformal_lat : "ndarray"
conformal latiude

Notes
-----
Equations from J. P. Snyder, "Map Projections - A Working Manual",
US Geological Survey Professional Paper 1395, US Government Printing
Office, Washington, DC, 1987, pp. 13-18.

"""
if use_numpy:
fun = vectorize(geodetic2conformal_point)
return fun(geodetic_lat, ell, deg)[()]
else:
return geodetic2conformal_point(geodetic_lat, ell, deg)

def geodetic2conformal_point(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
geodetic_lat, ell = sanitize(geodetic_lat, ell, deg)

e = ell.eccentricity
f1 = 1 - e * sin(geodetic_lat)
f2 = 1 + e * sin(geodetic_lat)
f3 = 1 - sin(geodetic_lat)
f4 = 1 + sin(geodetic_lat)

#  compute conformal latitudes with correction for points at +90
try:
conformal_lat = 2 * atan(sqrt((f4 / f3) * ((f1 / f2) ** e))) - (pi / 2)
except ZeroDivisionError:
conformal_lat = pi / 2

return degrees(conformal_lat) if deg else conformal_lat

# %% rectifying
def geodetic2rectifying(geodetic_lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
"""
converts from geodetic latitude to rectifying latitude

like Matlab map.geodesy.RectifyingLatitudeConverter.forward()

Parameters
----------
geodetic_lat : "ndarray"
geodetic latitude
ell : Ellipsoid, optional
reference ellipsoid (default WGS84)
deg : bool, optional

Returns
-------
rectifying_lat : "ndarray"
rectifying latiude

Notes
-----
Equations from J. P. Snyder, "Map Projections - A Working Manual",
US Geological Survey Professional Paper 1395, US Government Printing
Office, Washington, DC, 1987, pp. 13-18.

"""
geodetic_lat, ell = sanitize(geodetic_lat, ell, deg)

n = ell.thirdflattening
f1 = 3 * n / 2 - 9 * n ** 3 / 16
f2 = 15 * n ** 2 / 16 - 15 * n ** 4 / 32
f3 = 35 * n ** 3 / 48
f4 = 315 * n ** 4 / 512

rectifying_lat = (
geodetic_lat
- f1 * sin(2 * geodetic_lat)
+ f2 * sin(4 * geodetic_lat)
- f3 * sin(6 * geodetic_lat)
+ f4 * sin(8 * geodetic_lat)
)

return degrees(rectifying_lat) if deg else rectifying_lat

def rectifying2geodetic(rectifying_lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
"""
converts from rectifying latitude to geodetic latitude

like Matlab map.geodesy.RectifyingLatitudeConverter.inverse()

Parameters
----------
rectifying_lat : "ndarray"
latitude
ell : Ellipsoid, optional
reference ellipsoid (default WGS84)
deg : bool, optional

Returns
-------
geodetic_lat : "ndarray"
geodetic latiude

Notes
-----
Equations from J. P. Snyder, "Map Projections - A Working Manual",
US Geological Survey Professional Paper 1395, US Government Printing
Office, Washington, DC, 1987, pp. 13-18.
"""
rectifying_lat, ell = sanitize(rectifying_lat, ell, deg)

n = ell.thirdflattening
f1 = 3 * n / 2 - 27 * n ** 3 / 32
f2 = 21 * n ** 2 / 16 - 55 * n ** 4 / 32
f3 = 151 * n ** 3 / 96
f4 = 1097 * n ** 4 / 512

geodetic_lat = (
rectifying_lat
+ f1 * sin(2 * rectifying_lat)
+ f2 * sin(4 * rectifying_lat)
+ f3 * sin(6 * rectifying_lat)
+ f4 * sin(8 * rectifying_lat)
)

return degrees(geodetic_lat) if deg else geodetic_lat

# %% authalic
def geodetic2authalic(geodetic_lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
"""
converts from geodetic latitude to authalic latitude

like Matlab map.geodesy.AuthalicLatitudeConverter.forward()

Parameters
----------
geodetic_lat : "ndarray"
geodetic latitude
ell : Ellipsoid, optional
reference ellipsoid (default WGS84)
deg : bool, optional

Returns
-------
authalic_lat : "ndarray"
authalic latiude

Notes
-----
Equations from J. P. Snyder, "Map Projections - A Working Manual",
US Geological Survey Professional Paper 1395, US Government Printing
Office, Washington, DC, 1987, pp. 13-18.

"""
geodetic_lat, ell = sanitize(geodetic_lat, ell, deg)

e = ell.eccentricity
f1 = e ** 2 / 3 + 31 * e ** 4 / 180 + 59 * e ** 6 / 560
f2 = 17 * e ** 4 / 360 + 61 * e ** 6 / 1260
f3 = 383 * e ** 6 / 45360

authalic_lat = geodetic_lat - f1 * sin(2 * geodetic_lat) + f2 * sin(4 * geodetic_lat) - f3 * sin(6 * geodetic_lat)

return degrees(authalic_lat) if deg else authalic_lat

def authalic2geodetic(authalic_lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
"""
converts from authalic latitude to geodetic latitude

like Matlab map.geodesy.AuthalicLatitudeConverter.inverse()

Parameters
----------
authalic_lat : "ndarray"
latitude
ell : Ellipsoid, optional
reference ellipsoid (default WGS84)
deg : bool, optional

Returns
-------
geodetic_lat : "ndarray"
geodetic latiude

Notes
-----
Equations from J. P. Snyder, "Map Projections - A Working Manual",
US Geological Survey Professional Paper 1395, US Government Printing
Office, Washington, DC, 1987, pp. 13-18.
"""
authalic_lat, ell = sanitize(authalic_lat, ell, deg)
e = ell.eccentricity
f1 = e ** 2 / 3 + 31 * e ** 4 / 180 + 517 * e ** 6 / 5040
f2 = 23 * e ** 4 / 360 + 251 * e ** 6 / 3780
f3 = 761 * e ** 6 / 45360

geodetic_lat = authalic_lat + f1 * sin(2 * authalic_lat) + f2 * sin(4 * authalic_lat) + f3 * sin(6 * authalic_lat)

return degrees(geodetic_lat) if deg else geodetic_lat

# %% parametric
def geodetic2parametric(geodetic_lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
"""
converts from geodetic latitude to parametric latitude

like Matlab parametriclatitude()

Parameters
----------
geodetic_lat : "ndarray"
geodetic latitude
ell : Ellipsoid, optional
reference ellipsoid (default WGS84)
deg : bool, optional

Returns
-------
parametric_lat : "ndarray"
parametric latiude

Notes
-----
Equations from J. P. Snyder, "Map Projections - A Working Manual",
US Geological Survey Professional Paper 1395, US Government Printing
Office, Washington, DC, 1987, pp. 13-18.

"""
geodetic_lat, ell = sanitize(geodetic_lat, ell, deg)

parametric_lat = atan(sqrt(1 - (ell.eccentricity) ** 2) * tan(geodetic_lat))

return degrees(parametric_lat) if deg else parametric_lat

def parametric2geodetic(parametric_lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
"""
converts from parametric latitude to geodetic latitude

like Matlab geodeticLatitudeFromParametric()

Parameters
----------
parametric_lat : "ndarray"
latitude
ell : Ellipsoid, optional
reference ellipsoid (default WGS84)
deg : bool, optional

Returns
-------
geodetic_lat : "ndarray"
geodetic latiude

Notes
-----
Equations from J. P. Snyder, "Map Projections - A Working Manual",
US Geological Survey Professional Paper 1395, US Government Printing
Office, Washington, DC, 1987, pp. 13-18.
"""
parametric_lat, ell = sanitize(parametric_lat, ell, deg)

geodetic_lat = atan(tan(parametric_lat) / sqrt(1 - (ell.eccentricity) ** 2))

return degrees(geodetic_lat) if deg else geodetic_lat
```