Skip to content

Commit

Permalink
Revisions for Python 3 compatibility; many new features in pdstable.py
Browse files Browse the repository at this point in the history
  • Loading branch information
@markshowalter
markshowalter committed Apr 2, 2018
1 parent bc6dee6 commit 5f4aa89
Show file tree
Hide file tree
Showing 9 changed files with 950 additions and 437 deletions.
58 changes: 50 additions & 8 deletions gravity.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,9 @@
# - Redefined PLUTO_CHARON as PLUTO_CHARON_OLD
# - Redefined PLUTO_CHARON_AS_RINGS as PLUTO_CHARON
################################################################################

from __future__ import print_function

import numpy as np
import unittest
import warnings
Expand Down Expand Up @@ -106,14 +109,20 @@ def omega(self, a):
Gravity._jseries(self.omega_jn, self.r2/a2))
return np.sqrt(omega2)

def kappa(self, a):
"""Returns the radial oscillation frequency (radians/s) at semimajor
axis a."""
def kappa2(self, a):
"""Returns the square of the radial oscillation frequency (radians/s) at
semimajor axis a."""

a2 = a * a
kappa2 = self.gm/(a*a2) * (1. +
Gravity._jseries(self.kappa_jn, self.r2/a2))
return np.sqrt(kappa2)
return kappa2

def kappa(self, a):
"""Returns the radial oscillation frequency (radians/s) at semimajor
axis a."""

return np.sqrt(self.kappa2(a))

def nu(self, a):
"""Returns the vertical oscillation frequency (radians/s) at semimajor
Expand Down Expand Up @@ -263,7 +272,7 @@ def dcombo_da(self, a, factors):

return sum_values

def solve_a(self, freq, factors=(1,0,0), iters=5):
def solve_a(self, freq, factors=(1,0,0)):
"""Solves for the semimajor axis at which the frequency is equal to the
given combination of factors on omega, kappa and nu. Solution is via
Newton's method."""
Expand Down Expand Up @@ -302,12 +311,23 @@ def solve_a(self, freq, factors=(1,0,0), iters=5):
a = (self.gm * (term * self.r2 * self.r2 / freq)**2)**(1/11.)

# Iterate using Newton's method
for iter in range(iters):
da_prev_max = 1.e99
for iter in range(20):
# a step in Newton's method: x(i+1) = x(i) - f(xi) / fp(xi)
# our f(x) = self.combo() - freq
# fp(x) = self.dcombo()

a -= ((self.combo(a, factors) - freq) / self.dcombo_da(a, factors))
da = ((self.combo(a, factors) - freq) / self.dcombo_da(a, factors))
da_max = np.max(np.abs(da))
if da_max == 0.: break

a -= da

# If Newton's method stops converging, return what we've got
if iter > 4 and da_max >= da_prev_max:
break

da_prev_max = da_max

return a

Expand Down Expand Up @@ -977,7 +997,7 @@ def test_uncombo(self):

# Testing scalars in a loop...
tests = 100
planets = [JUPITER, SATURN, URANUS, NEPTUNE, PLUTO_CHARON]
planets = [JUPITER, SATURN, URANUS, NEPTUNE]
factors = [(1, 0, 0), (0, 1, 0), (0, 0, 1)]

for test in range(tests):
Expand All @@ -988,6 +1008,28 @@ def test_uncombo(self):
c = abs((b - a) / a)
self.assertTrue(c < ERROR_TOLERANCE)

# PLUTO_CHARON with factors (1,0,0) and (0,0,1)
for test in range(tests):
for obj in [PLUTO_CHARON]:
a = obj.rp * 10. ** (np.random.rand() * 2.)
for f in [(1,0,0),(0,0,1)]:
b = obj.solve_a(obj.combo(a, f), f)
c = abs((b - a) / a)
self.assertTrue(c < ERROR_TOLERANCE)

# PLUTO_CHARON with factors (0,1,0) can have duplicated values...
for test in range(tests):
for obj in [PLUTO_CHARON]:
a = obj.rp * 10. ** (np.random.rand() * 2.)
if obj.kappa2(a) < 0.: continue # this would raise RuntimeError

for f in [(0,1,0)]:
combo1 = obj.combo(a,f)
b = obj.solve_a(combo1, f)
combo2 = obj.combo(b,f)
c = abs((combo2 - combo1) / combo1)
self.assertTrue(c < ERROR_TOLERANCE)

# Testing a 100x100 array
for obj in planets:
a = obj.rp * 10. ** (np.random.rand(100,100) * 2.)
Expand Down
4 changes: 3 additions & 1 deletion interval.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,9 +62,11 @@ def __str__(self):

def __repr__(self): return str(self)

def __setitem__(self, (xlo, xhi), value):
def __setitem__(self, xlimits, value):
"""Overlays the given value within the limit s (xlo,xhi)"""

(xlo, xhi) = xlimits

# Ignore ranges completely outside the interval
if xhi <= self.xlist[0]: return
if xlo >= self.xlist[-1]: return
Expand Down
132 changes: 91 additions & 41 deletions julian.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1434,7 +1434,7 @@ def runTest(self):
# ISO format parsers
################################################################################

def day_from_iso(strings, validate=True):
def day_from_iso(strings, validate=True, strip=False):
"""Returns a day number based on a parsing of a date string in ISO format.
The format is strictly required to be either yyyy-mm-dd or yyyy-ddd.
Expand All @@ -1456,18 +1456,32 @@ def day_from_iso(strings, validate=True):

strings = np.array(strings).astype('S')

if validate:
if b" " in bytearray(strings):
raise ValueError("blank character in ISO date")
first_index = len(strings.shape) * (0,)
first = strings[first_index]

# Skip over leading and trailing blanks
if strip:
for k0 in range(len(first)):
if first[k0:k0+1] != b' ':
break
for k1 in range(len(first)):
if first[-k1-1] != b' ':
break
else:
k0 = 0
k1 = 0
if validate:
if b" " in bytearray(strings):
raise ValueError("blank character in ISO date")

# yyyy-mm-dd case:
if strings.dtype == np.dtype("|S10"):
dtype_dict = {"y": ("|S4",0),
"m": ("|S2",5),
"d": ("|S2",8)}
if strings.itemsize - k0 - k1 == 10:
dtype_dict = {"y": ("|S4", k0 + 0),
"m": ("|S2", k0 + 5),
"d": ("|S2", k0 + 8)}
if validate:
dtype_dict["dash1"] = ("|S1",4)
dtype_dict["dash2"] = ("|S1",7)
dtype_dict["dash1"] = ("|S1", k0 + 4)
dtype_dict["dash2"] = ("|S1", k0 + 7)

strings.dtype = np.dtype(dtype_dict)

Expand All @@ -1490,11 +1504,11 @@ def day_from_iso(strings, validate=True):
return day_from_ymd(y,m,d)

# yyyy-ddd case:
if strings.dtype == np.dtype("|S8"):
dtype_dict = {"y": ("|S4",0),
"d": ("|S3",5)}
if strings.itemsize - k0 == 8:
dtype_dict = {"y": ("|S4", k0 + 0),
"d": ("|S3", k0 + 5)}
if validate:
dtype_dict["dash"] = ("|S1",4)
dtype_dict["dash"] = ("|S1", k0 + 4)

strings.dtype = np.dtype(dtype_dict)

Expand All @@ -1517,7 +1531,7 @@ def day_from_iso(strings, validate=True):

########################################

def sec_from_iso(strings, validate=True):
def sec_from_iso(strings, validate=True, strip=False):
"""Returns a second value based on a parsing of a time string in ISO format.
The format is strictly required to be hh:mm:ss[.s...][Z].
Expand All @@ -1540,24 +1554,43 @@ def sec_from_iso(strings, validate=True):
# Convert to an array of strings
strings = np.array(strings).astype('S')

if validate:
merged = bytearray(strings)
if b" " in merged or b"-" in merged:
raise ValueError("blank character in ISO date")
first_index = len(strings.shape) * (0,)
first = strings[first_index]

# Skip over leading and trailing blanks
if strip:
for k0 in range(len(first)):
if first[k0:k0+1] != b' ':
break
for k1 in range(len(first)):
if first[-k1-1:-k1] != b' ':
break

lstring = len(first) - k1
else:
k0 = 0
k1 = 0
lstring = len(first)

if validate:
merged = bytearray(strings)
if b" " in merged or b"-" in merged:
raise ValueError("blank character in ISO time")

# Prepare a dictionary to define the string format
dtype_dict = {"h": ("|S2",0),
"m": ("|S2",3),
"s": ("|S2",6)}
dtype_dict = {"h": ("|S2", k0 + 0),
"m": ("|S2", k0 + 3),
"s": ("|S2", k0 + 6)}

if validate:
dtype_dict["colon1"] = ("|S1", 2)
dtype_dict["colon2"] = ("|S1", 5)
dtype_dict["colon1"] = ("|S1", k0 + 2)
dtype_dict["colon2"] = ("|S1", k0 + 5)

# Get the first string. Every subsequent string is assumed to match in
# format.
first = strings.ravel()[0]
lstring = len(first)
if k0 > 0:
dtype_dict["white0"] = ("|S" + str(k0), 0)

if k1 > 0:
dtype_dict["white1"] = ("|S" + str(k1), lstring)

# Check for a trailing "Z" to ignore
has_z = (first[-1:] == b"Z") # note first[-1] is an int in Python 3,
Expand All @@ -1569,15 +1602,16 @@ def sec_from_iso(strings, validate=True):
# Check for a period
has_dot = (lstring > 8)
if has_dot:
dtype_dict["dot"] = ("|S1", 8)
dtype_dict["dot"] = ("|S1", k0 + 8)

# Check for fractional seconds
lfrac = lstring - 9
lfrac = lstring - 9 - k0
has_frac = lfrac > 0
if has_frac:
dtype_dict["f"] = ("|S" + str(lfrac), 9)
dtype_dict["f"] = ("|S" + str(lfrac), k0 + 9)

# Extract hours, minutes, seconds
dtype = np.dtype(dtype_dict)
strings.dtype = np.dtype(dtype_dict)
h = strings["h"].astype("int")
m = strings["m"].astype("int")
Expand All @@ -1595,6 +1629,14 @@ def sec_from_iso(strings, validate=True):
np.any(strings["colon2"] != b":")):
raise ValueError("invalid ISO time punctuation")

if "white1" in dtype_dict:
if np.any(strings["white0"] != k0 * b" "):
raise ValueError("invalid ISO time punctuation")

if "white1" in dtype_dict:
if np.any(strings["white1"] != k1 * b" "):
raise ValueError("invalid ISO time punctuation")

if has_z:
if np.any(strings["z"] != b"Z"):
raise ValueError("invalid ISO time punctuation")
Expand All @@ -1621,7 +1663,7 @@ def sec_from_iso(strings, validate=True):

########################################

def day_sec_from_iso(strings, validate=True):
def day_sec_from_iso(strings, validate=True, strip=False):
"""Returns a day and second based on a parsing of the string in ISO
date-time format. The format is strictly enforced to be an ISO date plus an
ISO time, separated by a single space or a "T"."""
Expand All @@ -1639,30 +1681,38 @@ def day_sec_from_iso(strings, validate=True):

strings = np.array(strings).astype('S')

# Check for a T or blank separator
first = strings.ravel()[0]
first_index = len(strings.shape) * (0,)
first = strings[first_index]

# Check for a T or blank separator
csep = b"T"
isep = first.find(csep)

if isep == -1:
if strip:
for k0 in range(len(first)):
if first[k0:k0+1] != b' ':
break
else:
k0 = 0

csep = b" "
isep = first.find(csep)
isep = first.find(csep, k0)

# If no separator is found, assume it is just a date
if isep == -1:
return (day_from_iso(strings, validate), 0)
return (day_from_iso(strings, validate, strip), 0)

# Otherwise, parse the date and time separately
dtype_dict = {"day": ("|S" + str(isep), 0),
"sec": ("|S" + str(len(first) - isep - 1), isep+1)}
"sec": ("|S" + str(len(first) - isep - 1), isep + 1)}

if validate:
dtype_dict["sep"] = ("|S1", isep)

strings.dtype = np.dtype(dtype_dict)
day = day_from_iso(strings["day"], validate)
sec = sec_from_iso(strings["sec"], validate)
day = day_from_iso(strings["day"], validate, strip)
sec = sec_from_iso(strings["sec"], validate, strip)

if validate:
if (np.any(strings["sep"] != csep)):
Expand All @@ -1675,11 +1725,11 @@ def day_sec_from_iso(strings, validate=True):

########################################

def tai_from_iso(strings, validate=True):
def tai_from_iso(strings, validate=True, strip=False):
"""Returns the elapsed seconds TAI from January 1, 2000 given an ISO date
or date-time string. Works for scalars or arrays."""

(day, sec) = day_sec_from_iso(strings, validate)
(day, sec) = day_sec_from_iso(strings, validate, strip)
return tai_from_day(day) + sec

########################################
Expand Down
Loading

0 comments on commit 5f4aa89

Please sign in to comment.