diff --git a/CMakeLists.txt b/CMakeLists.txt index 776f386e4..34b05b21e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,3 +1,12 @@ +# 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. + # CMake project file for FOO ################################################## diff --git a/cmake/Modules/FindNETCDF.cmake b/cmake/Modules/FindNETCDF.cmake index f8e97dbca..2355ad900 100644 --- a/cmake/Modules/FindNETCDF.cmake +++ b/cmake/Modules/FindNETCDF.cmake @@ -1,3 +1,12 @@ +# 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. + # - Finds the NetCDF libraries find_path(NETCDF_INCLUDE_DIR NAMES netcdf.mod HINTS ENV NETCDF_FORTRAN_HOME) find_library(NETCDF_LIBRARY NAMES netcdf HINTS ENV NETCDF_FORTRAN_HOME) diff --git a/cmake/Modules/FindOpenMP_Fortran.cmake b/cmake/Modules/FindOpenMP_Fortran.cmake index bc440ae10..c8e0ca2b4 100644 --- a/cmake/Modules/FindOpenMP_Fortran.cmake +++ b/cmake/Modules/FindOpenMP_Fortran.cmake @@ -1,3 +1,12 @@ +# 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. + # - Finds OpenMP support # This module can be used to detect OpenMP support in a compiler. # If the compiler supports OpenMP, the flags required to compile with @@ -13,18 +22,6 @@ # OMP_NUM_PROCS - the max number of processors available to OpenMP #============================================================================= -# Copyright 2009 Kitware, Inc. -# Copyright 2008-2009 André Rigland Brodtkorb -# -# Distributed under the OSI-approved BSD License (the "License"); -# see accompanying file Copyright.txt for details. -# -# This software is distributed WITHOUT ANY WARRANTY; without even the -# implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. -# See the License for more information. -#============================================================================= -# (To distribute this file outside of CMake, substitute the full -# License text for the above reference.) INCLUDE (${CMAKE_ROOT}/Modules/FindPackageHandleStandardArgs.cmake) diff --git a/cmake/Modules/SetCompileFlag.cmake b/cmake/Modules/SetCompileFlag.cmake index 04ff3ffbd..1d110ae6d 100644 --- a/cmake/Modules/SetCompileFlag.cmake +++ b/cmake/Modules/SetCompileFlag.cmake @@ -1,3 +1,12 @@ +# 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. + ############################################################################# # Given a list of flags, this function will try each, one at a time, # and choose the first flag that works. If no flags work, then nothing diff --git a/cmake/Modules/SetFortranFlags.cmake b/cmake/Modules/SetFortranFlags.cmake index 926547bbb..e0b21862b 100644 --- a/cmake/Modules/SetFortranFlags.cmake +++ b/cmake/Modules/SetFortranFlags.cmake @@ -1,3 +1,12 @@ +# 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. + ###################################################### # Determine and set the Fortran compiler flags we want ###################################################### diff --git a/cmake/Modules/SetParallelizationLibrary.cmake b/cmake/Modules/SetParallelizationLibrary.cmake index 603d4299c..03ab970e6 100644 --- a/cmake/Modules/SetParallelizationLibrary.cmake +++ b/cmake/Modules/SetParallelizationLibrary.cmake @@ -1,3 +1,12 @@ +# 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. + # Turns on either OpenMP or MPI # If both are requested, the other is disabled # When one is turned on, the other is turned off diff --git a/cmake/Modules/SetUpNetCDF.cmake b/cmake/Modules/SetUpNetCDF.cmake index dbfe44c4b..48319ed9b 100644 --- a/cmake/Modules/SetUpNetCDF.cmake +++ b/cmake/Modules/SetUpNetCDF.cmake @@ -1,3 +1,12 @@ +# 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. + # Find NetCDF if not already found IF(NOT NETCDF_FOUND) ENABLE_LANGUAGE(C) # Some libraries need a C compiler to find diff --git a/distclean.cmake b/distclean.cmake index 8e24f9e49..00d2e454c 100644 --- a/distclean.cmake +++ b/distclean.cmake @@ -1,3 +1,12 @@ +# 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. + # This CMake script will delete build directories and files to bring the # package back to it's distribution state diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py index 7d7c20e14..78ef51be9 100644 --- a/python/swiftest/swiftest/init_cond.py +++ b/python/swiftest/swiftest/init_cond.py @@ -400,7 +400,7 @@ def vec2xr(param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, frame_str = info_vec_str.T if param['IN_TYPE'] == 'NETCDF_FLOAT': ftype=np.float32 - elif param['IN_TYPE'] == 'NETCDF_DOUBLE': + 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}) diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 176237b0a..4678b4784 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -999,7 +999,9 @@ def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,fram ------- A set of input files for a new Swiftest run """ - frame = select_active_from_frame(ds, param, framenum) + param_tmp = param.copy() + param_tmp['OUT_FORM'] = param['IN_FORM'] + frame = select_active_from_frame(ds, param_tmp, framenum) if in_type == "NETCDF_DOUBLE" or in_type == "NETCDF_FLOAT": # Convert strings back to byte form and save the NetCDF file @@ -1051,22 +1053,22 @@ def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,fram for i in pl.id: pli = pl.sel(id=i) if param['RHILL_PRESENT'] == 'YES': - print(pli['name'].values, pli['Gmass'].values, pli['rhill'].values, file=plfile) + print(pli['name'].values[0], pli['Gmass'].values[0], pli['rhill'].values[0], file=plfile) else: - print(pli['name'].values, pli['Gmass'].values, file=plfile) + print(pli['name'].values[0], pli['Gmass'].values[0], file=plfile) if param['CHK_CLOSE'] == 'YES': - print(pli['radius'].values, file=plfile) + print(pli['radius'].values[0], file=plfile) if param['IN_FORM'] == 'XV': - print(pli['xhx'].values, pli['xhy'].values, pli['xhz'].values, file=plfile) - print(pli['vhx'].values, pli['vhy'].values, pli['vhz'].values, file=plfile) + print(pli['xhx'].values[0], pli['xhy'].values[0], pli['xhz'].values[0], file=plfile) + print(pli['vhx'].values[0], pli['vhy'].values[0], pli['vhz'].values[0], file=plfile) elif param['IN_FORM'] == 'EL': - print(pli['a'].values, pli['e'].values, pli['inc'].values, file=plfile) - print(pli['capom'].values, pli['omega'].values, pli['capm'].values, file=plfile) + print(pli['a'].values[0], pli['e'].values[0], pli['inc'].values[0], file=plfile) + print(pli['capom'].values[0], pli['omega'].values[0], pli['capm'].values[0], file=plfile) else: print(f"{param['IN_FORM']} is not a valid input format type.") if param['ROTATION'] == 'YES': - print(pli['Ip1'].values, pli['Ip2'].values, pli['Ip3'].values, file=plfile) - print(pli['rotx'].values, pli['roty'].values, pli['rotz'].values, file=plfile) + print(pli['Ip1'].values[0], pli['Ip2'].values[0], pli['Ip3'].values[0], file=plfile) + print(pli['rotx'].values[0], pli['roty'].values[0], pli['rotz'].values[0], file=plfile) plfile.close() # TP file @@ -1074,105 +1076,16 @@ def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,fram print(tp.id.count().values, file=tpfile) for i in tp.id: tpi = tp.sel(id=i) - print(tpi['name'].values, file=tpfile) + print(tpi['name'].values[0], file=tpfile) if param['IN_FORM'] == 'XV': - print(tpi['xhx'].values, tpi['xhy'].values, tpi['xhz'].values, file=tpfile) - print(tpi['vhx'].values, tpi['vhy'].values, tpi['vhz'].values, file=tpfile) + print(tpi['xhx'].values[0], tpi['xhy'].values[0], tpi['xhz'].values[0], file=tpfile) + print(tpi['vhx'].values[0], tpi['vhy'].values[0], tpi['vhz'].values[0], file=tpfile) elif param['IN_FORM'] == 'EL': - print(tpi['a'].values, tpi['e'].values, tpi['inc'].values, file=tpfile) - print(tpi['capom'].values, tpi['omega'].values, tpi['capm'].values, file=tpfile) + print(tpi['a'].values[0], tpi['e'].values[0], tpi['inc'].values[0], file=tpfile) + print(tpi['capom'].values[0], tpi['omega'].values[0], tpi['capm'].values[0], file=tpfile) else: print(f"{param['IN_FORM']} is not a valid input format type.") tpfile.close() - elif in_type == 'REAL8': - # Now make Swiftest files - cbfile = FortranFile(param['CB_IN'], 'w') - cbfile.write_record(cbid) - cbfile.write_record(np.double(GMSun)) - cbfile.write_record(np.double(RSun)) - cbfile.write_record(np.double(J2)) - cbfile.write_record(np.double(J4)) - if param['ROTATION'] == 'YES': - cbfile.write_record(np.double(Ip1cb)) - cbfile.write_record(np.double(Ip2cb)) - cbfile.write_record(np.double(Ip3cb)) - cbfile.write_record(np.double(rotxcb)) - cbfile.write_record(np.double(rotycb)) - cbfile.write_record(np.double(rotzcb)) - - cbfile.close() - - plfile = FortranFile(param['PL_IN'], 'w') - npl = pl.id.count().values - plid = pl.id.values - if param['IN_FORM'] == 'XV': - v1 = pl['xhx'].values - v2 = pl['xhy'].values - v3 = pl['xhz'].values - v4 = pl['vhx'].values - v5 = pl['vhy'].values - v6 = pl['vhz'].values - elif param['IN_FORM'] == 'EL': - v1 = pl['a'].values - v2 = pl['e'].values - v3 = pl['inc'].values - v4 = pl['capom'].values - v5 = pl['omega'].values - v6 = pl['capm'].values - else: - print(f"{param['IN_FORM']} is not a valid input format type.") - Gmass = pl['Gmass'].values - if param['CHK_CLOSE'] == 'YES': - radius = pl['radius'].values - - plfile.write_record(npl) - plfile.write_record(plid) - plfile.write_record(v1) - plfile.write_record(v2) - plfile.write_record(v3) - plfile.write_record(v4) - plfile.write_record(v5) - plfile.write_record(v6) - plfile.write_record(Gmass) - if param['RHILL_PRESENT'] == 'YES': - plfile.write_record(pl['rhill'].values) - if param['CHK_CLOSE'] == 'YES': - plfile.write_record(radius) - if param['ROTATION'] == 'YES': - plfile.write_record(pl['Ip1'].values) - plfile.write_record(pl['Ip2'].values) - plfile.write_record(pl['Ip3'].values) - plfile.write_record(pl['rotx'].values) - plfile.write_record(pl['roty'].values) - plfile.write_record(pl['rotz'].values) - plfile.close() - tpfile = FortranFile(param['TP_IN'], 'w') - ntp = tp.id.count().values - tpid = tp.id.values - if param['IN_FORM'] == 'XV': - v1 = tp['xhx'].values - v2 = tp['xhy'].values - v3 = tp['xhz'].values - v4 = tp['vhx'].values - v5 = tp['vhy'].values - v6 = tp['vhz'].values - elif param['IN_FORM'] == 'EL': - v1 = tp['a'].values - v2 = tp['e'].values - v3 = tp['inc'].values - v4 = tp['capom'].values - v5 = tp['omega'].values - v6 = tp['capm'].values - else: - print(f"{param['IN_FORM']} is not a valid input format type.") - tpfile.write_record(ntp) - tpfile.write_record(tpid) - tpfile.write_record(v1) - tpfile.write_record(v2) - tpfile.write_record(v3) - tpfile.write_record(v4) - tpfile.write_record(v5) - tpfile.write_record(v6) else: print(f"{in_type} is an unknown file type") diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index cc4229a23..d1fa28ede 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -33,6 +33,9 @@ def __init__(self, codename="Swiftest", param_file="param.in", readbin=True, ver 'IN_FORM': "XV", 'IN_TYPE': "NETCDF_DOUBLE", 'NC_IN' : "init_cond.nc", + 'CB_IN' : "cb.in", + 'PL_IN' : "pl.in", + 'TP_IN' : "tp.in", 'ISTEP_OUT': "1", 'ISTEP_DUMP': "1", 'BIN_OUT': "bin.nc", @@ -191,6 +194,12 @@ def write_param(self, param_file, param=None): # Check to see if the parameter type matches the output type. If not, we need to convert codename = param['! VERSION'].split()[0] if codename == "Swifter" or codename == "Swiftest": + 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) io.write_labeled_param(param, param_file) elif codename == "Swift": io.write_swift_param(param, param_file) @@ -335,7 +344,7 @@ def save(self, param_file, framenum=-1, codename="Swiftest"): """ if codename == "Swiftest": - io.swiftest_xr2infile(ds=self.ds, param=self.param, framenum=framenum,infile_name=self.param['NC_IN']) + io.swiftest_xr2infile(ds=self.ds, param=self.param, in_type=self.param['IN_TYPE'], framenum=framenum,infile_name=self.param['NC_IN']) self.write_param(param_file) elif codename == "Swifter": if self.codename == "Swiftest": diff --git a/python/swiftest/swiftest/tool.py b/python/swiftest/swiftest/tool.py index 91f8c5a46..07e447ab4 100644 --- a/python/swiftest/swiftest/tool.py +++ b/python/swiftest/swiftest/tool.py @@ -100,4 +100,394 @@ def follow_swift(ds, ifol=None, nskp=None): except IOError: print(f"Error writing to follow.out") - return fol \ No newline at end of file + return fol + +def danby(M, ecc, accuracy=1e-14): + """ + Danby's method to solve Kepler's equation. + + Parameters + ---------- + M : float + the mean anomaly in radians + ecc : float + the eccentricity + accuracy : float + the relative accuracy to obtain a solution. Default is 1e-12. + + Returns + ---------- + E : float + the eccentric anomaly in radians + + References + __________ + Danby, J.M.A. 1988. Fundamentals of celestial mechanics. Richmond, Va., U.S.A., Willmann-Bell, 1988. 2nd ed. + Murray, C.D., Dermott, S.F., 1999. Solar system dynamics, New York, New York. ed, Cambridge University Press. + + """ + + def kepler_root(E, ecc, M, deriv): + """ + The Kepler equation root function. + + Parameters + ---------- + E : float + the eccentric anomaly in radians + ecc : float + the eccentricity + M : float + the mean anomaly in radians + deriv : int between 0 and 3 + The order of the derivative to compute + + Returns + ---------- + deriv = 0: E - e * np.sin(E) - M + deriv = 1: 1 - e * np.cos(E) + deriv = 2: e * np.sin(E) + deriv = 3: e * np.cos(E) + + Note: The function will return 0 when E is correct for a given e and M + """ + + if deriv == 0: + return E - ecc * np.sin(E) - M + elif deriv == 1: + return 1.0 - ecc * np.cos(E) + elif deriv == 2: + return ecc * np.sin(E) + elif deriv == 3: + return ecc * np.cos(E) + else: + print(f"deriv={deriv} is undefined") + return None + + def delta_ij(E, ecc, M, j): + """ + Danby's intermediate delta functions. This function is recursive for j>1. + + Parameters + ---------- + E : float + the eccentric anomaly in radians + ecc : float + the eccentricity + M : float + the mean anomaly in radians + j : int between 1 and 3 + The order of the delta function to compute + + Returns + ---------- + delta_ij value used in Danby's iterative Kepler equation solver + + """ + if j == 1: + return -kepler_root(E, ecc, M, 0) / kepler_root(E, ecc, M, 1) + elif j == 2: + return -kepler_root(E, ecc, M, 0) / (kepler_root(E, ecc, M, 1) + - delta_ij(E, ecc, M, 1) * kepler_root(E, ecc, M, 2) / 2.0) + elif j == 3: + return -kepler_root(E, ecc, M, 0) / (kepler_root(E, ecc, M, 1) + + delta_ij(E, ecc, M, 1) * kepler_root(E, ecc, M, 2) / 2.0 + + delta_ij(E, ecc, M, 2) ** 2 * kepler_root(E, ecc, M, + 3) / 6.0) + else: + print(f"j = {j} is not a valid input to the delta_ij function") + + # If this is a circular orbit, just return the mean anomaly as the eccentric anomaly + if ecc < np.finfo(np.float64).tiny: # Prevent floating point exception for 0 eccentricity orbits + return M + + # Initial guess + k = 0.85 + E = M + np.sign(np.sin(M)) * k * ecc + MAXLOOPS = 50 # Maximum number of times to iterate before we give up + for i in range(MAXLOOPS): + Enew = E + delta_ij(E, ecc, M, 3) + if np.abs((Enew - E) / E) < accuracy: + return Enew + E = Enew + + raise RuntimeError("The danby function did not converge on a solution.") + + + +def el2xv_one(mu, a, ecc, inc, Omega, omega, M): + """ + Compute osculating orbital elements from relative Cartesian position and velocity + All angular measures are returned in radians + If inclination < TINY, longitude of the ascending node is arbitrarily set to 0 + + If eccentricity < sqrt(TINY), argument of pericenter is arbitrarily set to 0 + + ALGORITHM: See Fitzpatrick "Principles of Cel. Mech." + + Adapted from Martin Duncan's el2xv.f + + Parameters + ---------- + mu : float + Central body gravitational constant + a : float + semimajor axis + ecc : float + eccentricity + inc : float + inclination (degrees) + Omega : float + longitude of ascending node (degrees) + omega : float + argument of periapsis (degrees) + M : flaot + Mean anomaly (degrees) + + Returns + ---------- + rvec, vvec : tuple of float vectors + + rvec : (3) float vector + Cartesian position vector + vvec : (3) float vector + Cartesian velocity vector + + """ + + if ecc < 0.0: + print("Error in el2xv! Eccentricity cannot be negative. Setting it to 0") + e = 0.0 + iorbit_type = "ellipse" + else: + e = ecc + em1 = e - 1. + if np.abs(em1) < 2 * np.finfo(float).eps: + iorbit_type = "parabola" + elif e > 1.0: + iorbit_type = "hyperbola" + else: + iorbit_type = "ellipse" + + def scget(angle): + """ + Efficiently compute the sine and cosine of an input angle + Input angle must be in radians + + Adapted from David E. Kaufmann's Swifter routine: orbel_scget.f90 + Adapted from Hal Levison's Swift routine orbel_scget.f + + Parameters + ---------- + angle : input angle + + Returns + ------- + sx : sin of angle + cx : cos of angle + + """ + TWOPI = 2 * np.pi + nper = int(angle / TWOPI) + x = angle - nper * TWOPI + if x < 0.0: + x += TWOPI + + sx = np.sin(x) + cx = np.sqrt(1.0 - sx ** 2) + if x > np.pi / 2.0 and x < 3 * np.pi / 2.0: + cx = -cx + return sx, cx + + sip, cip = scget(np.deg2rad(omega)) + so, co = scget(np.deg2rad(Omega)) + si, ci = scget(np.deg2rad(inc)) + + d11 = cip * co - sip * so * ci + d12 = cip * so + sip * co * ci + d13 = sip * si + d21 = -sip * co - cip * so * ci + d22 = -sip * so + cip * co * ci + d23 = cip * si + + # Get the other quantities depending on orbit type + if iorbit_type == "ellipse": + E = danby(np.deg2rad(M),e) + scap, ccap = scget(E) + sqe = np.sqrt(1. - e**2) + sqgma = np.sqrt(mu * a) + xfac1 = a * (ccap - e) + xfac2 = a * sqe * scap + ri = 1. / (a * (1. - e * ccap)) + vfac1 = -ri * sqgma * scap + vfac2 = ri * sqgma * sqe * ccap + else: + print(f"Orbit type: {iorbit_type} not yet implemented") + xfac1 = 0.0 + xfac2 = 0.0 + vfac1 = 0.0 + vfac2 = 0.0 + + rvec = np.array([d11 * xfac1 + d21 * xfac2, + d12 * xfac1 + d22 * xfac2, + d13 * xfac1 + d23 * xfac2]) + vvec = np.array([d11 * vfac1 + d21 * vfac2, + d12 * vfac1 + d22 * vfac2, + d13 * vfac1 + d23 * vfac2]) + + return rvec, vvec + + +def el2xv_vec(mu, a, ecc, inc, Omega, omega, M): + """ + + Vectorized call to el2xv_one + Parameters + ---------- + mu : float + Central body gravitational constant + a : (n) float array + semimajor axis + ecc : (n) float array + eccentricity + inc : (n) float array + inclination (degrees) + Omega : (n) float array + longitude of ascending node (degrees) + omega : (n) float array + argument of periapsis (degrees) + M : (n) flaot array + Mean anomaly (degrees) + + Returns + ---------- + rvec, vvec : tuple of float vectors + + rvec : (n,3) float rray + Cartesian position vector + vvec : (n,3) float array + Cartesian velocity vector + """ + vecfunc = np.vectorize(el2xv_one, signature='(),(),(),(),(),(),()->(3),(3)') + return vecfunc(mu, a, ecc, inc, Omega, omega, M) + +def xv2el_one(mu,rvec,vvec): + """ + Converts from cartesian position and velocity values to orbital elements + + Parameters + ---------- + mu : float + Central body gravitational constant + rvec : (3) float array + Cartesian position vector + vvec : (3) float array + Cartesian velocity vector + + Returns + ---------- + a, ecc, inc, Omega, omega, M, varpi, f, lam : tuple of floats + + a : float + semimajor axis + ecc : float + eccentricity + inc : float + inclination (degrees) + Omega : float + longitude of ascending node (degrees) + omega : float + argument of periapsis (degrees) + M : float + mean anomaly (degrees) + varpi : flaot + longitude of periapsis (degrees) + f : float + true anomaly (degrees) + lam : float + mean longitude (degrees) + + """ + + rmag = np.linalg.norm(rvec) + vmag2 = np.vdot(vvec,vvec) + h = np.cross(rvec,vvec) + hmag = np.linalg.norm(h) + + rdot = np.sign(np.vdot(rvec,vvec)) * np.sqrt(vmag2 - (hmag / rmag)**2) + + a = 1.0/(2.0 / rmag - vmag2/mu) + ecc = np.sqrt(1 - hmag**2 / (mu * a)) + inc = np.arccos(h[2]/hmag) + + goodinc = np.abs(inc) > np.finfo(np.float64).tiny + sO = np.where(goodinc, np.sign(h[2]) * h[0] / (hmag * np.sin(inc)),0.0) + cO = np.where(goodinc, -np.sign(h[2]) * h[1] / (hmag * np.sin(inc)),1.0) + + Omega = np.arctan2(sO, cO) + + sof = np.where(goodinc,rvec[2] / (rmag * np.sin(inc)), rvec[1]/rmag) + cof = np.where(goodinc,(rvec[0] / rmag + np.sin(Omega) * sof * np.cos(inc)) / np.cos(Omega), rvec[0]/rmag) + + of = np.arctan2(sof,cof) + + goodecc = ecc > np.finfo(np.float64).tiny + sf = np.where(goodecc,a * (1.0 - ecc**2) * rdot / (hmag * ecc), 0.0) + cf = np.where(goodecc,(1.0 / ecc) * (a * (1.0 - ecc**2) / rmag - 1.0), 1.0 ) + + f = np.arctan2(sf, cf) + + omega = of - f + + varpi = Omega + omega + + # Compute eccentric anomaly & mean anomaly in order to get mean longitude + E = np.where(ecc > np.finfo(np.float64).tiny, np.arccos(-(rmag - a) / (a * ecc)), 0) + if np.sign(np.vdot(rvec, vvec)) < 0.0: + E = 2 * np.pi - E + + M = E - ecc * np.sin(E) + lam = M + varpi + + return a, ecc, np.rad2deg(inc), np.rad2deg(Omega), np.rad2deg(omega), np.rad2deg(M), np.rad2deg(varpi), np.rad2deg(f), np.rad2deg(lam) + +def xv2el_vec(mu, rvec, vvec): + """ + Vectorized call to xv2el_one. + + Parameters + ---------- + mu : float + Central body gravitational constant + rvec : (n,3) float array + Cartesian position vector + vvec : (n,3) float array + Cartesian velocity vector + + Returns + ---------- + a, ecc, inc, Omega, omega, M, varpi, f, lam : tuple of float arrays + + a : (n) float array + semimajor axis + ecc : (n) float array + eccentricity + inc : (n) float array + inclination (degrees) + Omega : (n) float array + longitude of ascending node (degrees) + omega : (n) float array + argument of periapsis (degrees) + M : (n) float array + mean anomaly (degrees) + varpi : (n) flaot array + longitude of periapsis (degrees) + f : (n) float array + true anomaly (degrees) + lam : (n) float array + mean longitude (degrees) + + """ + + vecfunc = np.vectorize(xv2el_one, signature='(),(3),(3)->(),(),(),(),(),(),(),(),()') + return vecfunc(mu, rvec, vvec) \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5ea24cc0d..25411a4dc 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,3 +1,12 @@ +# 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. + ######################################## # Set up how to compile the source files ########################################