Merge pull request #233 from ArnaudLalague/ionosphere_radius

Update plot_ionosphere.py
fmihpc · Jun 11, 2024 · be5bcea · be5bcea
2 parents 5cb9e3b + 324ad91
commit be5bcea
Showing 1 changed file with 12 additions and 3 deletions.
diff --git a/pyPlots/plot_ionosphere.py b/pyPlots/plot_ionosphere.py
@@ -95,7 +95,7 @@ def plot_ionosphere(filename=None,
 
     '''
 
-    IONOSPHERE_RADIUS=6471e3 # R_E + 100 km
+    #IONOSPHERE_RADIUS=6471e3 # R_E + 100 km
     # Verify the location of this watermark image
     watermarkimage=os.path.join(os.path.dirname(__file__), 'logo_color.png')
     watermarkimageblack=os.path.join(os.path.dirname(__file__), 'logo_black.png')
@@ -260,6 +260,15 @@ def plot_ionosphere(filename=None,
 
     # Read ionosphere mesh node coordinates
     coords = f.get_ionosphere_node_coords()
+
+    # Read the ionosphere radius from the VLSV file
+    try:
+        ionosphere_radius = f.read_parameter("ionosphere_radius")
+    except:
+        # Recalculate the radius of the ionospheric grid
+        ionosphere_radius = np.max(np.linalg.norm(coords, axis=-1))
+
+
     # Read ionosphere mesh connectivity
     elements = f.get_ionosphere_element_corners()
     # Read selected variable valus
@@ -282,9 +291,9 @@ def plot_ionosphere(filename=None,
        cb_title_use = cb_title_use + "\,["+datamap_unit_latex+"]"
 
     if viewdir > 0:
-        r=np.degrees(np.arccos(coords[:,2]/IONOSPHERE_RADIUS))
+        r=np.degrees(np.arccos(coords[:,2]/ionosphere_radius))
     else:
-        r=np.degrees(np.arccos(-coords[:,2]/IONOSPHERE_RADIUS))
+        r=np.degrees(np.arccos(-coords[:,2]/ionosphere_radius))
     theta=np.arctan2(coords[:,1],coords[:,0])
 
     # Project nodes and elements into view plane