diff --git a/Documentation/sphinx/magnetopause2d.rst b/Documentation/sphinx/magnetopause2d.rst deleted file mode 100644 index 8e1717f1..00000000 --- a/Documentation/sphinx/magnetopause2d.rst +++ /dev/null @@ -1,5 +0,0 @@ -magnetopause2d --------------- - -.. automodule:: magnetopause2d - :members: \ No newline at end of file diff --git a/Documentation/sphinx/magnetopause3d.rst b/Documentation/sphinx/magnetopause3d.rst deleted file mode 100644 index 7d90752f..00000000 --- a/Documentation/sphinx/magnetopause3d.rst +++ /dev/null @@ -1,5 +0,0 @@ -magnetopause3d --------------- - -.. automodule:: magnetopause3d - :members: \ No newline at end of file diff --git a/Documentation/sphinx/plot_streamline_magnetopause_2d.rst b/Documentation/sphinx/plot_streamline_magnetopause_2d.rst new file mode 100644 index 00000000..38e7d4fc --- /dev/null +++ b/Documentation/sphinx/plot_streamline_magnetopause_2d.rst @@ -0,0 +1,6 @@ +plot_streamline_magnetopause_2d +------------------------------- + +.. automodule:: plot_streamline_magnetopause_2d + :members: + \ No newline at end of file diff --git a/Documentation/sphinx/plot_streamline_magnetopause_3d.rst b/Documentation/sphinx/plot_streamline_magnetopause_3d.rst new file mode 100644 index 00000000..f1794ffb --- /dev/null +++ b/Documentation/sphinx/plot_streamline_magnetopause_3d.rst @@ -0,0 +1,6 @@ +plot_streamline_magnetopause_3d +------------------------------- + +.. automodule:: plot_streamline_magnetopause_3d + :members: + \ No newline at end of file diff --git a/Documentation/sphinx/scripts.rst b/Documentation/sphinx/scripts.rst index 080c0f90..dffdadaf 100644 --- a/Documentation/sphinx/scripts.rst +++ b/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 diff --git a/analysator/pyCalculations/calculations.py b/analysator/pyCalculations/calculations.py index 7f5e940f..d9ea5bfd 100644 --- a/analysator/pyCalculations/calculations.py +++ b/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 diff --git a/analysator/pyCalculations/fieldtracer.py b/analysator/pyCalculations/fieldtracer.py index 7825423e..195c7d05 100644 --- a/analysator/pyCalculations/fieldtracer.py +++ b/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 diff --git a/analysator/pyCalculations/magnetopause_sw_streamline_2d.py b/analysator/pyCalculations/magnetopause_sw_streamline_2d.py new file mode 100644 index 00000000..4cf2a3b2 --- /dev/null +++ b/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 diff --git a/analysator/pyCalculations/magnetopause_sw_streamline_3d.py b/analysator/pyCalculations/magnetopause_sw_streamline_3d.py new file mode 100644 index 00000000..b955dadf --- /dev/null +++ b/analysator/pyCalculations/magnetopause_sw_streamline_3d.py @@ -0,0 +1,316 @@ +''' +Finds the magnetopause position by tracing streamlines of the plasma flow for three-dimensional Vlasiator runs. +''' + +import numpy as np +import analysator as pt + +def cartesian_to_polar(cartesian_coords): # for segments of plane + """Converts cartesian coordinates to polar (for the segments of the yz-planes). + + :param cartesian_coords: (y,z) coordinates as list or array + :returns: the polar coordinates r, phi (angle in degrees) + """ + y,z = cartesian_coords[0], cartesian_coords[1] + r = np.sqrt(z**2 + y**2) + phi = np.arctan2(z, y) + phi = np.rad2deg(phi) #angle in degrees + if phi < 0: phi = phi+360 + return(r, phi) + +def polar_to_cartesian(r, phi): + """Converts polar coordinates of the yz-plane to cartesian coordinates. + + :param r: radius of the segment (distance from the x-axis) + :param phi: the angle coordinate in degrees + :returns: y, z -coordinates in cartesian system + """ + phi = np.deg2rad(phi) + y = r * np.cos(phi) + z = r * np.sin(phi) + return(y, z) + + +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 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] + yp = arr[:,1][::-1] + zp = arr[:,2][::-1] + + #interpolate z from xz-plane + z_points = np.interp(x_points, xp, zp, left=np.NaN, right=np.NaN) + #interpolate y from xy-plane + y_points = np.interp(x_points, xp, yp, left=np.NaN, right=np.NaN) + + + return np.array([x_points, y_points, z_points]) + + + + +def make_surface(coords): + '''Defines a surface constructed of input coordinates as triangles. + + :param coords: points that make the surface + :returns: list of verts and vert indices of surface triangles as numpy arrays + + input coordinates must be in form [...[c11, c21, c31, ... cn1],[c12, c22, c32, ... cn2],... + where cij = [xij, yij, zij], i marks sector, j marks yz-plane (x_coord) index + + Three points make a triangle, triangles make the surface. + For every two planes next to each other: + - take every other point from plane1, every other from plane2 (in order!) + - from list of points: every three points closest to each other make a surface + + Example: + plane 1: [v1, v2, v3, v4] + plane 2: [v5, v6, v7, v8] + + -> list: [v1, v5, v2, v6, v3,...] + -> triangles: + v1 v5 v2 + v5 v2 v6 + v2 v6 v3 + . + . + . + + ''' + verts = [] #points + faces = [] #surface triangles + + slices_in_plane = len(coords[0]) + planes = len(coords) + + #get points + for plane in coords: + for vert in plane: + verts.append(vert) + + #get triangle (face) indices + #Let's consider the area between the first two planes + #ind1, ind2, ind3 for triangle indices + ind1 = 0 #starts from plane 1 + ind2 = slices_in_plane #starts from plane 2 + ind3 = 1 #starts from plane 1 + first_triangles = [] + + while len(first_triangles) < (2*slices_in_plane): + first_triangles.append([ind1, ind2, ind3]) + ind1 = ind2 + ind2 = ind3 + if (ind3 == (slices_in_plane*2)-1): #max index, go over to plane 1 first index + ind3 = 0 + elif (ind3 == 0): #last round, go to plane 2 first index + ind3 = slices_in_plane + else: + ind3 = ind1 + 1 + + first_triangles = np.array(first_triangles) + + #Now the rest of the triangles are just the first triangles + (index of area * slices_in_plane) + for area_index in range(planes-1): + next_triangles = [x + slices_in_plane*area_index for x in first_triangles] + faces.extend(next_triangles) + + # From last triangles remove every other triangle + # (a single subsolar point -> last triangles are actual triangles instead of rectangles sliced in two) + removed = 0 + for i in range(len(faces)-slices_in_plane*2, len(faces)): + if i%2!=0: + faces.pop(i-removed) + removed += 1 + + # From last triangles remove every other triangle + # (a single subsolar point -> last triangles are actual triangles instead of rectangles sliced in two) + # Also fix the last triangles so that they only point to one subsolar point and have normals towards outside + subsolar_index = int(len(verts)-slices_in_plane) + + for i,triangle in enumerate(reversed(faces)): + if i > (slices_in_plane): # faces not in last plane (we're going backwards) + break + + faces[len(faces)-i-1] = np.clip(triangle, a_min=0, a_max=subsolar_index) + + # this would remove duplicate subsolar points from vertices but makes 2d slicing harder + #verts = verts[:int(len(verts)-slices_in_plane+1)] + + # Change every other face triangle (except for last slice triangles) normal direction so that all face the same way (hopefully) + faces = np.array(faces) + for i in range(len(faces)-int(slices_in_plane)): + if i%2!=0: + faces[i,1], faces[i,2] = faces[i,2], faces[i,1] + + return np.array(verts), faces + + +def make_streamlines(vlsvfile, streamline_seeds=None, seeds_n=25, seeds_x0=20*6371000, seeds_range=[-5*6371000, 5*6371000], dl=2e6, iterations=200): + """Traces streamlines of velocity field from outside the magnetosphere to magnetotail. + + :param vlsvfile: directory and file name of .vlsv data file to use for VlsvReader + + :kword streamline_seeds: optional streamline starting points in numpy array + :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 coordinates (both y- and z-directions use the same range) + :kword dl: streamline iteration step length in m + :kword iterations: int, number of iteration steps + + :returns: streamlines as numpy array + """ + + f = pt.vlsvfile.VlsvReader(file_name=vlsvfile) + + # Create streamline starting points if needed + if streamline_seeds == None: + streamline_seeds = np.zeros([seeds_n**2, 3]) + + t = np.linspace(seeds_range[0], seeds_range[1], seeds_n) + k = 0 + for i in t: + for j in t: + streamline_seeds[k] = [seeds_x0, i, j] + k = k+1 + + # Trace the streamlines + streams = pt.calculations.static_field_tracer_3d( + vlsvReader=f, + seed_coords=streamline_seeds, + max_iterations=iterations, + dx=dl, + direction='+', + grid_var='vg_v') + + return streams + + +def make_magnetopause(streams, end_x=-15*6371000, x_point_n=50, sector_n=36): + """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 + :kword sector_n: integer, how many sectors the magnetopause will be divided in on each yz-plane + + :returns: the magnetopause position as coordinate points in numpy array, form [...[c11, c21, c31, ... cn1],[c12, c22, c32, ... cn2],... + where cij = [xij, yij, zij], i marks sector, j marks yz-plane (x-coordinate) index + """ + + RE = 6371000 + + #streams = streams*(1/RE) # streamlines in rE + streampoints = np.reshape(streams, (streams.shape[0]*streams.shape[1], 3)) #all the points in one array) + + ## find the subsolar dayside point in the x-axis + ## do this by finding a streamline point on positive x axis closest to the Earth + # streampoints closer than ~1 rE to positive x-axis: + x_axis_points = streampoints[(streampoints[:,1]-RE) & (streampoints[:,2]>-RE) & (streampoints[:,0]>0) & (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 + #dx = (subsolar_x-end_x)/x_point_n + next_from_subsolar_x = subsolar_x-1e3 # start making the magnetopause from a point slightly inwards from subsolar point + x_point_n = x_point_n-1 + x_points = np.linspace(next_from_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(streams), 2)) # new array for keeping interpolated streamlines in form new_streampoints[x_point, streamline, y and z -coordinates] + i=0 # streamline + for stream in streams: + interpolated_streamline = interpolate(stream, x_points) + + if type(interpolated_streamline) is np.ndarray: # don't use 'discarded' streamlines, see function interpolate() + for j in range(0, len(x_points)): + y,z = interpolated_streamline[1,j], interpolated_streamline[2,j] + new_streampoints[j, i,:] = np.array([y,z]) + i += 1 + + + + ## create a list of streamline points in polar coordinates (for every x_point) + polar_coords = np.zeros_like(new_streampoints) + for i in range(0,new_streampoints.shape[0]): + for j in range(0,new_streampoints.shape[1]): + polar_coords[i,j,:] = cartesian_to_polar(new_streampoints[i,j]) + + + ## now start making the magnetopause + ## in each x_point, divide the plane into sectors and look for the closest streamline to x-axis in the sector + + ## if given sector number isn't divisible by 4, make it so because we want to have magnetopause points at exactly y=0 and z=0 for 2d slices of the whole thing + while sector_n%4 != 0: + sector_n +=1 + + sector_width = 360/sector_n + magnetopause = np.zeros((len(x_points), sector_n, 3)) + + for i,x_point in enumerate(x_points): #loop over every chosen x-axis point + # divide the yz-plane into sectors + for j, mean_sector_angle in enumerate(np.arange(0, 360, sector_width)): + min_angle = mean_sector_angle-sector_width/2 + max_angle = mean_sector_angle+sector_width/2 + + # find points that are in the sector + if mean_sector_angle == 0: # special case as the first sector needs streamlines around phi=0 + min_angle = min_angle+360 + # divide into phi<360 and phi>0 + sector1 = polar_coords[i, (polar_coords[i,:,1] <= 360)*(polar_coords[i,:,1] > min_angle)] + sector2 = polar_coords[i, (polar_coords[i,:,1] <= max_angle)*(polar_coords[i,:,1] >= 0)] + sector_points = np.concatenate((sector1, sector2)) + + else: + sector_points = polar_coords[i, (polar_coords[i,:,1] <= max_angle)*(polar_coords[i,:,1] > min_angle)] + + # discard 'points' with r=0 and check that there's at least one streamline point in the sector + sector_points = sector_points[sector_points[:,0] != 0.0] + if sector_points.size == 0: + raise ValueError('No streamlines found in the sector, x_i=',i) + + # find the points closest to the x-axis + closest_point_radius = sector_points[sector_points[:,0].argmin(), 0] # smallest radius + + # return to cartesian coordinates and save as a magnetopause point at the middle of the sector + y,z = polar_to_cartesian(closest_point_radius, mean_sector_angle) + magnetopause[i,j,:] = [x_point, y, z] + + + # make a tip point for the magnetopause for prettier 3d plots + tip = np.array([subsolar_x, 0, 0]) + tips = np.tile(tip, (magnetopause.shape[1],1)) + magnetopause = np.vstack(([tips], magnetopause)) + + return magnetopause + + +def find_magnetopause_sw_streamline_3d(vlsvfile, streamline_seeds=None, seeds_n=25, seeds_x0=20*6371000, seeds_range=[-5*6371000, 5*6371000], dl=2e6, iterations=200, end_x=-15*6371000, x_point_n=50, sector_n=36): + """Finds the magnetopause position by tracing streamlines of the velocity field. + + :param vlsvfile: path to .vlsv bulk file to use for VlsvReader + :kword streamline_seeds: optional streamline starting points in numpy array + :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 coordinates (both y- and z-directions use the same range) + :kword dl: streamline iteration step length in m + :kword iterations: int, number of iteration steps + :kword end_x: tail end x-coordinate (how far along the 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 + :kword sector_n: integer, how many sectors the magnetopause will be divided in on each yz-plane + + :returns: vertices, faces of the magnetopause triangle mesh as numpy arrays + """ + + streams = make_streamlines(vlsvfile, streamline_seeds, seeds_n, seeds_x0, seeds_range, dl, iterations) + magnetopause = make_magnetopause(streams, end_x, x_point_n, sector_n) + vertices, faces = make_surface(magnetopause) + + return vertices, faces diff --git a/pyproject.toml b/pyproject.toml index ceeb7f0d..a1c26867 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,6 +26,7 @@ dependencies = [ "matplotlib", "packaging", "scikit-image", + "yt" ] [project.optional-dependencies] none = [ diff --git a/scripts/magnetopause2d.py b/scripts/magnetopause2d.py deleted file mode 100644 index 5bf4e5df..00000000 --- a/scripts/magnetopause2d.py +++ /dev/null @@ -1,185 +0,0 @@ -''' -Finds the magnetopause position by tracing steamines of the plasma flow for two-dimensional Vlasiator runs. Needs the yt package. -''' - -import numpy as np -import analysator as pt -import plot_colormap -import yt -from yt.visualization.api import Streamlines - - - - -def interpolate(streamline, x_points): - - 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(vlsvFileName): - ## make streamlines - boxre = [0,0] - - # bulk file - f = pt.vlsvfile.VlsvReader(file_name=vlsvFileName) - - #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_m =[xmin,xmax,ymin,ymax,zmin,zmax] - sizes = np.array([xsize,ysize,zsize]) - - def to_Re(m): #meters to Re - return m/6371000 - - simext = list(map(to_Re, simext_m)) - 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 streamline seeds (starting points for streamlines) - seedN = 200 # number of streamlines wanted - streamline_seeds = np.array([[20, 0 ,i] for i in np.linspace(-4, 4, seedN)]) - - #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, lenght of lines - streamlines_pos = yt.visualization.api.Streamlines(yt_dataset, streamline_seeds, - "Vx", "Vy", "Vz", length=40, direction=1) - - #where the magic happens - streamlines_pos.integrate_through_volume() - return np.array(streamlines_pos.streamlines) - -def make_magnetopause(streams): - - streampoints = np.reshape(streams, (streams.shape[0]*streams.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[np.floor(streampoints[:,2])==0] - x_axis_points[x_axis_points<0] = 800 - 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.arange(subsolar_x, -10, -0.2) - - ## interpolate more exact points for streamlines at exery x_point - new_streampoints = np.zeros((len(x_points), len(streams), 1)) # new array for keeping interpolated streamlines in form streamlines_new[x_point, streamline, z-coordinate] - i=0 - for stream in streams: - interpolated_streamline = interpolate(stream, x_points) - for j in range(0, len(x_points)): - new_streampoints[j, i,:] = interpolated_streamline[1,j] - i += 1 - - - ## 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 polar_to_cartesian(r, phi): - phi = np.deg2rad(phi) - y = r * np.cos(phi) - z = r * np.sin(phi) - return(y, z) - - - -def main(): - - ## get bulk data - run = 'BFD' - num = '2000' - - fileLocation="/wrk-vakka/group/spacephysics/vlasiator/2D/"+run+"/bulk/" - fileN = "bulk.000"+num+".vlsv" - - - ## STREAMLINES - streams = make_streamlines(fileLocation+fileN) - - ## MAGNETOPAUSE - magnetopause = make_magnetopause(streams) - - - ## PLOTTING - outdir="" - - # plot the magnetopause (and streamlines if needed) on top of colormap - def external_plot(ax,XmeshXY=None, YmeshXY=None, pass_maps=None): - plot_magnetopause = True - plot_streamlines = False - - if plot_streamlines: - for stream in streams: - ax.plot(stream[:,0], stream[:,2], color='paleturquoise', alpha=0.2) - - if plot_magnetopause: - ax.plot(magnetopause[:,0], magnetopause[:,1], color='cyan', linewidth=1.0) - - - plot_colormap.plot_colormap( - filename=fileLocation+fileN, - outputdir=outdir, - nooverwrite=None, - var="rho", #var="beta_star", - boxre = [-21,21,-21,21], - title=None, - draw=None, usesci=True, Earth=True, - external = external_plot, - run=run, - colormap='inferno', - ) - - - -if __name__ == "__main__": - main() diff --git a/scripts/magnetopause3d.py b/scripts/magnetopause3d.py deleted file mode 100644 index 906880d7..00000000 --- a/scripts/magnetopause3d.py +++ /dev/null @@ -1,376 +0,0 @@ -''' -Finds the magnetopause position by tracing steamines of the plasma flow for three-dimensional Vlasiator runs. Needs the yt package. -''' -from pyCalculations import ids3d -import matplotlib.pyplot as plt -import numpy as np -import analysator as pt -import plot_colormap3dslice -import yt -import math -from mpl_toolkits import mplot3d -from yt.visualization.api import Streamlines - - - - -def to_Re(m): #meters to Re - return m/6371000 - - -def cartesian_to_polar(cartesian_coords): # for segments of plane - y,z = cartesian_coords[0], cartesian_coords[1] - r = np.sqrt(z**2 + y**2) - phi = np.arctan2(z, y) - phi = np.rad2deg(phi) #angle in degrees - if phi < 0: phi = phi+360 - return(r, phi) - -def polar_to_cartesian(r, phi): - phi = np.deg2rad(phi) - y = r * np.cos(phi) - z = r * np.sin(phi) - return(y, z) - - -def interpolate(streamline, x_points): - arr = np.array(streamline) - - # set arrays for interpolation - xp = arr[:,0][::-1] - yp = arr[:,1][::-1] - zp = arr[:,2][::-1] - - #interpolate z from xz-plane - z_points = np.interp(x_points, xp, zp, left=np.NaN, right=np.NaN) - #interpolate y from xy-plane - y_points = np.interp(x_points, xp, yp, left=np.NaN, right=np.NaN) - - - return np.array([x_points, y_points, z_points]) - - - - -def make_surface(coords): - - ''' - Defines surface constructed of input coordinates as triangles - Returns list of verts and vert indices of surface triangles - - coordinates must be in form [...[c11, c21, c31, ... cn1],[c12, c22, c32, ... cn2],... - where cij = [xij, yij, zij], i marks slice, j marks yz-plane (x_coord) index - - How it works: - Three points make a triangle, triangles make the surface. - For every two planes next to each other: - - take every other point from plane1, every other from plane2 (in order!) - - from list of points: every three points closest to each other make a surface - - Example: - plane 1: [v1, v2, v3, v4] - plane 2: [v5, v6, v7, v8] - - -> list: [v1, v5, v2, v6, v3,...] - -> triangles: - v1 v5 v2 - v5 v2 v6 - v2 v6 v3 - . - . - . - - ''' - verts = [] #points - faces = [] #surface triangles - - slices_in_plane = len(coords[0]) - planes = len(coords) - - #get points - for plane in coords: - for vert in plane: - verts.append(vert) - - #get triangle (face) indices - #Let's consider the area between the first two planes - #ind1, ind2, ind3 for triangle indices - ind1 = 0 #starts from plane 1 - ind2 = slices_in_plane #starts from plane 2 - ind3 = 1 #starts from plane 1 - first_triangles = [] - - while len(first_triangles) < (2*slices_in_plane): - first_triangles.append([ind1, ind2, ind3]) - ind1 = ind2 - ind2 = ind3 - if (ind3 == (slices_in_plane*2)-1): #max index, go over to plane 1 first index - ind3 = 0 - elif (ind3 == 0): #last round, go to plane 2 first index - ind3 = slices_in_plane - else: - ind3 = ind1 + 1 - - first_triangles = np.array(first_triangles) - - #Now the rest of the triangles are just the first triangles + (index of area * slices_in_plane) - for area_index in range(planes-1): - next_triangles = [x + slices_in_plane*area_index for x in first_triangles] - faces.extend(next_triangles) - - return verts, faces - - - -def make_streamlines(vlsvFileName): - ## make streamlines - boxre = [0] - - # bulk file - f = pt.vlsvfile.VlsvReader(file_name=vlsvFileName) - - #get box coordinates from data - [xmin, ymin, zmin, xmax, ymax, zmax] = f.get_spatial_mesh_extent() - [xsize, ysize, zsize] = f.get_spatial_mesh_size() - [xsizefs, ysizefs, zsizefs] = f.get_fsgrid_mesh_size() - simext_m =[xmin,xmax,ymin,ymax,zmin,zmax] - sizes = np.array([xsize,ysize,zsize]) - sizesfs = np.array([xsizefs,ysizefs,zsizefs]) - - simext = list(map(to_Re, simext_m)) - boxcoords=list(simext) - - cellids = f.read_variable("CellID") - indexids = cellids.argsort() - cellids = cellids[indexids] - - reflevel = ids3d.refinement_level(xsize, ysize, zsize, cellids[-1]) - - - #Read the data from vlsv-file - V = f.read_variable("vg_v") - - #from m to Re - V = np.array([[to_Re(i) for i in j ]for j in V]) - - - V = V[indexids] - - if np.ndim(V)==1: - shape = None - elif np.ndim(V)==2: # vector variable - shape = V.shape[1] - elif np.ndim(V)==3: # tensor variable - shape = (V.shape[1], V.shape[2]) - - - Vdpoints = ids3d.idmesh3d2(cellids, V, reflevel, xsize, ysize, zsize, shape) - - - Vxs = Vdpoints[:,:,:,0] - Vys = Vdpoints[:,:,:,1] - Vzs = Vdpoints[:,:,:,2] - - data=dict(Vx=Vxs,Vy=Vys,Vz=Vzs) - - #Create streamline seeds (starting points for streamlines) - seedN = 50 #seeds per row, final seed count will be seedN*seedN ! - streamline_seeds = np.zeros([seedN**2, 3]) - - #range: np.arange(from, to, step) - t = np.linspace(-5, 5, seedN) - k = 0 - for i in t: - for j in t: - streamline_seeds[k] = [20, i, j] - k = k+1 - - - #dataset in yt-form - yt_dataset = yt.load_uniform_grid( - data, - sizesfs, - bbox=np.array([[boxcoords[0], boxcoords[1]], - [boxcoords[2],boxcoords[3]], - [boxcoords[4],boxcoords[5]]])) - - - # data, seeds, dictionary positions, length of streamlines - streamlines_pos = yt.visualization.api.Streamlines(yt_dataset, streamline_seeds, - "Vx", "Vy", "Vz", length=50, direction=1) - - # make the streamlines - streamlines_pos.integrate_through_volume() - - return np.array(streamlines_pos.streamlines) - - -def make_magnetopause(streams): - streampoints = np.reshape(streams, (streams.shape[0]*streams.shape[1], 3)) #all the points in one array - - ## find the subsolar dayside point in the x-axis - ## do this by finding a stremline point on positive x axis closest to the Earth - x_axis_points = streampoints[(np.floor(streampoints[:,1])==0) & (np.floor(streampoints[:,2])==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.arange(subsolar_x, -15, -0.5) - - ## interpolate more exact points for streamlines at exery x_point - new_streampoints = np.zeros((len(x_points), len(streams), 2)) # new array for keeping interpolated streamlines in form new_streampoints[x_point, streamline, y and z -coordinates] - i=0 # streamline - for stream in streams: - interpolated_streamline = interpolate(stream, x_points) - - if type(interpolated_streamline) is np.ndarray: # don't use 'discarded' streamlines, see function interpolate() - for j in range(0, len(x_points)): - y,z = interpolated_streamline[1,j], interpolated_streamline[2,j] - new_streampoints[j, i,:] = np.array([y,z]) - i += 1 - - - - ## create a list of streamline points in polar coordinates (for every x_point) - polar_coords = np.zeros_like(new_streampoints) - for i in range(0,new_streampoints.shape[0]): - for j in range(0,new_streampoints.shape[1]): - polar_coords[i,j,:] = cartesian_to_polar(new_streampoints[i,j]) - - - ## now start making the magnetopause - ## in each x_point, divide the plane into sectors and look for the closest streamline to x-axis in the sector - sector_n = 36 - - ## if given sector number isn't divisible by 4, make it so because we want to have magnetopause points at exactly y=0 and z=0 for 2d slices of the whole thing - while sector_n%4 != 0: - sector_n +=1 - - sector_width = 360/sector_n - magnetopause = np.zeros((len(x_points), sector_n, 3)) - - for i,x_point in enumerate(x_points): #loop over every chosen x-axis point - # divide the yz-plane into sectors - for j, mean_sector_angle in enumerate(np.arange(0, 360, sector_width)): - min_angle = mean_sector_angle-sector_width/2 - max_angle = mean_sector_angle+sector_width/2 - - # find points that are in the sector - if mean_sector_angle == 0: # special case as the first sector needs streamlines around phi=0 - min_angle = min_angle+360 - # divide into phi<360 and phi>0 - sector1 = polar_coords[i, (polar_coords[i,:,1] <= 360)*(polar_coords[i,:,1] > min_angle)] - sector2 = polar_coords[i, (polar_coords[i,:,1] <= max_angle)*(polar_coords[i,:,1] >= 0)] - sector_points = np.concatenate((sector1, sector2)) - - else: - sector_points = polar_coords[i, (polar_coords[i,:,1] <= max_angle)*(polar_coords[i,:,1] > min_angle)] - - # discard 'points' with r=0 and check that there's at least one streamline point in the sector - sector_points = sector_points[sector_points[:,0] != 0.0] - if sector_points.size == 0: - raise ValueError('No streamlines found in the sector') - - # find the points closest to the x-axis - closest_point_radius = sector_points[sector_points[:,0].argmin(), 0] # smallest radius - - # return to cartesian coordinates and save as a magnetopause point at the middle of the sector - y,z = polar_to_cartesian(closest_point_radius, mean_sector_angle) - magnetopause[i,j,:] = [x_point, y, z] - - - # make a tip point for the magnetopause for prettier 3d plots - tip = np.array([subsolar_x, 0, 0]) - tips = np.tile(tip, (magnetopause.shape[1],1)) - magnetopause = np.vstack(([tips], magnetopause)) - - return magnetopause - - - - -def main(): - - ## get bulk data - run = 'EGI' - fileLocation="/wrk-vakka/group/spacephysics/vlasiator/3D/"+run+"/bulk/" - fileN = "bulk5.0000070.vlsv" - - ## STREAMLINES - streams = make_streamlines(fileLocation+fileN) - ## MAGNETOPAUSE - magnetopause = make_magnetopause(streams) - - - ## PLOTTING - outdir="" - - ## take separate arrays for different 2d slice plots - slices = magnetopause.shape[1] - quarter_slice = int(slices/4) - # xy plane: z=0 - xy_slice = np.concatenate((magnetopause[:,0][::-1], magnetopause[:,2*quarter_slice])) - # xz plane: y=0 - xz_slice = np.concatenate((magnetopause[:,quarter_slice][::-1], magnetopause[:,3*quarter_slice])) - - - #2D plots - # analysator 3dcolormapslice y=0 - if True: - def external_plot(ax,XmeshXY=None, YmeshXY=None, pass_maps=None, requestvariables=False): - if requestvariables==True: - return ['vg_v'] - ax.plot(xz_slice[:,0], xz_slice[:,2], color='limegreen', linewidth=1.5) - - - plot_colormap3dslice.plot_colormap3dslice( - filename=fileLocation+fileN, - outputdir=outdir, - run=run, - nooverwrite=None, - boxre = [-21,21,-21,21], - title=None, - draw=None, usesci=True, Earth=True, - external = external_plot, - colormap='inferno', - normal='y' - ) - - # analysator 3dcolormapslice z=0 - if True: - def external_plot(ax,XmeshXY=None, YmeshXY=None, pass_maps=None, requestvariables=False): - if requestvariables==True: - return ['vg_v'] - ax.plot(xy_slice[:,0], xy_slice[:,1], color='limegreen', linewidth=1.5) - - plot_colormap3dslice.plot_colormap3dslice( - filename=fileLocation+fileN, - outputdir=outdir, - run=run, - nooverwrite=None, - boxre = [-21,21,-21,21], - title=None, - draw=None, usesci=True, Earth=True, - external = external_plot, - colormap='inferno', - normal='z' - ) - - # 3D plot - # matplotlib 3d surface plot for single image - if False: - verts, faces = make_surface(magnetopause) - verts = np.array(verts) - - fig = plt.figure() - ax2 = fig.add_subplot(projection='3d') - ax2.plot_trisurf(verts[:, 0], verts[:,1], faces, verts[:, 2], linewidth=0.2, antialiased=True) - ax2.view_init(azim=-60, elev=5) - ax2.set_xlabel('X') - ax2.set_ylabel('Y') - ax2.set_zlabel('Z') - fig.tight_layout() - plt.savefig(outdir+run+'_3d_magnetopause.png') - - -if __name__ == "__main__": - main() diff --git a/scripts/plot_streamline_magnetopause_2d.py b/scripts/plot_streamline_magnetopause_2d.py new file mode 100644 index 00000000..4c17e8b0 --- /dev/null +++ b/scripts/plot_streamline_magnetopause_2d.py @@ -0,0 +1,41 @@ +""" +An example of using find_magnetopause_sw_streamline_2d() and plotting the output over colormap. +""" + +import analysator as pt + +def main(): + + run = 'BFD' # for output file name + file_name="/wrk-vakka/group/spacephysics/vlasiator/2D/BFD/bulk/bulk.0002000.vlsv" + + # magnetopause points + magnetopause = pt.calculations.find_magnetopause_sw_streamline_2d( + file_name, + 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) + + magnetopause = magnetopause*(1/6371000) # in RE + + def external_plot(ax,XmeshXY=None, YmeshXY=None, pass_maps=None): + ax.plot(magnetopause[:,0], magnetopause[:,1], color='cyan', linewidth=1.0) + + pt.plot.plot_colormap( + filename=file_name, + var="rho", + boxre = [-21,21,-21,21], + Earth=True, + external = external_plot, + run=run, + colormap='inferno', + ) + + +if __name__ == "__main__": + main() + diff --git a/scripts/plot_streamline_magnetopause_3d.py b/scripts/plot_streamline_magnetopause_3d.py new file mode 100644 index 00000000..8c661b12 --- /dev/null +++ b/scripts/plot_streamline_magnetopause_3d.py @@ -0,0 +1,91 @@ +""" +An example of using find_magnetopause_sw_streamline_3d() and plotting the output. +""" + +import numpy as np +import matplotlib.pyplot as plt +import analysator as pt + +def main(): + ## get bulk data + run = 'EGI' # for plot output file name + file_name="/wrk-vakka/group/spacephysics/vlasiator/3D/EGI/bulk/bulk5.0000070.vlsv" + + ## magnetopause + x_point_n = 50 # number needed for 2d slice plots, default is 50 + + vertices, faces = pt.calculations.find_magnetopause_sw_streamline_3d( + file_name, + streamline_seeds=None, + seeds_n=25, + seeds_x0=20*6371000, + seeds_range=[-5*6371000, 5*6371000], + dl=2e6, + iterations=200, + end_x=-15*6371000, + x_point_n=x_point_n, + sector_n=36) + + + outdir="" # output plot save directory + + ### 3D surface plot ### + + fig = plt.figure() + ax = fig.add_subplot(projection='3d') + ax.plot_trisurf(vertices[:, 0], vertices[:,1], faces, vertices[:,2], linewidth=0.2) + plt.savefig(outdir+run+'_magnetopause_3D_surface.png') + plt.close() + + + ### 2D slice plots ### + + vertices = vertices/6371000 # to RE + magnetopause = np.array(np.split(vertices, x_point_n)) # magnetopause points grouped by x-axis coordinate + + ## take separate arrays for different 2d slice plots + + quarter_slice = int(np.shape(magnetopause)[1]/4) + # xy plane: z=0 + xy_slice = np.concatenate((magnetopause[:,0][::-1], magnetopause[:,2*quarter_slice])) + # xz plane: y=0 + xz_slice = np.concatenate((magnetopause[:,quarter_slice][::-1], magnetopause[:,3*quarter_slice])) + + + # analysator 3dcolormapslice y=0 + + def external_plot(ax,XmeshXY=None, YmeshXY=None, pass_maps=None, requestvariables=False): + if requestvariables==True: + return ['vg_v'] + ax.plot(xz_slice[:,0], xz_slice[:,2], color='limegreen', linewidth=1.5) + + pt.plot.plot_colormap3dslice( + filename=file_name, + run=run, + outputdir=outdir, + boxre = [-21,21,-21,21], + external = external_plot, + colormap='inferno', + normal='y' + ) + + # analysator 3dcolormapslice z=0 + + def external_plot(ax,XmeshXY=None, YmeshXY=None, pass_maps=None, requestvariables=False): + if requestvariables==True: + return ['vg_v'] + ax.plot(xy_slice[:,0], xy_slice[:,1], color='limegreen', linewidth=1.5) + + pt.plot.plot_colormap3dslice( + filename=file_name, + outputdir=outdir, + run=run, + boxre = [-21,21,-21,21], + external = external_plot, + colormap='inferno', + normal='z') + + + +if __name__ == "__main__": + main()