From a91b9a465a7ccc9b2efc3fa7575a972cb8d1f847 Mon Sep 17 00:00:00 2001 From: jreimi Date: Wed, 20 Nov 2024 13:56:03 +0200 Subject: [PATCH 1/2] New versions of magnetopause scripts and a small fix for 3dcolormapslice plot external function --- pyPlots/plot_colormap3dslice.py | 2 + scripts/magnetopause2d.py | 603 +++++---------------------- scripts/magnetopause3d.py | 710 +++++++++++--------------------- 3 files changed, 336 insertions(+), 979 deletions(-) diff --git a/pyPlots/plot_colormap3dslice.py b/pyPlots/plot_colormap3dslice.py index 77716b5b..d30d1e29 100644 --- a/pyPlots/plot_colormap3dslice.py +++ b/pyPlots/plot_colormap3dslice.py @@ -1328,6 +1328,8 @@ def exprMA_cust(exprmaps, requestvariables=False): # Optional external additional plotting routine overlayed on color plot # Uses the same pass_maps variable as expressions if external: + if 'pass_maps' not in vars(): + pass_maps = [] #extresult=external(ax1, XmeshXY,YmeshXY, pass_maps) if not axes: extresult=external(ax1, XmeshCentres,YmeshCentres, pass_maps) diff --git a/scripts/magnetopause2d.py b/scripts/magnetopause2d.py index a204deb6..a226cd7f 100644 --- a/scripts/magnetopause2d.py +++ b/scripts/magnetopause2d.py @@ -1,575 +1,180 @@ -import pytools as pt -import plot_colormap #gives error but doesn't work without -import sys, os, socket +import matplotlib.pylab as plt import numpy as np -from operator import attrgetter - +import pytools as pt +import plot_colormap import yt -import matplotlib.pylab as plt from yt.visualization.api import Streamlines -from matplotlib.collections import LineCollection -from matplotlib.colors import LinearSegmentedColormap -import logging - - -''' -running: --takes two (or zero) commad line arguments: how many inner streamlines to ignore and from how many next to count the magnetopause position - -Some issues: --Nose part unstable and slows the whole thing down, also requires lots of streamlines to work --Bulk file change is hard and must be done by hand --Everything else must also be done by changing the code by hand -''' - - - -class Point: - ''' - a 2D point - ''' - - def __init__(self, x, z): - self.x = x - self.z = z - - #radius and angle in regard to origo - self.r, self.phi = np.array(cartesian_to_polar(x,z)) - - #radius to x-axis - self.r_x = np.sqrt(z**2) - - #radius to z-axis - self.r_z = np.sqrt(x**2) - - - -class Slice: #slice of a plane - - def __init__(self, plane, constant, pos): - self.plane = plane # 'YZ' or 'XY' - self.constant = constant # x or z -coordinate depending on the plane - self.pos = pos # above or below the plane, above = 0, below = 1 - self.points = [] #contains points in slice - - - def add_point(self, p): - self.points.append(p) - - - def get_mean_point(self, ignore, count): - #get mean of points in slice s ignoring first 'ignore' points - #and counting the next 'count' slices - - #if self.plane == 'XY': - # ignore = 1 - # count = 2 - - if len(self.points) < (ignore+count): #if there are not enough points in slice just ignore it - return 0 - - - #sort points in slice by their radius - if self.plane == 'YZ': - pnts = sorted(self.points, key=attrgetter('r_x')) - elif self.plane == 'XY': - pnts = sorted(self.points, key=attrgetter('r_z')) - - - m = [] #list to save the radius of points that count - i = 0 #index of point under consideration - for p in pnts: - if i < ignore: #too close - i = i + 1 - continue - elif i < (ignore+count): #points that count - if self.plane == 'YZ': - m.append(p.r_x) - elif self.plane == 'XY': - m.append(p.r_z) - - else: #went over - break - i = i + 1 - - #get mean of the collected radii - mean_r = np.mean(m) - - #below or above - if self.pos == 0: r = mean_r - if self.pos == 1: r = -mean_r - - #make (x,z)-coordinates of the mean point - if self.plane == 'YZ': - coord_mean = (self.constant, r) - if self.plane == 'XY': - coord_mean = (r, self.constant) - return coord_mean - - -class Plane: - - def __init__(self, plane, constant): - self.plane = plane # 'XY' or 'YZ' - self.constant = constant # z or x -coordiante depending on the plane - self.slices = [] - - if self.plane == 'YZ': - self.slices.append(Slice(self.plane, self.constant, 0)) #above - self.slices.append(Slice(self.plane, self.constant, 1)) #below - elif self.plane == 'XY': - self.slices.append(Slice(self.plane, self.constant, 0)) #just above - - - def add_point(self, p): # add a point to the right slice in plane - - if self.plane == 'YZ': - if p.phi > 0 and p.phi <= 180: - ind = 0 - else: - ind = 1 - elif self.plane == 'XY': - ind = 0 - - sl = self.slices[ind] - sl.add_point(p) - - - - -def get_bulk(): #fetch bulk file name - - #BCQ: (bulks from 0001500 to 0002500) - run = 'BCQ' - num = 1550 - num = str(num) +def interpolate(streamline, x_points): - fileLocation="/wrk/group/spacephysics/vlasiator/2D/"+run+"/bulk/" - fileN = "bulk.000"+num+".vlsv" - - f = fileLocation+fileN - return f, run - - - -def to_Re(m): #meters to Re - Re = 6371000 - return m/Re - -def cartesian_to_polar(x, z): #cartesian coordinates to polar - r = np.sqrt(x**2 + z**2) - phi = np.degrees(np.arctan2(z, x)) - - - return(r, phi) + arr = np.array(streamline) -def polar_to_cartesian(r, phi): #polar coordinates to cartesian - phi = np.radians(phi) + # set arrays for interpolation + xp = arr[:,0][::-1] + zp = arr[:,2][::-1] - z = r * np.sin(phi) - x = r * np.cos(phi) + #interpolate z coordinates + z_points = np.interp(x_points, xp, zp, left=np.NaN, right=np.NaN) - return(x, z) + return np.array([x_points, z_points]) -def make_streamlines(): - boxre=[0,0] +def make_streamlines(vlsvFileName): + ## make streamlines + boxre = [0,0] # bulk file - name = get_bulk()[0] - f = pt.vlsvfile.VlsvReader(file_name=name) + 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=[xmin,xmax,ymin,ymax,zmin,zmax] - sizes=np.array([xsize,ysize,zsize]) - - simext = list(map(to_Re, simext)) - - #set box coordnates - if len(boxre)==4: - boxcoords=boxre - else: - boxcoords=list(simext) - - - # If box extents were provided manually, truncate to simulation extents - boxcoords[0] = max(boxcoords[0],simext[0]) #xmax - boxcoords[1] = min(boxcoords[1],simext[1]) #xmin - boxcoords[2] = max(boxcoords[2],simext[2]) #ymax - boxcoords[3] = min(boxcoords[3],simext[3]) #ymin - boxcoords[4] = max(boxcoords[4],simext[4]) #zmax - boxcoords[5] = min(boxcoords[5],simext[5]) #zmin - + 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) - streamline_seeds = [] - - #range: np.arange(from, to, step) - for i in np.arange(-10.0, 10.0, 0.01): # 2000 seeds is nice, 400 is on the low end - streamline_seeds.append([20, 0, i]) - - - streamline_seeds = np.array(streamline_seeds) + #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.frontends.stream.load_uniform_grid( + yt_dataset = yt.load_uniform_grid( data, sizes, bbox=np.array([[boxcoords[0], boxcoords[1]], - [boxcoords[2],boxcoords[3]], - [boxcoords[4],boxcoords[5]]]), - periodicity=[True, True, True]) # Has to be forced... + [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=60, direction=1) + "Vx", "Vy", "Vz", length=40, direction=1) #where the magic happens streamlines_pos.integrate_through_volume() - - - return streamlines_pos #returns yt streamlines object - - -def get_poynting_flux(coordinate_points): - - #input to meters - Re = 6371000 - coordinate_points = [[coord*Re for coord in arr] for arr in coordinate_points] - - #data file - name = get_bulk()[0] - f = pt.vlsvfile.VlsvReader(file_name=name) - - #get middle points of next-to-next coordinates - middle_points = [] - normals = [] - for i,j in zip(np.arange(0, len(coordinate_points)-1), np.arange(1, len(coordinate_points))): - A = coordinate_points[i] - B = coordinate_points[j] - - middle = np.mean([A, B], axis=0) - middle_points.append(middle) - - # normal vectors of lines between next-to-next coordinate points: - n = [-(B[2]-A[2]), 0, (B[0]-A[0])] #the outwards-facing surface vector - normal_vector = n/np.linalg.norm(n) #normalization - normals.append(normal_vector) - - - - middle = [[coord/Re for coord in arr] for arr in middle_points] - - - #get interpolated poynting vectors at magnetopause coordinate middle points - S = f.read_interpolated_variable("poynting", middle_points) - - #Pounting flux (S.n) - Poynting_flux = [] - Poynting_flux_normal =[] - for s, n in zip(S, normals): - P_flux = np.dot(s,n) - Poynting_flux.append(P_flux) - Poynting_flux_normal.append(P_flux*n) - - - - - return Poynting_flux, Poynting_flux_normal, middle - + return np.array(streamlines_pos.streamlines) +def make_magnetopause(streams): -def interpolate(YTstream, points, axis): - ''' - :kword YTstream: A single streamline in yt Streamline -form - :kword points: Array of some axis points where coordinates of streamline passing a plane are wanted - :kword axis: x or z depending on what axis the points are on - ''' - - streamline = YTstream.positions.to_value() - arr = np.array(streamline) + streampoints = np.reshape(streams, (streams.shape[0]*streams.shape[1], 3)) #all the points in one array - #maybe needs to be sorted??? - xp = arr[:,0] #x-coordinates of streamline - #yp = arr[:,1] #y-coordinates - zp = arr[:,2] #z-coordinates - - #reverse arrays so interpolation works - xp = xp[::-1] - zp = zp[::-1] - - valid = True - #interpolate missing coordinate at points - if axis == 'x': - x_points = points - z_points = np.interp(points, xp, zp, left=float('inf'), right=float('inf')) - for z in z_points: - if z == float('inf'): #if z-coordinate is inf at some point scrap the whole streamline from data - valid = False - - - #gather interpolated points to new stream coordinates - stream_points = [x_points, z_points] - - return stream_points, valid - - elif axis == 'z': - x_points = [] - z_points = [] - - for z0 in points: - x_coords = [] - for i in range(0, len(xp)-1): - if (zp[i] < z0 and zp[i+1] > z0) or (zp[i] > z0 and zp[i+1] < z0): - x_coords.append(xp[i] + (z0-zp[i])*(xp[i+1]-xp[i])/(zp[i+1]-zp[i])) - - if len(x_coords) > 0: - clean_x_coords = [x for x in x_coords if x > 5] - if len(clean_x_coords) > 0: - x_points.append(min(clean_x_coords)) - else: - x_points.append(float('inf')) - else: - x_points.append(float('inf')) - - stream_points = [x_points, points] - return stream_points, True + ## 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]) - else: - logging.info('axis must be x or z!') - exit() - - -def get_magnetopause(ignore, count): + ## define points in the x axis where to find magnetopause points on the yz-plane + x_points = np.arange(subsolar_x, -10, -0.2) - streams = make_streamlines() + ## 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 - #the tail of the magnetopause + ## start making the magnetopause + ## in each x_point, find the closest streamline to x-axis in the positive and negative z-axis - #wanted points in x-axis (YZ-plane places) - x_points = np.array(np.arange(-20, 7, 0.2)) + pos_z_mpause = np.zeros((len(x_points), 2)) + neg_z_mpause = np.zeros((len(x_points), 2)) - #planes - YZplanes = [] - for x in x_points: - YZplanes.append(Plane('YZ', x)) - - #go through all streamlines one by one - for s in np.arange(0, len(streams.streamlines)): - - YTstream = streams.path(s) + for i, x_point in enumerate(x_points): + pos = new_streampoints[i, new_streampoints[i,:] > 0] + neg = new_streampoints[i, new_streampoints[i,:] < 0] - #interpolate z-values of x_points - stream_points, valid = interpolate(YTstream, x_points, 'x') + 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') - if valid: - #add interpolated points to yz-planes - for i in range(0, len(x_points)): #for every YZ-plane - - if stream_points[0][i] != x_points[i]: #check that values for x match - logging.info("Something's very wrong") - exit() - - p = Point(stream_points[0][i], stream_points[1][i]) - YZplanes[i].add_point(p) - - - #list to save the magnetopause body coordinates by slice (above/below x-axis) - pause_coords = [[],[]] - - #for every yz-plane - for p in YZplanes: + # 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()]] - #for every slice get mean point - for i in range(0, len(p.slices)): - s = p.slices[i] #slice + magnetopause = np.concatenate((pos_z_mpause[::-1], np.array([[subsolar_x, 0]]), neg_z_mpause)) - mean_point = s.get_mean_point(ignore, count) - - if mean_point != 0: - pause_coords[i].append(mean_point) - - #get the ends of the body lines for the nose - above_end = pause_coords[0][-1] - below_end = pause_coords[i][-1] - - - - - #the nose + return magnetopause - #XY-plane places - dz=0.2 - z_points = np.array(np.arange(below_end[1]+dz, above_end[1]-dz, dz)) - XYplanes = [] - for z in z_points: - XYplanes.append(Plane('XY', z)) - - - #go through all streamlines one by one again - for s in np.arange(0, len(streams.streamlines)): - YTstream = streams.path(s) - - #linear interpolation of the points in the nose - tip_points, valid = interpolate(YTstream, z_points, 'z') - - if valid: - #add interpolated points to XY-planes - for i in range(0, len(z_points)): - if tip_points[0][i] != float('inf'): - p = Point(tip_points[0][i], tip_points[1][i]) - XYplanes[i].add_point(p) - - #construct the nose - - nose = [] - - for p in XYplanes: #for every plane - for s in p.slices: #for every slice +def polar_to_cartesian(r, phi): + phi = np.deg2rad(phi) + y = r * np.cos(phi) + z = r * np.sin(phi) + return(y, z) - mean_point = s.get_mean_point(ignore, count) - if mean_point != 0: - nose.append(mean_point) - nose = np.array(nose) +def main(): - - #make a continous line of the magnetopause - #structure: above, nose, below - coords = [] - for i in 0,1: - sl = pause_coords[i] #body coordinates - - if len(sl) == 0: #if slice is empty (shouldn't be) - continue - - if i == 1: #make the nose in second loop - nose = nose[::-1] - for x, z in zip(nose[:,0], nose[:,1]): - coords.append([x, 0, z]) - - #also reverse the below-line so it starts from where the nose ended - sl = sl[::-1] - - - sl = np.array(sl) - for x, z in zip(sl[:,0], sl[:,1]): - coords.append([x, 0, z]) + ## get bulk data + run = 'BFD' + num = '2000' + fileLocation="/wrk-vakka/group/spacephysics/vlasiator/2D/"+run+"/bulk/" + fileN = "bulk.000"+num+".vlsv" - return (streams, coords) #returns streamlines, array of magnetopause positions by sloce and by point, poynting vector ends in 2d-array - + ## STREAMLINES + streams = make_streamlines(fileLocation+fileN) -def plot_all(ax, XmeshXY=None, YmeshXY=None, pass_maps=None): #external function for plot_colormap - - #what to plot, boolean - streamlines = False #streamlines - magnetopause = True #magnetopause line - poynting_flux_vectors = False # S.n * normal vectors - poynting_flux = False #poynting flux as colours - - - #command line arguments (how many inner lines to ignore and count) - if len(sys.argv) < 2: - ignore, count = 3, 5 - else: - ignore = int(sys.argv[1]) - count = int(sys.argv[2]) - - - streams, coords = get_magnetopause(ignore, count) + ## MAGNETOPAUSE + magnetopause = make_magnetopause(streams) - if streamlines: - for s in np.arange(0, len(streams.streamlines)): - stream_pos = streams.streamlines[s] - ax.plot(stream_pos[:,0], stream_pos[:,2], color='#C4C4C4', lw=0.1) - - - - if magnetopause: - coords = np.array(coords) - ax.plot(coords[:,0], coords[:,2], 'k', lw = 0.6) - - - - if poynting_flux_vectors: - poynting, vectors, middle_coords = get_poynting_flux(coords) - - coords = np.array(middle_coords) - p = np.array(vectors) - p = p*10 - - ax.quiver(coords[:,0], coords[:,2], p[:,0], p[:,2], width=0.001) + ## 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 poynting_flux: - poynting = get_poynting_flux(coords)[0] - - coords = np.array(coords) - p = np.array(poynting)*100 - - maxval = abs(max(p, key=abs)) - - x = coords[:,0] - z = coords[:,2] - - points = np.array([x, z]).T.reshape(-1,1,2) - segments = np.concatenate([points[:-1],points[1:]], axis=1) - - lc = LineCollection(segments, cmap='PiYG', linewidth=1.5, norm=plt.Normalize(vmax=maxval, vmin=-maxval)) - lc.set_array(p) - plt.gca().add_collection(lc) - + if plot_magnetopause: + ax.plot(magnetopause[:,0], magnetopause[:,1], color='cyan', linewidth=1.0) -def main(): - - inputfile, run_name = get_bulk() - plot_colormap.plot_colormap( - filename=inputfile, - #step=j, + filename=fileLocation+fileN, + outputdir=outdir, nooverwrite=None, - var="rho", - boxre = [-25,20,-20,20], + var="rho", #var="beta_star", + boxre = [-21,21,-21,21], title=None, - draw=None, usesci=True, Earth=None, - external = plot_all, - run=run_name, - #vectors='B', vectordensity=500, vectorsize=0.5 - #streamlines="B", streamlinedensity=1.5, streamlinethick=0.7, - ) - - -if __name__ == "__main__": - main() - + draw=None, usesci=True, Earth=True, + external = external_plot, + run=run, + colormap='inferno', + ) +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/scripts/magnetopause3d.py b/scripts/magnetopause3d.py index 5461a7d4..936a0253 100644 --- a/scripts/magnetopause3d.py +++ b/scripts/magnetopause3d.py @@ -1,424 +1,51 @@ - -import os -import socket -import sys -from operator import attrgetter - +from pyCalculations import ids3d import matplotlib.pylab as plt import numpy as np import pytools as pt +import plot_colormap3dslice import yt +import math from mpl_toolkits import mplot3d -from skimage import measure from yt.visualization.api import Streamlines -import ids3d -import logging - -''' -Finds the magnetopause in 3D-run and plots it with matplotlib -Accepts 2 (or 0) command line arguments: 1) how many streamlines to ignore, 2) how many streamlines to count after the ignored ones - -output: -saves png-image to the current folder -''' - - - - -arcDeg = 10 # 360/arcDeg must be an integer. Best at 10, ok at 5, 10: 36 slices, 5: 72 slices - - -class Point: #a 3D-point - - def __init__(self, x, y, z): - self.x = x - self.y = y - self.z = z - - - self.x, self.r, self.phi = np.array(cartesian_to_polar(x,y,z)) - - def print(self): - print(self.x, self.y, self.z) - - -class YZSlice: #a slice in yz-plane - - def __init__(self, x, max_deg): - self.x = x - self.max_deg = max_deg - self.min_deg = max_deg-arcDeg - self.points = [] #contains point in slice - - def get_min_deg(self): - return self.min_deg - - def get_max_deg(self): - return self.max_deg - - def add_point(self, p): - p_phi = p.phi - if p_phi == 0: p_phi = 360 - - if (p_phi <= self.max_deg) and (p_phi >= self.min_deg): - self.points.append(p) - else: - logging.info("point has wrong polar angle! (x, max_deg, p_phi):", self.x, " ", self.max_deg, " ", p_phi) - exit() - - - def get_mean_point(self, ignore, count): - #get mean of points in slice s ignoring first 'ignore' points - #and counting the average of the next 'count' points - - if len(self.points) < (ignore+count): #if there are no enough points in slice just ignore it - return 0 - - pnts = sorted(self.points, key=attrgetter('r')) #sort points by radius to the x-axis - - m = [] - i = 0 - for p in pnts: - if i < ignore: - i = i + 1 - continue - elif i < (ignore+count): - m.append(p.r) - else: - break - i = i + 1 - - - mean_r = np.mean(m) - - #make (x,y,z)-coordinates of the mean point - coord_mean = (polar_to_cartesian(self.x, mean_r, (self.max_deg-(arcDeg/2)))) - return coord_mean - -class YZPlane: #an yz-plane at a certain x - - def __init__(self, x): - self.x = x - self.slices = [] - - for i in range (0, int(360/arcDeg)): #initialize slices list - s = YZSlice(self.x, (i*arcDeg+arcDeg)) - self.slices.append(s) - - - def add_point(self, p): - - if (p.x == self.x): - ind = get_slice_index(p.phi) - sl = self.slices[ind] - sl.add_point(p) - else: - logging.info("point is in the wrong yz-plane!") - - - -def get_bulk(): #fetch bulk file name - - fileLocation="/wrk/group/spacephysics/vlasiator/3D/EGE/bulk/" - fileN = "bulk.0002193.vlsv" - - return (fileLocation+fileN) - - -def get_slice_index(deg): #aux function for indexing arcDeg-degree slices - - if deg == 0: #goes into last index slice - return int((360/arcDeg)-1) - - if deg < 0: - deg = 360 + deg - - ind = 0 - j = 0 - for i in range (arcDeg, 360+arcDeg, arcDeg): - if (deg <= i) and (deg > j): - return ind - j = i - ind = ind + 1 - - logging.info("Nope, degrees: " + str(deg) + " ind: " + str(ind)) #something's wrong - exit() def to_Re(m): #meters to Re - Re = 6371000 - return m/Re - + return m/6371000 -#Coordinate transforms -def cartesian_to_polar(x, y, z): #cartesian y,z coordinates to polar - x = x - r = np.sqrt(y**2 + z**2) - phi_0 = np.arctan2(z, y) - phi = np.degrees(phi_0) #angle in degrees +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(x, r, phi) - -def polar_to_cartesian(x, r, phi): #polar r, phi coordinates to cartesian - x = x - phi_0 = np.radians(phi) - y = r * np.cos(phi_0) - z = r * np.sin(phi_0) - return(x, y, z) - - + 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 get_poynting_vectors(coordinate_points): #TODO - ''' - input: coordinate points in form [[x0,y0,z0],[x1,y1,z1],...] - returns Poynting vectors at input points - DOESN'T WORK - ''' - - #coordinates to m - Re = 6371000 - coordinate_points = [[coord*Re for coord in arr] for arr in coordinate_points] - - - mu_0 = 1.25663706212e-06 #vacuum permeability, unit N - - name = get_bulk() - f = pt.vlsvfile.VlsvReader(file_name=name) - - - #The non-working part, no interpolation for fsgrid variables yet - B = f.read_interpolated_variable("fg_b", coordinate_points) - E = f.read_interpolated_variable("fg_e", coordinate_points) - - #Poyntnig vectors pseudocode - P = (1/mu_0) * np.cross(E,B) - - return P - - - - -def make_streamlines(): - - #boxre=[-30, 30, -30, 30, -30, 30] - boxre = [0] - - # bulk file - name = get_bulk() - f = pt.vlsvfile.VlsvReader(file_name=name) - - #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)) - - - #set box coordnates - if len(boxre)==6: - boxcoords=boxre - else: - boxcoords=list(simext) - - - # If box extents were provided manually, truncate to simulation extents - boxcoords[0] = max(boxcoords[0],simext[0]) #xmax - boxcoords[1] = min(boxcoords[1],simext[1]) #xmin - boxcoords[2] = max(boxcoords[2],simext[2]) #ymax - boxcoords[3] = min(boxcoords[3],simext[3]) #ymin - boxcoords[4] = max(boxcoords[4],simext[4]) #zmax - boxcoords[5] = min(boxcoords[5],simext[5]) #zmin - - - 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]) - #logging.info(V) - - V = V[indexids] - - Vdpoints = ids3d.idmesh3d2(cellids, V, reflevel, xsize, ysize, zsize, 3) - - - Vxs = Vdpoints[:,:,:,0] - Vys = Vdpoints[:,:,:,1] - Vzs = Vdpoints[:,:,:,2] - - data=dict(Vx=Vxs,Vy=Vys,Vz=Vzs) - - - #Create streamline seeds (starting points for streamlines) - streamline_seeds = [] - - #range: np.arange(from, to, step) - for i in np.arange(-15.0, 15.0, 0.5): - for j in np. arange(-15.0, 15.0, 0.5): - streamline_seeds.append([20, i, j]) - - - streamline_seeds = np.array(streamline_seeds) - - #dataset in yt-form - yt_dataset = yt.frontends.stream.load_uniform_grid( - data, - sizesfs, - 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=60, direction=1) - - #where the magic happens - streamlines_pos.integrate_through_volume() - - #get rid of lines that don't work - #for s in np.arange(0, len(streamlines_pos.streamlines)): - # logging.info(streamlines_pos.streamlines[s]) - # stream_pos = streamlines_pos.streamlines[s] - # stream_pos = stream_pos[(np.all(stream_pos != 0.0) | (np.linalg.norm(stream_pos, axis = 1) > 4.7)) & (stream_pos[0,2] != 0.0)] - - - return streamlines_pos #returns yt streamlines object - - -def interpolate(YTstream, x_points): - ''' - :kword YTstream: A single streamline in yt Streamline -form - :kword x_points: Array of x-axis points where coordinates of streamline passing YZ-plane are wanted - ''' - - streamline = YTstream.positions.to_value() - +def interpolate(streamline, x_points): arr = np.array(streamline) - #maybe needs to be sorted??? - xp = arr[:,0] #x-coordinates of streamline - yp = arr[:,1] #y-coordinates - zp = arr[:,2] #z-coordinates - - #reverse arrays so interpolation works - xp = xp[::-1] - yp = yp[::-1] - zp = zp[::-1] - + # 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=float('inf'), right=float('inf')) - + 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=float('inf'), right=float('inf')) - - #if z or y -coordinate is inf at some point scrap the whole streamline from data - for y, z in zip(y_points, z_points): - if z == float('inf') or y == float('inf'): - return None - - return [x_points, y_points, z_points] - - -def get_magnetopause(streams, ignore, count): - - #wanted points in x-axis (YZ-plane places) - x_points = np.array(np.arange(8, -30, -0.5)) - - #planes - planes = [] - for x in x_points: - planes.append(YZPlane(x)) - - #get coordinates of every stream at given x points (goes through all streamlines on by one) - for s in np.arange(0, len(streams.streamlines)): - YTstream = streams.path(s) - - stream_points = interpolate(YTstream, x_points) - if stream_points is None: #discard ugly streamlines (those that end too soon for example) - continue - - #add interpolated points to yz-planes - for i in range(0, len(x_points)): - if stream_points[0][i] != x_points[i]: - logging.info("Something's very wrong") - exit() - p = Point(stream_points[0][i], stream_points[1][i], stream_points[2][i]) - planes[i].add_point(p) - - - - #Now we (should) have every interpolated Point saved - - #where to save the magnetopause coordinates by slice - pause_coords_by_slice = [] - for i in range(0, int(360/arcDeg)): - pause_coords_by_slice.append([]) - - #where to save the magnetopause coordinates by x-value - pause_coords_by_x = [] - - #for every yz-plane - for p in planes: - p_ind = len(pause_coords_by_x) - pause_coords_by_x.append([]) + y_points = np.interp(x_points, xp, yp, left=np.NaN, right=np.NaN) - #for every slice get mean point - for i in range(0, len(p.slices)): - s = p.slices[i] - mean_point = s.get_mean_point(ignore, count) - if mean_point != 0: - pause_coords_by_slice[i].append(mean_point) - pause_coords_by_x[p_ind].append(mean_point) - pause_coords_by_slice = np.array(pause_coords_by_slice) + return np.array([x_points, y_points, z_points]) - return (pause_coords_by_slice, pause_coords_by_x) #returns streamlines and array of magnetopause positions -def plot_all(ax, XmeshXY=None, YmeshXY=None, pass_maps=None): #external function for plot_colormap - #command line args - if len(sys.argv) != 2: - ignore, count = 3, 3 - else: - ignore = int(sys.argv[1]) - count = int(sys.argv[2]) - - streams, magnetopause_coords = get_magnetopause(ignore, count) - - #plot streamlines - for s in np.arange(0, len(streams.streamlines)): - stream_pos = streams.streamlines[s] - ax.plot(stream_pos[:,0], stream_pos[:,2], color='#C4C4C4', lw=0.3) - - #plot magnetopause - #for every slice - for sl in magnetopause_coords: - if not sl: #if slice is empty(works on nested lists?) - continue - - sl = np.array(sl) - ax.plot3D(sl[:,0], sl[:1], sl[:,2], 'k', lw = 0.4) def make_surface(coords): @@ -426,8 +53,8 @@ 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 + 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. @@ -455,16 +82,12 @@ def make_surface(coords): 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 @@ -484,134 +107,261 @@ def make_surface(coords): 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) - #maybe could be done better? for area_index in range(planes-1): next_triangles = [x + slices_in_plane*area_index for x in first_triangles] faces.extend(next_triangles) - #faces = np.array(faces, dtype=np.int32, order='C') - #verts = np.array(verts) - #logging.info(faces.max()) - #logging.info(len(verts[:,0])) + return verts, faces + + + +def make_streamlines(vlsvFileName): + ## make streamlines + boxre = [0] + + # bulk file + f = pt.vlsvfile.VlsvReader(file_name=vlsvFileName) - return verts, faces #verts is a list of coordinates, faces an array of arrays including indices of verts that make the triangles + #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] -def triangle(A,B,C): - #A, B, C are 3D points defining the triangle - #NOTE: if input points are straight from make_surface(), every other normal will be outwards and every other inwards! + reflevel = ids3d.refinement_level(xsize, ysize, zsize, cellids[-1]) - #Centre point of the triangle (centroid) - c_x = (A[0]+B[0]+C[0])/3 - c_y = (A[1]+B[1]+C[1])/3 - c_z = (A[2]+B[2]+C[2])/3 - centre = [c_x, c_y, c_z] + #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]) - #Area of the triangle: - #vetors with A as starting point - AB = [(B[0]-A[0]), (B[1]-A[1]), (B[2]-A[2])] - AC = [(C[0]-A[0]), (C[1]-A[1]), (C[2]-A[2])] - #Cross product (non-unit length normal vector) - cross = np.cross(AB, AC) - area = np.linalg.norm(cross)/2 + V = V[indexids] - #unit-length normal vector - normal = cross / np.linalg.norm(cross) + 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]) - # vector area = area * normal vector - vector_area = area * normal - #check that area vector points inwards (TODO) + Vdpoints = ids3d.idmesh3d2(cellids, V, reflevel, xsize, ysize, zsize, shape) - #Poynting vector at the centre (not working yet) - P = get_poynting_vectors(centre) - flow = np.dot(P, vector_area) + Vxs = Vdpoints[:,:,:,0] + Vys = Vdpoints[:,:,:,1] + Vzs = Vdpoints[:,:,:,2] - return flow + 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 -def main(): - ''' - command line arguments: how many inner streamlines to ignore and how many to count in - ''' + #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]]])) + - #command line args - if len(sys.argv) != 2: - ignore, count = 3, 3 - else: - ignore = int(sys.argv[1]) - count = int(sys.argv[2]) + # data, seeds, dictionary positions, length of streamlines + streamlines_pos = yt.visualization.api.Streamlines(yt_dataset, streamline_seeds, + "Vx", "Vy", "Vz", length=50, direction=1) - #What to plot: - #[streamlines, magnetopause, poynting vectors, surface] - plotting = [False, False, False, True] + # make the streamlines + streamlines_pos.integrate_through_volume() + return np.array(streamlines_pos.streamlines) - streams = make_streamlines() - pause_by_slice, pause_by_x = get_magnetopause(streams, ignore, count) - ax = plt.axes(projection='3d', xlabel='X (Re)', ylabel='Y (Re)', zlabel='Z (Re)') +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 - if plotting[0]: #plot streamlines - for s in np.arange(0, len(streams.streamlines)): - stream_pos = streams.streamlines[s] - ax.scatter(stream_pos[:,0], stream_pos[:,1], stream_pos[:,2]) - if plotting[1]: #plot magnetopause by coordinates + ## 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]) - coords = [] - for sl in pause_by_slice: - if len(sl) == 0: continue #if slice is empty(works on nested lists?) - sl = np.array(sl) - ax.plot3D(sl[:,0], sl[:,1], sl[:,2]) - coords.append([sl[:,0], sl[:,1], sl[:,2]]) + ## 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 - for plane in pause_by_x: - plane = np.array(plane) - ax.plot3D(plane[:,0], plane[:,1], plane[:,2]) + ## 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)) - if plotting[2]: #plot Poynting vetors (not working!) - raise NotImplementedError("Poynting vector plotting not implemented!") - # coords = np.array(coords) - # poynting = np.array(poynting) + 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 - # ax.quiver(coords[:,0], coords[:,1], coords[:,2], P[:,0], P[:,1], P[:,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)] - if plotting[3]: #plot surface - logging.info('Making the surface...') - verts, faces = make_surface(pause_by_x) - verts = np.array(verts) + # 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') - ax.plot_trisurf(verts[:, 0], verts[:,1], faces, verts[:, 2], - linewidth=0.2, antialiased=True) + # 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] - plt.savefig('magnetopause3d.png') - #save 2d projection - #ax.view_init(azim=45, elev=0) - #plt.savefig('magnetopause3dprojection.png') + # 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 - logging.info('Ready!') +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): + 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): + 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 True: + 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() + main() \ No newline at end of file From 07f496a47a11441e11e90c3932845202ff2e526d Mon Sep 17 00:00:00 2001 From: jreimi Date: Thu, 28 Nov 2024 14:31:51 +0200 Subject: [PATCH 2/2] Reverted changes to plot_colormap3dslice and added short docstirngs to magnetopause scripts --- pyPlots/plot_colormap3dslice.py | 2 -- scripts/magnetopause2d.py | 7 ++++++- scripts/magnetopause3d.py | 19 ++++++++++++++----- 3 files changed, 20 insertions(+), 8 deletions(-) diff --git a/pyPlots/plot_colormap3dslice.py b/pyPlots/plot_colormap3dslice.py index d30d1e29..77716b5b 100644 --- a/pyPlots/plot_colormap3dslice.py +++ b/pyPlots/plot_colormap3dslice.py @@ -1328,8 +1328,6 @@ def exprMA_cust(exprmaps, requestvariables=False): # Optional external additional plotting routine overlayed on color plot # Uses the same pass_maps variable as expressions if external: - if 'pass_maps' not in vars(): - pass_maps = [] #extresult=external(ax1, XmeshXY,YmeshXY, pass_maps) if not axes: extresult=external(ax1, XmeshCentres,YmeshCentres, pass_maps) diff --git a/scripts/magnetopause2d.py b/scripts/magnetopause2d.py index a226cd7f..910cda71 100644 --- a/scripts/magnetopause2d.py +++ b/scripts/magnetopause2d.py @@ -6,6 +6,11 @@ from yt.visualization.api import Streamlines +''' +Finds the magnetopause position by tracing steamines of the plasma flow for two-dimensional Vlasiator runs +Needs the yt package +''' + def interpolate(streamline, x_points): arr = np.array(streamline) @@ -177,4 +182,4 @@ def external_plot(ax,XmeshXY=None, YmeshXY=None, pass_maps=None): if __name__ == "__main__": - main() \ No newline at end of file + main() diff --git a/scripts/magnetopause3d.py b/scripts/magnetopause3d.py index 936a0253..550c4461 100644 --- a/scripts/magnetopause3d.py +++ b/scripts/magnetopause3d.py @@ -8,6 +8,11 @@ from mpl_toolkits import mplot3d from yt.visualization.api import Streamlines +''' +Finds the magnetopause position by tracing steamines of the plasma flow for three-dimensional Vlasiator runs +Needs the yt package +''' + def to_Re(m): #meters to Re return m/6371000 @@ -311,8 +316,10 @@ def main(): #2D plots # analysator 3dcolormapslice y=0 if True: - def external_plot(ax,XmeshXY=None, YmeshXY=None, pass_maps=None): - ax.plot(xz_slice[:,0], xz_slice[:,2], color='limegreen', linewidth=1.5) + 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( @@ -330,7 +337,9 @@ def external_plot(ax,XmeshXY=None, YmeshXY=None, pass_maps=None): # analysator 3dcolormapslice z=0 if True: - def external_plot(ax,XmeshXY=None, YmeshXY=None, pass_maps=None): + 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( @@ -348,7 +357,7 @@ def external_plot(ax,XmeshXY=None, YmeshXY=None, pass_maps=None): # 3D plot # matplotlib 3d surface plot for single image - if True: + if False: verts, faces = make_surface(magnetopause) verts = np.array(verts) @@ -364,4 +373,4 @@ def external_plot(ax,XmeshXY=None, YmeshXY=None, pass_maps=None): if __name__ == "__main__": - main() \ No newline at end of file + main()