;https://arm-doe.github.io/pyart/_modules/pyart/core/transforms.html#cartesian_to_geographic_aeqd
;def cartesian_to_geographic_aeqd(x,y,lon0,lat0,r=6370997


pro cartesian_to_geographic_aeqd,x,y,lon_0,lat_0,lon_deg,lat_deg

;  Constant
R=6370997.0


;Azimuthal equidistant Cartesian to geographic coordinate transform.

;Transform a set of Cartesian/Cartographic coordinates (x, y) to
;geographic coordinate system (lat, lon) using a azimuthal equidistant
;map projection [1]_.

;\\=floor

;lat = \\arcsin(\\cos(c) * \\sin(lat_0) + (y * \\sin(c) * \\cos(lat_0) / \\rho))

;lon = lon_0 + \\arctan2(x * \\sin(c),\\rho * \\cos(lat_0) * \\cos(c) - y * \\sin(lat_0) * \\sin(c))

;\\rho = \\sqrt(x^2 + y^2)

;c = \\rho / R

;Where x, y are the Cartesian position from the center of projection;
;lat, lon the corresponding latitude and longitude; lat_0, lon_0 are the
;latitude and longitude of the center of the projection; R is the radius of
;the earth (defaults to ~6371 km). lon is adjusted to be between -180 and
;180.

;Parameters
;----------
;x, y : array-like
;Cartesian coordinates in the same units as R, typically meters.
;lon_0, lat_0 : float
;Longitude and latitude, in degrees, of the center of the projection.
;R : float, optional
;Earth radius in the same units as x and y. The default value is in
;units of meters.

;Returns
;-------
;lon, lat : array
;Longitude and latitude of Cartesian coordinates in degrees.

;References
;----------
;.. [1] Snyder, J. P. Map Projections--A Working Manual. U. S. Geological
;Survey Professional Paper 1395, 1987, pp. 191-202.

;x & y need no modification
;x = np.atleast_1d(np.asarray(x))
;y = np.atleast_1d(np.asarray(y))

;lat_0_rad = np.deg2rad(lat_0)
;lon_0_rad = np.deg2rad(lon_0)
lat_0_rad=lat_0*!dtor
lon_0_rad=lon_0*!dtor


rho = sqrt(x*x + y*y)
c = rho / R

;with warnings.catch_warnings():
;# division by zero may occur here but is properly addressed below so
;# the warnings can be ignored
;warnings.simplefilter("ignore", RuntimeWarning)

lat_rad = asin(cos(c) * sin(lat_0_rad) + y * sin(c) * cos(lat_0_rad) / rho)
lat_deg = lat_rad*!radeg
;# fix cases where the distance from the center of the projection is zero
idx=where(rho eq 0,count)
if count gt 0 then lat_deg[idx]=lat_0

x1 = x * sin(c)
x2 = rho*cos(lat_0_rad)*cos(c) - y*sin(lat_0_rad)*sin(c)
lon_rad = lon_0_rad + atan(x1, x2)
lon_deg = lon_rad*!radeg

;# Longitudes should be from -180 to 180 degrees
idx=where(lon_deg gt 180,c)
if c gt 0 then lon_deg[idx]=lon_deg[idx]-360.
idx=where(lon_deg lt -180,c)
if c gt 0 then lon_deg[idx]=lon_deg[idx]+360.

return
end