Index: /issm/trunk-jpl/src/m/coordsystems/ll2xy.py
===================================================================
--- /issm/trunk-jpl/src/m/coordsystems/ll2xy.py	(revision 17587)
+++ /issm/trunk-jpl/src/m/coordsystems/ll2xy.py	(revision 17587)
@@ -0,0 +1,62 @@
+import numpy as npy 
+
+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 'Info: creating coordinates in polar stereographic (Std Latitude: 70N Meridian: 45)'
+	else: 
+		delta = central_meridian
+		slat = standard_parallel
+		print 'Info: creating coordinates in 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 = npy.sqrt(ex2)
+	
+	latitude = npy.abs(lat) * npy.pi/180.
+	longitude = (lon + delta) * npy.pi/180.
+	
+	# compute X and Y in grid coordinates.
+	T = npy.tan(npy.pi/4-latitude/2) / ((1-ex*npy.sin(latitude))/(1+ex*npy.sin(latitude)))**(ex/2)
+	
+	if (90 - slat) <  1.e-5:
+		rho = 2.*re*T/npy.sqrt((1.+ex)**(1.+ex)*(1.-ex)**(1.-ex))
+	else:
+		sl  = slat*npy.pi/180.
+		tc  = npy.tan(npy.pi/4.-sl/2.)/((1.-ex*npy.sin(sl))/(1.+ex*npy.sin(sl)))**(ex/2.)
+		mc  = npy.cos(sl)/npy.sqrt(1.0-ex2*(npy.sin(sl)**2))
+		rho = re*mc*T/tc
+	
+	y = -rho * sgn * npy.cos(sgn*longitude)
+	x =  rho * sgn * npy.sin(sgn*longitude)
+
+	cnt1=npy.nonzero(latitude>= npy.pi/2.)[0]
+	
+	if cnt1:
+		x[cnt1,0] = 0.0
+		y[cnt1,0] = 0.0
+	return x,y