From 0517a48214abc9babf45af98d12d0b94c42fecc2 Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Thu, 4 Sep 2025 16:17:22 +0300 Subject: [PATCH 1/7] prototype to ditch Delaunay alpha-hulls when we actually have usable connectivity information --- scripts/magnetopause.py | 87 +++++++++++++++++++++++++++++++++++------ scripts/regions.py | 30 ++++++++++++-- 2 files changed, 100 insertions(+), 17 deletions(-) diff --git a/scripts/magnetopause.py b/scripts/magnetopause.py index 9f8dd66ea..1ea51a7c9 100644 --- a/scripts/magnetopause.py +++ b/scripts/magnetopause.py @@ -67,7 +67,7 @@ def write_SDF_to_file(SDF, datafilen, outfilen): -def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds=None, return_surface=True, return_SDF=True, SDF_points=None, Delaunay_alpha=None, beta_star_range=[0.4, 0.5], method_args={}): # TODO: separate streamline suface and vtkDelaunay3d surface in streamline method +def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds=None, return_surface=True, return_SDF=True, SDF_points=None, Delaunay_alpha=None, beta_star_range=[0.0, 0.5], method_args={}): # TODO: separate streamline suface and vtkDelaunay3d surface in streamline method """Finds the magnetopause using the specified method. Surface is constructed using vtk's Delaunay3d triangulation which results in a convex hull if no Delaunay_alpha is given. Returns vtk.vtkDataSetSurfaceFilter object and/or signed distances (negative -> inside magnetopause) (=SDF) to all cells Note that using alpha for Delaunay might make SDF different from expected inside the magnetosphere, especially if surface is constructed with points not everywhere in the magnetosphere (e.g. beta* 0.4-0.5) or if simulation grid size is larger than alpha @@ -119,19 +119,80 @@ def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds= elif method == "beta_star_with_connectivity": # magnetopause from beta_star, with connectivity if possible betastar_region = regions.treshold_mask(f.read_variable("vg_beta_star"), beta_star_range) - try: - connectivity_region = regions.treshold_mask(f.read_variable("vg_connection"), 0) # closed-closed magnetic field lines - magnetosphere_proper = np.where((connectivity_region==1) | (betastar_region==1), 1, 0) - contour_coords = f.get_cell_coordinates(cellids[magnetosphere_proper==1]) - np.save("pointcloud.npy", contour_coords) - except: - logging.warning("using field line connectivity for magnetosphere did not work, using only beta*") - #condition_dict = {"beta_star": [0.5, 0.6]} # FIC: [0.4, 0.5]) # EGE: [0.9, 1.0]) # max 0.6 in FHA to not take flyaways from outside magnetopause - mpause_flags = np.where(betastar_region==1, 1, 0) - contour_coords = f.get_cell_coordinates(cellids[mpause_flags!=0]) + if True: + vtkreader = pt.vlsvfile.VlsvVtkReader() + vtkreader.SetFileName(f.file_name) + vtkreader.Update() + # vars =vtkreader.findVariablesFromVlsv(getReducers=False) + # add here more/other datareducer outputs for downstream use if needed + vars = ["proton/vg_beta_star", "vg_connection","cellid"] + + for var in [v for v in vars if ("vg_" in v.lower()) or (v.lower() == "cellid")]: + vtkreader.addArrayFromVlsv(var) + + vtkreader.Modified() + vtkreader.Update() + dataport = vtkreader.GetOutputPort() + + dual = vtk.vtkHyperTreeGridToDualGrid() + dual.SetInputConnection(dataport) + vtkreader.Update() + dual.Update() + + data = dual.GetOutputDataObject(0) + renamer = vtk.vtkArrayRename() + renamer.SetPointArrayName("proton/vg_beta_star", "vg_beta_star") + renamer.SetInputData(data) + renamer.Update() + + data = renamer.GetOutputDataObject(0) + # print(data) + + newarr = vtk.vtkDoubleArray() + newarr.SetName("BL") + newarr.SetNumberOfComponents(1) # This will be a scalar results + newarr.SetNumberOfTuples(data.GetPointData().GetNumberOfTuples()) # Set how many values to allocate + data.GetPointData().AddArray(newarr) # add the array to the dataset + + + # Setting up the actual calulcation + calcBL = vtk.vtkArrayCalculator() + calcBL.SetInputData(data) # This data we ingest + # We need to specify which arrays we will use to derive the value + calcBL.AddScalarArrayName('vg_beta_star') + calcBL.AddScalarArrayName('vg_connection') + calcBL.SetReplaceInvalidValues(True) + calcBL.SetReplacementValue(np.nan) + # The actual function to calculate + calcBL.SetFunction("vg_beta_star*min(vg_connection,1)") + # Set result array name + calcBL.SetResultArrayName("BL") + calcBL.Update() + data = calcBL.GetOutputDataObject(0) + + threshold0 = vtk.vtkThreshold() + threshold0.SetInputArrayToProcess(0,0,0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "BL") + threshold0.SetInputData(data) + threshold0.SetLowerThreshold(beta_star_range[0]) + threshold0.SetUpperThreshold(beta_star_range[1]) + threshold0.Update() + + vtkSurface, SDF = regions.vtkSDF(query_points, threshold0.GetOutputDataObject(0)) + + else: + try: + connectivity_region = regions.treshold_mask(f.read_variable("vg_connection"), 0) # closed-closed magnetic field lines + magnetosphere_proper = np.where((connectivity_region==1) | (betastar_region==1), 1, 0) + contour_coords = f.get_cell_coordinates(cellids[magnetosphere_proper==1]) + np.save("pointcloud.npy", contour_coords) + except: + logging.warning("using field line connectivity for magnetosphere did not work, using only beta*") + #condition_dict = {"beta_star": [0.5, 0.6]} # FIC: [0.4, 0.5]) # EGE: [0.9, 1.0]) # max 0.6 in FHA to not take flyaways from outside magnetopause + mpause_flags = np.where(betastar_region==1, 1, 0) + contour_coords = f.get_cell_coordinates(cellids[mpause_flags!=0]) - # make a convex hull surface with vtk's Delaunay - vtkSurface, SDF = regions.vtkDelaunay3d_SDF(query_points, contour_coords, Delaunay_alpha) + # make a convex hull surface with vtk's Delaunay + vtkSurface, SDF = regions.vtkDelaunay3d_SDF(query_points, contour_coords, Delaunay_alpha) #elif method == "beta_star_with_fieldlines": # either incredibly slow or does not work, don't use without fixing #TODO proprer measure of actual field line backwall point # # magnetopause from beta_star, with field lines connecting to back wall if possible diff --git a/scripts/regions.py b/scripts/regions.py index 7b7895f77..807da4569 100644 --- a/scripts/regions.py +++ b/scripts/regions.py @@ -67,6 +67,27 @@ def vtkDelaunay3d_SDF(query_points, coordinates, alpha=None): return surface, convexhull_sdf +def vtkSDF(query_points, dualgrid): + ''' Obtain the SDF in relation to a waterproof volumetric grid "dualgrid" + ''' + + # print(dualgrid) + surface = vtk.vtkDataSetSurfaceFilter() + surface.SetInputData(dualgrid) + surface.Update() + + # print(surface.GetOutput()) + + + implicitPolyDataDistance = vtk.vtkImplicitPolyDataDistance() + implicitPolyDataDistance.SetInput(surface.GetOutput()) + import sys + convexhull_sdf = np.zeros(len(query_points)) + for i,coord in enumerate(query_points): + convexhull_sdf[i] = implicitPolyDataDistance.EvaluateFunction(coord) + # sys.exit() + + return surface.GetOutputDataObject(0), convexhull_sdf def treshold_mask(data_array, value): @@ -347,7 +368,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo if magnetopause_kwargs: __, magnetopause_SDF = magnetopause.magnetopause(datafile, **magnetopause_kwargs) else: - __, magnetopause_SDF = magnetopause.magnetopause(datafile, method="beta_star_with_connectivity", Delaunay_alpha=None) # default magnetopause: beta*+ B connectivity convex hull + __, magnetopause_SDF = magnetopause.magnetopause(datafile, method="beta_star_with_connectivity", Delaunay_alpha=2e6, ) # default magnetopause: beta*+ B connectivity convex hull write_flags(writer, magnetopause_SDF, 'SDF_magnetopause') write_flags(writer, np.where(np.abs(magnetopause_SDF) < 5e6, 1, 0), "flag_magnetopause") @@ -495,10 +516,11 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo def main(): - datafile = "/wrk-vakka/group/spacephysics/vlasiator/3D/EGE/bulk/bulk.0002000.vlsv" - outfilen = "EGE_regions_t2000.vlsv" + fileid = 1000 + datafile = "/wrk-vakka/group/spacephysics/vlasiator/3D/FID/bulk1/bulk1.{:07d}.vlsv".format(fileid) + outfilen = "/wrk-vakka/group/spacephysics/vlasiator/3D/FID/postprocessing/prototyping/FID_mpause_{:07d}.vlsv".format(fileid) - RegionFlags(datafile, outfilen, regions=["all"]) + RegionFlags(datafile, outfilen, regions=["magnetopause"]) if __name__ == "__main__": From e844c9ef6c18f9341153bcb0ea6dd15585f24f60 Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Thu, 4 Sep 2025 17:12:29 +0300 Subject: [PATCH 2/7] vg_-prefix to outputs --- scripts/regions.py | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/scripts/regions.py b/scripts/regions.py index 807da4569..d6da1fd4e 100644 --- a/scripts/regions.py +++ b/scripts/regions.py @@ -369,8 +369,8 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo __, magnetopause_SDF = magnetopause.magnetopause(datafile, **magnetopause_kwargs) else: __, magnetopause_SDF = magnetopause.magnetopause(datafile, method="beta_star_with_connectivity", Delaunay_alpha=2e6, ) # default magnetopause: beta*+ B connectivity convex hull - write_flags(writer, magnetopause_SDF, 'SDF_magnetopause') - write_flags(writer, np.where(np.abs(magnetopause_SDF) < 5e6, 1, 0), "flag_magnetopause") + write_flags(writer, magnetopause_SDF, 'vg_SDF_magnetopause') + write_flags(writer, np.where(np.abs(magnetopause_SDF) < 5e6, 1, 0), "vg_flag_magnetopause") # save some magnetopause values for later magnetopause_density = np.mean(variables["density"][np.abs(magnetopause_SDF) < 5e6]) @@ -382,7 +382,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo # magnetosphere from magnetopause SDF if "magnetosphere" in regions: magnetosphere = np.where(magnetopause_SDF<0, 1, 0) - write_flags(writer, magnetosphere, 'flag_magnetosphere') + write_flags(writer, magnetosphere, 'vg_flag_magnetosphere') ## BOW SHOCK ## #TODO: similar kwargs system as magnetopause? @@ -392,8 +392,8 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo bowshock = bowshock_SDF(f, variables, all_points, own_condition_dict=region_conditions["bowshock"]) else: bowshock = bowshock_SDF(f, variables, all_points) # default upstream rho method, might fail with foreshock - write_flags(writer, bowshock, 'SDF_bowshock') - write_flags(writer, np.where(np.abs(bowshock) < 5e6, 1, 0), "flag_bowshock") + write_flags(writer, bowshock, 'vg_SDF_bowshock') + write_flags(writer, np.where(np.abs(bowshock) < 5e6, 1, 0), "vg_flag_bowshock") # magnetosphere+magnetosheath -area inside_bowshock = np.where(bowshock<0, 1, 0) @@ -402,7 +402,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo if "magnetosheath" in regions: # magnetosheath from bow shock-magnetosphere difference magnetosheath_flags = np.where((inside_bowshock & 1-magnetosphere), 1, 0) - write_flags(writer, magnetosheath_flags, 'flag_magnetosheath') + write_flags(writer, magnetosheath_flags, 'vg_flag_magnetosheath') # save magnetosheath density and temperature for further use #magnetosheath_density = np.mean(variables["density"][magnetosheath_flags == 1]) @@ -413,7 +413,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo ## UPSTREAM ## if "upstream" in regions: # upstream from !bowshock - write_flags(writer, 1-inside_bowshock, 'flag_upstream') + write_flags(writer, 1-inside_bowshock, 'vg_flag_upstream') #write_flags(writer, inside_bowshock, 'flag_inside_bowshock') @@ -444,7 +444,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo } cusp_flags = make_region_flags(variables, cusp_conditions, flag_type=region_flag_type, mask=mask_inMagnetosphere) - write_flags(writer, cusp_flags, 'flag_cusps', mask_inMagnetosphere) + write_flags(writer, cusp_flags, 'vg_flag_cusps', mask_inMagnetosphere) # magnetotail lobes @@ -474,11 +474,11 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo } lobe_S_flags = make_region_flags(variables, lobe_S_conditions, flag_type=region_flag_type) - write_flags(writer, lobe_N_flags, 'flag_lobe_N') - write_flags(writer, lobe_S_flags, 'flag_lobe_S') + write_flags(writer, lobe_N_flags, 'vg_flag_lobe_N') + write_flags(writer, lobe_S_flags, 'vg_flag_lobe_S') lobes_flags = make_region_flags(variables, lobes_conditions, flag_type=region_flag_type) - write_flags(writer, lobes_flags, 'flag_lobes') + write_flags(writer, lobes_flags, 'vg_flag_lobes') # lobe density from median densities? @@ -498,7 +498,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo } central_plasma_sheet_flags = make_region_flags(variables, central_plasma_sheet_conditions,flag_type=region_flag_type, mask=mask_inMagnetosphere) - write_flags(writer, central_plasma_sheet_flags, 'flag_central_plasma_sheet', mask_inMagnetosphere) + write_flags(writer, central_plasma_sheet_flags, 'vg_flag_central_plasma_sheet', mask_inMagnetosphere) ## Other boundary layers, PSBL sometimes works @@ -516,9 +516,14 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo def main(): - fileid = 1000 + import sys + + try: + fileid = int(sys.argv[1]) + except Exception as e: + print("Need fileid") datafile = "/wrk-vakka/group/spacephysics/vlasiator/3D/FID/bulk1/bulk1.{:07d}.vlsv".format(fileid) - outfilen = "/wrk-vakka/group/spacephysics/vlasiator/3D/FID/postprocessing/prototyping/FID_mpause_{:07d}.vlsv".format(fileid) + outfilen = "/wrk-vakka/group/spacephysics/vlasiator/3D/FID/postprocessing/magnetopause_sdfs/FID_mpause_{:07d}.vlsv".format(fileid) RegionFlags(datafile, outfilen, regions=["magnetopause"]) From 43bb0126bc6215955abde4fda78bc0bf0e0e2c7c Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Wed, 17 Sep 2025 15:32:31 +0300 Subject: [PATCH 3/7] Add variable cache injection function --- analysator/vlsvfile/vlsvreader.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/analysator/vlsvfile/vlsvreader.py b/analysator/vlsvfile/vlsvreader.py index 915bbc185..d0543edd8 100644 --- a/analysator/vlsvfile/vlsvreader.py +++ b/analysator/vlsvfile/vlsvreader.py @@ -2491,6 +2491,19 @@ def read_variable_to_cache(self, name, operator="pass"): # self.__read_fileindex_for_cellid() return self.__variable_cache[(name,operator)] + def add_cached_variable(self, data, name, operator="pass"): + ''' Add a variable to cache, for the whole grid and after applying + operator. Works by injecting a (variable,operator) key and the data to cache. + + :param data: data array - make sure it is of the correct shape and size. + :param name: Name of the variable (or datareducer) + :param operator: Datareduction operator. "pass" does no operation on data. + + ''' + + # add data to dict, use a tuple of (name,operator) as the key [tuples are immutable and hashable] + self.__variable_cache[(name,operator)] = data + def read_variable(self, name, cellids=-1,operator="pass"): ''' Read variables from the open vlsv file. Arguments: From 238f4013faff1e6481a0213a080df78f8346d880 Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Thu, 18 Sep 2025 13:46:13 +0300 Subject: [PATCH 4/7] Region VTK pipeline cleanup, temp variable cache injection for vlsvreader and reader-based init for VlsvVtkInterface --- analysator/vlsvfile/vlsvvtkinterface.py | 18 +++++++--- scripts/magnetopause.py | 44 +++++++------------------ 2 files changed, 24 insertions(+), 38 deletions(-) diff --git a/analysator/vlsvfile/vlsvvtkinterface.py b/analysator/vlsvfile/vlsvvtkinterface.py index e0e529534..2e2d010b4 100644 --- a/analysator/vlsvfile/vlsvvtkinterface.py +++ b/analysator/vlsvfile/vlsvvtkinterface.py @@ -78,13 +78,21 @@ def __init__(self): def SetFileName(self, filename): if filename != self.__FileName: - self.Modified() - self.__FileName = filename - if self.__FileName is not None: - self.__reader = pt.vlsvfile.VlsvReader(self.__FileName) - fn = os.path.basename(self.__FileName) + if filename is not None: + if self.__reader is None: + self.__reader = pt.vlsvfile.VlsvReader(self.__FileName) + else: + raise ValueError("Tried to change an existing reader ("+self.__reader.file_name+") to " + filename) + self.Modified() + self.__FileName = filename self.__metafile = os.path.join(self.__reader.get_cache_folder(),"vlsvvtkcache.pkl") + def SetReader(self, reader): + self.__reader = reader + self.Modified() + self.__FileName = reader.file_name + self.__metafile = os.path.join(self.__reader.get_cache_folder(),"vlsvvtkcache.pkl") + def GetFileName(self): return self.__FileName diff --git a/scripts/magnetopause.py b/scripts/magnetopause.py index 1ea51a7c9..b466b53af 100644 --- a/scripts/magnetopause.py +++ b/scripts/magnetopause.py @@ -118,14 +118,19 @@ def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds= elif method == "beta_star_with_connectivity": # magnetopause from beta_star, with connectivity if possible - betastar_region = regions.treshold_mask(f.read_variable("vg_beta_star"), beta_star_range) + vg_beta_star =f.read_variable("vg_beta_star") + vg_conn = f.read_variable("vg_connection") + vg_classifier = vg_beta_star*np.min(vg_conn,1) # if closed-closed fieldlines, set var to zero if True: + vtkreader = pt.vlsvfile.VlsvVtkReader() - vtkreader.SetFileName(f.file_name) + vtkreader.SetReader(f) + f.add_cached_variable(vg_classifier, "vg_betastar_classifier") vtkreader.Update() + # vars =vtkreader.findVariablesFromVlsv(getReducers=False) # add here more/other datareducer outputs for downstream use if needed - vars = ["proton/vg_beta_star", "vg_connection","cellid"] + vars = ["vg_betastar_classifier", "cellid"] for var in [v for v in vars if ("vg_" in v.lower()) or (v.lower() == "cellid")]: vtkreader.addArrayFromVlsv(var) @@ -140,38 +145,9 @@ def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds= dual.Update() data = dual.GetOutputDataObject(0) - renamer = vtk.vtkArrayRename() - renamer.SetPointArrayName("proton/vg_beta_star", "vg_beta_star") - renamer.SetInputData(data) - renamer.Update() - - data = renamer.GetOutputDataObject(0) - # print(data) - - newarr = vtk.vtkDoubleArray() - newarr.SetName("BL") - newarr.SetNumberOfComponents(1) # This will be a scalar results - newarr.SetNumberOfTuples(data.GetPointData().GetNumberOfTuples()) # Set how many values to allocate - data.GetPointData().AddArray(newarr) # add the array to the dataset - - - # Setting up the actual calulcation - calcBL = vtk.vtkArrayCalculator() - calcBL.SetInputData(data) # This data we ingest - # We need to specify which arrays we will use to derive the value - calcBL.AddScalarArrayName('vg_beta_star') - calcBL.AddScalarArrayName('vg_connection') - calcBL.SetReplaceInvalidValues(True) - calcBL.SetReplacementValue(np.nan) - # The actual function to calculate - calcBL.SetFunction("vg_beta_star*min(vg_connection,1)") - # Set result array name - calcBL.SetResultArrayName("BL") - calcBL.Update() - data = calcBL.GetOutputDataObject(0) threshold0 = vtk.vtkThreshold() - threshold0.SetInputArrayToProcess(0,0,0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "BL") + threshold0.SetInputArrayToProcess(0,0,0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "vg_betastar_classifier") threshold0.SetInputData(data) threshold0.SetLowerThreshold(beta_star_range[0]) threshold0.SetUpperThreshold(beta_star_range[1]) @@ -180,6 +156,8 @@ def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds= vtkSurface, SDF = regions.vtkSDF(query_points, threshold0.GetOutputDataObject(0)) else: + betastar_region = regions.treshold_mask(vg_classifier, beta_star_range) + try: connectivity_region = regions.treshold_mask(f.read_variable("vg_connection"), 0) # closed-closed magnetic field lines magnetosphere_proper = np.where((connectivity_region==1) | (betastar_region==1), 1, 0) From ac845c2a5626f8ccac0fcd5d993f4d69046e45e8 Mon Sep 17 00:00:00 2001 From: jreimi Date: Mon, 8 Jun 2026 10:44:37 +0300 Subject: [PATCH 5/7] Added a quick fix for vtk to file writer func to accept both connection and dataobject, bug and typo fix --- scripts/magnetopause.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/scripts/magnetopause.py b/scripts/magnetopause.py index b466b53af..1c77f8b2e 100644 --- a/scripts/magnetopause.py +++ b/scripts/magnetopause.py @@ -51,9 +51,16 @@ def write_vtk_surface_to_file(vtkSurface, outfilen): - writer = vtk.vtkXMLPolyDataWriter() - writer.SetInputConnection(vtkSurface.GetOutputPort()) + + writer = vtk.vtkXMLPolyDataWriter() writer.SetFileName(outfilen) + + try: + writer.SetInputConnection(vtkSurface.GetOutputPort()) + except: + writer.SetInputData(vtkSurface) + + writer.Write() logging.info("wrote ", outfilen) @@ -79,7 +86,7 @@ def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds= :kwarg return_SDF: True/False, return array of distances in m to SDF_points in point input order, negative distance inside the surface :kwarg SDF_points: optionally give array of own points to calculate signed distances to. If not given, distances will be to cell centres in the order of f.read_variable("CellID") output :kwarg Delaunay_alpha: alpha (float) to give to vtkDelaunay3d, None -> convex hull, alpha=__: surface egdes longer than __ will be excluded - :kwarg beta_star_range: [min, max] treshold rage to use with methods "beta_star" and "beta_star_with_connectivity" + :kwarg beta_star_range: [min, max] treshold range to use with methods "beta_star" and "beta_star_with_connectivity" :kwarg method_args: dict of keyword arguments to be passed down to external functions (for streamlines and shue) :returns: vtkDataSetSurfaceFilter object of convex hull or alpha shape if return_surface=True, signed distance field of convex hull or alpha shape of magnetopause if return_SDF=True """ @@ -120,7 +127,7 @@ def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds= # magnetopause from beta_star, with connectivity if possible vg_beta_star =f.read_variable("vg_beta_star") vg_conn = f.read_variable("vg_connection") - vg_classifier = vg_beta_star*np.min(vg_conn,1) # if closed-closed fieldlines, set var to zero + vg_classifier = vg_beta_star*np.minimum(vg_conn,1) # if closed-closed fieldlines, set var to zero if True: vtkreader = pt.vlsvfile.VlsvVtkReader() From 6e4c0b4c474998279309f2d24eb3f993eb874515 Mon Sep 17 00:00:00 2001 From: jreimi Date: Mon, 15 Jun 2026 17:14:19 +0300 Subject: [PATCH 6/7] magnetopause streamline creation condition fix --- analysator/calculations/magnetopause_sw_streamline_3d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysator/calculations/magnetopause_sw_streamline_3d.py b/analysator/calculations/magnetopause_sw_streamline_3d.py index 1cab1e051..8483db09a 100644 --- a/analysator/calculations/magnetopause_sw_streamline_3d.py +++ b/analysator/calculations/magnetopause_sw_streamline_3d.py @@ -245,7 +245,7 @@ def make_streamlines(vlsvfile, streamline_seeds=None, seeds_n=25, seeds_x0=20*63 f = pt.vlsvfile.VlsvReader(file_name=vlsvfile) # Create streamline starting points if needed - if streamline_seeds == None: + if not np.any(streamline_seeds): streamline_seeds = np.zeros([seeds_n**2, 3]) t = np.linspace(seeds_range[0], seeds_range[1], seeds_n) From a2542fb52c747163e6ea9cf42fb95bbd7e313c27 Mon Sep 17 00:00:00 2001 From: Jonas Suni <38424612+JonasSuni@users.noreply.github.com> Date: Wed, 17 Jun 2026 14:58:11 +0300 Subject: [PATCH 7/7] Fixed logic for dayside and nightside x_points --- analysator/calculations/magnetopause_sw_streamline_3d.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/analysator/calculations/magnetopause_sw_streamline_3d.py b/analysator/calculations/magnetopause_sw_streamline_3d.py index 8483db09a..8c46162ec 100644 --- a/analysator/calculations/magnetopause_sw_streamline_3d.py +++ b/analysator/calculations/magnetopause_sw_streamline_3d.py @@ -301,7 +301,7 @@ def make_magnetopause(streams, end_x=-15*6371000, x_point_n=50, sector_n=36, ign subsolar_x = np.partition(x_axis_points[:,0], ignore)[ignore] # take the nth point as subsolar point # divide the x point numbers between x > 0 (radial) an x < 0 (yz-planes) by ratio - dayside_x_point_n = int((subsolar_x/np.abs(end_x))*x_point_n) + dayside_x_point_n = int((subsolar_x/np.abs(subsolar_x - end_x))*x_point_n) ### dayside magnetopause ### # for x > 0, look for magnetopause radially @@ -355,11 +355,13 @@ def grid_mid_point(theta_idx, phi_idx): dayside_magnetopause[ring_idx] = ring_points - ### x < 0 magnetopause ### # rest: look for magnetopause in yz-planes ## define points in the x axis where to find magnetopause points on the yz-plane - x_points = np.linspace(0.0, end_x, x_point_n-dayside_x_point_n) + if end_x < 0: + x_points = np.linspace(0.0, end_x, int((np.abs(end_x)/np.abs(subsolar_x - end_x))*x_point_n)) + else: + x_points = np.array([],dtype=float) ## 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]