From a8a2423c5372655ad4301277f333db8509a92dd8 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 10 Nov 2022 11:31:58 -0500 Subject: [PATCH 1/7] Started the process of restructuring the methods for adding bodies. It's partially done, so none of the methods work quite yet. I also added a new method for setting an "ephemeris_date" instance variable that is used when pulling bodies down from JPL Horizons. --- python/swiftest/swiftest/init_cond.py | 4 +- python/swiftest/swiftest/simulation_class.py | 518 ++++++++++++------- 2 files changed, 333 insertions(+), 189 deletions(-) diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py index bf2720a3c..641fa2b68 100644 --- a/python/swiftest/swiftest/init_cond.py +++ b/python/swiftest/swiftest/init_cond.py @@ -18,7 +18,7 @@ from datetime import date import xarray as xr -def solar_system_horizons(plname, idval, param, ephemerides_start_date, ds): +def solar_system_horizons(plname, idval, param, ephemerides_start_date): """ Initializes a Swiftest dataset containing the major planets of the Solar System at a particular data from JPL/Horizons @@ -28,8 +28,6 @@ def solar_system_horizons(plname, idval, param, ephemerides_start_date, ds): Swiftest paramuration parameters. This method uses the unit conversion factors to convert from JPL's AU-day system into the system specified in the param file ephemerides_start_date : string Date to use when obtaining the ephemerides in the format YYYY-MM-DD. - ds : xarray Dataset - Dataset to append Returns ------- diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index f3c786af1..0d26ebdf0 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -14,7 +14,7 @@ from swiftest import init_cond from swiftest import tool from swiftest import constants -from datetime import date +import datetime import xarray as xr import numpy as np import os @@ -26,13 +26,15 @@ Any ) + class Simulation: """ This is a class that defines the basic Swift/Swifter/Swiftest simulation object """ + def __init__(self, codename: Literal["Swiftest", "Swifter", "Swift"] = "Swiftest", - param_file: os.PathLike | str ="param.in", + param_file: os.PathLike | str = "param.in", read_param: bool = False, t0: float = 0.0, tstart: float = 0.0, @@ -44,9 +46,10 @@ def __init__(self, init_cond_file_type: Literal["NETCDF_DOUBLE", "NETCDF_FLOAT", "ASCII"] = "NETCDF_DOUBLE", init_cond_file_name: str | os.PathLike | Dict[str, str] | Dict[str, os.PathLike] | None = None, init_cond_format: Literal["EL", "XV"] = "EL", - output_file_type: Literal["NETCDF_DOUBLE","NETCDF_FLOAT","REAL4","REAL8","XDR4","XDR8"] = "NETCDF_DOUBLE", + output_file_type: Literal[ + "NETCDF_DOUBLE", "NETCDF_FLOAT", "REAL4", "REAL8", "XDR4", "XDR8"] = "NETCDF_DOUBLE", output_file_name: os.PathLike | str | None = None, - output_format: Literal["XV","XVEL"] = "XVEL", + output_format: Literal["XV", "XVEL"] = "XVEL", read_old_output_file: bool = False, MU: str = "MSUN", DU: str = "AU", @@ -68,10 +71,11 @@ def __init__(self, big_discard: bool = False, rhill_present: bool = False, restart: bool = False, - interaction_loops: Literal["TRIANGULAR","FLAT","ADAPTIVE"] = "TRIANGULAR", - encounter_check_loops: Literal["TRIANGULAR","SORTSWEEP","ADAPTIVE"] = "TRIANGULAR", + interaction_loops: Literal["TRIANGULAR", "FLAT", "ADAPTIVE"] = "TRIANGULAR", + encounter_check_loops: Literal["TRIANGULAR", "SORTSWEEP", "ADAPTIVE"] = "TRIANGULAR", + ephemeris_date: str = "MBCL", verbose: bool = True - ): + ): """ Parameters @@ -229,6 +233,7 @@ def __init__(self, If set to True, then more information is printed by Simulation methods as they are executed. Setting to False suppresses most messages other than errors. """ + self.ds = xr.Dataset() self.param = { '! VERSION': f"Swiftest parameter input", @@ -240,6 +245,8 @@ def __init__(self, self.verbose = verbose self.restart = restart + # Width of the column in the printed name of the parameter in parameter getters + self._getter_column_width = '32' # If the parameter file is in a different location than the current working directory, we will need # to use it to properly open bin files self.sim_dir = os.path.dirname(os.path.realpath(param_file)) @@ -250,37 +257,37 @@ def __init__(self, print(f"{param_file} not found.") self.set_parameter(t0=t0, - tstart=tstart, - tstop=tstop, - dt=dt, - tstep_out=tstep_out, - istep_out=istep_out, - istep_dump=istep_dump, - rmin=rmin, rmax=rmax, - MU=MU, DU=DU, TU=TU, - MU2KG=MU2KG, DU2M=DU2M, TU2S=TU2S, - MU_name=MU_name, DU_name=DU_name, TU_name=TU_name, - recompute_unit_values=False, - init_cond_file_type=init_cond_file_type, - init_cond_file_name=init_cond_file_name, - init_cond_format=init_cond_format, - output_file_type=output_file_type, - output_file_name=output_file_name, - output_format=output_format, - close_encounter_check=close_encounter_check, - general_relativity=general_relativity, - fragmentation=fragmentation, - rotation=rotation, - compute_conservation_values=compute_conservation_values, - extra_force=extra_force, - big_discard=big_discard, - rhill_present=rhill_present, - restart=restart, - verbose = False) - + tstart=tstart, + tstop=tstop, + dt=dt, + tstep_out=tstep_out, + istep_out=istep_out, + istep_dump=istep_dump, + rmin=rmin, rmax=rmax, + MU=MU, DU=DU, TU=TU, + MU2KG=MU2KG, DU2M=DU2M, TU2S=TU2S, + MU_name=MU_name, DU_name=DU_name, TU_name=TU_name, + recompute_unit_values=False, + init_cond_file_type=init_cond_file_type, + init_cond_file_name=init_cond_file_name, + init_cond_format=init_cond_format, + output_file_type=output_file_type, + output_file_name=output_file_name, + output_format=output_format, + close_encounter_check=close_encounter_check, + general_relativity=general_relativity, + fragmentation=fragmentation, + rotation=rotation, + compute_conservation_values=compute_conservation_values, + extra_force=extra_force, + big_discard=big_discard, + rhill_present=rhill_present, + restart=restart, + ephemeris_date=ephemeris_date, + verbose=False) if read_old_output_file: - binpath = os.path.join(self.sim_dir,self.param['BIN_OUT']) + binpath = os.path.join(self.sim_dir, self.param['BIN_OUT']) if os.path.exists(binpath): self.bin2xr() else: @@ -334,10 +341,10 @@ def set_simulation_time(self, tstart: float | None = None, tstop: float | None = None, dt: float | None = None, - istep_out : int | None = None, - tstep_out : float | None = None, - istep_dump : int | None = None, - verbose : bool | None = None, + istep_out: int | None = None, + tstep_out: float | None = None, + istep_dump: int | None = None, + verbose: bool | None = None, **kwargs: Any ): """ @@ -375,7 +382,6 @@ def set_simulation_time(self, update_list = [] - if t0 is None: t0 = self.param.pop("T0", None) if t0 is None: @@ -414,7 +420,7 @@ def set_simulation_time(self, if dt is not None and tstop is not None: if dt > (tstop - tstart): print("Error! dt must be smaller than tstop-tstart") - print(f"Setting dt = {tstop-tstart} instead of {dt}") + print(f"Setting dt = {tstop - tstart} instead of {dt}") dt = tstop - tstart if dt is not None: @@ -435,7 +441,7 @@ def set_simulation_time(self, update_list.append("istep_out") if istep_dump is None: - istep_dump = self.param.pop("ISTEP_DUMP",None) + istep_dump = self.param.pop("ISTEP_DUMP", None) if istep_dump is None: istep_dump = istep_out @@ -443,7 +449,7 @@ def set_simulation_time(self, self.param['ISTEP_DUMP'] = istep_dump update_list.append("istep_dump") - time_dict = self.get_simulation_time(update_list,verbose=verbose) + time_dict = self.get_simulation_time(update_list, verbose=verbose) return time_dict @@ -481,13 +487,13 @@ def get_simulation_time(self, arg_list: str | List[str] | None = None, verbose: "istep_dump": "ISTEP_DUMP", } - units = {"t0" : self.TU_name, - "tstart" : self.TU_name, - "tstop" : self.TU_name, - "dt" : self.TU_name, - "tstep_out" : self.TU_name, - "istep_out" : "", - "istep_dump" : ""} + units = {"t0": self.TU_name, + "tstart": self.TU_name, + "tstop": self.TU_name, + "dt": self.TU_name, + "tstep_out": self.TU_name, + "istep_out": "", + "istep_dump": ""} tstep_out = None if arg_list is None or "tstep_out" in arg_list or "istep_out" in arg_list: @@ -505,11 +511,11 @@ def get_simulation_time(self, arg_list: str | List[str] | None = None, verbose: for arg in valid_arg: key = valid_var[arg] if key in time_dict: - print(f"{arg:<32} {time_dict[key]} {units[arg]}") + print(f"{arg:<{self._getter_column_width}} {time_dict[key]} {units[arg]}") else: - print(f"{arg:<32} NOT SET") + print(f"{arg:<{self._getter_column_width}} NOT SET") if tstep_out is not None: - print(f"{'tstep_out':<32} {tstep_out} {units['tstep_out']}") + print(f"{'tstep_out':<{self._getter_column_width}} {tstep_out} {units['tstep_out']}") return time_dict @@ -525,6 +531,8 @@ def set_parameter(self, **kwargs): param : A dictionary of all Simulation parameters that changed """ + + # Setters returning parameter dictionary values param_dict = self.set_unit_system(**kwargs) param_dict.update(self.set_distance_range(**kwargs)) param_dict.update(self.set_feature(**kwargs)) @@ -532,6 +540,9 @@ def set_parameter(self, **kwargs): param_dict.update(self.set_output_files(**kwargs)) param_dict.update(self.set_simulation_time(**kwargs)) + # Non-returning setters + self.set_ephemeris_date(**kwargs) + return param_dict def get_parameter(self, **kwargs): @@ -546,6 +557,8 @@ def get_parameter(self, **kwargs): param : A dictionary of all Simulation parameters requested """ + + # Getters returning parameter dictionary values param_dict = self.get_simulation_time(**kwargs) param_dict.update(self.get_init_cond_files(**kwargs)) param_dict.update(self.get_output_files(**kwargs)) @@ -553,6 +566,10 @@ def get_parameter(self, **kwargs): param_dict.update(self.get_unit_system(**kwargs)) param_dict.update(self.get_feature(**kwargs)) + # Non-returning getters + if not bool(kwargs) or "ephemeris_date" in kwargs: + self.get_ephemeris_date(**kwargs) + return param_dict def set_feature(self, @@ -560,7 +577,7 @@ def set_feature(self, general_relativity: bool | None = None, fragmentation: bool | None = None, rotation: bool | None = None, - compute_conservation_values: bool | None=None, + compute_conservation_values: bool | None = None, extra_force: bool | None = None, big_discard: bool | None = None, rhill_present: bool | None = None, @@ -570,7 +587,7 @@ def set_feature(self, encounter_check_loops: Literal["TRIANGULAR", "SORTSWEEP", "ADAPTIVE"] | None = None, verbose: bool | None = None, **kwargs: Any - ): + ): """ Turns on or off various features of a simulation. @@ -683,7 +700,7 @@ def set_feature(self, update_list.append("restart") if interaction_loops is not None: - valid_vals = ["TRIANGULAR","FLAT","ADAPTIVE"] + valid_vals = ["TRIANGULAR", "FLAT", "ADAPTIVE"] if interaction_loops not in valid_vals: print(f"{interaction_loops} is not a valid option for interaction loops.") print(f"Must be one of {valid_vals}") @@ -705,7 +722,6 @@ def set_feature(self, feature_dict = self.get_feature(update_list, verbose) return feature_dict - def get_feature(self, arg_list: str | List[str] | None = None, verbose: bool | None = None, **kwargs: Any): """ @@ -738,9 +754,9 @@ def get_feature(self, arg_list: str | List[str] | None = None, verbose: bool | N "extra_force": "EXTRA_FORCE", "big_discard": "BIG_DISCARD", "rhill_present": "RHILL_PRESENT", - "restart" : "RESTART", - "interaction_loops" : "INTERACTION_LOOPS", - "encounter_check_loops" : "ENCOUNTER_CHECK" + "restart": "RESTART", + "interaction_loops": "INTERACTION_LOOPS", + "encounter_check_loops": "ENCOUNTER_CHECK" } valid_arg, feature_dict = self._get_valid_arg_list(arg_list, valid_var) @@ -751,13 +767,14 @@ def get_feature(self, arg_list: str | List[str] | None = None, verbose: bool | N if verbose: for arg in valid_arg: key = valid_var[arg] - print(f"{arg:<32} {feature_dict[key]}") + print(f"{arg:<{self._getter_column_width}} {feature_dict[key]}") return feature_dict def set_init_cond_files(self, init_cond_file_type: Literal["NETCDF_DOUBLE", "NETCDF_FLOAT", "ASCII"] | None = None, - init_cond_file_name: str | os.PathLike | Dict[str, str] | Dict[str, os.PathLike] | None = None, + init_cond_file_name: str | os.PathLike | Dict[str, str] | Dict[ + str, os.PathLike] | None = None, init_cond_format: Literal["EL", "XV"] | None = None, verbose: bool | None = None, **kwargs: Any @@ -838,13 +855,12 @@ def ascii_file_input_error_msg(codename): init_cond_keys = ["PL", "TP"] if init_cond_file_type != "ASCII": print(f"{init_cond_file_type} is not supported by {self.codename}. Using ASCII instead") - init_cond_file_type="ASCII" + init_cond_file_type = "ASCII" if init_cond_format != "XV": print(f"{init_cond_format} is not supported by {self.codename}. Using XV instead") init_cond_format = "XV" - - valid_formats={"EL", "XV"} + valid_formats = {"EL", "XV"} if init_cond_format not in valid_formats: print(f"{init_cond_format} is not a valid input format") else: @@ -883,12 +899,10 @@ def ascii_file_input_error_msg(codename): else: self.param["NC_IN"] = init_cond_file_name - - init_cond_file_dict = self.get_init_cond_files(update_list,verbose) + init_cond_file_dict = self.get_init_cond_files(update_list, verbose) return init_cond_file_dict - def get_init_cond_files(self, arg_list: str | List[str] | None = None, verbose: bool | None = None, **kwargs): """ @@ -917,10 +931,10 @@ def get_init_cond_files(self, arg_list: str | List[str] | None = None, verbose: valid_var = {"init_cond_file_type": "IN_TYPE", "init_cond_format": "IN_FORM", - "init_cond_file_name" : "NC_IN", - "init_cond_file_name['CB']" : "CB_IN", - "init_cond_file_name['PL']" : "PL_IN", - "init_cond_file_name['TP']" : "TP_IN", + "init_cond_file_name": "NC_IN", + "init_cond_file_name['CB']": "CB_IN", + "init_cond_file_name['PL']": "PL_IN", + "init_cond_file_name['TP']": "TP_IN", } three_file_args = ["init_cond_file_name['CB']", @@ -958,14 +972,13 @@ def get_init_cond_files(self, arg_list: str | List[str] | None = None, verbose: if verbose: for arg in valid_arg: key = valid_var[arg] - print(f"{arg:<32} {init_cond_file_dict[key]}") + print(f"{arg:<{self._getter_column_width}} {init_cond_file_dict[key]}") return init_cond_file_dict - - def set_output_files(self, - output_file_type: Literal[ "NETCDF_DOUBLE", "NETCDF_FLOAT", "REAL4", "REAL8", "XDR4", "XDR8"] | None = None, + output_file_type: Literal[ + "NETCDF_DOUBLE", "NETCDF_FLOAT", "REAL4", "REAL8", "XDR4", "XDR8"] | None = None, output_file_name: os.PathLike | str | None = None, output_format: Literal["XV", "XVEL"] | None = None, verbose: bool | None = None, @@ -1021,7 +1034,7 @@ def set_output_files(self, output_file_type = self.param.pop("OUT_TYPE", None) if output_file_type is None: output_file_type = "REAL8" - elif output_file_type not in ["REAL4","REAL8","XDR4","XDR8"]: + elif output_file_type not in ["REAL4", "REAL8", "XDR4", "XDR8"]: print(f"{output_file_type} is not compatible with Swifter. Setting to REAL8") output_file_type = "REAL8" elif self.codename == "Swift": @@ -1056,7 +1069,6 @@ def set_output_files(self, return output_file_dict - def get_output_files(self, arg_list: str | List[str] | None = None, verbose: bool | None = None, **kwargs): """ @@ -1096,11 +1108,10 @@ def get_output_files(self, arg_list: str | List[str] | None = None, verbose: boo if verbose: for arg in valid_arg: key = valid_var[arg] - print(f"{arg:<32} {output_file_dict[key]}") + print(f"{arg:<{self._getter_column_width}} {output_file_dict[key]}") return output_file_dict - def set_unit_system(self, MU: str | None = None, DU: str | None = None, @@ -1194,7 +1205,7 @@ def set_unit_system(self, update_list.append("TU") if MU2KG is not None or MU is not None: - MU2KG_old = self.param.pop('MU2KG',None) + MU2KG_old = self.param.pop('MU2KG', None) if MU2KG is not None: self.param['MU2KG'] = MU2KG self.MU_name = None @@ -1217,7 +1228,7 @@ def set_unit_system(self, self.MU_name = "MSun" if DU2M is not None or DU is not None: - DU2M_old = self.param.pop('DU2M',None) + DU2M_old = self.param.pop('DU2M', None) if DU2M is not None: self.param['DU2M'] = DU2M self.DU_name = None @@ -1240,7 +1251,7 @@ def set_unit_system(self, self.DU_name = "AU" if TU2S is not None or TU is not None: - TU2S_old = self.param.pop('TU2S',None) + TU2S_old = self.param.pop('TU2S', None) if TU2S is not None: self.param['TU2S'] = TU2S self.TU_name = None @@ -1259,7 +1270,6 @@ def set_unit_system(self, self.param['TU2S'] = constants.YR2S self.TU_name = "y" - if MU_name is not None: self.MU_name = MU_name if DU_name is not None: @@ -1268,7 +1278,7 @@ def set_unit_system(self, self.TU_name = TU_name self.VU_name = f"{self.DU_name}/{self.TU_name}" - self.GU = constants.GC * self.param['TU2S']**2 * self.param['MU2KG'] / self.param['DU2M']**3 + self.GU = constants.GC * self.param['TU2S'] ** 2 * self.param['MU2KG'] / self.param['DU2M'] ** 3 if recompute_unit_values: self.update_param_units(MU2KG_old, DU2M_old, TU2S_old) @@ -1301,10 +1311,10 @@ def get_unit_system(self, arg_list: str | List[str] | None = None, verbose: bool """ valid_var = { - "MU": "MU2KG", - "DU": "DU2M", - "TU": "TU2S", - } + "MU": "MU2KG", + "DU": "DU2M", + "TU": "TU2S", + } if self.MU_name is None: MU_name = "mass unit" @@ -1320,15 +1330,15 @@ def get_unit_system(self, arg_list: str | List[str] | None = None, verbose: bool TU_name = self.TU_name units1 = { - "MU" : MU_name, - "DU" : DU_name, - "TU" : TU_name - } + "MU": MU_name, + "DU": DU_name, + "TU": TU_name + } units2 = { - "MU" : f"kg / {MU_name}", - "DU" : f"m / {DU_name}", - "TU" : f"s / {TU_name}" - } + "MU": f"kg / {MU_name}", + "DU": f"m / {DU_name}", + "TU": f"s / {TU_name}" + } valid_arg, unit_dict = self._get_valid_arg_list(arg_list, valid_var) @@ -1338,11 +1348,12 @@ def get_unit_system(self, arg_list: str | List[str] | None = None, verbose: bool if verbose: for arg in valid_arg: key = valid_var[arg] - print(f"{arg}: {units1[arg]:<28} {unit_dict[key]} {units2[arg]}") + col_width = str(int(self._getter_column_width) - 4) + print(f"{arg}: {units1[arg]:<{col_width}} {unit_dict[key]} {units2[arg]}") return unit_dict - def update_param_units(self,MU2KG_old,DU2M_old,TU2S_old): + def update_param_units(self, MU2KG_old, DU2M_old, TU2S_old): """ Updates the values of parameters that have units when the units have changed. @@ -1359,8 +1370,8 @@ def update_param_units(self,MU2KG_old,DU2M_old,TU2S_old): """ mass_keys = ['GMTINY', 'MIN_GMFRAG'] - distance_keys = ['CHK_QMIN','CHK_RMIN','CHK_RMAX', 'CHK_EJECT'] - time_keys = ['T0','TSTOP','DT'] + distance_keys = ['CHK_QMIN', 'CHK_RMIN', 'CHK_RMAX', 'CHK_EJECT'] + time_keys = ['T0', 'TSTOP', 'DT'] if MU2KG_old is not None: for k in mass_keys: @@ -1383,11 +1394,10 @@ def update_param_units(self,MU2KG_old,DU2M_old,TU2S_old): if TU2S_old is not None: for k in time_keys: if k in self.param: - self.param[k] *= TU2S_old / self.param['TU2S'] + self.param[k] *= TU2S_old / self.param['TU2S'] return - def set_distance_range(self, rmin: float | None = None, rmax: float | None = None, @@ -1414,7 +1424,7 @@ def set_distance_range(self, """ update_list = [] - CHK_QMIN_RANGE = self.param.pop('CHK_QMIN_RANGE',None) + CHK_QMIN_RANGE = self.param.pop('CHK_QMIN_RANGE', None) if CHK_QMIN_RANGE is None: CHK_QMIN_RANGE = [-1, -1] else: @@ -1430,13 +1440,12 @@ def set_distance_range(self, CHK_QMIN_RANGE[1] = rmax update_list.append("rmax") - self.param['CHK_QMIN_RANGE'] =f"{CHK_QMIN_RANGE[0]} {CHK_QMIN_RANGE[1]}" + self.param['CHK_QMIN_RANGE'] = f"{CHK_QMIN_RANGE[0]} {CHK_QMIN_RANGE[1]}" - range_dict = self.get_distance_range(update_list,verbose=verbose) + range_dict = self.get_distance_range(update_list, verbose=verbose) return range_dict - def get_distance_range(self, arg_list: str | List[str] | None = None, verbose: bool | None = None, **kwargs: Any): """ @@ -1448,6 +1457,8 @@ def get_distance_range(self, arg_list: str | List[str] | None = None, verbose: b arg_list: str | List[str], optional A single string or list of strings containing the names of the features to extract. Default is all of: ["rmin", "rmax"] + verbose: bool, optional + If passed, it will override the Simulation object's verbose flag **kwargs A dictionary of additional keyword argument. This allows this method to be called by the more general get_parameter method, which takes all possible Simulation parameters as arguments, so these are ignored. @@ -1462,13 +1473,13 @@ def get_distance_range(self, arg_list: str | List[str] | None = None, verbose: b valid_var = {"rmin": "CHK_RMIN", "rmax": "CHK_RMAX", "qmin": "CHK_QMIN", - "qminR" : "CHK_QMIN_RANGE" + "qminR": "CHK_QMIN_RANGE" } - units = {"rmin" : self.DU_name, - "rmax" : self.DU_name, - "qmin" : self.DU_name, - "qminR" : self.DU_name, + units = {"rmin": self.DU_name, + "rmax": self.DU_name, + "qmin": self.DU_name, + "qminR": self.DU_name, } if type(arg_list) is str: @@ -1487,66 +1498,203 @@ def get_distance_range(self, arg_list: str | List[str] | None = None, verbose: b if verbose: if "rmin" in valid_arg: key = valid_var["rmin"] - print(f"{'rmin':<32} {range_dict[key]} {units['rmin']}") + print(f"{'rmin':<{self._getter_column_width}} {range_dict[key]} {units['rmin']}") if "rmax" in valid_arg: key = valid_var["rmax"] - print(f"{'rmax':<32} {range_dict[key]} {units['rmax']}") + print(f"{'rmax':<{self._getter_column_width}} {range_dict[key]} {units['rmax']}") return range_dict - - def add(self, plname, date=date.today().isoformat(), idval=None): + def add_solar_system_body(self, + name: str | List[str] | None = None, + id: int | List[int] | None = None, + date: str | None = None, + origin_type: str = "initial_conditions", + source: str = "HORIZONS"): """ - Adds a solar system body to an existing simulation DataSet. + Adds a solar system body to an existing simulation Dataset from the JPL Horizons ephemeris service. Parameters ---------- - plname : string - Name of planet to add (e.g. "Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune" - date : string - Date to use when obtaining the ephemerides in the format YYYY-MM-DD. Defaults to "today" + name : str | List[str], optional + Add solar system body by name. + Currently bodies from the following list will result in fully-massive bodies (they include mass, radius, + and rotation parameters). "Sun" (added as a central body), "Mercury", "Venus", "Earth", "Mars", + "Jupiter", "Saturn", "Uranus", "Neptune", "Pluto" + + Bodies not on this list will be added as test particles, but additional properties can be added later if + desired. + id : int | List[int], optional + Add solar system body by id number. + date : str, optional + ISO-formatted date sto use when obtaining the ephemerides in the format YYYY-MM-DD. Defaults to value + set by `set_ephemeris_date`. + origin_type : str, default "initial_conditions" + The string that will be added to the `origin_type` variable for all bodies added to the list + source : str, default "Horizons" + The source of the ephemerides. + >*Note.* Currently only the JPL Horizons ephemeris is implemented, so this is ignored. Returns ------- - self.ds : xarray dataset + ds : Xarray dataset with body or bodies added. """ - #self.ds = init_cond.solar_system_horizons(plname, idval, self.param, date, self.ds) - self.addp(*init_cond.solar_system_horizons(plname, idval, self.param, date, self.ds)) + + if self.ephemeris_date is None: + self.set_ephemeris_date() + + if date is None: + date = self.ephemeris_date + try: + datetime.datetime.fromisoformat(date) + except: + print(f"{date} is not a valid date format. Must be 'YYYY-MM-DD'. Setting to {self.ephemeris_date}") + date = self.ephemeris_date + + if source.upper() != "HORIZONS": + print("Currently only the JPL Horizons ephemeris service is supported") + + if id is not None and name is not None: + print("Warning! Requesting both id and name could lead to duplicate bodies.") + dsnew = [] + if name is not None: + if type(name) is str: + name = [name] + + if origin_type is None: + origin_type = ['initial_conditions'] * len(name) + + for n in name: + dsnew.append(self.addp(*init_cond.solar_system_horizons(n, self.param, date))) + + + if id is not None: + if type(id) is str: + id = [id] + + if origin_type is None: + origin_type = ['initial_conditions'] * len(id) + + for i in id: + dsnew.append(self.addp(*init_cond.solar_system_horizons(i, self.param, date))) + + return - - - def addp(self, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, rhill=None, Ip1=None, Ip2=None, Ip3=None, rotx=None, roty=None, rotz=None, J2=None,J4=None,t=None): + + + def set_ephemeris_date(self, + ephemeris_date: str | None = None, + verbose: bool | None = None, + **kwargs: Any): + """ + + Parameters + ---------- + ephemeris_date : str, optional + Date to use when obtaining the ephemerides. + Valid options are "today", "MBCL", or date in the format YYYY-MM-DD. + verbose: bool, optional + If passed, it will override the Simulation object's verbose flag + **kwargs + A dictionary of additional keyword argument. This allows this method to be called by the more general + set_parameter method, which takes all possible Simulation parameters as arguments, so these are ignored. + + Returns + ------- + Sets the `ephemeris_date` instance variable. + + """ + + if ephemeris_date is None: + return + + + # The default value is Prof. Minton's Brimley/Cocoon line crossing date (aka MBCL) + minton_bday = datetime.date.fromisoformat('1976-08-05') + brimley_cocoon_line = datetime.timedelta(days=18530) + minton_bcl = (minton_bday + brimley_cocoon_line).isoformat() + + if ephemeris_date is None or ephemeris_date.upper() == "MBCL": + ephemeris_date = minton_bcl + elif ephemeris_date.upper() == "TODAY": + ephemeris_date = datetime.date.today().isoformat() + else: + try: + datetime.datetime.fromisoformat(ephemeris_date) + except: + valid_date_args = ['"MBCL"', '"TODAY"', '"YYYY-MM-DD"'] + print(f"{ephemeris_date} is not a valid format. Valid options include:", ', '.join(valid_date_args)) + print("Using MBCL for date.") + ephemeris_date = minton_bcl + + self.ephemeris_date = ephemeris_date + + ephemeris_date = self.get_ephemeris_date(verbose=verbose) + + return ephemeris_date + + def get_ephemeris_date(self, verbose: bool | None = None, **kwargs: Any): + """ + + Prints the current value of the ephemeris date + + Parameters + ---------- + verbose: bool, optional + If passed, it will override the Simulation object's verbose flag + **kwargs + A dictionary of additional keyword argument. This allows this method to be called by the more general + get_parameter method, which takes all possible Simulation parameters as arguments, so these are ignored. + + Returns + ------- + ephemeris_date: str + The ISO-formatted date string for the ephemeris computation + + """ + if self.ephemeris_date is None: + print(f"ephemeris_date is not set") + return + + if verbose is None: + verbose = self.verbose + if verbose: + print(f"{'ephemeris_date':<{self._getter_column_width}} {self.ephemeris_date}") + + return self.ephemeris_date + + def add_body(self, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, rhill=None, Ip1=None, Ip2=None, + Ip3=None, rotx=None, roty=None, rotz=None, J2=None, J4=None, t=None): """ Adds a body (test particle or massive body) to the internal DataSet given a set up 6 vectors (orbital elements or cartesian state vectors, depending on the value of self.param). Input all angles in degress Parameters ---------- - idvals : integer - Array of body index values. - v1 : float - xh for param['IN_FORM'] == "XV"; a for param['IN_FORM'] == "EL" - v2 : float - yh for param['IN_FORM'] == "XV"; e for param['IN_FORM'] == "EL" - v3 : float - zh for param['IN_FORM'] == "XV"; inc for param['IN_FORM'] == "EL" - v4 : float - vhxh for param['IN_FORM'] == "XV"; capom for param['IN_FORM'] == "EL" - v5 : float - vhyh for param['IN_FORM'] == "XV"; omega for param['IN_FORM'] == "EL" - v6 : float - vhzh for param['IN_FORM'] == "XV"; capm for param['IN_FORM'] == "EL" - Gmass : float - Optional: Array of G*mass values if these are massive bodies - radius : float - Optional: Array radius values if these are massive bodies - rhill : float - Optional: Array rhill values if these are massive bodies - Ip1,y,z : float - Optional: Principal axes moments of inertia - rotx,y,z: float - Optional: Rotation rate vector components - t : float - Optional: Time at start of simulation + v1 : float + xh for param['IN_FORM'] == "XV"; a for param['IN_FORM'] == "EL" + v2 : float + yh for param['IN_FORM'] == "XV"; e for param['IN_FORM'] == "EL" + v3 : float + zh for param['IN_FORM'] == "XV"; inc for param['IN_FORM'] == "EL" + v4 : float + vhxh for param['IN_FORM'] == "XV"; capom for param['IN_FORM'] == "EL" + v5 : float + vhyh for param['IN_FORM'] == "XV"; omega for param['IN_FORM'] == "EL" + v6 : float + vhzh for param['IN_FORM'] == "XV"; capm for param['IN_FORM'] == "EL" + Gmass : float + Optional: Array of G*mass values if these are massive bodies + radius : float + Optional: Array radius values if these are massive bodies + rhill : float + Optional: Array rhill values if these are massive bodies + Ip1,y,z : float + Optional: Principal axes moments of inertia + rotx,y,z: float + Optional: Rotation rate vector components + t : float + Optional: Time at start of simulation + Returns ------- self.ds : xarray dataset @@ -1554,20 +1702,20 @@ def addp(self, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, rh if t is None: t = self.param['T0'] - dsnew = init_cond.vec2xr(self.param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl, Rpl, rhill, Ip1, Ip2, Ip3, rotx, roty, rotz, J2, J4, t) + dsnew = init_cond.vec2xr(self.param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl, Rpl, rhill, Ip1, Ip2, Ip3, + rotx, roty, rotz, J2, J4, t) if dsnew is not None: self.ds = xr.combine_by_coords([self.ds, dsnew]) self.ds['ntp'] = self.ds['id'].where(np.isnan(self.ds['Gmass'])).count(dim="id") self.ds['npl'] = self.ds['id'].where(np.invert(np.isnan(self.ds['Gmass']))).count(dim="id") - 1 if self.param['OUT_TYPE'] == "NETCDF_DOUBLE": - self.ds = io.fix_types(self.ds,ftype=np.float64) + self.ds = io.fix_types(self.ds, ftype=np.float64) elif self.param['OUT_TYPE'] == "NETCDF_FLOAT": - self.ds = io.fix_types(self.ds,ftype=np.float32) + self.ds = io.fix_types(self.ds, ftype=np.float32) return - - + def read_param(self, param_file, codename="Swiftest", verbose=True): """ Reads in a param.in file and determines whether it is a Swift/Swifter/Swiftest parameter file. @@ -1596,8 +1744,7 @@ def read_param(self, param_file, codename="Swiftest", verbose=True): print(f'{codename} is not a recognized code name. Valid options are "Swiftest", "Swifter", or "Swift".') self.codename = "Unknown" return - - + def write_param(self, param_file, param=None): """ Writes to a param.in file and determines whether the output format needs to be converted between Swift/Swifter/Swiftest. @@ -1618,18 +1765,19 @@ def write_param(self, param_file, param=None): if param['IN_TYPE'] == "ASCII": param.pop("NC_IN", None) else: - param.pop("CB_IN",None) - param.pop("PL_IN",None) - param.pop("TP_IN",None) + param.pop("CB_IN", None) + param.pop("PL_IN", None) + param.pop("TP_IN", None) io.write_labeled_param(param, param_file) elif codename == "Swift": io.write_swift_param(param, param_file) else: - print('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".') + print( + 'Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".') return - - - def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", tpname="tp.swiftest.in", cbname="cb.swiftest.in", conversion_questions={}): + + def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", tpname="tp.swiftest.in", + cbname="cb.swiftest.in", conversion_questions={}): """ Converts simulation input files from one format to another (Swift, Swifter, or Swiftest). @@ -1675,14 +1823,13 @@ def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", t goodconversion = False else: goodconversion = False - + if goodconversion: self.write_param(param_file) else: print(f"Conversion from {self.codename} to {newcodename} is not supported.") return oldparam - - + def bin2xr(self): """ Converts simulation output files from a flat binary file to a xarray dataset. @@ -1699,7 +1846,7 @@ def bin2xr(self): # This is done to handle cases where the method is called from a different working directory than the simulation # results param_tmp = self.param.copy() - param_tmp['BIN_OUT'] = os.path.join(self.dir_path,self.param['BIN_OUT']) + param_tmp['BIN_OUT'] = os.path.join(self.dir_path, self.param['BIN_OUT']) if self.codename == "Swiftest": self.ds = io.swiftest2xr(param_tmp, verbose=self.verbose) if self.verbose: print('Swiftest simulation data stored as xarray DataSet .ds') @@ -1709,10 +1856,10 @@ def bin2xr(self): elif self.codename == "Swift": print("Reading Swift simulation data is not implemented yet") else: - print('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".') + print( + 'Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".') return - - + def follow(self, codestyle="Swifter"): """ An implementation of the Swift tool_follow algorithm. Under development. Currently only for Swift simulations. @@ -1731,10 +1878,10 @@ def follow(self, codestyle="Swifter"): if codestyle == "Swift": try: with open('follow.in', 'r') as f: - line = f.readline() # Parameter file (ignored because bin2xr already takes care of it - line = f.readline() # PL file (ignored) - line = f.readline() # TP file (ignored) - line = f.readline() # ifol + line = f.readline() # Parameter file (ignored because bin2xr already takes care of it + line = f.readline() # PL file (ignored) + line = f.readline() # TP file (ignored) + line = f.readline() # ifol i_list = [i for i in line.split(" ") if i.strip()] ifol = int(i_list[0]) line = f.readline() # nskp @@ -1747,11 +1894,10 @@ def follow(self, codestyle="Swifter"): fol = tool.follow_swift(self.ds, ifol=ifol, nskp=nskp) else: fol = None - + if self.verbose: print('follow.out written') return fol - - + def save(self, param_file, framenum=-1, codename="Swiftest"): """ Saves an xarray dataset to a set of input files. @@ -1785,7 +1931,8 @@ def save(self, param_file, framenum=-1, codename="Swiftest"): return - def initial_conditions_from_bin(self, framenum=-1, new_param=None, new_param_file="param.new.in", new_initial_conditions_file="bin_in.nc", restart=False, codename="Swiftest"): + def initial_conditions_from_bin(self, framenum=-1, new_param=None, new_param_file="param.new.in", + new_initial_conditions_file="bin_in.nc", restart=False, codename="Swiftest"): """ Generates a set of input files from a old output file. @@ -1809,7 +1956,6 @@ def initial_conditions_from_bin(self, framenum=-1, new_param=None, new_param_fil frame : NetCDF dataset """ - if codename != "Swiftest": self.save(new_param_file, framenum, codename) return @@ -1830,7 +1976,7 @@ def initial_conditions_from_bin(self, framenum=-1, new_param=None, new_param_fil if self.param['BIN_OUT'] != new_param['BIN_OUT'] and restart: print(f"Restart run with new output file. Copying {self.param['BIN_OUT']} to {new_param['BIN_OUT']}") - shutil.copy2(self.param['BIN_OUT'],new_param['BIN_OUT']) + shutil.copy2(self.param['BIN_OUT'], new_param['BIN_OUT']) new_param['IN_FORM'] = 'XV' if restart: @@ -1842,7 +1988,7 @@ def initial_conditions_from_bin(self, framenum=-1, new_param=None, new_param_fil new_param.pop('TP_IN', None) new_param.pop('CB_IN', None) print(f"Extracting data from dataset at time frame number {framenum} and saving it to {new_param['NC_IN']}") - frame = io.swiftest_xr2infile(self.ds, self.param, infile_name=new_param['NC_IN'],framenum=framenum) + frame = io.swiftest_xr2infile(self.ds, self.param, infile_name=new_param['NC_IN'], framenum=framenum) print(f"Saving parameter configuration file to {new_param_file}") self.write_param(new_param_file, param=new_param) From d07c603f832e78237b2680cc865560502f0032d5 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 10 Nov 2022 18:10:06 -0500 Subject: [PATCH 2/7] Made improvements to the solar system and body add methods. Still not quite finished, but getting close. --- python/swiftest/swiftest/init_cond.py | 222 ++++++-------- python/swiftest/swiftest/simulation_class.py | 294 +++++++++++++------ 2 files changed, 306 insertions(+), 210 deletions(-) diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py index 641fa2b68..c90f1d59b 100644 --- a/python/swiftest/swiftest/init_cond.py +++ b/python/swiftest/swiftest/init_cond.py @@ -8,17 +8,26 @@ You should have received a copy of the GNU General Public License along with Swiftest. If not, see: https://www.gnu.org/licenses. """ +from __future__ import annotations import swiftest import numpy as np +import numpy.typing as npt from astroquery.jplhorizons import Horizons import astropy.units as u from astropy.coordinates import SkyCoord import datetime -from datetime import date import xarray as xr - -def solar_system_horizons(plname, idval, param, ephemerides_start_date): +from typing import ( + Literal, + Dict, + List, + Any +) +def solar_system_horizons(plname: str, + param: Dict, + ephemerides_start_date: str, + idval: int | None = None): """ Initializes a Swiftest dataset containing the major planets of the Solar System at a particular data from JPL/Horizons @@ -118,10 +127,10 @@ def solar_system_horizons(plname, idval, param, ephemerides_start_date): THIRDLONG = np.longdouble(1.0) / np.longdouble(3.0) # Central body value vectors - GMcb = np.array([swiftest.GMSun * param['TU2S'] ** 2 / param['DU2M'] ** 3]) - Rcb = np.array([swiftest.RSun / param['DU2M']]) - J2RP2 = np.array([swiftest.J2Sun * (swiftest.RSun / param['DU2M']) ** 2]) - J4RP4 = np.array([swiftest.J4Sun * (swiftest.RSun / param['DU2M']) ** 4]) + GMcb = swiftest.GMSun * param['TU2S'] ** 2 / param['DU2M'] ** 3 + Rcb = swiftest.RSun / param['DU2M'] + J2RP2 = swiftest.J2Sun * (swiftest.RSun / param['DU2M']) ** 2 + J4RP4 = swiftest.J4Sun * (swiftest.RSun / param['DU2M']) ** 4 solarpole = SkyCoord(ra=286.13 * u.degree, dec=63.87 * u.degree) solarrot = planetrot['Sun'] * param['TU2S'] @@ -145,12 +154,12 @@ def solar_system_horizons(plname, idval, param, ephemerides_start_date): J2 = J2RP2 J4 = J4RP4 if param['ROTATION']: - Ip1 = [Ipsun[0]] - Ip2 = [Ipsun[1]] - Ip3 = [Ipsun[2]] - rotx = [rotcb.x] - roty = [rotcb.y] - rotz = [rotcb.z] + Ip1 = Ipsun[0] + Ip2 = Ipsun[1] + Ip3 = Ipsun[2] + rotx = rotcb.x.value + roty = rotcb.y.value + rotz = rotcb.z.value else: Ip1 = None Ip2 = None @@ -168,12 +177,6 @@ def solar_system_horizons(plname, idval, param, ephemerides_start_date): ephemerides_end_date = tend.isoformat() ephemerides_step = '1d' - v1 = [] - v2 = [] - v3 = [] - v4 = [] - v5 = [] - v6 = [] J2 = None J4 = None @@ -183,42 +186,33 @@ def solar_system_horizons(plname, idval, param, ephemerides_start_date): 'step': ephemerides_step}) if param['IN_FORM'] == 'XV': - v1.append(pldata[plname].vectors()['x'][0] * DCONV) - v2.append(pldata[plname].vectors()['y'][0] * DCONV) - v3.append(pldata[plname].vectors()['z'][0] * DCONV) - v4.append(pldata[plname].vectors()['vx'][0] * VCONV) - v5.append(pldata[plname].vectors()['vy'][0] * VCONV) - v6.append(pldata[plname].vectors()['vz'][0] * VCONV) + v1 = pldata[plname].vectors()['x'][0] * DCONV + v2 = pldata[plname].vectors()['y'][0] * DCONV + v3 = pldata[plname].vectors()['z'][0] * DCONV + v4 = pldata[plname].vectors()['vx'][0] * VCONV + v5 = pldata[plname].vectors()['vy'][0] * VCONV + v6 = pldata[plname].vectors()['vz'][0] * VCONV elif param['IN_FORM'] == 'EL': - v1.append(pldata[plname].elements()['a'][0] * DCONV) - v2.append(pldata[plname].elements()['e'][0]) - v3.append(pldata[plname].elements()['incl'][0]) - v4.append(pldata[plname].elements()['Omega'][0]) - v5.append(pldata[plname].elements()['w'][0]) - v6.append(pldata[plname].elements()['M'][0]) + v1 = pldata[plname].elements()['a'][0] * DCONV + v2 = pldata[plname].elements()['e'][0] + v3 = pldata[plname].elements()['incl'][0] + v4 = pldata[plname].elements()['Omega'][0] + v5 = pldata[plname].elements()['w'][0] + v6 = pldata[plname].elements()['M'][0] if ispl: - GMpl = [] - GMpl.append(GMcb[0] / MSun_over_Mpl[plname]) + GMpl = GMcb / MSun_over_Mpl[plname] if param['CHK_CLOSE']: - Rpl = [] - Rpl.append(planetradius[plname] * DCONV) + Rpl = planetradius[plname] * DCONV else: Rpl = None # Generate planet value vectors if (param['RHILL_PRESENT']): - rhill = [] - rhill.append(pldata[plname].elements()['a'][0] * DCONV * (3 * MSun_over_Mpl[plname]) ** (-THIRDLONG)) + rhill = pldata[plname].elements()['a'][0] * DCONV * (3 * MSun_over_Mpl[plname]) ** (-THIRDLONG) else: rhill = None if (param['ROTATION']): - Ip1 = [] - Ip2 = [] - Ip3 = [] - rotx = [] - roty = [] - rotz = [] RA = pldata[plname].ephemerides()['NPole_RA'][0] DEC = pldata[plname].ephemerides()['NPole_DEC'][0] @@ -226,12 +220,12 @@ def solar_system_horizons(plname, idval, param, ephemerides_start_date): rotrate = planetrot[plname] * param['TU2S'] rot = rotpole.cartesian * rotrate Ip = np.array([0.0, 0.0, planetIpz[plname]]) - Ip1.append(Ip[0]) - Ip2.append(Ip[1]) - Ip3.append(Ip[2]) - rotx.append(rot.x) - roty.append(rot.y) - rotz.append(rot.z) + Ip1 = Ip[0] + Ip2 = Ip[1] + Ip3 = Ip[2] + rotx = rot.x.value + roty = rot.y.value + rotz = rot.z.value else: Ip1 = None Ip2 = None @@ -243,13 +237,33 @@ def solar_system_horizons(plname, idval, param, ephemerides_start_date): GMpl = None if idval is None: - plid = np.array([planetid[plname]], dtype=int) + plid = planetid[plname] else: - plid = np.array([idval], dtype=int) + plid = idval - return plid,[plname],v1,v2,v3,v4,v5,v6,GMpl,Rpl,rhill,Ip1,Ip2,Ip3,rotx,roty,rotz,J2,J4 + return plname,v1,v2,v3,v4,v5,v6,idval,GMpl,Rpl,rhill,Ip1,Ip2,Ip3,rotx,roty,rotz,J2,J4 -def vec2xr(param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, rhill=None, Ip1=None, Ip2=None, Ip3=None, rotx=None, roty=None, rotz=None, J2=None, J4=None,t=0.0): +def vec2xr(param: Dict, + namevals: npt.NDArray[np.str_], + v1: npt.NDArray[np.float_], + v2: npt.NDArray[np.float_], + v3: npt.NDArray[np.float_], + v4: npt.NDArray[np.float_], + v5: npt.NDArray[np.float_], + v6: npt.NDArray[np.float_], + idvals: npt.NDArray[np.int_], + GMpl: npt.NDArray[np.float_] | None=None, + Rpl: npt.NDArray[np.float_] | None=None, + rhill: npt.NDArray[np.float_] | None=None, + Ip1: npt.NDArray[np.float_] | None=None, + Ip2: npt.NDArray[np.float_] | None=None, + Ip3: npt.NDArray[np.float_] | None=None, + rotx: npt.NDArray[np.float_] | None=None, + roty: npt.NDArray[np.float_] | None=None, + rotz: npt.NDArray[np.float_] | None=None, + J2: npt.NDArray[np.float_] | None=None, + J4: npt.NDArray[np.float_] | None=None, + t: float=0.0): """ Converts and stores the variables of all bodies in an xarray dataset. @@ -297,11 +311,6 @@ def vec2xr(param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, ------- ds : xarray dataset """ - if v1 is None: # This is the central body - iscb = True - else: - iscb = False - if param['ROTATION']: if Ip1 is None: Ip1 = np.full_like(v1, 0.4) @@ -318,11 +327,18 @@ def vec2xr(param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, dims = ['time', 'id', 'vec'] infodims = ['id', 'vec'] - if not iscb and GMpl is not None: + + # The central body is always given id 0 + icb = idvals == 0 + iscb = any(icb) + + if GMpl is not None: ispl = True + ipl = ~np.isnan(GMpl) + itp = np.isnan(GMpl) else: ispl = False - + if ispl and param['CHK_CLOSE'] and Rpl is None: print("Massive bodies need a radius value.") return None @@ -337,83 +353,33 @@ def vec2xr(param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, param['OUT_FORM'] = old_out_form vec_str = np.vstack([namevals]) label_str = ["name"] + particle_type = np.empty_like(namevals) if iscb: label_float = clab.copy() - vec_float = np.vstack([GMpl,Rpl,J2,J4]) + vec_float = np.vstack([GMpl[icb],Rpl[icb],J2[icb],J4[icb]]) + if param['ROTATION']: + vec_float = np.vstack([vec_float, Ip1[icb], Ip2[icb], Ip3[icb], rotx[icb], roty[icb], rotz[icb]]) + particle_type[icb] = "Central Body" + # vec_float = np.vstack([v1, v2, v3, v4, v5, v6]) + if ispl: + label_float = plab.copy() + vec_float = np.vstack([vec_float, GMpl]) + if param['CHK_CLOSE']: + vec_float = np.vstack([vec_float, Rpl]) + if param['RHILL_PRESENT']: + vec_float = np.vstack([vec_float, rhill]) if param['ROTATION']: vec_float = np.vstack([vec_float, Ip1, Ip2, Ip3, rotx, roty, rotz]) - particle_type = "Central Body" + particle_type[ipl] = np.repeat("Massive Body",idvals.size) else: - vec_float = np.vstack([v1, v2, v3, v4, v5, v6]) - if ispl: - label_float = plab.copy() - vec_float = np.vstack([vec_float, GMpl]) - if param['CHK_CLOSE']: - vec_float = np.vstack([vec_float, Rpl]) - if param['RHILL_PRESENT']: - vec_float = np.vstack([vec_float, rhill]) - if param['ROTATION']: - vec_float = np.vstack([vec_float, Ip1, Ip2, Ip3, rotx, roty, rotz]) - particle_type = np.repeat("Massive Body",idvals.size) - else: - label_float = tlab.copy() - particle_type = np.repeat("Test Particle",idvals.size) - origin_type = np.repeat("User Added Body",idvals.size) - origin_time = np.full_like(v1,t) - collision_id = np.full_like(idvals,0) - origin_xhx = v1 - origin_xhy = v2 - origin_xhz = v3 - origin_vhx = v4 - origin_vhy = v5 - origin_vhz = v6 - discard_time = np.full_like(v1,-1.0) - status = np.repeat("ACTIVE",idvals.size) - discard_xhx = np.zeros_like(v1) - discard_xhy = np.zeros_like(v1) - discard_xhz = np.zeros_like(v1) - discard_vhx = np.zeros_like(v1) - discard_vhy = np.zeros_like(v1) - discard_vhz = np.zeros_like(v1) - discard_body_id = np.full_like(idvals,-1) - info_vec_float = np.vstack([ - origin_time, - origin_xhx, - origin_xhy, - origin_xhz, - origin_vhx, - origin_vhy, - origin_vhz, - discard_time, - discard_xhx, - discard_xhy, - discard_xhz, - discard_vhx, - discard_vhy, - discard_vhz]) - info_vec_int = np.vstack([collision_id, discard_body_id]) - info_vec_str = np.vstack([particle_type, origin_type, status]) - frame_float = info_vec_float.T - frame_int = info_vec_int.T - frame_str = info_vec_str.T - if param['IN_TYPE'] == 'NETCDF_FLOAT': - ftype=np.float32 - elif param['IN_TYPE'] == 'NETCDF_DOUBLE' or param['IN_TYPE'] == 'ASCII': - ftype=np.float64 - da_float = xr.DataArray(frame_float, dims=infodims, coords={'id': idvals, 'vec': infolab_float}).astype(ftype) - da_int = xr.DataArray(frame_int, dims=infodims, coords={'id': idvals, 'vec': infolab_int}) - da_str = xr.DataArray(frame_str, dims=infodims, coords={'id': idvals, 'vec': infolab_str}) - ds_float = da_float.to_dataset(dim="vec") - ds_int = da_int.to_dataset(dim="vec") - ds_str = da_str.to_dataset(dim="vec") - info_ds = xr.combine_by_coords([ds_float, ds_int, ds_str]) - + label_float = tlab.copy() + particle_type[itp] = np.repeat("Test Particle",idvals.size) frame_float = np.expand_dims(vec_float.T, axis=0) frame_str = vec_str.T - da_float = xr.DataArray(frame_float, dims=dims, coords={'time': [t], 'id': idvals, 'vec': label_float}).astype(ftype) + da_float = xr.DataArray(frame_float, dims=dims, coords={'time': [t], 'id': idvals, 'vec': label_float}) da_str= xr.DataArray(frame_str, dims=infodims, coords={'id': idvals, 'vec': label_str}) ds_float = da_float.to_dataset(dim="vec") ds_str = da_str.to_dataset(dim="vec") - ds = xr.combine_by_coords([ds_float, ds_str,info_ds]) + ds = xr.combine_by_coords([ds_float, ds_str]) return ds \ No newline at end of file diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 0d26ebdf0..01c3b8414 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -14,10 +14,11 @@ from swiftest import init_cond from swiftest import tool from swiftest import constants +import os import datetime import xarray as xr import numpy as np -import os +import numpy.typing as npt import shutil from typing import ( Literal, @@ -1506,39 +1507,58 @@ def get_distance_range(self, arg_list: str | List[str] | None = None, verbose: b return range_dict def add_solar_system_body(self, - name: str | List[str] | None = None, - id: int | List[int] | None = None, + name: str | List[str], + ephemeris_id: int | List[int] | None = None, date: str | None = None, - origin_type: str = "initial_conditions", source: str = "HORIZONS"): """ Adds a solar system body to an existing simulation Dataset from the JPL Horizons ephemeris service. - + + The following are name/ephemeris_id pairs that are currently known to Swiftest, and therefore have + physical properties that can be used to make massive bodies. + + Sun : 0 + Mercury : 1 + Venus : 2 + Earth : 3 + Mars : 4 + Jupiter : 5 + Saturn : 6 + Uranus : 7 + Neptune : 8 + Pluto : 9 + Parameters ---------- - name : str | List[str], optional - Add solar system body by name. - Currently bodies from the following list will result in fully-massive bodies (they include mass, radius, - and rotation parameters). "Sun" (added as a central body), "Mercury", "Venus", "Earth", "Mars", - "Jupiter", "Saturn", "Uranus", "Neptune", "Pluto" - - Bodies not on this list will be added as test particles, but additional properties can be added later if - desired. - id : int | List[int], optional - Add solar system body by id number. - date : str, optional - ISO-formatted date sto use when obtaining the ephemerides in the format YYYY-MM-DD. Defaults to value - set by `set_ephemeris_date`. - origin_type : str, default "initial_conditions" - The string that will be added to the `origin_type` variable for all bodies added to the list - source : str, default "Horizons" - The source of the ephemerides. - >*Note.* Currently only the JPL Horizons ephemeris is implemented, so this is ignored. + name : str | List[str] + Add solar system body by name. + Bodies not on this list will be added as test particles, but additional properties can be added later if + desired. + ephemeris_id : int | List[int], optional but must be the same length as `name` if passed. + Use id if the body you wish to add is recognized by Swiftest. In that case, the id is passed to the + ephemeris service and the name is used. The body specified by `id` supercedes that given by `name`. + date : str, optional + ISO-formatted date sto use when obtaining the ephemerides in the format YYYY-MM-DD. Defaults to value + set by `set_ephemeris_date`. + source : str, default "Horizons" + The source of the ephemerides. + >*Note.* Currently only the JPL Horizons ephemeris is implemented, so this is ignored. Returns ------- ds : Xarray dataset with body or bodies added. """ + if type(name) is str: + name = [name] + if ephemeris_id is not None: + if type(ephemeris_id) is int: + ephemeris_id = [ephemeris_id] + if len(ephemeris_id) != len(name): + print(f"Error! The length of ephemeris_id ({len(ephemeris_id)}) does not match the length of name ({len(name)})") + return None + else: + ephemeris_id = [None] * len(name) + if self.ephemeris_date is None: self.set_ephemeris_date() @@ -1553,32 +1573,59 @@ def add_solar_system_body(self, if source.upper() != "HORIZONS": print("Currently only the JPL Horizons ephemeris service is supported") - if id is not None and name is not None: - print("Warning! Requesting both id and name could lead to duplicate bodies.") - dsnew = [] - if name is not None: - if type(name) is str: - name = [name] - - if origin_type is None: - origin_type = ['initial_conditions'] * len(name) - - for n in name: - dsnew.append(self.addp(*init_cond.solar_system_horizons(n, self.param, date))) - - - if id is not None: - if type(id) is str: - id = [id] - - if origin_type is None: - origin_type = ['initial_conditions'] * len(id) - - for i in id: - dsnew.append(self.addp(*init_cond.solar_system_horizons(i, self.param, date))) - - - return + body_list = [] + for i,n in enumerate(name): + body_list.append(init_cond.solar_system_horizons(n, self.param, date, idval=ephemeris_id[i])) + + #Convert the list receieved from the solar_system_horizons output and turn it into arguments to vec2xr + name,v1,v2,v3,v4,v5,v6,ephemeris_id,GMpl,Rpl,rhill,Ip1,Ip2,Ip3,rotx,roty,rotz,J2,J4 = tuple(np.squeeze(np.hsplit(np.array(body_list),19))) + + v1 = v1.astype(np.float64) + v2 = v2.astype(np.float64) + v3 = v3.astype(np.float64) + v4 = v4.astype(np.float64) + v5 = v5.astype(np.float64) + v6 = v6.astype(np.float64) + ephemeris_id = ephemeris_id.astype(int) + GMpl = GMpl.astype(np.float64) + Rpl = Rpl.astype(np.float64) + rhill = rhill.astype(np.float64) + Ip1 = Ip1.astype(np.float64) + Ip2 = Ip2.astype(np.float64) + Ip3 = Ip3.astype(np.float64) + rotx = rotx.astype(np.float64) + roty = roty.astype(np.float64) + rotz = rotz.astype(np.float64) + J2 = J2.astype(np.float64) + J4 = J4.astype(np.float64) + + if all(np.isnan(GMpl)): + GMpl = None + if all(np.isnan(Rpl)): + Rpl = None + if all(np.isnan(rhill)): + rhill = None + if all(np.isnan(Ip1)): + Ip1 = None + if all(np.isnan(Ip2)): + Ip2 = None + if all(np.isnan(Ip3)): + Ip3 = None + if all(np.isnan(rotx)): + rotx = None + if all(np.isnan(roty)): + roty = None + if all(np.isnan(rotz)): + rotz = None + if all(np.isnan(J2)): + J2 = None + if all(np.isnan(J4)): + J4 = None + + + dsnew = init_cond.vec2xr(self.param,name,v1,v2,v3,v4,v5,v6,ephemeris_id,GMpl,Rpl,rhill,Ip1,Ip2,Ip3,rotx,roty,rotz,J2,J4) + + return body_list def set_ephemeris_date(self, @@ -1662,48 +1709,131 @@ def get_ephemeris_date(self, verbose: bool | None = None, **kwargs: Any): return self.ephemeris_date - def add_body(self, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, rhill=None, Ip1=None, Ip2=None, - Ip3=None, rotx=None, roty=None, rotz=None, J2=None, J4=None, t=None): + def add_body(self, + name: str | List[str] | npt.NDArray[np.str_], + v1: float | List[float] | npt.NDArray[np.float_], + v2: float | List[float] | npt.NDArray[np.float_], + v3: float | List[float] | npt.NDArray[np.float_], + v4: float | List[float] | npt.NDArray[np.float_], + v5: float | List[float] | npt.NDArray[np.float_], + v6: float | List[float] | npt.NDArray[np.float_], + idvals: int | list[int] | npt.NDArray[np.int_] | None=None, + GMpl: float | List[float] | npt.NDArray[np.float_] | None=None, + Rpl: float | List[float] | npt.NDArray[np.float_] | None=None, + rhill: float | List[float] | npt.NDArray[np.float_] | None=None, + Ip1: float | List[float] | npt.NDArray[np.float_] | None=None, + Ip2: float | List[float] | npt.NDArray[np.float_] | None=None, + Ip3: float | List[float] | npt.NDArray[np.float_] | None=None, + rotx: float | List[float] | npt.NDArray[np.float_] | None=None, + roty: float | List[float] | npt.NDArray[np.float_] | None=None, + rotz: float | List[float] | npt.NDArray[np.float_] | None=None, + J2: float | List[float] | npt.NDArray[np.float_] | None=None, + J4: float | List[float] | npt.NDArray[np.float_] | None=None): """ Adds a body (test particle or massive body) to the internal DataSet given a set up 6 vectors (orbital elements - or cartesian state vectors, depending on the value of self.param). Input all angles in degress + or cartesian state vectors, depending on the value of self.param). Input all angles in degress. + + This method will update self.ds with the new body or bodies added to the existing Dataset. Parameters ---------- - v1 : float - xh for param['IN_FORM'] == "XV"; a for param['IN_FORM'] == "EL" - v2 : float - yh for param['IN_FORM'] == "XV"; e for param['IN_FORM'] == "EL" - v3 : float - zh for param['IN_FORM'] == "XV"; inc for param['IN_FORM'] == "EL" - v4 : float - vhxh for param['IN_FORM'] == "XV"; capom for param['IN_FORM'] == "EL" - v5 : float - vhyh for param['IN_FORM'] == "XV"; omega for param['IN_FORM'] == "EL" - v6 : float - vhzh for param['IN_FORM'] == "XV"; capm for param['IN_FORM'] == "EL" - Gmass : float - Optional: Array of G*mass values if these are massive bodies - radius : float - Optional: Array radius values if these are massive bodies - rhill : float - Optional: Array rhill values if these are massive bodies - Ip1,y,z : float - Optional: Principal axes moments of inertia - rotx,y,z: float - Optional: Rotation rate vector components - t : float - Optional: Time at start of simulation + name : str or array-like of str + Name or names of + v1 : float or array-like of float + xhx for param['IN_FORM'] == "XV"; a for param['IN_FORM'] == "EL" + v2 : float or array-like of float + xhy for param['IN_FORM'] == "XV"; e for param['IN_FORM'] == "EL" + v3 : float or array-like of float + xhz for param['IN_FORM'] == "XV"; inc for param['IN_FORM'] == "EL" + v4 : float or array-like of float + vhx for param['IN_FORM'] == "XV"; capom for param['IN_FORM'] == "EL" + v5 : float or array-like of float + vhy for param['IN_FORM'] == "XV"; omega for param['IN_FORM'] == "EL" + v6 : float or array-like of float + vhz for param['IN_FORM'] == "XV"; capm for param['IN_FORM'] == "EL" + idvals : int or array-like of int, optional + Unique id values. If not passed, this will be computed based on the pre-existing Dataset ids. + Gmass : float or array-like of float, optional + G*mass values if these are massive bodies + radius : float or array-like of float, optional + Radius values if these are massive bodies + rhill : float, optional + Hill's radius values if these are massive bodies + Ip1,y,z : float, optional + Principal axes moments of inertia these are massive bodies with rotation enabled + rotx,y,z: float, optional + Rotation rate vector components if these are massive bodies with rotation enabled Returns ------- - self.ds : xarray dataset + ds : Xarray Dataset + Dasaset containing the body or bodies that were added + """ - if t is None: - t = self.param['T0'] - dsnew = init_cond.vec2xr(self.param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl, Rpl, rhill, Ip1, Ip2, Ip3, - rotx, roty, rotz, J2, J4, t) + #convert all inputs to numpy arrays + def input_to_array(val,t,n=None): + if t == "f": + t = np.float64 + elif t == "i": + t = np.int64 + elif t == "s": + t = np.str + if val is None: + return None + elif np.isscalar(val): + val = np.array([val],dtype=t) + elif type(val) is list: + val = np.array(val,dtype=t) + + if n is None: + return val, len(val) + else: + if n != len(val): + raise ValueError(f"Error! Mismatched array lengths in add_body. Got {len(val)} when expecting {n}") + return val + + + name,nbodies = input_to_array(name,"s") + v1 = input_to_array(v1,"f",nbodies) + v2 = input_to_array(v2,"f",nbodies) + v3 = input_to_array(v3,"f",nbodies) + v4 = input_to_array(v4,"f",nbodies) + v5 = input_to_array(v5,"f",nbodies) + v6 = input_to_array(v6,"f",nbodies) + idvals = input_to_array(idvals,"i",nbodies) + GMpl = input_to_array(GMpl,"f",nbodies) + rhill = input_to_array(rhill,"f",nbodies) + Rpl = input_to_array(Rpl,"f",nbodies) + Ip1 = input_to_array(Ip1,"f",nbodies) + Ip2 = input_to_array(Ip2,"f",nbodies) + Ip3 = input_to_array(Ip3,"f",nbodies) + rotx = input_to_array(rotx,"f",nbodies) + roty = input_to_array(roty,"f",nbodies) + rotz = input_to_array(rotz,"f",nbodies) + J2 = input_to_array(J2,"f",nbodies) + J4 = input_to_array(J4,"f",nbodies) + + if len(self.ds) == 0: + maxid = -1 + else: + maxid = self.ds.id.max().values[()] + + if idvals is None: + idvals = np.arange(start=maxid+1,stop=maxid+1+nbodies,dtype=int) + + if len(self.ds) > 0: + dup_id = np.in1d(idvals,self.ds.id) + if any(dup_id): + raise ValueError(f"Duplicate ids detected: ", *idvals[dup_id]) + + t = self.param['TSTART'] + + dsnew = init_cond.vec2xr(self.param, idvals, name, v1, v2, v3, v4, v5, v6, + GMpl=GMpl, Rpl=Rpl, rhill=rhill, + Ip1=Ip1, Ip2=Ip2, Ip3=Ip3, + rotx=rotx, roty=roty, rotz=rotz, + J2=J2, J4=J4, t=t) if dsnew is not None: self.ds = xr.combine_by_coords([self.ds, dsnew]) self.ds['ntp'] = self.ds['id'].where(np.isnan(self.ds['Gmass'])).count(dim="id") @@ -1714,7 +1844,7 @@ def add_body(self, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, rhill= elif self.param['OUT_TYPE'] == "NETCDF_FLOAT": self.ds = io.fix_types(self.ds, ftype=np.float32) - return + return dsnew def read_param(self, param_file, codename="Swiftest", verbose=True): """ From 3c4c2b638ea3cc606695145644a67d30a9337a2b Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 10 Nov 2022 19:56:39 -0500 Subject: [PATCH 3/7] Finished add_solar_system_body method --- python/swiftest/swiftest/init_cond.py | 65 ++++++++++++-------- python/swiftest/swiftest/simulation_class.py | 4 +- 2 files changed, 42 insertions(+), 27 deletions(-) diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py index c90f1d59b..a5bf89c2b 100644 --- a/python/swiftest/swiftest/init_cond.py +++ b/python/swiftest/swiftest/init_cond.py @@ -329,15 +329,21 @@ def vec2xr(param: Dict, infodims = ['id', 'vec'] # The central body is always given id 0 - icb = idvals == 0 - iscb = any(icb) if GMpl is not None: - ispl = True - ipl = ~np.isnan(GMpl) - itp = np.isnan(GMpl) + icb = (~np.isnan(GMpl)) & (idvals == 0) + ipl = (~np.isnan(GMpl)) & (idvals != 0) + itp = (np.isnan(GMpl)) & (idvals != 0) + iscb = any(icb) + ispl = any(ipl) + istp = any(itp) else: + icb = np.full_like(idvals,False) + ipl = np.full_like(idvals,False) + itp = idvals != 0 + iscb = False ispl = False + istp = any(itp) if ispl and param['CHK_CLOSE'] and Rpl is None: print("Massive bodies need a radius value.") @@ -351,35 +357,42 @@ def vec2xr(param: Dict, param['OUT_FORM'] = param['IN_FORM'] clab, plab, tlab, infolab_float, infolab_int, infolab_str = swiftest.io.make_swiftest_labels(param) param['OUT_FORM'] = old_out_form - vec_str = np.vstack([namevals]) - label_str = ["name"] particle_type = np.empty_like(namevals) + vec = np.vstack([v1,v2,v3,v4,v5,v6]) + if iscb: - label_float = clab.copy() - vec_float = np.vstack([GMpl[icb],Rpl[icb],J2[icb],J4[icb]]) + lab_cb = clab.copy() + vec_cb = np.vstack([GMpl[icb],Rpl[icb],J2[icb],J4[icb]]) if param['ROTATION']: - vec_float = np.vstack([vec_float, Ip1[icb], Ip2[icb], Ip3[icb], rotx[icb], roty[icb], rotz[icb]]) + vec_cb = np.vstack([vec_cb, Ip1[icb], Ip2[icb], Ip3[icb], rotx[icb], roty[icb], rotz[icb]]) particle_type[icb] = "Central Body" - # vec_float = np.vstack([v1, v2, v3, v4, v5, v6]) + vec_cb = np.expand_dims(vec_cb.T,axis=0) # Make way for the time dimension! + ds_cb = xr.DataArray(vec_cb, dims=dims, coords={'time': [t], 'id': idvals[icb], 'vec': lab_cb}).to_dataset(dim='vec') + else: + ds_cb = xr.Dataset() if ispl: - label_float = plab.copy() - vec_float = np.vstack([vec_float, GMpl]) + lab_pl = plab.copy() + vec_pl = np.vstack([vec[:,ipl], GMpl[ipl]]) if param['CHK_CLOSE']: - vec_float = np.vstack([vec_float, Rpl]) + vec_pl = np.vstack([vec_pl, Rpl[ipl]]) if param['RHILL_PRESENT']: - vec_float = np.vstack([vec_float, rhill]) + vec_pl = np.vstack([vec_pl, rhill[ipl]]) if param['ROTATION']: - vec_float = np.vstack([vec_float, Ip1, Ip2, Ip3, rotx, roty, rotz]) - particle_type[ipl] = np.repeat("Massive Body",idvals.size) + vec_pl = np.vstack([vec_pl, Ip1[ipl], Ip2[ipl], Ip3[ipl], rotx[ipl], roty[ipl], rotz[ipl]]) + particle_type[ipl] = np.repeat("Massive Body",idvals[ipl].size) + vec_pl = np.expand_dims(vec_pl.T,axis=0) # Make way for the time dimension! + ds_pl = xr.DataArray(vec_pl, dims=dims, coords={'time': [t], 'id': idvals[ipl], 'vec': lab_pl}).to_dataset(dim='vec') else: - label_float = tlab.copy() - particle_type[itp] = np.repeat("Test Particle",idvals.size) - frame_float = np.expand_dims(vec_float.T, axis=0) - frame_str = vec_str.T - da_float = xr.DataArray(frame_float, dims=dims, coords={'time': [t], 'id': idvals, 'vec': label_float}) - da_str= xr.DataArray(frame_str, dims=infodims, coords={'id': idvals, 'vec': label_str}) - ds_float = da_float.to_dataset(dim="vec") - ds_str = da_str.to_dataset(dim="vec") - ds = xr.combine_by_coords([ds_float, ds_str]) + ds_pl = xr.Dataset() + if istp: + lab_tp = tlab.copy() + vec_tp = np.expand_dims(vec[:,itp].T,axis=0) # Make way for the time dimension! + ds_tp = xr.DataArray(vec_tp, dims=dims, coords={'time': [t], 'id': idvals[itp], 'vec': lab_tp}).to_dataset(dim='vec') + particle_type[itp] = np.repeat("Test Particle",idvals[itp].size) + else: + ds_tp = xr.Dataset() + + ds_info = xr.DataArray(np.vstack([namevals,particle_type]).T, dims=infodims, coords={'id': idvals, 'vec' : ["name", "particle_type"]}).to_dataset(dim='vec') + ds = xr.combine_by_coords([ds_cb, ds_pl, ds_tp, ds_info]) return ds \ No newline at end of file diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 01c3b8414..9ecbc1e51 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -1625,7 +1625,9 @@ def add_solar_system_body(self, dsnew = init_cond.vec2xr(self.param,name,v1,v2,v3,v4,v5,v6,ephemeris_id,GMpl,Rpl,rhill,Ip1,Ip2,Ip3,rotx,roty,rotz,J2,J4) - return body_list + self.ds = xr.combine_by_coords([self.ds,dsnew]) + + return dsnew def set_ephemeris_date(self, From d166719e3b44a2a05fcbe96e8f23088395439512 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 10 Nov 2022 20:33:35 -0500 Subject: [PATCH 4/7] Finished the new add_body and add_solar_system_body methods --- python/swiftest/swiftest/init_cond.py | 13 +++-- python/swiftest/swiftest/simulation_class.py | 52 +++++++++++++++++--- 2 files changed, 53 insertions(+), 12 deletions(-) diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py index a5bf89c2b..ea28bcb8c 100644 --- a/python/swiftest/swiftest/init_cond.py +++ b/python/swiftest/swiftest/init_cond.py @@ -369,7 +369,7 @@ def vec2xr(param: Dict, vec_cb = np.expand_dims(vec_cb.T,axis=0) # Make way for the time dimension! ds_cb = xr.DataArray(vec_cb, dims=dims, coords={'time': [t], 'id': idvals[icb], 'vec': lab_cb}).to_dataset(dim='vec') else: - ds_cb = xr.Dataset() + ds_cb = None if ispl: lab_pl = plab.copy() vec_pl = np.vstack([vec[:,ipl], GMpl[ipl]]) @@ -383,16 +383,21 @@ def vec2xr(param: Dict, vec_pl = np.expand_dims(vec_pl.T,axis=0) # Make way for the time dimension! ds_pl = xr.DataArray(vec_pl, dims=dims, coords={'time': [t], 'id': idvals[ipl], 'vec': lab_pl}).to_dataset(dim='vec') else: - ds_pl = xr.Dataset() + ds_pl = None if istp: lab_tp = tlab.copy() vec_tp = np.expand_dims(vec[:,itp].T,axis=0) # Make way for the time dimension! ds_tp = xr.DataArray(vec_tp, dims=dims, coords={'time': [t], 'id': idvals[itp], 'vec': lab_tp}).to_dataset(dim='vec') particle_type[itp] = np.repeat("Test Particle",idvals[itp].size) else: - ds_tp = xr.Dataset() + ds_tp = None ds_info = xr.DataArray(np.vstack([namevals,particle_type]).T, dims=infodims, coords={'id': idvals, 'vec' : ["name", "particle_type"]}).to_dataset(dim='vec') - ds = xr.combine_by_coords([ds_cb, ds_pl, ds_tp, ds_info]) + ds = [d for d in [ds_cb, ds_pl, ds_tp] if d is not None] + if len(ds) > 1: + ds = xr.combine_by_coords(ds) + else: + ds = ds[0] + ds = xr.merge([ds_info,ds]) return ds \ No newline at end of file diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 9ecbc1e51..8e1ac24e0 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -1622,10 +1622,15 @@ def add_solar_system_body(self, if all(np.isnan(J4)): J4 = None + t = self.param['TSTART'] - dsnew = init_cond.vec2xr(self.param,name,v1,v2,v3,v4,v5,v6,ephemeris_id,GMpl,Rpl,rhill,Ip1,Ip2,Ip3,rotx,roty,rotz,J2,J4) + dsnew = init_cond.vec2xr(self.param,name,v1,v2,v3,v4,v5,v6,ephemeris_id, + GMpl=GMpl, Rpl=Rpl, rhill=rhill, + Ip1=Ip1, Ip2=Ip2, Ip3=Ip3, + rotx=rotx, roty=roty, rotz=rotz, + J2=J2, J4=J4, t=t) - self.ds = xr.combine_by_coords([self.ds,dsnew]) + dsnew = self._combine_and_fix_dsnew(dsnew) return dsnew @@ -1831,19 +1836,50 @@ def input_to_array(val,t,n=None): t = self.param['TSTART'] - dsnew = init_cond.vec2xr(self.param, idvals, name, v1, v2, v3, v4, v5, v6, + dsnew = init_cond.vec2xr(self.param, name, v1, v2, v3, v4, v5, v6, idvals, GMpl=GMpl, Rpl=Rpl, rhill=rhill, Ip1=Ip1, Ip2=Ip2, Ip3=Ip3, rotx=rotx, roty=roty, rotz=rotz, - J2=J2, J4=J4, t=t) - if dsnew is not None: - self.ds = xr.combine_by_coords([self.ds, dsnew]) - self.ds['ntp'] = self.ds['id'].where(np.isnan(self.ds['Gmass'])).count(dim="id") - self.ds['npl'] = self.ds['id'].where(np.invert(np.isnan(self.ds['Gmass']))).count(dim="id") - 1 + J2=J2, J4=J4,t=t) + + dsnew = self._combine_and_fix_dsnew(dsnew) + + return dsnew + + def _combine_and_fix_dsnew(self,dsnew): + """ + Combines the new Dataset with the old one. Also computes the values of ntp and npl and sets the proper types. + Parameters + ---------- + dsnew : xarray Dataset + Dataset with new bodies + + Returns + ------- + dsnew : xarray Dataset + Updated Dataset with ntp, npl values and types fixed. + + """ + + self.ds = xr.combine_by_coords([self.ds, dsnew]) + + def get_nvals(ds): + if "Gmass" in dsnew: + ds['ntp'] = ds['id'].where(np.isnan(ds['Gmass'])).count(dim="id") + ds['npl'] = ds['id'].where(np.invert(np.isnan(ds['Gmass']))).count(dim="id") - 1 + else: + ds['ntp'] = ds['id'].count(dim="id") + ds['npl'] = xr.full_like(ds['ntp'],0) + return ds + + dsnew = get_nvals(dsnew) + self.ds = get_nvals(self.ds) if self.param['OUT_TYPE'] == "NETCDF_DOUBLE": + dsnew = io.fix_types(dsnew, ftype=np.float64) self.ds = io.fix_types(self.ds, ftype=np.float64) elif self.param['OUT_TYPE'] == "NETCDF_FLOAT": + dsnew = io.fix_types(dsnew, ftype=np.float32) self.ds = io.fix_types(self.ds, ftype=np.float32) return dsnew From 4bd9b6a30f099249a5ae37f0753090bf65d114f0 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 10 Nov 2022 20:37:38 -0500 Subject: [PATCH 5/7] Fixed typo that gave wrong npl and ntp values --- python/swiftest/swiftest/simulation_class.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 8e1ac24e0..d98a4ec1b 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -1864,9 +1864,9 @@ def _combine_and_fix_dsnew(self,dsnew): self.ds = xr.combine_by_coords([self.ds, dsnew]) def get_nvals(ds): - if "Gmass" in dsnew: + if "Gmass" in ds: ds['ntp'] = ds['id'].where(np.isnan(ds['Gmass'])).count(dim="id") - ds['npl'] = ds['id'].where(np.invert(np.isnan(ds['Gmass']))).count(dim="id") - 1 + ds['npl'] = ds['id'].where(~(np.isnan(ds['Gmass']))).count(dim="id") - 1 else: ds['ntp'] = ds['id'].count(dim="id") ds['npl'] = xr.full_like(ds['ntp'],0) From ff6da0f6188a8222e774060ebdd814a6a15ba531 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 10 Nov 2022 20:38:49 -0500 Subject: [PATCH 6/7] Updated example for the Basic Simulation --- examples/Basic_Simulation/cb.in | 9 - .../Basic_Simulation/initial_conditions.py | 34 ++-- examples/Basic_Simulation/param.in | 37 ++++ examples/Basic_Simulation/pl.in | 163 +++++++++--------- examples/Basic_Simulation/tp.in | 49 +++--- 5 files changed, 148 insertions(+), 144 deletions(-) create mode 100644 examples/Basic_Simulation/param.in diff --git a/examples/Basic_Simulation/cb.in b/examples/Basic_Simulation/cb.in index 29c8a7fca..b2cb85c35 100644 --- a/examples/Basic_Simulation/cb.in +++ b/examples/Basic_Simulation/cb.in @@ -1,12 +1,3 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - Sun 39.476926408897626 0.004650467260962157 diff --git a/examples/Basic_Simulation/initial_conditions.py b/examples/Basic_Simulation/initial_conditions.py index c9d0823af..f12bf6e07 100644 --- a/examples/Basic_Simulation/initial_conditions.py +++ b/examples/Basic_Simulation/initial_conditions.py @@ -30,31 +30,20 @@ from numpy.random import default_rng # Initialize the simulation object as a variable -sim = swiftest.Simulation(init_cond_file_type="ASCII") - +sim = swiftest.Simulation() sim.set_simulation_time(tstart=0.0, tstop=10.0, dt=0.005, tstep_out=1.0) - # Add parameter attributes to the simulation object sim.param['GMTINY'] = 1e-6 sim.param['MIN_GMFRAG'] = 1e-9 # Add the modern planets and the Sun using the JPL Horizons Database -sim.add("Sun", idval=0, date="2022-08-08") -sim.add("Mercury", idval=1, date="2022-08-08") -sim.add("Venus", idval=2, date="2022-08-08") -sim.add("Earth", idval=3, date="2022-08-08") -sim.add("Mars", idval=4, date="2022-08-08") -sim.add("Jupiter", idval=5, date="2022-08-08") -sim.add("Saturn", idval=6, date="2022-08-08") -sim.add("Uranus", idval=7, date="2022-08-08") -sim.add("Neptune", idval=8, date="2022-08-08") +sim.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune","Pluto"]) # Add 5 user-defined massive bodies npl = 5 density_pl = 3000.0 / (sim.param['MU2KG'] / sim.param['DU2M'] ** 3) -id_pl = np.array([9, 10, 11, 12, 13]) -name_pl = np.array(["MassiveBody_01", "MassiveBody_02", "MassiveBody_03", "MassiveBody_04", "MassiveBody_05"]) +name_pl = ["MassiveBody_01", "MassiveBody_02", "MassiveBody_03", "MassiveBody_04", "MassiveBody_05"] a_pl = default_rng().uniform(0.3, 1.5, npl) e_pl = default_rng().uniform(0.0, 0.3, npl) inc_pl = default_rng().uniform(0.0, 90, npl) @@ -64,19 +53,18 @@ GM_pl = (np.array([6e23, 8e23, 1e24, 3e24, 5e24]) / sim.param['MU2KG']) * sim.GU R_pl = np.full(npl, (3 * (GM_pl / sim.GU) / (4 * np.pi * density_pl)) ** (1.0 / 3.0)) Rh_pl = a_pl * ((GM_pl) / (3 * sim.GU)) ** (1.0 / 3.0) -Ip1_pl = np.array([0.4, 0.4, 0.4, 0.4, 0.4]) -Ip2_pl = np.array([0.4, 0.4, 0.4, 0.4, 0.4]) -Ip3_pl = np.array([0.4, 0.4, 0.4, 0.4, 0.4]) -rotx_pl = np.array([0.0, 0.0, 0.0, 0.0, 0.0]) -roty_pl = np.array([0.0, 0.0, 0.0, 0.0, 0.0]) -rotz_pl = np.array([0.0, 0.0, 0.0, 0.0, 0.0]) +Ip1_pl = [0.4, 0.4, 0.4, 0.4, 0.4] +Ip2_pl = [0.4, 0.4, 0.4, 0.4, 0.4] +Ip3_pl = [0.4, 0.4, 0.4, 0.4, 0.4] +rotx_pl = [0.0, 0.0, 0.0, 0.0, 0.0] +roty_pl = [0.0, 0.0, 0.0, 0.0, 0.0] +rotz_pl = [0.0, 0.0, 0.0, 0.0, 0.0] -sim.addp(id_pl, name_pl, a_pl, e_pl, inc_pl, capom_pl, omega_pl, capm_pl, GMpl=GM_pl, Rpl=R_pl, rhill=Rh_pl, Ip1=Ip1_pl, Ip2=Ip2_pl, Ip3=Ip3_pl, rotx=rotx_pl, roty=roty_pl, rotz=rotz_pl) +sim.add_body(name_pl, a_pl, e_pl, inc_pl, capom_pl, omega_pl, capm_pl, GMpl=GM_pl, Rpl=R_pl, rhill=Rh_pl, Ip1=Ip1_pl, Ip2=Ip2_pl, Ip3=Ip3_pl, rotx=rotx_pl, roty=roty_pl, rotz=rotz_pl) # Add 10 user-defined test particles ntp = 10 -id_tp = np.array([14, 15, 16, 17, 18, 19, 20, 21, 22, 23]) name_tp = np.array(["TestParticle_01", "TestParticle_02", "TestParticle_03", "TestParticle_04", "TestParticle_05", "TestParticle_06", "TestParticle_07", "TestParticle_08", "TestParticle_09", "TestParticle_10"]) a_tp = default_rng().uniform(0.3, 1.5, ntp) e_tp = default_rng().uniform(0.0, 0.3, ntp) @@ -85,7 +73,7 @@ omega_tp = default_rng().uniform(0.0, 360.0, ntp) capm_tp = default_rng().uniform(0.0, 360.0, ntp) -sim.addp(id_tp, name_tp, a_tp, e_tp, inc_tp, capom_tp, omega_tp, capm_tp) +sim.add_body(name_tp, a_tp, e_tp, inc_tp, capom_tp, omega_tp, capm_tp) # Save everything to a set of initial conditions files sim.save('param.in') diff --git a/examples/Basic_Simulation/param.in b/examples/Basic_Simulation/param.in new file mode 100644 index 000000000..47498726c --- /dev/null +++ b/examples/Basic_Simulation/param.in @@ -0,0 +1,37 @@ +! VERSION Swiftest parameter input +T0 0.0 +TSTOP 10.0 +DT 0.005 +ISTEP_OUT 200 +ISTEP_DUMP 200 +OUT_FORM XVEL +OUT_TYPE NETCDF_DOUBLE +OUT_STAT REPLACE +IN_TYPE NETCDF_DOUBLE +BIN_OUT bin.nc +CHK_QMIN 0.004650467260962157 +CHK_RMIN 0.004650467260962157 +CHK_RMAX 10000.0 +CHK_EJECT 10000.0 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 0.004650467260962157 10000.0 +MU2KG 1.988409870698051e+30 +TU2S 31557600.0 +DU2M 149597870700.0 +RESTART NO +INTERACTION_LOOPS TRIANGULAR +ENCOUNTER_CHECK TRIANGULAR +CHK_CLOSE YES +GR YES +FRAGMENTATION YES +ROTATION YES +ENERGY NO +EXTRA_FORCE NO +BIG_DISCARD NO +RHILL_PRESENT NO +TIDES NO +IN_FORM EL +NC_IN init_cond.nc +TSTART 0.0 +GMTINY 1e-06 +MIN_GMFRAG 1e-09 diff --git a/examples/Basic_Simulation/pl.in b/examples/Basic_Simulation/pl.in index 3f18aca0e..215b91f1d 100644 --- a/examples/Basic_Simulation/pl.in +++ b/examples/Basic_Simulation/pl.in @@ -1,88 +1,85 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -13 -Mercury 6.553709809565314146e-06 0.0014751262621647182575 -1.6306381826061645943e-05 -0.38709864823618972407 0.20562513690973019398 7.0036250456070128223 -48.30204974520415817 29.18823342267911869 114.9720047697775982 -0.0 0.0 0.34599999999999997424 -3.5734889863322150192 -18.38008501876561206 34.361513668512199956 -Venus 9.6633133995815381836e-05 0.006759085242739840662 -4.0453784346544178454e-05 -0.72332603580811538624 0.0067486079191121668003 3.394406261149808035 -76.61738837179606776 54.777760273590551776 315.37095689555837907 -0.0 0.0 0.4000000000000000222 -0.17650282045605921225 -3.6612475825356215592 8.702866268072763821 -Earth 0.000120026935827952456416 0.010044668314295209318 -4.25875607065040958e-05 -0.9999943732711822353 0.016707309394717231171 0.0028984767436730000077 -174.05498211951208987 289.04709044403989537 213.07530468023790604 -0.0 0.0 0.33069999999999999396 -5.002093202481912218 0.055213529850334066125 2301.2110537292529557 -Mars 1.2739802010675941808e-05 0.007246950762048707243 -2.265740805092889601e-05 -1.5237812078483019551 0.0935087708803710449 1.8479353068000929916 -49.489305419773351957 286.70300191753761965 24.878418068365981242 -0.0 0.0 0.3644000000000000017 -997.9357048213454125 -909.4072592492943007 1783.4501726537997323 -Jupiter 0.03769225108898567778 0.3552222491747608486 -0.00046732617030490929307 -5.2028063728088866924 0.048395118271449058533 1.3035670146561249005 -100.516498130230701236 273.44233262595901124 346.26538105843917492 -0.0 0.0 0.27560000000000001164 --80.96619889339111482 -2388.0060524649362916 5008.7314931237953832 -Saturn 0.01128589982009127331 0.43757948578866074266 -0.00038925687730393611812 -9.580020069168169172 0.053193613750490406633 2.4864365613724639381 -113.597044717589099605 335.10179422401358806 237.66485199561481068 -0.0 0.0 0.22000000000000000111 -441.93538182505989814 378.5284220382117538 5135.9110455622733884 -Uranus 0.001723658947826773068 0.4705353566089944894 -0.00016953449859497231466 -19.272143108769419939 0.043779687288749750962 0.7707536154556786645 -74.077748995180698444 93.42271392662131291 242.37685081109759722 -0.0 0.0 0.23000000000000000999 --677.3000258209181323 -3008.109907190578637 -836.301326618569835 -Neptune 0.0020336100526728302882 0.78184929587893868845 -0.000164587904124493665 -30.305539399096510067 0.014544938874222059638 1.7686697746048700708 -131.73604731224671127 249.9779420269553043 332.54824537252648042 -0.0 0.0 0.23000000000000000999 -1231.1804455066093229 -2178.0887091151860042 2329.6411363603121418 -MassiveBody_01 1.1912109366578087428e-05 0.0016092923734511708263 -2.425055692051244981e-05 -0.3460404950890429432 0.2093906512182220625 0.11109012870384793459 -114.31328763688792094 347.82259114762894114 96.0534391561842682 -0.4000000000000000222 0.4000000000000000222 0.4000000000000000222 +14 +Mercury 6.553709809565314e-06 +1.6306381826061646e-05 +0.3870985843095394 0.2056234010897001 7.003302508001384 +48.29611837378607 29.20442403952454 338.3394874682879 +0.0 0.0 0.346 +3.573018077015318 -18.380317085416586 34.36143850429876 +Venus 9.663313399581537e-05 +4.0453784346544176e-05 +0.7233297579736101 0.006717605698865438 3.394439273342282 +76.60235891771119 54.96037946082961 200.4789339550648 +0.0 0.0 0.4 +0.17650282045605922 -3.6612475825356214 8.702866268072764 +Earth 0.00012002693582795245 +4.25875607065041e-05 +0.9999904874543223 0.01671400376545858 0.003637862608863003 +175.025172600231 287.9619628812575 114.3482934042427 +0.0 0.0 0.3307 +6.157239621449141 0.057224133705898524 2301.2082528350797 +Mars 1.2739802010675942e-05 +2.2657408050928896e-05 +1.523711925589535 0.09344151133508208 1.847441673557901 +49.4728572124747 286.7379771285891 209.3396773477138 +0.0 0.0 0.3644 +997.9224351226384 -909.5549030011778 1783.3823046046184 +Jupiter 0.037692251088985676 +0.0004673261703049093 +5.2027278008516 0.04824497711637968 1.303631134570075 +100.5192588433081 273.5898402882514 129.5536700659942 +0.0 0.0 0.2756 +-80.9864396731672 -2388.0246092955053 5008.722318533006 +Saturn 0.011285899820091273 +0.00038925687730393614 +9.532011952667288 0.05486329870433341 2.487906363280301 +113.6305781676206 339.5467356402391 290.8995806568904 +0.0 0.0 0.22 +441.95954822014636 378.52638822638795 5135.909115928892 +Uranus 0.001723658947826773 +0.00016953449859497232 +19.24498838290236 0.04796174942301296 0.7730102596086205 +74.0125809801658 93.59554912280227 262.8658637277515 +0.0 0.0 0.23 +-677.3000258209181 -3008.1099071905787 -836.3013266185699 +Neptune 0.0020336100526728304 +0.00016458790412449367 +30.03895991152209 0.008955570138096731 1.771119354296142 +131.8221159748827 284.4748429674216 308.4513720536233 +0.0 0.0 0.23 +1232.224106980634 -2177.3040821077648 2329.8227878119233 +Pluto 2.924216771029454e-07 +7.943294877391593e-06 +39.36791814672583 0.2487178537481577 17.1705505990969 +110.3314332962701 113.0826635900664 55.11416408345664 +0.0 0.0 0.4 +-243.59404988903637 261.28663002814227 -38.57352022187049 +MassiveBody_01 1.1912109366578089e-05 +2.425055692051245e-05 +0.7452368716298337 0.0633011418780484 0.11151363780595558 +203.01823417718037 284.9353898127118 266.79344592519305 +0.4 0.4 0.4 0.0 0.0 0.0 -MassiveBody_02 1.5882812488770779849e-05 0.0016802531895603555184 -2.6691191565570073646e-05 -0.32826188156947710972 0.27866488696288682636 77.21223337306255985 -251.99014895640269174 53.772702227560969845 165.6085387284213084 -0.4000000000000000222 0.4000000000000000222 0.4000000000000000222 +MassiveBody_02 1.5882812488770783e-05 +2.6691191565570074e-05 +0.7671203280602826 0.10149029388964753 53.46706656938751 +61.74738152068808 68.20565593722856 271.3352706902475 +0.4 0.4 0.4 0.0 0.0 0.0 -MassiveBody_03 1.9853515610963475658e-05 0.004324919007577881884 -2.8752214513575297366e-05 -0.7843689028022314824 0.06176128116947356833 78.74144231136708072 -286.63765100951468412 347.55933571120488068 266.36960496595537506 -0.4000000000000000222 0.4000000000000000222 0.4000000000000000222 +MassiveBody_03 1.985351561096348e-05 +2.8752214513575297e-05 +1.4698824276418418 0.13621250684495437 23.635498327264845 +38.071905339231236 283.134612455057 250.67457601578352 +0.4 0.4 0.4 0.0 0.0 0.0 -MassiveBody_04 5.9560546832890430362e-05 0.0056765613035874530265 -4.146786902759040254e-05 -0.7138176832994267418 0.28016098557400787028 22.725690778108500467 -203.41845532080247949 219.74297850728484605 14.730732982803269593 -0.4000000000000000222 0.4000000000000000222 0.4000000000000000222 +MassiveBody_04 5.9560546832890444e-05 +4.14678690275904e-05 +0.9741731590760117 0.1519713326893784 20.51335588582416 +350.53780805624825 304.05941264938997 142.62713592644738 +0.4 0.4 0.4 0.0 0.0 0.0 -MassiveBody_05 9.9267578054817388455e-05 0.010382929139161686458 -4.916559523190238318e-05 -1.101215402063684401 0.076651567404004070094 52.41961577462824806 -142.90070862650665617 293.70542448390904156 318.5666754758643151 -0.4000000000000000222 0.4000000000000000222 0.4000000000000000222 +MassiveBody_05 9.926757805481742e-05 +4.916559523190238e-05 +0.6633500122731075 0.13550930562584215 10.680323453316653 +328.45970148163457 49.93948991697533 316.2109831817007 +0.4 0.4 0.4 0.0 0.0 0.0 diff --git a/examples/Basic_Simulation/tp.in b/examples/Basic_Simulation/tp.in index 162f8c0af..9e3c6b9bc 100644 --- a/examples/Basic_Simulation/tp.in +++ b/examples/Basic_Simulation/tp.in @@ -1,40 +1,31 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - 10 TestParticle_01 -1.0288790290699558749 0.088535376967516773994 39.679062610233010844 -51.107327480220099858 20.90477961082723013 183.92936958121950397 +0.5697767200549144 0.11834332057809717 80.87229166892182 +222.5373287698991 111.76944690551065 270.3754252090885 TestParticle_02 -0.54831541717225318333 0.20048853000030755767 32.784952066779297297 -66.59436607119042151 196.2609764918545352 76.40623742948525887 +1.0092273242276735 0.05399051523964523 82.35601343884242 +117.598050656921 116.26893959032198 208.13367102381633 TestParticle_03 -0.9815578696795972391 0.25717438840188583393 38.959194313128776344 -158.31201212846775661 93.86793546512863884 126.96040079919036714 +0.7590606700221183 0.027936243833100036 88.12926503625694 +188.5414923780248 169.44578516513008 282.4333664924941 TestParticle_04 -0.70510786000431613374 0.068740260615181125736 39.981235917453354034 -26.668314027440779057 12.507089902982141183 53.606668142734086757 +0.9972940668463142 0.2569338531459909 22.75744046640894 +199.388648818177 152.8772634782716 307.54746324357103 TestParticle_05 -1.336737159289705712 0.25044351028750111432 74.189544264066626056 -313.1647935522670423 152.39951433565357775 322.52271715518395467 +0.341620427452812 0.052524493325270837 18.839916665085173 +327.2732232808312 334.4396604858765 359.4553699087638 TestParticle_06 -0.64018428687513928566 0.1478972702874425671 73.39555666508663023 -74.03379498826825511 124.22185942531125136 106.36095293685497154 +1.1670239341669624 0.07778574626824612 22.716135789279228 +196.30730270195494 154.1907148406917 159.6860762761786 TestParticle_07 -1.2738161048947760356 0.17129911220230967239 78.723790909435408025 -111.971184037421835455 315.6294407604535195 108.74538288631850946 +1.0954005977667731 0.10128096358705894 9.878980088486013 +101.85525076444998 124.8686145955563 254.7757898831765 TestParticle_08 -0.5196197141250760154 0.26607581023277782073 79.221465395061358095 -240.07768052067768849 177.4793327047061382 189.85775920180086018 +0.39125485873865873 0.17327367123519463 22.276414655828646 +105.99598258123197 235.2803528505754 348.706620365582 TestParticle_09 -0.7160151851892424535 0.29584187625470087513 58.1214116614385361 -24.22606869850931588 338.4678256522021229 136.32070882113120547 +0.48720268138208256 0.29707768397862033 39.717820434869466 +307.9106242922398 221.91704762104942 333.56626647648994 TestParticle_10 -1.0782431797945351004 0.11096856177737297877 54.729767236739554903 -280.85362091874827684 116.780853088052467115 286.99855363222661708 +0.6458552755486036 0.07522818043309405 45.225887389492826 +74.40348173261458 92.77221707157337 260.30279498657796 From 0be34922757a59eb9d0f20166bac6b00097acf27 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 10 Nov 2022 20:41:02 -0500 Subject: [PATCH 7/7] switched the test particle names in the Basic_Simulation example to be a list rather than a numpy array just to make it cleaner (it accepts either) --- examples/Basic_Simulation/initial_conditions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/Basic_Simulation/initial_conditions.py b/examples/Basic_Simulation/initial_conditions.py index f12bf6e07..1edc13879 100644 --- a/examples/Basic_Simulation/initial_conditions.py +++ b/examples/Basic_Simulation/initial_conditions.py @@ -65,7 +65,7 @@ # Add 10 user-defined test particles ntp = 10 -name_tp = np.array(["TestParticle_01", "TestParticle_02", "TestParticle_03", "TestParticle_04", "TestParticle_05", "TestParticle_06", "TestParticle_07", "TestParticle_08", "TestParticle_09", "TestParticle_10"]) +name_tp = ["TestParticle_01", "TestParticle_02", "TestParticle_03", "TestParticle_04", "TestParticle_05", "TestParticle_06", "TestParticle_07", "TestParticle_08", "TestParticle_09", "TestParticle_10"] a_tp = default_rng().uniform(0.3, 1.5, ntp) e_tp = default_rng().uniform(0.0, 0.3, ntp) inc_tp = default_rng().uniform(0.0, 90, ntp)