Magnetopause script updates

Documentation/sphinx/plot_streamline_magnetopause_2d.rst

@@ -0,0 +1,6 @@
+plot_streamline_magnetopause_2d
+-------------------------------
+
+.. automodule:: plot_streamline_magnetopause_2d
+	:members:
+

Documentation/sphinx/plot_streamline_magnetopause_3d.rst

@@ -0,0 +1,6 @@
+plot_streamline_magnetopause_3d
+-------------------------------
+
+.. automodule:: plot_streamline_magnetopause_3d
+	:members:
+

Documentation/sphinx/scripts.rst

@@ -21,21 +21,21 @@ gics
 
 ------------
 
-magnetopause2d
---------------
-:doc:`magnetopause2d`
+plot_streamline_magnetopause_2d
+-------------------------------
+:doc:`plot_streamline_magnetopause_2d`
 
-.. automodule:: magnetopause2d
-    :no-index:
+.. automodule:: plot_streamline_magnetopause_2d
+	:no-index:
 
 ------------
 
-magnetopause3d
---------------
-:doc:`magnetopause3d`
+plot_streamline_magnetopause_3d
+-------------------------------
+:doc:`plot_streamline_magnetopause_3d`
 
-.. automodule:: magnetopause3d
-    :no-index:
+.. automodule:: plot_streamline_magnetopause_3d
+	:no-index:
 
 ------------
 
@@ -81,8 +81,8 @@ tsyganenko
 
    biot_savart
    gics
-   magnetopause2d
-   magnetopause3d
+   plot_streamline_magnetopause_2d
+   plot_streamline_magnetopause_3d
    obliqueshock
    obliqueshock_nif
    shue

analysator/pyCalculations/calculations.py

@@ -61,3 +61,5 @@
 from non_maxwellianity import epsilon_M
 from null_lines import LMN_null_lines_FOTE
 from interpolator_amr import AMRInterpolator, supported_amr_interpolators
+from magnetopause_sw_streamline_2d import find_magnetopause_sw_streamline_2d
+from magnetopause_sw_streamline_3d import find_magnetopause_sw_streamline_3d

analysator/pyCalculations/fieldtracer.py

@@ -51,7 +51,7 @@ def static_field_tracer( vlsvReader, x0, max_iterations, dx, direction='+', bvar
    ''' Field tracer in a static frame
 
        :param vlsvReader:         An open vlsv file
-       :param x:                  Starting point for the field trace
+       :param x0:                  Starting point for the field trace
        :param max_iterations:     The maximum amount of iteractions before the algorithm stops
        :param dx:                 One iteration step length
        :param direction:          '+' or '-' or '+-' Follow field in the plus direction or minus direction
@@ -61,6 +61,10 @@ def static_field_tracer( vlsvReader, x0, max_iterations, dx, direction='+', bvar
        :returns:                  List of coordinates
    '''
 
+   # Check the starting point
+   if (np.shape(x0) != (3,)): 
+      raise ValueError("Starting point should be a single [x,y,z] coordinate in 1-dimensional list or array")
+
    if(bvar != 'B'):
      warnings.warn("User defined tracing variable detected. fg, volumetric variable results may not work as intended, use face-values instead.")
 
@@ -117,7 +121,7 @@ def static_field_tracer( vlsvReader, x0, max_iterations, dx, direction='+', bvar
       x = np.arange(mins[0], maxs[0], dcell[0]) + 0.5*dcell[0]
       y = np.arange(mins[1], maxs[1], dcell[1]) + 0.5*dcell[1]
       z = np.arange(mins[2], maxs[2], dcell[2]) + 0.5*dcell[2]
-      coordinates = np.array([x,y,z])
+      coordinates = [x,y,z]
       # Debug:
       if( len(x) != sizes[0] ):
          logging.info("SIZE WRONG: " + str(len(x)) + " " + str(sizes[0]))
@@ -142,7 +146,7 @@ def static_field_tracer( vlsvReader, x0, max_iterations, dx, direction='+', bvar
       x = np.arange(mins[0], maxs[0], dcell[0]) + 0.5*dcell[0]
       y = np.arange(mins[1], maxs[1], dcell[1]) + 0.5*dcell[1]
       z = np.arange(mins[2], maxs[2], dcell[2]) + 0.5*dcell[2]
-      coordinates = np.array([x,y,z])
+      coordinates = [x,y,z]
       # Debug:
       if( len(x) != sizes[0] ):
          logging.info("SIZE WRONG: " + str(len(x)) + " " + str(sizes[0]))
@@ -359,6 +363,10 @@ def static_field_tracer_3d( vlsvReader, seed_coords, max_iterations, dx, directi
    # Standardize input: (N,3) np.array
    if type(seed_coords) != np.ndarray:
       raise TypeError("Please give a numpy array.")
+
+   # Transform a single (3,) coordinate to (1,3) if needed
+   if np.shape(seed_coords)==(3,):
+      seed_coords = np.array([seed_coords])
 
    # Cache and read variables:
    vg = None

analysator/pyCalculations/magnetopause_sw_streamline_2d.py

@@ -0,0 +1,162 @@
+'''
+Finds the magnetopause position by tracing streamlines of the plasma flow for two-dimensional Vlasiator runs. Needs the yt package.
+'''
+
+import numpy as np
+import analysator as pt
+import yt
+
+
+def interpolate(streamline, x_points):
+    """Interpolates a single streamline for make_magnetopause().
+
+        :param streamline: a single streamline to be interpolated
+        :param x_points: points in the x-axis to use for interpolation
+        :returns: the streamline as numpy array of x,z coordinate points where the x-axis coordinates are the points given to the function
+    """
+
+    arr = np.array(streamline)
+
+    # set arrays for interpolation
+    xp = arr[:,0][::-1]
+    zp = arr[:,2][::-1]
+
+    #interpolate z coordinates
+    z_points = np.interp(x_points, xp, zp, left=np.NaN, right=np.NaN)
+
+    return np.array([x_points, z_points])
+
+
+def make_streamlines(vlsvfile, streamline_seeds=None,seeds_n=200, seeds_x0=20*6371000, seeds_range=[-5*6371000, 5*6371000], streamline_length=40*6371000):
+    """Traces streamlines of velocity field using the yt package.
+
+        :param vlsvfile: directory and file name of .vlsv data file to use for VlsvReader
+        :kword streamline_seeds: optional streamline starting points in numpy array (coordinates in meters including the y-coordinate 0.0)
+        :kword seeds_n: instead of streamline_seeds provide a number of streamlines to be traced 
+        :kword seeds_x0: instead of streamline_seeds provide an x-coordinate for streamline starting points 
+        :kword seeds_range: instead of streamline_seeds provide [min, max] range to use for streamline starting point z-coordinates
+        :kword streamline_length: streamline length
+
+        :returns: streamlines as numpy array
+    """
+
+    # bulk file
+    f = pt.vlsvfile.VlsvReader(file_name=vlsvfile)
+
+    # get box coordinates from data
+    [xmin, ymin, zmin, xmax, ymax, zmax] = f.get_spatial_mesh_extent()
+    [xsize, ysize, zsize] = f.get_spatial_mesh_size()
+    simext =[xmin,xmax,ymin,ymax,zmin,zmax]
+    sizes = np.array([xsize,ysize,zsize])
+    boxcoords=list(simext)
+
+    cellids = f.read_variable("CellID")
+
+    #Read the data from vlsv-file
+    Vx = f.read_variable("v", operator="x")
+    Vz = f.read_variable("v", operator="z")
+
+
+    #Re-shape variable data
+    Vxs=Vx[np.argsort(cellids)].reshape(f.get_spatial_mesh_size(), order="F")
+    Vys = np.zeros_like(Vxs)
+    Vzs=Vz[np.argsort(cellids)].reshape(f.get_spatial_mesh_size(), order="F")
+
+    data=dict(Vx=Vxs,Vy=Vys,Vz=Vzs)
+
+    #Create starting points for streamlines if they are not given
+    if streamline_seeds == None:
+        streamline_seeds = np.array([[seeds_x0, 0 ,i] for i in  np.linspace(seeds_range[0], seeds_range[1], seeds_n)])
+
+    #streamline_seeds = np.array(streamline_seeds)
+    #dataset in yt-form
+    yt_dataset = yt.load_uniform_grid(
+        data,
+        sizes,
+        bbox=np.array([[boxcoords[0], boxcoords[1]],
+                    [boxcoords[2],boxcoords[3]],
+                    [boxcoords[4],boxcoords[5]]]))
+
+
+    #data, seeds, dictionary positions, step size
+    streamlines = yt.visualization.api.Streamlines(yt_dataset, streamline_seeds,
+                                                    "Vx", "Vy", "Vz", length=streamline_length, direction=1)
+
+    #trace the streamlines with yt
+    streamlines.integrate_through_volume()
+    # return streamline positions
+    return np.array(streamlines.streamlines)
+
+
+def make_magnetopause(streamlines, end_x=-15*6371000, x_point_n=50):
+    """Finds the mangetopause location based on streamlines.
+
+        :param streams: streamlines (coordinates in m)
+        :kword end_x: tail end x-coordinate (how far along the negative x-axis the magnetopause is calculated)
+        :kword x_point_n: integer, how many x-axis points the magnetopause will be divided in between the subsolar point and tail
+
+        :returns:   the magnetopause position as coordinate points in numpy array
+    """
+
+    RE = 6371000
+
+    streampoints = np.reshape(streamlines, (streamlines.shape[0]*streamlines.shape[1], 3)) #all the points in one array
+
+    ## find the subsolar dayside point in the positive x-axis
+    ## do this by finding a stremline point on positive x axis closest to the Earth
+    x_axis_points = streampoints[(streampoints[:,2]< RE) & (streampoints[:,2]> -RE) & (streampoints[:,0]> 0)]
+    subsolar_x =np.min(x_axis_points[:,0])
+
+    ## define points in the x axis where to find magnetopause points on the yz-plane
+    x_points = np.linspace(subsolar_x, end_x, x_point_n)
+
+    ## interpolate more exact points for streamlines at exery x_point
+    new_streampoints = np.zeros((len(x_points), len(streamlines), 1)) # new array for keeping interpolated streamlines in form streamlines_new[x_point, streamline, z-coordinate] 
+
+    for i,stream in enumerate(streamlines):
+        interpolated_streamline = interpolate(stream, x_points)
+        for j in range(0, len(x_points)):
+            new_streampoints[j, i,:] = interpolated_streamline[1,j]
+
+
+    ## start making the magnetopause
+    ## in each x_point, find the closest streamline to x-axis in the positive and negative z-axis
+
+    pos_z_mpause = np.zeros((len(x_points), 2))
+    neg_z_mpause =  np.zeros((len(x_points), 2))
+
+    for i, x_point in enumerate(x_points):
+        pos = new_streampoints[i, new_streampoints[i,:] > 0]
+        neg = new_streampoints[i, new_streampoints[i,:] < 0]
+
+        if (pos.size == 0) or (neg.size == 0):
+            raise ValueError('No streamlines found for x axis point, try adding streamlines or checking the x_points')
+
+        # find points closest to x-axis and save found points
+        pos_z_mpause[i] = [x_point, pos[pos.argmin()]]
+        neg_z_mpause[i] = [x_point, neg[neg.argmax()]]
+
+    magnetopause = np.concatenate((pos_z_mpause[::-1], np.array([[subsolar_x, 0]]),  neg_z_mpause))
+
+    return magnetopause
+
+
+def find_magnetopause_sw_streamline_2d(vlsvfile, streamline_seeds=None, seeds_n=200, seeds_x0=20*6371000, seeds_range=[-5*6371000, 5*6371000], streamline_length=45*6371000, end_x=-15*6371000, x_point_n=50):
+    """Finds the magnetopause position by tracing streamlines of the velocity field for 2d runs.
+
+        :param vlsvfile: directory and file name of .vlsv data file to use for VlsvReader
+        :kword streamline_seeds: optional streamline starting points in numpy array (coordinates in meters including the y-coordinate 0.0)
+        :kword seeds_n: instead of streamline_seeds provide a number of streamlines to be traced 
+        :kword seeds_x0: instead of streamline_seeds provide an x-coordinate for streamline starting points 
+        :kword seeds_range: instead of streamline_seeds provide [min, max] range to use for streamline starting point z-coordinates
+        :kword streamline_length: streamline length for tracing
+        :kword end_x: tail end x-coordinate (how far along the negative x-axis the magnetopause is calculated)
+        :kword x_point_n: integer, how many x-axis points the magnetopause will be divided in between the subsolar point and tail
+
+        :returns:   the magnetopause position as coordinate points in numpy array
+    """
+
+    streamlines = make_streamlines(vlsvfile, streamline_seeds, seeds_n, seeds_x0, seeds_range, streamline_length)
+    magnetopause = make_magnetopause(streamlines, end_x, x_point_n)
+
+    return magnetopause