Changeset 807 for codes/icosagcm/devel/Python/dynamico/maps.py

Timestamp:

02/08/19 18:53:34 (5 years ago)

Author:

jisesh

Message:

devel/Python : polar projection + Williamson test case 2 on LAM

File:

• codes/icosagcm/devel/Python/dynamico/maps.py (modified) (2 diffs)

codes/icosagcm/devel/Python/dynamico/maps.py

@@ -14 +14 @@
     def map(self, x, y):
         """x,y are coordinates in the reference domain."""
-        return self.xc + self.scale*x, self.yc + self.scale*y 
+        return self.xc + self.scale*x, self.yc + self.scale*y
     def map_factor(self, lon,lat): return 0.*lon + self.scale
+    def map_angle(self, lon, lat): return 0.*lon
 
 class SphereMap:
@@ -27 +28 @@
     def map(self, lon,lat): return lon, lat
     def map_factor(self, lon,lat): return 0.*lon + self.radius
+    def map_angle(self, lon, lat): return 0.*lon
 
 # https://en.wikipedia.org/wiki/Mercator_projection
+# https://en.wikipedia.org/wiki/Stereographic_projection
+class PolarStereoMap:
+    def __init__(self, radius,Omega, delta, phic=0., lonc=0.):
+        self.a, self.delta, self.Omega = radius, delta, Omega
+        self.phic, self.lonc = phic, lonc # latitude of domain center
+    def coriolis(self, lon, lat):
+        return 2.*self.Omega*np.sin(lat)
+    def map_factor(self,X,Y):
+        lon, lat, factor, angle = self.compute_map_factor_angle(X,Y)
+        return factor
+    def map_angle(self,X,Y):
+        lon, lat, factor, angle = self.compute_map_factor_angle(X,Y)
+        return angle
+    def map(self,X,Y):
+        lon, lat, factor, angle = self.compute_map_factor_angle(X,Y)
+        return lon, lat
+    def compute_map_factor_angle(self,X,Y):
+        """ X, Y : reference coordinates on a mesh where cells have unit size"""
+        alpha = .5*self.delta/self.a
+        X, Y = alpha*X, alpha*Y # scale X and Y to obtain desired resolution delta at domain center
+        # now X and Y are as in wikipedia
+        # we let (x,y) => (-y,x) to align X axis with equator when phic=0.
+        factor =  1./(1.+X*X+Y*Y)
+        y = 2.*X*factor
+        x = -2.*Y*factor
+        z = (1.-X*X-Y*Y)*factor
+        # compute d/dX to track mesh orientation
+        f2 = factor**2
+        dy = 2.*(1.-X*X+Y*Y)*f2
+        dx = 4.*X*Y*f2
+        dz = -4.*X*f2
+        # TODO : rotate x,y,z about polar axis
+        x,y,z = self.RotationPhi(x,y,z)
+        dx,dy,dz = self.RotationPhi(dx,dy,dz)
+        cos_angle = (dy*x-dx*y)/np.sqrt((dx**2+dy**2+dz**2)*(x**2+y**2))
+        sin_angle = (dz*(1.-z**2)-z*(dy*y+dx*x))/np.sqrt((dx**2+dy**2+dz**2)*(x**2+y**2))
+        angle = np.arctan2(dz*(1.-z**2)-z*(dy*y+dx*x), dy*x-dx*y)
+        err = cos_angle-np.cos(angle)
+        print 'Consistency check :', err.min(), err.max()
+        # TODO : rotate x,y,z so that North Pole gets to desired latitude
+        lat = np.arcsin(z)
+        lon = np.arctan2(y,x) + self.lonc
+        return lon, lat, self.delta*factor, angle
+    def RotationPhi(self,x2,y2,z2):
+        cos, sin = np.cos(self.phic), np.sin(self.phic)
+        x3 = sin * x2 + cos*z2
+        z3 = -cos * x2 + sin*z2
+        return x3,y2,z3