aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorVolker Hoffmann <volker@cheleb.net>2014-09-04 15:54:38 +0200
committerVolker Hoffmann <volker@cheleb.net>2014-09-04 15:54:38 +0200
commitba7a98c3c0d588c5670e685ab45cf95b400e793c (patch)
tree342981d2346c0a08eaf0eacf954139a6ef27fe8b
parent6b471a6c8f8e2a72c64c0d48532b205b069b5d6f (diff)
option to use quad_precision on kepler elements, also store precision and version field for snapshot
-rw-r--r--Commons/Loaders.py16
-rw-r--r--Commons/Structs.py25
-rw-r--r--Reduce/reduce_gengaout.py19
-rw-r--r--Reduce/reduce_ssascii.py19
4 files changed, 58 insertions, 21 deletions
diff --git a/Commons/Loaders.py b/Commons/Loaders.py
index 1aa422e..3f5eb67 100644
--- a/Commons/Loaders.py
+++ b/Commons/Loaders.py
@@ -9,20 +9,22 @@ from Structs import Snapshot, Particle
import numpy as np
class Loader():
- def __init__(self, nstep, fname, ellipses):
+ def __init__(self, nstep, fname, ellipses, precision):
self.nstep = nstep
self.fname = fname
self.snapshot = Snapshot()
self.snapshot.nstep = nstep
self.snapshot.ellipses = ellipses
+ self.snapshot.precision = precision
class GengaIC(Loader):
pass
class GengaOut(Loader):
- def __init__(self, nstep=1, ellipses=False, run_name="gasrun"):
+ def __init__(self, nstep=1, ellipses=False, run_name="gasrun", \
+ precision="double"):
fname = "Out%s_%012d.dat" % (run_name, nstep)
- Loader.__init__(self, nstep, fname, ellipses)
+ Loader.__init__(self, nstep, fname, ellipses, precision)
def load(self):
with open(self.fname, 'r') as f:
@@ -35,7 +37,7 @@ class GengaOut(Loader):
if first:
self.snapshot.tout = float(line[0]) # yr
first = False
- particle = Particle()
+ particle = Particle(self.snapshot.precision)
particle.id = float(line[1]) # -
particle.m = float(line[2]) # Msun
particle.r = float(line[3]) # AU
@@ -52,9 +54,9 @@ class GengaOut(Loader):
self.snapshot.particles = particles
class SSAscii(Loader):
- def __init__(self, nstep=1, ellipses=False):
+ def __init__(self, nstep=1, ellipses=False, precision="double"):
fname = "Out.%012d.dat" % nstep
- Loader.__init__(self, nstep, fname, ellipses)
+ Loader.__init__(self, nstep, fname, ellipses, precision)
def load(self):
with open(self.fname, 'r') as f:
@@ -67,7 +69,7 @@ class SSAscii(Loader):
if first:
self.snapshot.tout = float(line[0]) / 2. / np.pi # TU -> yr
first = False
- particle = Particle()
+ particle = Particle(self.snapshot.precision)
particle.id = float(line[1]) # -
particle.m = float(line[2]) # Msun
particle.r = float(line[3]) # AU
diff --git a/Commons/Structs.py b/Commons/Structs.py
index 5cae38b..ada3608 100644
--- a/Commons/Structs.py
+++ b/Commons/Structs.py
@@ -19,7 +19,7 @@ class Particle():
NB: For G=M=1, 1 Year = 2 Pi (From Kepler's Third Law).
"""
- def __init__(self):
+ def __init__(self, precision='double'):
# General
self.id = None
# Cartesian
@@ -32,14 +32,25 @@ class Particle():
self.m, self.r = None, None
# Ellipse
self.xell = None; self.yell = None; self.zell = None
+ # Set Precision
+ if precision == 'double' or precision == 'quad':
+ self.precision = precision
+ else:
+ raise Exception("Invalid Kepler Elements Precision Requested.")
+ # File Version
+ self.version = 2
def cart2kep(self):
+ # Set Precision
+ if self.precision == 'quad':
+ x = np.array([self.x, self.y, self.z], dtype=np.float128)
+ v = np.array([self.vx, self.vy, self.vz], dtype=np.float128)
+ elif self.precision == 'double':
+ x = np.array([self.x, self.y, self.z], dtype=np.float64)
+ v = np.array([self.vx, self.vy, self.vz], dtype=np.float64)
+ # Compute Kepler Elements
self.a, self.ecc, self.inc, self.Omega, self.omega, self.M0 = \
- cart2kep(np.array([self.x, self.y, self.z], dtype=np.float128), \
- np.array([self.vx, self.vy, self.vz], \
- dtype=np.float128), \
- self.m, \
- central_mass=1.0)
+ cart2kep(x, v, self.m, central_mass=1.0)
def kep2cart(self):
x, v = \
kep2cart(self.a, self.ecc, self.inc, \
@@ -71,3 +82,5 @@ class Snapshot():
self.nstep = None
self.nparticles = None
self.ellipses = None
+ self.version = 2
+ self.precision = 'double'
diff --git a/Reduce/reduce_gengaout.py b/Reduce/reduce_gengaout.py
index 2e85d10..60db83f 100644
--- a/Reduce/reduce_gengaout.py
+++ b/Reduce/reduce_gengaout.py
@@ -10,11 +10,16 @@ parser.add_argument("--run_name", default='gasrun', \
help='Name of Simulation Run.')
parser.add_argument("--ellipses", action='store_true', \
help='Compute & Store Orbit Ellipses for Particles.')
-group = parser.add_mutually_exclusive_group(required=True)
-group.add_argument('--all', action='store_true', \
+group1 = parser.add_mutually_exclusive_group(required=True)
+group1.add_argument('--all', action='store_true', \
help="Reduce Full Set of Snapshots.")
-group.add_argument('--custom', type=int, nargs='+', \
+group1.add_argument('--custom', type=int, nargs='+', \
help="Plot Custom Snapshot.")
+group2 = parser.add_mutually_exclusive_group(required=True)
+group2.add_argument('--quad', action='store_true', \
+ help="Store Kepler Elements in Quad Precision.")
+group2.add_argument('--double', action='store_true', \
+ help="Store Kepler Elements in Double Precision.")
args = parser.parse_args()
# Sanity Check
@@ -23,6 +28,12 @@ if args.custom:
print "!! Output set must be defined by three numbers."
sys.exit()
+# Precision?
+if args.double:
+ precision = 'double'
+if args.quad:
+ precision = 'quad'
+
# Full Set
if args.all:
nsteps = []
@@ -45,7 +56,7 @@ print "// Starting -- %s UTC" % strftime("%H:%M:%S", gmtime())
for istep, nstep in enumerate(nsteps):
print "// (%s UTC) Processing Snapshot %012d/%012d" % \
(strftime("%H:%M:%S", gmtime()), nstep, nsteps[-1])
- loader = Loaders.GengaOut(nstep, args.ellipses, args.run_name)
+ loader = Loaders.GengaOut(nstep, args.ellipses, args.run_name, precision)
loader.load()
np.savez('Snapshot_%012d.npz' % nstep, snapshot=loader.snapshot)
print "// Done -- %s UTC" % strftime("%H:%M:%S", gmtime())
diff --git a/Reduce/reduce_ssascii.py b/Reduce/reduce_ssascii.py
index f461dfd..7cffdae 100644
--- a/Reduce/reduce_ssascii.py
+++ b/Reduce/reduce_ssascii.py
@@ -10,11 +10,16 @@ parser.add_argument("--run_name", default='Out', \
help='Name of Simulation Run.')
parser.add_argument("--ellipses", action='store_true', \
help='Compute & Store Orbit Ellipses for Particles.')
-group = parser.add_mutually_exclusive_group(required=True)
-group.add_argument('--all', action='store_true', \
+group1 = parser.add_mutually_exclusive_group(required=True)
+group1.add_argument('--all', action='store_true', \
help="Reduce Full Set of Snapshots.")
-group.add_argument('--custom', type=int, nargs='+', \
+group1.add_argument('--custom', type=int, nargs='+', \
help="Plot Custom Snapshot.")
+group2 = parser.add_mutually_exclusive_group(required=True)
+group2.add_argument('--quad', action='store_true', \
+ help="Store Kepler Elements in Quad Precision.")
+group2.add_argument('--double', action='store_true', \
+ help="Store Kepler Elements in Double Precision.")
args = parser.parse_args()
# Sanity Check
@@ -23,6 +28,12 @@ if args.custom:
print "!! Output set must be defined by three numbers."
sys.exit()
+# Precision?
+if args.double:
+ precision = 'double'
+if args.quad:
+ precision = 'quad'
+
# Full Set
if args.all:
nsteps = []
@@ -45,7 +56,7 @@ print "// Starting -- %s UTC" % strftime("%H:%M:%S", gmtime())
for istep, nstep in enumerate(nsteps):
print "// (%s UTC) Processing Snapshot %012d/%012d" % \
(strftime("%H:%M:%S", gmtime()), nstep, nsteps[-1])
- loader = Loaders.SSAscii(nstep, args.ellipses)
+ loader = Loaders.SSAscii(nstep, args.ellipses, precision)
loader.load()
np.savez('Snapshot_%012d.npz' % nstep, snapshot=loader.snapshot)
print "// Done -- %s UTC" % strftime("%H:%M:%S", gmtime())