From 4d533ab4645fa75d9ffa8b4a72464415b4c3c9fc Mon Sep 17 00:00:00 2001
From: Noe Brucy <noe.brucy@cea.fr>
Date: Thu, 11 Feb 2021 09:06:25 +0100
Subject: [PATCH] [postprocessor] [bug] take dx and dy into account in
 gradients

---
 postprocessor.py | 28 ++++++++++++++++++----------
 1 file changed, 18 insertions(+), 10 deletions(-)

diff --git a/postprocessor.py b/postprocessor.py
index 8009518..0783270 100644
--- a/postprocessor.py
+++ b/postprocessor.py
@@ -902,17 +902,25 @@ class PostProcessor(HDF5Container):
             dy = (im_extent[3] - im_extent[2]) / map_size
 
             # Physical coordinates of the center of the cells
-            x = np.linspace(im_extent[0], im_extent[1], map_size) + 0.5 * dx
-            y = np.linspace(im_extent[2], im_extent[3], map_size) + 0.5 * dy
+            x = (
+                np.linspace(im_extent[0], im_extent[1], map_size)
+                + 0.5 * dx
+                - pos_star[0]
+            )
+            y = (
+                np.linspace(im_extent[2], im_extent[3], map_size)
+                + 0.5 * dy
+                - pos_star[1]
+            )
 
             xx, yy = np.meshgrid(x, y)
-            rr = np.sqrt((xx - pos_star[0]) ** 2 + (yy - pos_star[1]) ** 2)
+            rr = np.sqrt(xx ** 2 + yy ** 2)
 
             # Compute kappa**2 = (2omega/R)*d(R**2*omega)/dR
-            Gx, Gy = np.gradient(rr ** 2 * omega)
+            Gy, Gx = np.gradient(rr ** 2 * omega, dy, dx, edge_order=2)
             Gr = (xx * Gx + yy * Gy) / rr
             kappa_square = (2 * omega / rr) * Gr
-            kappa = np.square(kappa_square)
+            kappa = np.sqrt(kappa_square)
 
         map_Q = (np.sqrt(cs2) * kappa) / (np.pi * self.G * coldens)
 
@@ -1303,22 +1311,22 @@ class PostProcessor(HDF5Container):
         dy = (im_extent[3] - im_extent[2]) / map_size
 
         # Physical coordinates of the center of the cells
-        x = np.linspace(im_extent[0], im_extent[1], map_size) + 0.5 * dx
-        y = np.linspace(im_extent[2], im_extent[3], map_size) + 0.5 * dy
+        x = np.linspace(im_extent[0], im_extent[1], map_size) + 0.5 * dx - pos_star[0]
+        y = np.linspace(im_extent[2], im_extent[3], map_size) + 0.5 * dy - pos_star[1]
 
         xx, yy = np.meshgrid(x, y)
-        rr = np.sqrt((xx - pos_star[0]) ** 2 + (yy - pos_star[1]) ** 2)
+        rr = np.sqrt(xx ** 2 + yy ** 2)
 
         # Rotational support
         R = vphi ** 2 / rr
 
         # Equilibrium
-        Gvrx, Gvry = np.gradient(vr)
+        Gvry, Gvrx = np.gradient(vr, dy, dx, edge_order=2)
         gradvr = (xx * Gvrx + yy * Gvry) / rr
         dvr = gradvr + vr * gradvr  # Complete derivative
 
         # Thermal support
-        GPx, GPy = np.gradient(Pz)
+        GPy, GPx = np.gradient(Pz, dy, dx, edge_order=2)
         gradPr = (xx * GPx + yy * GPy) / rr
         fP = gradPr / rho