From 19f9e7d0b735f38b5b2e0f5d21fe9b4fb2b79838 Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 30 Jun 2021 15:12:41 -0400 Subject: [PATCH] Added follow method to swiftest that recreates the functionality of tool_follow.f (Swift version only so far) --- python/swiftest/swiftest/init_cond.py | 175 ++++++++++ python/swiftest/swiftest/io.py | 180 +---------- python/swiftest/swiftest/simulation_class.py | 133 ++++---- python/swiftest/swiftest/tool.py | 85 +++++ .../tests/bin2xr/swiftest/bin2xr_swiftest.py | 5 + .../tests/bin2xr/swiftest/cb.swiftest.in | 4 + .../tests/bin2xr/swiftest/param.swiftest.in | 31 ++ .../tests/bin2xr/swiftest/pl.swiftest.in | 29 ++ .../tests/bin2xr/swiftest/tp.swiftest.in | 4 + python/swiftest/tests/pl.in | 28 -- python/swiftest/tests/swift.in | 3 - python/swiftest/tests/tp.in | 301 ------------------ 12 files changed, 419 insertions(+), 559 deletions(-) create mode 100644 python/swiftest/swiftest/init_cond.py create mode 100644 python/swiftest/swiftest/tool.py create mode 100644 python/swiftest/tests/bin2xr/swiftest/bin2xr_swiftest.py create mode 100644 python/swiftest/tests/bin2xr/swiftest/cb.swiftest.in create mode 100644 python/swiftest/tests/bin2xr/swiftest/param.swiftest.in create mode 100644 python/swiftest/tests/bin2xr/swiftest/pl.swiftest.in create mode 100644 python/swiftest/tests/bin2xr/swiftest/tp.swiftest.in delete mode 100644 python/swiftest/tests/pl.in delete mode 100644 python/swiftest/tests/swift.in delete mode 100644 python/swiftest/tests/tp.in diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py new file mode 100644 index 000000000..d54fb1162 --- /dev/null +++ b/python/swiftest/swiftest/init_cond.py @@ -0,0 +1,175 @@ +import swiftest +import numpy as np +from astroquery.jplhorizons import Horizons +import datetime +import xarray as xr + +def solar_system_pl(param, ephemerides_start_date): + """ + Initializes a Swiftest dataset containing the major planets of the Solar System at a particular data from JPL/Horizons + + Parameters + ---------- + param : dict + 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 + + Returns + ------- + xarray dataset + """ + # Planet ids + planetid = { + 'mercury': '1', + 'venus': '2', + 'earthmoon': '3', + 'mars': '4', + 'jupiter': '5', + 'saturn': '6', + 'uranus': '7', + 'neptune': '8', + 'plutocharon': '9' + } + + # Planet MSun/M ratio + MSun_over_Mpl = { + 'mercury': np.longdouble(6023600.0), + 'venus': np.longdouble(408523.71), + 'earthmoon': np.longdouble(328900.56), + 'mars': np.longdouble(3098708.), + 'jupiter': np.longdouble(1047.3486), + 'saturn': np.longdouble(3497.898), + 'uranus': np.longdouble(22902.98), + 'neptune': np.longdouble(19412.24), + 'plutocharon': np.longdouble(1.35e8) + } + + # Planet radii in AU + planetradius = { + 'mercury': np.longdouble(2439.4e3) / swiftest.AU2M, + 'venus': np.longdouble(6051.8e3) / swiftest.AU2M, + 'earthmoon': np.longdouble(6371.0084e3) / swiftest.AU2M, # Earth only for radius + 'mars': np.longdouble(3389.50e3) / swiftest.AU2M, + 'jupiter': np.longdouble(69911e3) / swiftest.AU2M, + 'saturn': np.longdouble(58232.0e3) / swiftest.AU2M, + 'uranus': np.longdouble(25362.e3) / swiftest.AU2M, + 'neptune': np.longdouble(24622.e3) / swiftest.AU2M, + 'plutocharon': np.longdouble(1188.3e3 / swiftest.AU2M) + } + + # Unit conversion factors + DCONV = swiftest.AU2M / param['DU2M'] + VCONV = (swiftest.AU2M / swiftest.JD2S) / (param['DU2M'] / param['TU2S']) + THIRDLONG = np.longdouble(1.0) / np.longdouble(3.0) + + # Central body value vectors + GMcb = np.array([swiftest.GMSunSI * 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]) + cbid = np.array([0]) + cvec = np.vstack([GMcb, Rcb, J2RP2, J4RP4]) + + # Horizons date time internal variables + tstart = datetime.date.fromisoformat(ephemerides_start_date) + tstep = datetime.timedelta(days=1) + tend = tstart + tstep + ephemerides_end_date = tend.isoformat() + ephemerides_step = '1d' + + pldata = {} + p1 = [] + p2 = [] + p3 = [] + p4 = [] + p5 = [] + p6 = [] + p7 = [] + p8 = [] + p9 = [] + p10 = [] + p11 = [] + p12 = [] + Rhill = [] + Rpl = [] + GMpl = [] + + # Fetch solar system ephemerides from Horizons + for key, val in planetid.items(): + pldata[key] = Horizons(id=val, id_type='majorbody', location='@sun', + epochs={'start': ephemerides_start_date, 'stop': ephemerides_end_date, + 'step': ephemerides_step}) + if param['OUT_FORM'] == 'XV': + p1.append(pldata[key].vectors()['x'][0] * DCONV) + p2.append(pldata[key].vectors()['y'][0] * DCONV) + p3.append(pldata[key].vectors()['z'][0] * DCONV) + p4.append(pldata[key].vectors()['vx'][0] * VCONV) + p5.append(pldata[key].vectors()['vy'][0] * VCONV) + p6.append(pldata[key].vectors()['vz'][0] * VCONV) + p7.append(pldata[key].elements()['a'][0] * DCONV) + p8.append(pldata[key].elements()['e'][0]) + p9.append(pldata[key].elements()['incl'][0] * np.pi / 180.0) + p10.append(pldata[key].elements()['Omega'][0] * np.pi / 180.0) + p11.append(pldata[key].elements()['w'][0] * np.pi / 180.0) + p12.append(pldata[key].elements()['M'][0] * np.pi / 180.0) + elif param['OUT_FORM'] == 'EL': + p1.append(pldata[key].elements()['a'][0] * DCONV) + p2.append(pldata[key].elements()['e'][0]) + p3.append(pldata[key].elements()['incl'][0] * np.pi / 180.0) + p4.append(pldata[key].elements()['Omega'][0] * np.pi / 180.0) + p5.append(pldata[key].elements()['w'][0] * np.pi / 180.0) + p6.append(pldata[key].elements()['M'][0] * np.pi / 180.0) + p7.append(pldata[key].vectors()['x'][0] * DCONV) + p8.append(pldata[key].vectors()['y'][0] * DCONV) + p9.append(pldata[key].vectors()['z'][0] * DCONV) + p10.append(pldata[key].vectors()['vx'][0] * VCONV) + p11.append(pldata[key].vectors()['vy'][0] * VCONV) + p12.append(pldata[key].vectors()['vz'][0] * VCONV) + Rhill.append(pldata[key].elements()['a'][0] * (3 * MSun_over_Mpl[key]) ** (-THIRDLONG)) + Rpl.append(planetradius[key] * DCONV) + GMpl.append(GMcb[0] / MSun_over_Mpl[key]) + # Generate planet value vectors + plid = np.fromiter(planetid.values(), dtype=int) + pvec = np.vstack([p1, p2, p3, p4, p5, p6, GMpl, Rpl, p7, p8, p9, p10, p11, p12, Rhill]) + + dims = ['time', 'id', 'vec'] + cb = [] + pl = [] + t = np.array([0.0]) + + clab, plab, tlab = swiftest.io.make_swiftest_labels(param) + if param['OUT_FORM'] == 'XV': + plab.append('a') + plab.append('e') + plab.append('inc') + plab.append('capom') + plab.append('omega') + plab.append('capm') + elif param['OUT_FORM'] == 'EL': + plab.append('px') + plab.append('py') + plab.append('pz') + plab.append('vx') + plab.append('vy') + plab.append('vz') + plab.append('Rhill') + + # Prepare frames by adding an extra axis for the time coordinate + cbframe = np.expand_dims(cvec.T, axis=0) + plframe = np.expand_dims(pvec.T, axis=0) + + # Create xarray DataArrays out of each body type + cbxr = xr.DataArray(cbframe, dims=dims, coords={'time': t, 'id': cbid, 'vec': clab}) + plxr = xr.DataArray(plframe, dims=dims, coords={'time': t, 'id': plid, 'vec': plab}) + + cb.append(cbxr) + pl.append(plxr) + + cbda = xr.concat(cb, dim='time') + plda = xr.concat(pl, dim='time') + + cbds = cbda.to_dataset(dim='vec') + plds = plda.to_dataset(dim='vec') + ds = xr.combine_by_coords([cbds, plds]) + return ds \ No newline at end of file diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 81d5758f6..14742c8a5 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -2,8 +2,6 @@ import numpy as np from scipy.io import FortranFile import xarray as xr -from astroquery.jplhorizons import Horizons -import datetime import sys import tempfile @@ -546,7 +544,6 @@ def swiftest_stream(f, param): npl, plid, pvec.T, plab, \ ntp, tpid, tvec.T, tlab - def swifter2xr(param): """ Converts a Swifter binary data file into an xarray DataSet. @@ -576,16 +573,19 @@ def swifter2xr(param): pl.append(plxr) tp.append(tpxr) + sys.stdout.write('\r' + f"Reading in time {t[0]:.3e}") + sys.stdout.flush() plda = xr.concat(pl, dim='time') tpda = xr.concat(tp, dim='time') plds = plda.to_dataset(dim='vec') tpds = tpda.to_dataset(dim='vec') + print('\nCreating Dataset') ds = xr.combine_by_coords([plds, tpds]) + print(f"Successfully converted {ds.sizes['time']} output frames.") return ds - def swiftest2xr(param): """ Converts a Swiftest binary data file into an xarray DataSet. @@ -621,6 +621,8 @@ def swiftest2xr(param): cb.append(cbxr) pl.append(plxr) tp.append(tpxr) + sys.stdout.write('\r' + f"Reading in time {t[0]:.3e}") + sys.stdout.flush() cbda = xr.concat(cb, dim='time') plda = xr.concat(pl, dim='time') @@ -629,179 +631,13 @@ def swiftest2xr(param): cbds = cbda.to_dataset(dim='vec') plds = plda.to_dataset(dim='vec') tpds = tpda.to_dataset(dim='vec') + print('\nCreating Dataset') ds = xr.combine_by_coords([cbds, plds, tpds]) + print(f"Successfully converted {ds.sizes['time']} output frames.") return ds -def solar_system_pl(param, ephemerides_start_date): - """ - Initializes a Swiftest dataset containing the major planets of the Solar System at a particular data from JPL/Horizons - - Parameters - ---------- - param : dict - 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 - - Returns - ------- - xarray dataset - """ - # Planet ids - planetid = { - 'mercury': '1', - 'venus': '2', - 'earthmoon': '3', - 'mars': '4', - 'jupiter': '5', - 'saturn': '6', - 'uranus': '7', - 'neptune': '8', - 'plutocharon': '9' - } - - # Planet MSun/M ratio - MSun_over_Mpl = { - 'mercury': np.longdouble(6023600.0), - 'venus': np.longdouble(408523.71), - 'earthmoon': np.longdouble(328900.56), - 'mars': np.longdouble(3098708.), - 'jupiter': np.longdouble(1047.3486), - 'saturn': np.longdouble(3497.898), - 'uranus': np.longdouble(22902.98), - 'neptune': np.longdouble(19412.24), - 'plutocharon': np.longdouble(1.35e8) - } - - # Planet radii in AU - planetradius = { - 'mercury': np.longdouble(2439.4e3) / swiftest.AU2M, - 'venus': np.longdouble(6051.8e3) / swiftest.AU2M, - 'earthmoon': np.longdouble(6371.0084e3) / swiftest.AU2M, # Earth only for radius - 'mars': np.longdouble(3389.50e3) / swiftest.AU2M, - 'jupiter': np.longdouble(69911e3) / swiftest.AU2M, - 'saturn': np.longdouble(58232.0e3) / swiftest.AU2M, - 'uranus': np.longdouble(25362.e3) / swiftest.AU2M, - 'neptune': np.longdouble(24622.e3) / swiftest.AU2M, - 'plutocharon': np.longdouble(1188.3e3 / swiftest.AU2M) - } - - # Unit conversion factors - DCONV = swiftest.AU2M / param['DU2M'] - VCONV = (swiftest.AU2M / swiftest.JD2S) / (param['DU2M'] / param['TU2S']) - THIRDLONG = np.longdouble(1.0) / np.longdouble(3.0) - - # Central body value vectors - GMcb = np.array([swiftest.GMSunSI * 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]) - cbid = np.array([0]) - cvec = np.vstack([GMcb, Rcb, J2RP2, J4RP4]) - - # Horizons date time internal variables - tstart = datetime.date.fromisoformat(ephemerides_start_date) - tstep = datetime.timedelta(days=1) - tend = tstart + tstep - ephemerides_end_date = tend.isoformat() - ephemerides_step = '1d' - - pldata = {} - p1 = [] - p2 = [] - p3 = [] - p4 = [] - p5 = [] - p6 = [] - p7 = [] - p8 = [] - p9 = [] - p10 = [] - p11 = [] - p12 = [] - Rhill = [] - Rpl = [] - GMpl = [] - - # Fetch solar system ephemerides from Horizons - for key, val in planetid.items(): - pldata[key] = Horizons(id=val, id_type='majorbody', location='@sun', - epochs={'start': ephemerides_start_date, 'stop': ephemerides_end_date, - 'step': ephemerides_step}) - if param['OUT_FORM'] == 'XV': - p1.append(pldata[key].vectors()['x'][0] * DCONV) - p2.append(pldata[key].vectors()['y'][0] * DCONV) - p3.append(pldata[key].vectors()['z'][0] * DCONV) - p4.append(pldata[key].vectors()['vx'][0] * VCONV) - p5.append(pldata[key].vectors()['vy'][0] * VCONV) - p6.append(pldata[key].vectors()['vz'][0] * VCONV) - p7.append(pldata[key].elements()['a'][0] * DCONV) - p8.append(pldata[key].elements()['e'][0]) - p9.append(pldata[key].elements()['incl'][0] * np.pi / 180.0) - p10.append(pldata[key].elements()['Omega'][0] * np.pi / 180.0) - p11.append(pldata[key].elements()['w'][0] * np.pi / 180.0) - p12.append(pldata[key].elements()['M'][0] * np.pi / 180.0) - elif param['OUT_FORM'] == 'EL': - p1.append(pldata[key].elements()['a'][0] * DCONV) - p2.append(pldata[key].elements()['e'][0]) - p3.append(pldata[key].elements()['incl'][0] * np.pi / 180.0) - p4.append(pldata[key].elements()['Omega'][0] * np.pi / 180.0) - p5.append(pldata[key].elements()['w'][0] * np.pi / 180.0) - p6.append(pldata[key].elements()['M'][0] * np.pi / 180.0) - p7.append(pldata[key].vectors()['x'][0] * DCONV) - p8.append(pldata[key].vectors()['y'][0] * DCONV) - p9.append(pldata[key].vectors()['z'][0] * DCONV) - p10.append(pldata[key].vectors()['vx'][0] * VCONV) - p11.append(pldata[key].vectors()['vy'][0] * VCONV) - p12.append(pldata[key].vectors()['vz'][0] * VCONV) - Rhill.append(pldata[key].elements()['a'][0] * (3 * MSun_over_Mpl[key]) ** (-THIRDLONG)) - Rpl.append(planetradius[key] * DCONV) - GMpl.append(GMcb[0] / MSun_over_Mpl[key]) - # Generate planet value vectors - plid = np.fromiter(planetid.values(), dtype=int) - pvec = np.vstack([p1, p2, p3, p4, p5, p6, GMpl, Rpl, p7, p8, p9, p10, p11, p12, Rhill]) - - dims = ['time', 'id', 'vec'] - cb = [] - pl = [] - t = np.array([0.0]) - - clab, plab, tlab = make_swiftest_labels(param) - if param['OUT_FORM'] == 'XV': - plab.append('a') - plab.append('e') - plab.append('inc') - plab.append('capom') - plab.append('omega') - plab.append('capm') - elif param['OUT_FORM'] == 'EL': - plab.append('px') - plab.append('py') - plab.append('pz') - plab.append('vx') - plab.append('vy') - plab.append('vz') - plab.append('Rhill') - - # Prepare frames by adding an extra axis for the time coordinate - cbframe = np.expand_dims(cvec.T, axis=0) - plframe = np.expand_dims(pvec.T, axis=0) - # Create xarray DataArrays out of each body type - cbxr = xr.DataArray(cbframe, dims=dims, coords={'time': t, 'id': cbid, 'vec': clab}) - plxr = xr.DataArray(plframe, dims=dims, coords={'time': t, 'id': plid, 'vec': plab}) - - cb.append(cbxr) - pl.append(plxr) - - cbda = xr.concat(cb, dim='time') - plda = xr.concat(pl, dim='time') - - cbds = cbda.to_dataset(dim='vec') - plds = plda.to_dataset(dim='vec') - ds = xr.combine_by_coords([cbds, plds]) - return ds def swiftest_xr2infile(ds, param, framenum=-1): diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 93b4c51cd..94f922dfe 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -1,84 +1,84 @@ from swiftest import io - +from swiftest import tool class Simulation: """ This is a class that define the basic Swift/Swifter/Swiftest simulation object """ - def __init__(self, codename="Swiftest", source=""): - if source == "": - self.param = { - '! VERSION': f"Swiftest parameter input", - 'T0': "0.0", - 'TSTOP': "0.0", - 'DT': "0.0", - 'PL_IN': "", - 'TP_IN': "", - 'CB_IN': "", - 'IN_TYPE': "REAL8", - 'ISTEP_OUT': "1", - 'ISTEP_DUMP': "1", - 'BIN_OUT': "bin.dat", - 'OUT_TYPE': 'REAL8', - 'OUT_FORM': "XV", - 'OUT_STAT': "REPLACE", - 'J2': "0.0", - 'J4': "0.0", - 'CHK_RMIN': "-1.0", - 'CHK_RMAX': "-1.0", - 'CHK_EJECT': "-1.0", - 'CHK_QMIN': "-1.0", - 'CHK_QMIN_COORD': "HELIO", - 'CHK_QMIN_RANGE': "", - 'ENC_OUT': "", - 'MTINY': "-1.0", - 'MU2KG': "-1.0", - 'TU2S': "-1.0", - 'DU2M': "-1.0", - 'GU': "-1.0", - 'EXTRA_FORCE': "NO", - 'BIG_DISCARD': "NO", - 'CHK_CLOSE': "NO", - 'FRAGMENTATION': "NO", - 'MTINY_SET': "NO", - 'ROTATION': "NO", - 'TIDES': "NO", - 'ENERGY': "NO", - 'GR': "NO", - 'YARKOVSKY': "NO", - 'YORP': "NO", - } - else: - self.read_param(source, codename) + def __init__(self, codename="Swiftest", param_file=""): + self.ds = None + self.param = { + '! VERSION': f"Swiftest parameter input", + 'T0': "0.0", + 'TSTOP': "0.0", + 'DT': "0.0", + 'PL_IN': "", + 'TP_IN': "", + 'CB_IN': "", + 'IN_TYPE': "REAL8", + 'ISTEP_OUT': "1", + 'ISTEP_DUMP': "1", + 'BIN_OUT': "bin.dat", + 'OUT_TYPE': 'REAL8', + 'OUT_FORM': "XV", + 'OUT_STAT': "REPLACE", + 'J2': "0.0", + 'J4': "0.0", + 'CHK_RMIN': "-1.0", + 'CHK_RMAX': "-1.0", + 'CHK_EJECT': "-1.0", + 'CHK_QMIN': "-1.0", + 'CHK_QMIN_COORD': "HELIO", + 'CHK_QMIN_RANGE': "", + 'ENC_OUT': "", + 'MTINY': "-1.0", + 'MU2KG': "-1.0", + 'TU2S': "-1.0", + 'DU2M': "-1.0", + 'GU': "-1.0", + 'EXTRA_FORCE': "NO", + 'BIG_DISCARD': "NO", + 'CHK_CLOSE': "NO", + 'FRAGMENTATION': "NO", + 'MTINY_SET': "NO", + 'ROTATION': "NO", + 'TIDES': "NO", + 'ENERGY': "NO", + 'GR': "NO", + 'YARKOVSKY': "NO", + 'YORP': "NO", + } + if param_file != "" : + self.read_param(param_file, codename) return - def read_param(self, param_file_name, codename="Swiftest"): + def read_param(self, param_file, codename="Swiftest"): if codename == "Swiftest": - self.param = io.read_swiftest_param(param_file_name) + self.param = io.read_swiftest_param(param_file) self.codename = "Swiftest" elif codename == "Swifter": - self.param = io.read_swifter_param(param_file_name) + self.param = io.read_swifter_param(param_file) self.codename = "Swifter" elif codename == "Swift": - self.param = io.read_swift_param(param_file_name) + self.param = io.read_swift_param(param_file) self.codename = "Swift" else: 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_name): + def write_param(self, param_file): # Check to see if the parameter type matches the output type. If not, we need to convert codename = self.param['! VERSION'].split()[0] if codename == "Swifter" or codename == "Swiftest": - io.write_labeled_param(self.param, param_file_name) + io.write_labeled_param(self.param, param_file) elif codename == "Swift": - io.write_swift_param(self.param, param_file_name) + io.write_swift_param(self.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".') return - def convert(self, param_file_name, 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 code type to another (Swift, Swifter, or Swiftest). Returns the old parameter configuration. """ @@ -104,7 +104,7 @@ def convert(self, param_file_name, newcodename="Swiftest", plname="pl.swiftest.i goodconversion = False if goodconversion: - self.write_param(param_file_name) + self.write_param(param_file) else: print(f"Conversion from {self.codename} to {newcodename} is not supported.") return oldparam @@ -123,3 +123,26 @@ def bin2xr(self): return + def follow(self, codestyle="Swifter"): + if self.ds is None: + self.bin2xr() + 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 + i_list = [i for i in line.split(" ") if i.strip()] + ifol = int(i_list[0]) + line = f.readline() # nskp + i_list = [i for i in line.split(" ") if i.strip()] + nskp = int(i_list[0]) + except IOError: + print('No follow.in file found') + ifol = None + nskp = None + fol = tool.follow_swift(self.ds, ifol=ifol, nskp=nskp) + + print('follow.out written') + return fol \ No newline at end of file diff --git a/python/swiftest/swiftest/tool.py b/python/swiftest/swiftest/tool.py new file mode 100644 index 000000000..fba725f4e --- /dev/null +++ b/python/swiftest/swiftest/tool.py @@ -0,0 +1,85 @@ +import swiftest +import numpy as np +import os +import glob +from pyslalib import slalib +import xarray as xr +""" +Functions that recreate the Swift/Swifter tool programs +""" + +def sla_dranrm(angle): + func = np.vectorize(slalib.sla_dranrm) + return xr.apply_ufunc(func, angle) + +def follow_swift(ds, ifol=None, nskp=None): + """ + Emulates the Swift version of tool_follow.f + + + Parameters + ---------- + ds : Dataset containing orbital elements + + Returns + ------- + fol : Dataset containing only the followed body with angles converted to degrees + + Generates a file called follow.out containing 10 columns (angles are all in degrees): + 1 2 3 4 5 6 7 8 9 10 + t,a,e,inc,capom,omega,capm,peri,apo,obar + + """ + fol = None + if ifol is None: + intxt = input(' Input the particle number to follow \n') + ifol = int(intxt) + print(f"Following particle {ifol}") + if ifol < 0: # Negative numbers are planets + fol = ds.where(np.invert(np.isnan(ds['Mass'])), drop=True) + fol = fol.where(np.invert(np.isnan(fol['a'])), drop=True) # Remove times where this body doesn't exist (but this also gets rid of the central body) + fol = fol.isel(id = -ifol - 2) # Take 1 off for 0-indexed arrays in Python, and take 1 more off because the central body is gone + elif ifol > 0: # Positive numbers are test particles + fol = ds.where(np.isnan(ds['Mass']), drop=True).drop_vars(['Mass', 'Radius']) + fol = fol.where(np.invert(np.isnan(fol['a'])), drop=True) + fol = fol.isel(id = ifol - 1) # Take 1 off for 0-indexed arrays in Python + + if nskp is None: + intxt = input('Input the print frequency\n') + nskp = int(intxt) + + + dr = 180.0 / np.pi + fol['obar'] = fol['capom'] + fol['omega'] + fol['obar'] = fol['obar'].fillna(0) + fol['obar'] = sla_dranrm(fol['obar']) + fol['obar'] = fol['obar'] * dr + fol['inc'] = fol['inc'] * dr + fol['capom'] = fol['capom'] * dr + fol['omega'] = fol['omega'] * dr + fol['capm'] = fol['capm'] * dr + fol['peri'] = fol['a'] * (1.0 - fol['e']) + fol['apo'] = fol['a'] * (1.0 + fol['e']) + + print('1 2 3 4 5 6 7 8 9 10') + print('t,a,e,inc,capom,omega,capm,peri,apo,obar') + + tslice = slice(None, None, nskp) + try: + with open('follow.out', 'w') as f: + for t in fol.isel(time=tslice).time: + a = fol['a'].sel(time=t).values + e = fol['e'].sel(time=t).values + inc = fol['inc'].sel(time=t).values + capom = fol['capom'].sel(time=t).values + omega = fol['omega'].sel(time=t).values + capm = fol['capm'].sel(time=t).values + peri = fol['peri'].sel(time=t).values + apo = fol['apo'].sel(time=t).values + obar = fol['obar'].sel(time=t).values + print(f"{t.values:15.7e} {a:10.4f} {e:7.5f} {inc:9.4f} {capom:9.4f} {omega:9.4f} {capm:9.4f} {peri:10.4f} {apo:10.4f} {obar:9.4f}", file=f) + + except IOError: + print(f"Error writing to follow.out") + + return fol \ No newline at end of file diff --git a/python/swiftest/tests/bin2xr/swiftest/bin2xr_swiftest.py b/python/swiftest/tests/bin2xr/swiftest/bin2xr_swiftest.py new file mode 100644 index 000000000..4437247be --- /dev/null +++ b/python/swiftest/tests/bin2xr/swiftest/bin2xr_swiftest.py @@ -0,0 +1,5 @@ +import swiftest +import swiftest.io as swio +param_file = "param.swiftest.in" +sim = swiftest.Simulation(source="param.swiftest.in") +ds = swio.swiftest2xr(sim.param) \ No newline at end of file diff --git a/python/swiftest/tests/bin2xr/swiftest/cb.swiftest.in b/python/swiftest/tests/bin2xr/swiftest/cb.swiftest.in new file mode 100644 index 000000000..e50255f4a --- /dev/null +++ b/python/swiftest/tests/bin2xr/swiftest/cb.swiftest.in @@ -0,0 +1,4 @@ +0.000295923385935 +0.00468 +0.0 +0.0 diff --git a/python/swiftest/tests/bin2xr/swiftest/param.swiftest.in b/python/swiftest/tests/bin2xr/swiftest/param.swiftest.in new file mode 100644 index 000000000..6a9e599f9 --- /dev/null +++ b/python/swiftest/tests/bin2xr/swiftest/param.swiftest.in @@ -0,0 +1,31 @@ +! VERSION Swiftest parameter file converted from Swift +BIG_DISCARD NO +BIN_OUT bin.dat +CB_IN cb.swiftest.in +CHK_CLOSE YES +CHK_QMIN 0.00468 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 4.68e-03 100.0 +CHK_RMAX 100.0 +CHK_RMIN 0.00468 +DT 5.0 +DU2M 149597870700.0 +ENC_OUT enc.dat +ENERGY NO +EXTRA_FORCE NO +FRAGMENTATION NO +GR NO +ISTEP_DUMP 7305000 +ISTEP_OUT 7305000 +MU2KG 1.988409870698051e+30 +OUT_FORM EL +OUT_STAT UNKNOWN +PL_IN pl.swiftest.in +ROTATION NO +T0 0.0 +TIDES NO +TP_IN tp.swiftest.in +TSTOP 365250000000.0 +TU2S 86400 +YARKOVSKY NO +YORP NO diff --git a/python/swiftest/tests/bin2xr/swiftest/pl.swiftest.in b/python/swiftest/tests/bin2xr/swiftest/pl.swiftest.in new file mode 100644 index 000000000..b7814ac98 --- /dev/null +++ b/python/swiftest/tests/bin2xr/swiftest/pl.swiftest.in @@ -0,0 +1,29 @@ +7 +2 4.91254745e-11 +1.63083872e-05 +-0.1407280799640445 -0.4439009577663329 -0.02334555971312333 +0.021170988258808896 -0.007097975438870806 -0.002522830951443755 +3 7.24345249e-10 +4.04544528e-05 +-0.7186302169039649 -0.02250380105571625 0.04117184137682461 +0.0005135327579269579 -0.02030614162239802 -0.0003071745100210849 +4 8.99701139e-10 +4.26352325e-05 +-0.1685504489908262 0.9687636586748941 -1.151849464023427e-06 +-0.01723001047964496 -0.003013273823615044 2.412522755326018e-08 +5 9.54953511e-11 +2.27075425e-05 +1.390361066044069 -0.02100972149003466 -0.03461801461346188 +0.0007479271121979245 0.01518629868805657 0.0002997532108600616 +6 2.8238813e-07 +0.000477894503 +4.003460074251473 2.93535365924021 -0.1018232730751003 +-0.004563473375353001 0.006446757832442719 7.545633314092875e-05 +7 8.4549815e-08 +0.000402866697 +6.408554725090048 6.568044145835263 -0.3691278858265656 +-0.004291119710689724 0.003891580681674363 0.0001028767891901977 +8 1.291492e-08 +0.000170851362 +14.43051751807564 -13.73565829091644 -0.238128468762112 +0.002678379897839438 0.002672442537311352 -2.477647055276606e-05 diff --git a/python/swiftest/tests/bin2xr/swiftest/tp.swiftest.in b/python/swiftest/tests/bin2xr/swiftest/tp.swiftest.in new file mode 100644 index 000000000..b3b8f423c --- /dev/null +++ b/python/swiftest/tests/bin2xr/swiftest/tp.swiftest.in @@ -0,0 +1,4 @@ +1 +9 +1.407578089380049 0.5926170247361682 -0.7567948854164774 +-0.004970782813745737 0.01294038841446245 3.261913169537027e-05 diff --git a/python/swiftest/tests/pl.in b/python/swiftest/tests/pl.in deleted file mode 100644 index 4d479d3cf..000000000 --- a/python/swiftest/tests/pl.in +++ /dev/null @@ -1,28 +0,0 @@ - 8 - 2.9592338592955439E-004 - 0.0000000000000000 0.0000000000000000 0.0000000000000000 - 0.0000000000000000 0.0000000000000000 0.0000000000000000 - 4.9127330156310911E-011 1.6308387199999999E-005 - -0.15273942296363005 -0.44077188683118762 -2.1987973899201478E-002 - 2.0924205227543746E-002 -7.6522013575130771E-003 -2.5455577585009612E-003 - 7.2437260968171072E-010 4.0454452799999999E-005 - -0.68052645343902751 -0.21255610959625498 3.6375571604128076E-002 - 5.7676478402592110E-003 -1.9654618148759011E-002 -6.0150804827979417E-004 - 8.9973512337453730E-010 4.2635232500000000E-005 - 0.21153912017632401 0.96910814953928759 -1.5976841060778844E-006 - -1.6944135164995145E-002 3.6911953785631699E-003 1.4525071481945407E-008 - 9.5498958252779674E-011 2.2707542499999999E-005 - 1.4232971230754801 0.21815483859431103 -3.0416153305576600E-002 - -1.9274688569975952E-003 1.4590108763912674E-002 3.5304139123317309E-004 - 2.8254526327676804E-007 4.7789450300000003E-004 - 4.0036874054921974 2.9350393270474617 -0.10182563723380476 - -4.5629876169643530E-003 6.4471059617920590E-003 7.5448856044019538E-005 - 8.4600347388504292E-008 4.0286669700000002E-004 - 6.4083018390909805 6.5682871940009111 -0.36911425745370574 - -4.2912688970106801E-003 3.8914184987695602E-003 1.0289743699759401E-004 - 1.2920737211033427E-008 1.7085136200000001E-004 - 14.430329793174279 -13.735823636348456 -0.23812543891764848 - 2.6784223946886604E-003 2.6724061291222238E-003 -2.4777620883141212E-005 - 1.5244164811971947E-008 1.6553711599999999E-004 - 16.808300160620913 -24.994308693440338 0.12729909149223279 - 2.5795430409598025E-003 1.7765162604947151E-003 -9.5908529499802122E-005 diff --git a/python/swiftest/tests/swift.in b/python/swiftest/tests/swift.in deleted file mode 100644 index 74b6b1061..000000000 --- a/python/swiftest/tests/swift.in +++ /dev/null @@ -1,3 +0,0 @@ -param.in -pl.in -tp.in diff --git a/python/swiftest/tests/tp.in b/python/swiftest/tests/tp.in deleted file mode 100644 index 1a843a151..000000000 --- a/python/swiftest/tests/tp.in +++ /dev/null @@ -1,301 +0,0 @@ - 50 - 1.75536839836414 -0.397301559901775 -2.886320950191011E-002 - 2.828957660824792E-003 1.250539366801959E-002 -8.833159624221525E-005 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.75462547108338 -0.397439253344038 -5.771762699017597E-002 - 2.826684041441025E-003 1.250497227754510E-002 -1.766362858232247E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.75338751027446 -0.397668695834012 -8.655446554537932E-002 - 2.822895444946565E-003 1.250427010264027E-002 -2.648871776719500E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.75165489340570 -0.397989817412071 -0.115364933939178 - 2.817593026527604E-003 1.250328735740620E-002 -3.530573675304868E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.74942814808290 -0.398402520292259 -0.144140261060723 - 2.810778400843813E-003 1.250202434110241E-002 -4.411200127075503E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.74670795159950 -0.398906678945642 -0.172871691292285 - 2.802453640651961E-003 1.250048143789167E-002 -5.290483179262686E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.74349513437761 -0.399502139462572 -0.201550453497502 - 2.792621287336390E-003 1.249865911879178E-002 -6.168154404172742E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.73979067424302 -0.400188720613739 -0.230167821489439 - 2.781284333388745E-003 1.249655793842830E-002 -7.043946749723006E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.73559569961003 -0.400966213259897 -0.258715078141622 - 2.768446232154610E-003 1.249417853684106E-002 -7.917593441112837E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.73091148830840 -0.401834380569283 -0.287183527684168 - 2.754110894243519E-003 1.249152163881874E-002 -8.788828357128555E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.72026103118157 -0.389355528881347 -2.828594532477477E-002 - 2.886114093303328E-003 1.275805340156797E-002 -9.011625318008804E-005 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.71953296244609 -0.389490468454826 -5.656327447617424E-002 - 2.883794537575375E-003 1.275762349729473E-002 -1.802050560751523E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.71831976085279 -0.389715322095103 -8.482337627316461E-002 - 2.879929396047131E-003 1.275690713560546E-002 -2.702389743052803E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.71662179632064 -0.390030021241745 -0.113057635311967 - 2.874519847244236E-003 1.275590453492696E-002 -3.601905600373064E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.71443958590329 -0.390434470064514 -0.141257455903944 - 2.867567538449044E-003 1.275461600056502E-002 -4.500324282485798E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.71177379334835 -0.390928545545054 -0.169414257543719 - 2.859074584296430E-003 1.275304192444406E-002 -5.397372422888266E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.70862523246946 -0.391512096851912 -0.197519444517652 - 2.849043577517016E-003 1.275118278709834E-002 -6.292776170557357E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.70499486153590 -0.392184946380362 -0.225564465162543 - 2.837477571062932E-003 1.274903915435912E-002 -7.186263077873903E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.70088378639370 -0.392946889173945 -0.253540776694444 - 2.824380088051367E-003 1.274661167919763E-002 -8.077560980102980E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.69629325931601 -0.393797693137531 -0.281439857258866 - 2.809755118102047E-003 1.274390110104623E-002 -8.966398379301735E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.68515366399900 -0.381409497860919 -2.770868114763943E-002 - 2.944472485597326E-003 1.301602639267184E-002 -9.193844020566386E-005 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.68444045380880 -0.381541683565614 -5.540892196217252E-002 - 2.942106027515987E-003 1.301558779556061E-002 -1.838488750704571E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.68325201143113 -0.381761948356195 -8.309228700094988E-002 - 2.938162731265411E-003 1.301485694874873E-002 -2.757033154802284E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.68158869923558 -0.382070225071419 -0.110750336684755 - 2.932643799201507E-003 1.301383407507970E-002 -3.674737585955427E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.67945102372369 -0.382466419836768 -0.138374650747165 - 2.925550911915339E-003 1.301251948603269E-002 -4.591322656575339E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.67683963509720 -0.382950412144466 -0.165956823795153 - 2.916886226800555E-003 1.301091358145700E-002 -5.506509472578220E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.67375533056131 -0.383522054241252 -0.193488435537802 - 2.906652389013839E-003 1.300901685160348E-002 -6.420018645562698E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.67019904882879 -0.384181172146986 -0.220961108835647 - 2.894852513239243E-003 1.300682987374474E-002 -7.331572218908740E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.66617187317737 -0.384927565087993 -0.248366475247266 - 2.881490193832806E-003 1.300435331405534E-002 -8.240892524600587E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.66167503032362 -0.385761005705779 -0.275696186833564 - 2.866569501085982E-003 1.300158792691927E-002 -9.147702574897361E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.65004629681643 -0.373463466840492 -2.713141697050409E-002 - 3.004107385165185E-003 1.327964217801076E-002 -9.380048499464011E-005 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.64934794517152 -0.373592898676403 -5.425456944817078E-002 - 3.001692998807830E-003 1.327919469791802E-002 -1.875723974515014E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.64818426200946 -0.373808574617286 -8.136119772873515E-002 - 2.997669838310936E-003 1.327844904914197E-002 -2.812871813881668E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.64655560215051 -0.374110428901093 -0.108443038057544 - 2.992039130381927E-003 1.327740545903941E-002 -3.749162668189549E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.64446246154409 -0.374498369609023 -0.135491845590387 - 2.984802589649166E-003 1.327606424539795E-002 -4.684311491365882E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.64190547684605 -0.374972278743879 -0.162499390046587 - 2.975962417200357E-003 1.327442581616511E-002 -5.618033740837786E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.63888542865316 -0.375532011630592 -0.189457426557952 - 2.965521311764990E-003 1.327249067152087E-002 -6.550044369703316E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.63540323612168 -0.376177397913609 -0.216357752508752 - 2.953482451109342E-003 1.327025940042943E-002 -7.480059791839970E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.63145995996104 -0.376908241002040 -0.243192173800089 - 2.939849502386553E-003 1.326773268255752E-002 -8.407796715574675E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.62705680133123 -0.377724318274027 -0.269952516408261 - 2.924626618324377E-003 1.326491128756776E-002 -9.332972543286958E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.61493892963385 -0.365517435820064 -2.655415279336874E-002 - 3.065098326561108E-003 1.354925233969722E-002 -9.570487094018532E-005 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.61425543653423 -0.365644113787191 -5.310021693416905E-002 - 3.062634922083575E-003 1.354879577463227E-002 -1.913805892481565E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.61311651258780 -0.365855200878378 -7.963010845652042E-002 - 3.058530081302113E-003 1.354803498730923E-002 -2.869980191831766E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.61152250506545 -0.366150632730767 -0.106135739430332 - 3.052785055829394E-003 1.354697020970083E-002 -3.825280107169432E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.60947389936449 -0.366530319381278 -0.132609040433608 - 3.045401595105050E-003 1.354560176604655E-002 -4.779414805269034E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.60697131859490 -0.366994145343291 -0.159041956298021 - 3.036381944904408E-003 1.354393007257628E-002 -5.732093966627254E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.60401552674501 -0.367541969019932 -0.185426417578102 - 3.025728858747958E-003 1.354195563962491E-002 -6.683026757172563E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.60060742341456 -0.368173623680233 -0.211754396181856 - 3.013445578918633E-003 1.353967906811411E-002 -7.631923833270372E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.59674804674471 -0.368888916916088 -0.238017872352911 - 2.999535847022010E-003 1.353710105150953E-002 -8.578496151713587E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.59243857233884 -0.369687630842276 -0.264208845982959 - 2.984003900096668E-003 1.353422237509990E-002 -9.522455377438729E-004 - 0 0 0 0 0 0 - 0 0 0 0 0 0 - 0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0