import numpy as np


def ll2xy(lat, lon, sgn=-1, central_meridian=0, standard_parallel=71):
    '''
    LL2XY - converts lat lon to polar stereographic

   Converts from geodetic latitude and longitude to Polar
   Stereographic (X, Y) coordinates for the polar regions.
   Author: Michael P. Schodlok, December 2003 (map2ll)

   Usage:
      x, y = ll2xy(lat, lon, sgn)
      x, y = ll2xy(lat, lon, sgn, central_meridian, standard_parallel)

 - sgn = Sign of latitude	1 : north latitude (default is mer = 45 lat = 70)
 						   -1 : south latitude (default is mer = 0  lat = 71)
    '''

    assert sgn == 1 or sgn == -1, 'error: sgn should be either 1 or -1'

    #Get central_meridian and standard_parallel depending on hemisphere
    if sgn == 1:
        delta = 45
        slat = 70
        print('        ll2xy: creating coordinates in north polar stereographic (Std Latitude: 70N Meridian: 45)')
    else:
        delta = central_meridian
        slat = standard_parallel
        print('        ll2xy: creating coordinates in south polar stereographic (Std Latitude: 71S Meridian: 0)')

    # Conversion constant from degrees to radians
    #cde = 57.29577951
    # Radius of the earth in meters
    re = 6378.273 * 10**3
    # Eccentricity of the Hughes ellipsoid squared
    ex2 = .006693883
    # Eccentricity of the Hughes ellipsoid
    ex = np.sqrt(ex2)

    latitude = np.abs(lat) * np.pi / 180.
    longitude = (lon + delta) * np.pi / 180.

    # compute X and Y in grid coordinates.
    T = np.tan(np.pi / 4 - latitude / 2) / ((1 - ex * np.sin(latitude)) / (1 + ex * np.sin(latitude)))**(ex / 2)

    if (90 - slat) < 1.e-5:
        rho = 2. * re * T / np.sqrt((1. + ex)**(1. + ex) * (1. - ex)**(1. - ex))
    else:
        sl = slat * np.pi / 180.
        tc = np.tan(np.pi / 4. - sl / 2.) / ((1. - ex * np.sin(sl)) / (1. + ex * np.sin(sl)))**(ex / 2.)
        mc = np.cos(sl) / np.sqrt(1.0 - ex2 * (np.sin(sl)**2))
        rho = re * mc * T / tc

    y = -rho * sgn * np.cos(sgn * longitude)
    x = rho * sgn * np.sin(sgn * longitude)

    cnt1 = np.nonzero(latitude >= np.pi / 2.)[0]

    if cnt1:
        x[cnt1, 0] = 0.0
        y[cnt1, 0] = 0.0
    return x, y
