source: issm/trunk/src/py3/coordsystems/xy2ll.py@ 20500

import numpy as npy
from math import pi

def xy2ll(x, y, sgn, *args):
        '''
        XY2LL - converts xy to lat long
        
        Converts Polar  Stereographic (X, Y) coordinates for the polar regions to
        latitude and longitude Stereographic (X, Y) coordinates for the polar
        regions.
        Author: Michael P. Schodlok, December 2003 (map2xy.m)
        
        Usage:
           [lat, lon] = xy2ll(x, y, sgn);
           [lat, lon] = xy2ll(x, y, 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)
        '''

        #Get central_meridian and standard_parallel depending on hemisphere
        if len(args) == 2:
                delta = args[0]
                slat  = args[1]
        elif len(args) == 0:
                if sgn == 1:
                        delta = 45. 
                        slat = 70.
                        print('         xy2ll: creating coordinates in north polar stereographic (Std Latitude: 70degN Meridian: 45deg)')
                elif sgn == -1:
                        delta = 0.  
                        slat = 71.
                        print('         xy2ll: creating coordinates in south polar stereographic (Std Latitude: 71degS Meridian: 0deg)')
                else:
                        raise ValueError('sgn should be either +1 or -1')
        else:
                raise Exception('bad usage: type "help(xy2ll)" for details')

        # if x,y passed as lists, convert to numpy arrays
        if type(x) != "numpy.ndarray":
                x=npy.array(x)
        if type(y) != "numpy.ndarray":
                y=npy.array(y)

        ## 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 = npy.sqrt(ex2)
        
        sl = slat*pi/180.
        rho = npy.sqrt(x**2 + y**2)
        cm = npy.cos(sl) / npy.sqrt(1.0 - ex2 * (npy.sin(sl)**2))
        T = npy.tan((pi/4.0) - (sl/2.0)) / ((1.0 - ex*npy.sin(sl)) / (1.0 + ex*npy.sin(sl)))**(ex / 2.0)
        
        if abs(slat-90.) < 1.e-5:
                T = rho*npy.sqrt((1. + ex)**(1. + ex) * (1. - ex)**(1. - ex)) / 2. / re
        else:
                T = rho * T / (re * cm)
        
        chi = (pi / 2.0) - 2.0 * npy.arctan(T)
        lat = chi + ((ex2 / 2.0) + (5.0 * ex2**2.0 / 24.0) + (ex2**3.0 / 12.0)) * \
                npy.sin(2 * chi) + ((7.0 * ex2**2.0 / 48.0) + (29.0 * ex2**3 / 240.0)) * \
                npy.sin(4.0 * chi) + (7.0 * ex2**3.0 / 120.0) * npy.sin(6.0 * chi) 
        
        lat = sgn * lat
        lon = npy.arctan2(sgn * x,-sgn * y)
        lon = sgn * lon
        
        res1 = npy.nonzero(rho <= 0.1)[0]
        if len(res1) > 0:
                lat[res1] = 90. * sgn
                lon[res1] = 0.0
        
        lon = lon * 180. / pi
        lat = lat * 180. / pi
        lon = lon - delta 

        return lat, lon