diff --git a/pytides/astro.py b/pytides/astro.py index b9d8ecf..9c5b730 100644 --- a/pytides/astro.py +++ b/pytides/astro.py @@ -8,201 +8,221 @@ # to make a significant difference to the resulting accuracy of harmonic # analysis. -#Convert a sexagesimal angle into decimal degrees -def s2d(degrees, arcmins = 0, arcsecs = 0, mas = 0, muas = 0): - return ( - degrees - + (arcmins / 60.0) - + (arcsecs / (60.0*60.0)) - + (mas / (60.0*60.0*1e3)) - + (muas / (60.0*60.0*1e6)) - ) - -#Evaluate a polynomial at argument + +def s2d(degrees, arcmins=0, arcsecs=0, mas=0, muas=0): + # Convert a sexagesimal angle into decimal degrees + return ( + degrees + + (arcmins / 60.0) + + (arcsecs / (60.0*60.0)) + + (mas / (60.0*60.0*1e3)) + + (muas / (60.0*60.0*1e6)) + ) + + def polynomial(coefficients, argument): - return sum([c * (argument ** i) for i,c in enumerate(coefficients)]) + # Evaluate a polynomial at argument + return sum([c * (argument ** i) for i, c in enumerate(coefficients)]) + -#Evaluate the first derivative of a polynomial at argument def d_polynomial(coefficients, argument): - return sum([c * i * (argument ** (i-1)) for i,c in enumerate(coefficients)]) + # Evaluate the first derivative of a polynomial at argument + return sum([c * i * (argument ** (i - 1)) + for i, c in enumerate(coefficients)]) + -#Meeus formula 11.1 def T(t): - return (JD(t) - 2451545.0)/36525 + # Meeus formula 11.1 + return (JD(t) - 2451545.0)/36525 + -#Meeus formula 7.1 def JD(t): - Y, M = t.year, t.month - D = ( - t.day - + t.hour / (24.0) - + t.minute / (24.0*60.0) - + t.second / (24.0*60.0*60.0) - + t.microsecond / (24.0 * 60.0 * 60.0 * 1e6) - ) - if M <= 2: - Y = Y - 1 - M = M + 12 - A = np.floor(Y / 100.0) - B = 2 - A + np.floor(A / 4.0) - return np.floor(365.25*(Y+4716)) + np.floor(30.6001*(M+1)) + D + B - 1524.5 - -#Meeus formula 21.3 + # Meeus formula 7.1 + Y, M = t.year, t.month + D = ( + t.day + + t.hour / (24.0) + + t.minute / (24.0*60.0) + + t.second / (24.0*60.0*60.0) + + t.microsecond / (24.0 * 60.0 * 60.0 * 1e6) + ) + if M <= 2: + Y = Y - 1 + M = M + 12 + A = np.floor(Y / 100.0) + B = 2 - A + np.floor(A / 4.0) + return np.floor(365.25*(Y+4716)) + np.floor(30.6001*(M+1)) + D + B - 1524.5 + +# Meeus formula 21.3 terrestrial_obliquity_coefficients = ( - s2d(23,26,21.448), - -s2d(0,0,4680.93), - -s2d(0,0,1.55), - s2d(0,0,1999.25), - -s2d(0,0,51.38), - -s2d(0,0,249.67), - -s2d(0,0,39.05), - s2d(0,0,7.12), - s2d(0,0,27.87), - s2d(0,0,5.79), - s2d(0,0,2.45) + s2d(23, 26, 21.448), + -s2d(0, 0, 4680.93), + -s2d(0, 0, 1.55), + s2d(0, 0, 1999.25), + -s2d(0, 0, 51.38), + -s2d(0, 0, 249.67), + -s2d(0, 0, 39.05), + s2d(0, 0, 7.12), + s2d(0, 0, 27.87), + s2d(0, 0, 5.79), + s2d(0, 0, 2.45) ) -#Adjust these coefficients for parameter T rather than U +# Adjust these coefficients for parameter T rather than U terrestrial_obliquity_coefficients = [ - c * (1e-2) ** i for i,c in enumerate(terrestrial_obliquity_coefficients) + c * (1e-2) ** i for i, c in enumerate(terrestrial_obliquity_coefficients) ] -#Not entirely sure about this interpretation, but this is the difference -#between Meeus formulae 24.2 and 24.3 and seems to work +# Not entirely sure about this interpretation, but this is the difference +# between Meeus formulae 24.2 and 24.3 and seems to work solar_perigee_coefficients = ( - 280.46645 - 357.52910, - 36000.76932 - 35999.05030, - 0.0003032 + 0.0001559, - 0.00000048 + 280.46645 - 357.52910, + 36000.76932 - 35999.05030, + 0.0003032 + 0.0001559, + 0.00000048 ) -#Meeus formula 24.2 +# Meeus formula 24.2 solar_longitude_coefficients = ( - 280.46645, - 36000.76983, - 0.0003032 + 280.46645, + 36000.76983, + 0.0003032 ) -#This value is taken from JPL Horizon and is essentially constant +# This value is taken from JPL Horizon and is essentially constant lunar_inclination_coefficients = ( - 5.145, + 5.145, ) -#Meeus formula 45.1 +# Meeus formula 45.1 lunar_longitude_coefficients = ( - 218.3164591, - 481267.88134236, - -0.0013268, - 1/538841.0 - -1/65194000.0 + 218.3164591, + 481267.88134236, + -0.0013268, + 1 / 538841.0 + -1 / 65194000.0 ) -#Meeus formula 45.7 +# Meeus formula 45.7 lunar_node_coefficients = ( - 125.0445550, - -1934.1361849, - 0.0020762, - 1/467410.0, - -1/60616000.0 + 125.0445550, + -1934.1361849, + 0.0020762, + 1/467410.0, + -1/60616000.0 ) -#Meeus, unnumbered formula directly preceded by 45.7 +# Meeus, unnumbered formula directly preceded by 45.7 lunar_perigee_coefficients = ( - 83.3532430, - 4069.0137111, - -0.0103238, - -1/80053.0, - 1/18999000.0 + 83.3532430, + 4069.0137111, + -0.0103238, + -1/80053.0, + 1/18999000.0 ) -#Now follow some useful auxiliary values, we won't need their speed. -#See notes on Table 6 in Schureman for I, nu, xi, nu', 2nu'' +# Now follow some useful auxiliary values, we won't need their speed. +# See notes on Table 6 in Schureman for I, nu, xi, nu', 2nu'' + + def _I(N, i, omega): - N, i, omega = d2r * N, d2r*i, d2r*omega - cosI = np.cos(i)*np.cos(omega)-np.sin(i)*np.sin(omega)*np.cos(N) - return r2d*np.arccos(cosI) + N, i, omega = d2r * N, d2r*i, d2r*omega + cosI = np.cos(i) * np.cos(omega) - np.sin(i) * np.sin(omega) * np.cos(N) + return r2d * np.arccos(cosI) + def _xi(N, i, omega): - N, i, omega = d2r * N, d2r*i, d2r*omega - e1 = np.cos(0.5*(omega-i))/np.cos(0.5*(omega+i)) * np.tan(0.5*N) - e2 = np.sin(0.5*(omega-i))/np.sin(0.5*(omega+i)) * np.tan(0.5*N) - e1, e2 = np.arctan(e1), np.arctan(e2) - e1, e2 = e1 - 0.5*N, e2 - 0.5*N - return -(e1 + e2)*r2d + N, i, omega = d2r * N, d2r*i, d2r*omega + e1 = (np.cos(0.5 * (omega - i)) + / np.cos(0.5 * (omega + i)) * np.tan(0.5 * N)) + e2 = (np.sin(0.5 * (omega - i)) + / np.sin(0.5 * (omega + i)) * np.tan(0.5 * N)) + e1, e2 = np.arctan(e1), np.arctan(e2) + e1, e2 = e1 - 0.5*N, e2 - 0.5 * N + return -(e1 + e2)*r2d + def _nu(N, i, omega): - N, i, omega = d2r * N, d2r*i, d2r*omega - e1 = np.cos(0.5*(omega-i))/np.cos(0.5*(omega+i)) * np.tan(0.5*N) - e2 = np.sin(0.5*(omega-i))/np.sin(0.5*(omega+i)) * np.tan(0.5*N) - e1, e2 = np.arctan(e1), np.arctan(e2) - e1, e2 = e1 - 0.5*N, e2 - 0.5*N - return (e1 - e2)*r2d - -#Schureman equation 224 -#Can we be more precise than B "the solar coefficient" = 0.1681? + N, i, omega = d2r * N, d2r*i, d2r*omega + e1 = (np.cos(0.5 * (omega - i)) + / np.cos(0.5 * (omega + i)) * np.tan(0.5 * N)) + e2 = (np.sin(0.5 * (omega - i)) + / np.sin(0.5 * (omega + i)) + * np.tan(0.5 * N)) + e1, e2 = np.arctan(e1), np.arctan(e2) + e1, e2 = e1 - 0.5 * N, e2 - 0.5 * N + return (e1 - e2) * r2d + + +# Can we be more precise than B "the solar coefficient" = 0.1681? def _nup(N, i, omega): - I = d2r * _I(N, i, omega) - nu = d2r * _nu(N, i, omega) - return r2d * np.arctan(np.sin(2*I)*np.sin(nu)/(np.sin(2*I)*np.cos(nu)+0.3347)) + # Schureman equation 224 + I = d2r * _I(N, i, omega) + nu = d2r * _nu(N, i, omega) + return (r2d * np.arctan( + np.sin(2 * I) * np.sin(nu) / (np.sin(2 * I) * np.cos(nu) + 0.3347)) + ) + -#Schureman equation 232 def _nupp(N, i, omega): - I = d2r * _I(N, i, omega) - nu = d2r * _nu(N, i, omega) - tan2nupp = (np.sin(I)**2*np.sin(2*nu))/(np.sin(I)**2*np.cos(2*nu)+0.0727) - return r2d * 0.5 * np.arctan(tan2nupp) + # Schureman equation 232 + I = d2r * _I(N, i, omega) + nu = d2r * _nu(N, i, omega) + tan2nupp = (np.sin(I)**2*np.sin(2*nu))/(np.sin(I)**2*np.cos(2*nu)+0.0727) + return r2d * 0.5 * np.arctan(tan2nupp) AstronomicalParameter = namedtuple('AstronomicalParameter', ['value', 'speed']) + def astro(t): - a = {} - #We can use polynomial fits from Meeus to obtain good approximations to - #some astronomical values (and therefore speeds). - polynomials = { - 's': lunar_longitude_coefficients, - 'h': solar_longitude_coefficients, - 'p': lunar_perigee_coefficients, - 'N': lunar_node_coefficients, - 'pp': solar_perigee_coefficients, - '90': (90.0,), - 'omega': terrestrial_obliquity_coefficients, - 'i': lunar_inclination_coefficients - } - #Polynomials are in T, that is Julian Centuries; we want our speeds to be - #in the more convenient unit of degrees per hour. - dT_dHour = 1 / (24 * 365.25 * 100) - for name, coefficients in polynomials.items(): - a[name] = AstronomicalParameter( - np.mod(polynomial(coefficients, T(t)), 360.0), - d_polynomial(coefficients, T(t)) * dT_dHour - ) - - #Some other parameters defined by Schureman which are dependent on the - #parameters N, i, omega for use in node factor calculations. We don't need - #their speeds. - args = list(each.value for each in [a['N'], a['i'], a['omega']]) - for name, function in { - 'I': _I, - 'xi': _xi, - 'nu': _nu, - 'nup': _nup, - 'nupp': _nupp - }.items(): - a[name] = AstronomicalParameter(np.mod(function(*args), 360.0), None) - - #We don't work directly with the T (hours) parameter, instead our spanning - #set for equilibrium arguments #is given by T+h-s, s, h, p, N, pp, 90. - #This is in line with convention. - hour = AstronomicalParameter((JD(t) - np.floor(JD(t))) * 360.0, 15.0) - a['T+h-s'] = AstronomicalParameter( - hour.value + a['h'].value - a['s'].value, - hour.speed + a['h'].speed - a['s'].speed - ) - #It is convenient to calculate Schureman's P here since several node - #factors need it, although it could be argued that these - #(along with I, xi, nu etc) belong somewhere else. - a['P'] = AstronomicalParameter( - np.mod(a['p'].value -a['xi'].value,360.0), - None - ) - return a + a = {} + # We can use polynomial fits from Meeus to obtain good approximations to + # some astronomical values (and therefore speeds). + polynomials = { + 's': lunar_longitude_coefficients, + 'h': solar_longitude_coefficients, + 'p': lunar_perigee_coefficients, + 'N': lunar_node_coefficients, + 'pp': solar_perigee_coefficients, + '90': (90.0,), + 'omega': terrestrial_obliquity_coefficients, + 'i': lunar_inclination_coefficients + } + # Polynomials are in T, that is Julian Centuries; we want our speeds to be + # in the more convenient unit of degrees per hour. + dT_dHour = 1 / (24 * 365.25 * 100) + for name, coefficients in polynomials.items(): + a[name] = AstronomicalParameter( + np.mod(polynomial(coefficients, T(t)), 360.0), + d_polynomial(coefficients, T(t)) * dT_dHour + ) + + # Some other parameters defined by Schureman which are dependent on the + # parameters N, i, omega for use in node factor calculations. We don't need + # their speeds. + args = list(each.value for each in [a['N'], a['i'], a['omega']]) + for name, function in { + 'I': _I, + 'xi': _xi, + 'nu': _nu, + 'nup': _nup, + 'nupp': _nupp + }.items(): + a[name] = AstronomicalParameter(np.mod(function(*args), 360.0), None) + + # We don't work directly with the T (hours) parameter, instead our spanning + # set for equilibrium arguments # is given by T+h-s, s, h, p, N, pp, 90. + # This is in line with convention. + hour = AstronomicalParameter((JD(t) - np.floor(JD(t))) * 360.0, 15.0) + a['T+h-s'] = AstronomicalParameter( + hour.value + a['h'].value - a['s'].value, + hour.speed + a['h'].speed - a['s'].speed + ) + # It is convenient to calculate Schureman's P here since several node + # factors need it, although it could be argued that these + # (along with I, xi, nu etc) belong somewhere else. + a['P'] = AstronomicalParameter( + np.mod(a['p'].value - a['xi'].value, 360.0), + None + ) + return a diff --git a/pytides/constituent.py b/pytides/constituent.py index 3ff8197..6c6230a 100644 --- a/pytides/constituent.py +++ b/pytides/constituent.py @@ -4,157 +4,166 @@ import numpy as np import nodal_corrections as nc -class BaseConstituent(object): - xdo_int = { - 'A': 1, 'B': 2, 'C': 3, 'D': 4, 'E': 5, 'F': 6, 'G': 7, 'H': 8, 'I': 9, - 'J': 10, 'K': 11, 'L': 12, 'M': 13, 'N': 14, 'O': 15, 'P': 16, 'Q': 17, - 'R': -8, 'S': -7, 'T': -6, 'U': -5, 'V': -4, 'W': -3, 'X': -2, 'Y': -1, - 'Z': 0 - } - - int_xdo = {v:k for k, v in xdo_int.items()} - - def __init__(self, name, xdo='', coefficients=[], u=nc.u_zero, f=nc.f_unity): - if xdo == '': - self.coefficients = np.array(coefficients) - else: - self.coefficients = np.array(self.xdo_to_coefficients(xdo)) - self.name = name - self.u = u - self.f = f - - def xdo_to_coefficients(self, xdo): - return [self.xdo_int[l.upper()] for l in xdo if l in string.ascii_letters] - - def coefficients_to_xdo(self, coefficients): - return ''.join([self.int_xdo[c] for c in cooefficients]) - - def V(self, astro): - return np.dot(self.coefficients, self.astro_values(astro)) - - def xdo(self): - return self.coefficients_to_xdo(self.coefficients) - - def speed(self, a): - return np.dot(self.coefficients, self.astro_speeds(a)) - - def astro_xdo(self, a): - return [a['T+h-s'], a['s'], a['h'], a['p'], a['N'], a['pp'], a['90']] - - def astro_speeds(self, a): - return np.array([each.speed for each in self.astro_xdo(a)]) - - def astro_values(self, a): - return np.array([each.value for each in self.astro_xdo(a)]) - - #Consider two out of phase constituents which travel at the same speed to - #be identical - def __eq__(self, c): - return np.all(self.coefficients[:-1] == c.coefficients[:-1]) - def __hash__(self): - return hash(tuple(self.coefficients[:-1])) - -class CompoundConstituent(BaseConstituent): - - def __init__(self, members = [], **kwargs): - self.members = members - - if 'u' not in kwargs: - kwargs['u'] = self.u - if 'f' not in kwargs: - kwargs['f'] = self.f - - super(CompoundConstituent,self).__init__(**kwargs) - - self.coefficients = reduce(op.add,[c.coefficients * n for (c,n) in members]) - - def speed(self, a): - return reduce(op.add, [n * c.speed(a) for (c,n) in self.members]) - - def V(self, a): - return reduce(op.add, [n * c.V(a) for (c,n) in self.members]) - - def u(self, a): - return reduce(op.add, [n * c.u(a) for (c,n) in self.members]) +class BaseConstituent(object): + xdo_int = { + 'A': 1, 'B': 2, 'C': 3, 'D': 4, 'E': 5, 'F': 6, 'G': 7, 'H': 8, 'I': 9, + 'J': 10, 'K': 11, 'L': 12, 'M': 13, 'N': 14, 'O': 15, 'P': 16, 'Q': 17, + 'R': -8, 'S': -7, 'T': -6, 'U': -5, 'V': -4, 'W': -3, 'X': -2, 'Y': -1, + 'Z': 0 + } - def f(self, a): - return reduce(op.mul, [c.f(a) ** abs(n) for (c,n) in self.members]) + int_xdo = {v: k for k, v in xdo_int.items()} -###### Base Constituents -#Long Term -_Z0 = BaseConstituent(name = 'Z0', xdo = 'Z ZZZ ZZZ', u = nc.u_zero, f = nc.f_unity) -_Sa = BaseConstituent(name = 'Sa', xdo = 'Z ZAZ ZZZ', u = nc.u_zero, f = nc.f_unity) -_Ssa = BaseConstituent(name = 'Ssa', xdo = 'Z ZBZ ZZZ', u = nc.u_zero, f = nc.f_unity) -_Mm = BaseConstituent(name = 'Mm', xdo = 'Z AZY ZZZ', u = nc.u_zero, f = nc.f_Mm) -_Mf = BaseConstituent(name = 'Mf', xdo = 'Z BZZ ZZZ', u = nc.u_Mf, f = nc.f_Mf) + def __init__(self, name, xdo='', coefficients=[], + u=nc.u_zero, f=nc.f_unity): + if xdo == '': + self.coefficients = np.array(coefficients) + else: + self.coefficients = np.array(self.xdo_to_coefficients(xdo)) + self.name = name + self.u = u + self.f = f -#Diurnals -_Q1 = BaseConstituent(name = 'Q1', xdo = 'A XZA ZZA', u = nc.u_O1, f = nc.f_O1) -_O1 = BaseConstituent(name = 'O1', xdo = 'A YZZ ZZA', u = nc.u_O1, f = nc.f_O1) -_K1 = BaseConstituent(name = 'K1', xdo = 'A AZZ ZZY', u = nc.u_K1, f = nc.f_K1) -_J1 = BaseConstituent(name = 'J1', xdo = 'A BZY ZZY', u = nc.u_J1, f = nc.f_J1) + def xdo_to_coefficients(self, xdo): + return [self.xdo_int[l.upper()] + for l in xdo + if l in string.ascii_letters] -#M1 is a tricky business for reasons of convention, rather than theory. The -#reasons for this are best summarised by Schureman paragraphs 126, 127 and in -#the comments found in congen_input.txt of xtides, so I won't go over all this -#again here. + def coefficients_to_xdo(self, coefficients): + return ''.join([self.int_xdo[c] for c in cooefficients]) -_M1 = BaseConstituent(name = 'M1', xdo = 'A ZZZ ZZA', u = nc.u_M1, f = nc.f_M1) -_P1 = BaseConstituent(name = 'P1', xdo = 'A AXZ ZZA', u = nc.u_zero, f = nc.f_unity) -_S1 = BaseConstituent(name = 'S1', xdo = 'A AYZ ZZZ', u = nc.u_zero, f = nc.f_unity) -_OO1 = BaseConstituent(name = 'OO1', xdo = 'A CZZ ZZY', u = nc.u_OO1, f = nc.f_OO1) + def V(self, astro): + return np.dot(self.coefficients, self.astro_values(astro)) -#Semi-Diurnals -_2N2 = BaseConstituent(name = '2N2', xdo = 'B XZB ZZZ', u = nc.u_M2, f = nc.f_M2) -_N2 = BaseConstituent(name = 'N2', xdo = 'B YZA ZZZ', u = nc.u_M2, f = nc.f_M2) -_nu2 = BaseConstituent(name = 'nu2', xdo = 'B YBY ZZZ', u = nc.u_M2, f = nc.f_M2) -_M2 = BaseConstituent(name = 'M2', xdo = 'B ZZZ ZZZ', u = nc.u_M2, f = nc.f_M2) -_lambda2 = BaseConstituent(name = 'lambda2', xdo = 'B AXA ZZB', u = nc.u_M2, f = nc.f_M2) -_L2 = BaseConstituent(name = 'L2', xdo = 'B AZY ZZB', u = nc.u_L2, f = nc.f_L2) -_T2 = BaseConstituent(name = 'T2', xdo = 'B BWZ ZAZ', u = nc.u_zero, f = nc.f_unity) -_S2 = BaseConstituent(name = 'S2', xdo = 'B BXZ ZZZ', u = nc.u_zero, f = nc.f_unity) -_R2 = BaseConstituent(name = 'R2', xdo = 'B BYZ ZYB', u = nc.u_zero, f = nc.f_unity) -_K2 = BaseConstituent(name = 'K2', xdo = 'B BZZ ZZZ', u = nc.u_K2, f = nc.f_K2) + def xdo(self): + return self.coefficients_to_xdo(self.coefficients) -#Third-Diurnals -_M3 = BaseConstituent(name = 'M3', xdo = 'C ZZZ ZZZ', u = lambda a: nc.u_Modd(a,3), f = lambda a: nc.f_Modd(a,3)) + def speed(self, a): + return np.dot(self.coefficients, self.astro_speeds(a)) -###### Compound Constituents -#Long Term -_MSF = CompoundConstituent(name = 'MSF', members = [(_S2, 1), (_M2, -1)]) + def astro_xdo(self, a): + return [a['T+h-s'], a['s'], a['h'], a['p'], a['N'], a['pp'], a['90']] -#Diurnal -_2Q1 = CompoundConstituent(name = '2Q1', members = [(_N2, 1), (_J1, -1)]) -_rho1 = CompoundConstituent(name = 'rho1', members = [(_nu2, 1), (_K1, -1)]) + def astro_speeds(self, a): + return np.array([each.speed for each in self.astro_xdo(a)]) -#Semi-Diurnal + def astro_values(self, a): + return np.array([each.value for each in self.astro_xdo(a)]) -_mu2 = CompoundConstituent(name = 'mu2', members = [(_M2, 2), (_S2, -1)]) #2MS2 -_2SM2 = CompoundConstituent(name = '2SM2', members = [(_S2, 2), (_M2, -1)]) + # Consider two out of phase constituents which travel at the same speed to + # be identical + def __eq__(self, c): + return np.all(self.coefficients[:-1] == c.coefficients[:-1]) -#Third-Diurnal -_2MK3 = CompoundConstituent(name = '2MK3', members = [(_M2, 1), (_O1, 1)]) -_MK3 = CompoundConstituent(name = 'MK3', members = [(_M2, 1), (_K1, 1)]) + def __hash__(self): + return hash(tuple(self.coefficients[:-1])) -#Quarter-Diurnal -_MN4 = CompoundConstituent(name = 'MN4', members = [(_M2, 1), (_N2, 1)]) -_M4 = CompoundConstituent(name = 'M4', members = [(_M2, 2)]) -_MS4 = CompoundConstituent(name = 'MS4', members = [(_M2, 1), (_S2, 1)]) -_S4 = CompoundConstituent(name = 'S4', members = [(_S2, 2)]) -#Sixth-Diurnal -_M6 = CompoundConstituent(name = 'M6', members = [(_M2, 3)]) -_S6 = CompoundConstituent(name = 'S6', members = [(_S2, 3)]) +class CompoundConstituent(BaseConstituent): -#Eighth-Diurnals -_M8 = CompoundConstituent(name = 'M8', members = [(_M2, 4)]) + def __init__(self, members=[], **kwargs): + self.members = members + + if 'u' not in kwargs: + kwargs['u'] = self.u + if 'f' not in kwargs: + kwargs['f'] = self.f + + super(CompoundConstituent, self).__init__(**kwargs) + + self.coefficients = reduce( + op.add, + [c.coefficients * n for (c, n) in members]) + + def speed(self, a): + return reduce(op.add, [n * c.speed(a) for (c, n) in self.members]) + + def V(self, a): + return reduce(op.add, [n * c.V(a) for (c, n) in self.members]) + + def u(self, a): + return reduce(op.add, [n * c.u(a) for (c, n) in self.members]) + + def f(self, a): + return reduce(op.mul, [c.f(a) ** abs(n) for (c, n) in self.members]) + +# ## ## # Base Constituents +# Long Term +_Z0 = BaseConstituent(name='Z0', xdo='Z ZZZ ZZZ', u=nc.u_zero, f=nc.f_unity) +_Sa = BaseConstituent(name='Sa', xdo='Z ZAZ ZZZ', u=nc.u_zero, f=nc.f_unity) +_Ssa = BaseConstituent(name='Ssa', xdo='Z ZBZ ZZZ', u=nc.u_zero, f=nc.f_unity) +_Mm = BaseConstituent(name='Mm', xdo='Z AZY ZZZ', u=nc.u_zero, f=nc.f_Mm) +_Mf = BaseConstituent(name='Mf', xdo='Z BZZ ZZZ', u=nc.u_Mf, f=nc.f_Mf) + +# Diurnals +_Q1 = BaseConstituent(name='Q1', xdo='A XZA ZZA', u=nc.u_O1, f=nc.f_O1) +_O1 = BaseConstituent(name='O1', xdo='A YZZ ZZA', u=nc.u_O1, f=nc.f_O1) +_K1 = BaseConstituent(name='K1', xdo='A AZZ ZZY', u=nc.u_K1, f=nc.f_K1) +_J1 = BaseConstituent(name='J1', xdo='A BZY ZZY', u=nc.u_J1, f=nc.f_J1) + +# M1 is a tricky business for reasons of convention, rather than theory. The +# reasons for this are best summarised by Schureman paragraphs 126, 127 and in +# the comments found in congen_input.txt of xtides, so I won't go over all this +# again here. + +_M1 = BaseConstituent(name='M1', xdo='A ZZZ ZZA', u=nc.u_M1, f=nc.f_M1) +_P1 = BaseConstituent(name='P1', xdo='A AXZ ZZA', u=nc.u_zero, f=nc.f_unity) +_S1 = BaseConstituent(name='S1', xdo='A AYZ ZZZ', u=nc.u_zero, f=nc.f_unity) +_OO1 = BaseConstituent(name='OO1', xdo='A CZZ ZZY', u=nc.u_OO1, f=nc.f_OO1) + +# Semi-Diurnals +_2N2 = BaseConstituent(name='2N2', xdo='B XZB ZZZ', u=nc.u_M2, f=nc.f_M2) +_N2 = BaseConstituent(name='N2', xdo='B YZA ZZZ', u=nc.u_M2, f=nc.f_M2) +_nu2 = BaseConstituent(name='nu2', xdo='B YBY ZZZ', u=nc.u_M2, f=nc.f_M2) +_M2 = BaseConstituent(name='M2', xdo='B ZZZ ZZZ', u=nc.u_M2, f=nc.f_M2) +_lambda2 = BaseConstituent(name='lambda2', xdo='B AXA ZZB', + u=nc.u_M2, f=nc.f_M2) +_L2 = BaseConstituent(name='L2', xdo='B AZY ZZB', u=nc.u_L2, f=nc.f_L2) +_T2 = BaseConstituent(name='T2', xdo='B BWZ ZAZ', u=nc.u_zero, f=nc.f_unity) +_S2 = BaseConstituent(name='S2', xdo='B BXZ ZZZ', u=nc.u_zero, f=nc.f_unity) +_R2 = BaseConstituent(name='R2', xdo='B BYZ ZYB', u=nc.u_zero, f=nc.f_unity) +_K2 = BaseConstituent(name='K2', xdo='B BZZ ZZZ', u=nc.u_K2, f=nc.f_K2) + +# Third-Diurnals +_M3 = BaseConstituent(name='M3', xdo='C ZZZ ZZZ', + u=lambda a: nc.u_Modd(a, 3), + f=lambda a: nc.f_Modd(a, 3)) + +# ## ## # Compound Constituents +# Long Term +_MSF = CompoundConstituent(name='MSF', members=[(_S2, 1), (_M2, -1)]) + +# Diurnal +_2Q1 = CompoundConstituent(name='2Q1', members=[(_N2, 1), (_J1, -1)]) +_rho1 = CompoundConstituent(name='rho1', members=[(_nu2, 1), (_K1, -1)]) + +# Semi-Diurnal + +_mu2 = CompoundConstituent(name='mu2', members=[(_M2, 2), (_S2, -1)]) # 2MS2 +_2SM2 = CompoundConstituent(name='2SM2', members=[(_S2, 2), (_M2, -1)]) + +# Third-Diurnal +_2MK3 = CompoundConstituent(name='2MK3', members=[(_M2, 1), (_O1, 1)]) +_MK3 = CompoundConstituent(name='MK3', members=[(_M2, 1), (_K1, 1)]) + +# Quarter-Diurnal +_MN4 = CompoundConstituent(name='MN4', members=[(_M2, 1), (_N2, 1)]) +_M4 = CompoundConstituent(name='M4', members=[(_M2, 2)]) +_MS4 = CompoundConstituent(name='MS4', members=[(_M2, 1), (_S2, 1)]) +_S4 = CompoundConstituent(name='S4', members=[(_S2, 2)]) + +# Sixth-Diurnal +_M6 = CompoundConstituent(name='M6', members=[(_M2, 3)]) +_S6 = CompoundConstituent(name='S6', members=[(_S2, 3)]) + +# Eighth-Diurnals +_M8 = CompoundConstituent(name='M8', members=[(_M2, 4)]) noaa = [ - _M2, _S2, _N2, _K1, _M4, _O1, _M6, _MK3, _S4, _MN4, _nu2, _S6, _mu2, - _2N2, _OO1, _lambda2, _S1, _M1, _J1, _Mm, _Ssa, _Sa, _MSF, _Mf, - _rho1, _Q1, _T2, _R2, _2Q1, _P1, _2SM2, _M3, _L2, _2MK3, _K2, - _M8, _MS4 + _M2, _S2, _N2, _K1, _M4, _O1, _M6, _MK3, _S4, _MN4, _nu2, _S6, _mu2, + _2N2, _OO1, _lambda2, _S1, _M1, _J1, _Mm, _Ssa, _Sa, _MSF, _Mf, + _rho1, _Q1, _T2, _R2, _2Q1, _P1, _2SM2, _M3, _L2, _2MK3, _K2, + _M8, _MS4 ] - diff --git a/pytides/nodal_corrections.py b/pytides/nodal_corrections.py index 170f0d6..4114d81 100644 --- a/pytides/nodal_corrections.py +++ b/pytides/nodal_corrections.py @@ -2,140 +2,171 @@ import numpy as np d2r, r2d = np.pi/180.0, 180.0/np.pi -#The following functions take a dictionary of astronomical values (in degrees) -#and return dimensionless scale factors for constituent amplitudes. +# The following functions take a dictionary of astronomical values (in degrees) +# and return dimensionless scale factors for constituent amplitudes. + def f_unity(a): - return 1.0 + return 1.0 + -#Schureman equations 73, 65 def f_Mm(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - mean = (2/3.0 - np.sin(omega)**2)*(1 - 3/2.0 * np.sin(i)**2) - return (2/3.0 - np.sin(I)**2) / mean + # Schureman equations 73, 65 + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + mean = (2/3.0 - np.sin(omega)**2)*(1 - 3/2.0 * np.sin(i)**2) + return (2/3.0 - np.sin(I)**2) / mean + -#Schureman equations 74, 66 def f_Mf(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - mean = np.sin(omega)**2 * np.cos(0.5*i)**4 - return np.sin(I)**2 / mean + # Schureman equations 74, 66 + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + mean = np.sin(omega)**2 * np.cos(0.5*i)**4 + return np.sin(I)**2 / mean + -#Schureman equations 75, 67 def f_O1(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - mean = np.sin(omega) * np.cos(0.5*omega)**2 * np.cos(0.5*i)**4 - return (np.sin(I) * np.cos(0.5*I)**2) / mean + # Schureman equations 75, 67 + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + mean = np.sin(omega) * np.cos(0.5*omega)**2 * np.cos(0.5*i)**4 + return (np.sin(I) * np.cos(0.5*I)**2) / mean + -#Schureman equations 76, 68 def f_J1(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - mean = np.sin(2*omega) * (1-3/2.0 * np.sin(i)**2) - return np.sin(2*I) / mean + # Schureman equations 76, 68 + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + mean = np.sin(2*omega) * (1-3/2.0 * np.sin(i)**2) + return np.sin(2*I) / mean + -#Schureman equations 77, 69 def f_OO1(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - mean = np.sin(omega) * np.sin(0.5*omega)**2 * np.cos(0.5*i)**4 - return np.sin(I) * np.sin(0.5*I)**2 / mean + # Schureman equations 77, 69 + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + mean = np.sin(omega) * np.sin(0.5*omega)**2 * np.cos(0.5*i)**4 + return np.sin(I) * np.sin(0.5*I)**2 / mean + -#Schureman equations 78, 70 def f_M2(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - mean = np.cos(0.5*omega)**4 * np.cos(0.5*i)**4 - return np.cos(0.5*I)**4 / mean - -#Schureman equations 227, 226, 68 -#Should probably eventually include the derivations of the magic numbers (0.5023 etc). + # Schureman equations 78, 70 + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + mean = np.cos(0.5*omega)**4 * np.cos(0.5*i)**4 + return np.cos(0.5*I)**4 / mean + + def f_K1(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - nu = d2r*a['nu'].value - sin2Icosnu_mean = np.sin(2*omega) * (1-3/2.0 * np.sin(i)**2) - mean = 0.5023*sin2Icosnu_mean + 0.1681 - return (0.2523*np.sin(2*I)**2 + 0.1689*np.sin(2*I)*np.cos(nu)+0.0283)**(0.5) / mean - -#Schureman equations 215, 213, 204 -#It can be (and has been) confirmed that the exponent for R_a reads 1/2 via Schureman Table 7 + # Schureman equations 227, 226, 68 + # Should probably eventually include the derivations of the magic numbers + # (0.5023 etc). + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + nu = d2r*a['nu'].value + sin2Icosnu_mean = np.sin(2*omega) * (1-3/2.0 * np.sin(i)**2) + mean = 0.5023*sin2Icosnu_mean + 0.1681 + return (0.2523 * np.sin(2 * I) ** 2 + 0.1689 * np.sin(2 * I) * np.cos(nu) + + 0.0283) ** (0.5) / mean + + def f_L2(a): - P = d2r*a['P'].value - I = d2r*a['I'].value - R_a_inv = (1 - 12*np.tan(0.5*I)**2 * np.cos(2*P)+36*np.tan(0.5*I)**4)**(0.5) - return f_M2(a) * R_a_inv + # Schureman equations 215, 213, 204 + # It can be (and has been) confirmed that the exponent for R_a reads 1/2 + # via Schureman Table 7 + P = d2r*a['P'].value + I = d2r*a['I'].value + R_a_inv = (1 - 12 * np.tan(0.5 * I) ** 2 * np.cos(2 * P) + 36 + * np.tan(0.5 * I) ** 4) ** (0.5) + return f_M2(a) * R_a_inv + -#Schureman equations 235, 234, 71 -#Again, magic numbers def f_K2(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - nu = d2r*a['nu'].value - sinsqIcos2nu_mean = np.sin(omega)**2 * (1-3/2.0 * np.sin(i)**2) - mean = 0.5023*sinsqIcos2nu_mean + 0.0365 - return (0.2533*np.sin(I)**4 + 0.0367*np.sin(I)**2 *np.cos(2*nu)+0.0013)**(0.5) / mean - -#Schureman equations 206, 207, 195 + # Schureman equations 235, 234, 71 + # Again, magic numbers + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + nu = d2r*a['nu'].value + sinsqIcos2nu_mean = np.sin(omega)**2 * (1-3/2.0 * np.sin(i)**2) + mean = 0.5023 * sinsqIcos2nu_mean + 0.0365 + return (0.2533 * (np.sin(I) ** 4) + 0.0367 * np.sin(I) ** 2 + * np.cos(2 * nu) + 0.0013) ** (0.5) / mean + + def f_M1(a): - P = d2r*a['P'].value - I = d2r*a['I'].value - Q_a_inv = (0.25 + 1.5*np.cos(I)*np.cos(2*P)*np.cos(0.5*I)**(-0.5) + 2.25*np.cos(I)**2 * np.cos(0.5*I)**(-4))**(0.5) - return f_O1(a) * Q_a_inv + # Schureman equations 206, 207, 195 + P = d2r*a['P'].value + I = d2r*a['I'].value + Q_a_inv = (0.25 + 1.5 * np.cos(I) * np.cos(2 * P) + * np.cos(0.5 * I) ** (-0.5) + 2.25 * np.cos(I) ** 2 + * np.cos(0.5 * I) ** (-4)) ** (0.5) + return f_O1(a) * Q_a_inv + -#See e.g. Schureman equation 149 def f_Modd(a, n): - return f_M2(a) ** (n / 2.0) + # See e.g. Schureman equation 149 + return f_M2(a) ** (n / 2.0) + +# Node factors u, see Table 2 of Schureman. -#Node factors u, see Table 2 of Schureman. def u_zero(a): - return 0.0 + return 0.0 + def u_Mf(a): - return -2.0 * a['xi'].value + return -2.0 * a['xi'].value + def u_O1(a): - return 2.0 * a['xi'].value - a['nu'].value + return 2.0 * a['xi'].value - a['nu'].value + def u_J1(a): - return -a['nu'].value + return -a['nu'].value + def u_OO1(a): - return -2.0 * a['xi'].value - a['nu'].value + return -2.0 * a['xi'].value - a['nu'].value + def u_M2(a): - return 2.0 * a['xi'].value - 2.0 * a['nu'].value + return 2.0 * a['xi'].value - 2.0 * a['nu'].value + def u_K1(a): - return -a['nup'].value + return -a['nup'].value + -#Schureman 214 def u_L2(a): - I = d2r*a['I'].value - P = d2r*a['P'].value - R = r2d*np.arctan(np.sin(2*P)/(1/6.0 * np.tan(0.5*I) **(-2) -np.cos(2*P))) - return 2.0 * a['xi'].value - 2.0 * a['nu'].value - R + # Schureman 214 + I = d2r * a['I'].value + P = d2r * a['P'].value + R = r2d * np.arctan( + np.sin(2 * P) / (1 / 6.0 * np.tan(0.5 * I) ** (-2) - np.cos(2 * P))) + return 2.0 * a['xi'].value - 2.0 * a['nu'].value - R + def u_K2(a): - return -2.0 * a['nupp'].value + return -2.0 * a['nupp'].value + -#Schureman 202 def u_M1(a): - I = d2r*a['I'].value - P = d2r*a['P'].value - Q = r2d*np.arctan((5*np.cos(I)-1)/(7*np.cos(I)+1)*np.tan(P)) - return a['xi'].value - a['nu'].value + Q + # Schureman 202 + I = d2r*a['I'].value + P = d2r*a['P'].value + Q = r2d*np.arctan((5*np.cos(I)-1)/(7*np.cos(I)+1)*np.tan(P)) + return a['xi'].value - a['nu'].value + Q + def u_Modd(a, n): - return n/2.0 * u_M2(a) + return n/2.0 * u_M2(a) diff --git a/pytides/tide.py b/pytides/tide.py index bc7089c..34fb8a0 100644 --- a/pytides/tide.py +++ b/pytides/tide.py @@ -2,10 +2,10 @@ from collections import OrderedDict, Iterable from itertools import takewhile, count try: - from itertools import izip, ifilter -except ImportError: #Python3 - izip = zip - ifilter = filter + from itertools import izip, ifilter +except ImportError: # Python3 + izip = zip + ifilter = filter from datetime import datetime, timedelta import numpy as np from scipy.optimize import leastsq, fsolve @@ -14,393 +14,435 @@ d2r, r2d = np.pi/180.0, 180.0/np.pi + class Tide(object): - dtype = np.dtype([ - ('constituent', object), - ('amplitude', float), - ('phase', float)]) - - def __init__( - self, - constituents = None, - amplitudes = None, - phases = None, - model = None, - radians = False - ): - """ - Initialise a tidal model. Provide constituents, amplitudes and phases OR a model. - Arguments: - constituents -- list of constituents used in the model. - amplitudes -- list of amplitudes corresponding to constituents - phases -- list of phases corresponding to constituents - model -- an ndarray of type Tide.dtype representing the constituents, amplitudes and phases. - radians -- boolean representing whether phases are in radians (default False) - """ - if None not in [constituents, amplitudes, phases]: - if len(constituents) == len(amplitudes) == len(phases): - model = np.zeros(len(phases), dtype=Tide.dtype) - model['constituent'] = np.array(constituents) - model['amplitude'] = np.array(amplitudes) - model['phase'] = np.array(phases) - else: - raise ValueError("Constituents, amplitudes and phases should all be arrays of equal length.") - elif model is not None: - if not model.dtype == Tide.dtype: - raise ValueError("Model must be a numpy array with dtype == Tide.dtype") - else: - raise ValueError("Must be initialised with constituents, amplitudes and phases; or a model.") - if radians: - model['phase'] = r2d*model['phase'] - self.model = model[:] - self.normalize() - - def prepare(self, *args, **kwargs): - return Tide._prepare(self.model['constituent'], *args, **kwargs) - - @staticmethod - def _prepare(constituents, t0, t = None, radians = True): - """ - Return constituent speed and equilibrium argument at a given time, and constituent node factors at given times. - Arguments: - constituents -- list of constituents to prepare - t0 -- time at which to evaluate speed and equilibrium argument for each constituent - t -- list of times at which to evaluate node factors for each constituent (default: t0) - radians -- whether to return the angular arguments in radians or degrees (default: True) - """ - #The equilibrium argument is constant and taken at the beginning of the - #time series (t0). The speed of the equilibrium argument changes very - #slowly, so again we take it to be constant over any length of data. The - #node factors change more rapidly. - if isinstance(t0, Iterable): - t0 = t0[0] - if t is None: - t = [t0] - if not isinstance(t, Iterable): - t = [t] - a0 = astro(t0) - a = [astro(t_i) for t_i in t] - - #For convenience give u, V0 (but not speed!) in [0, 360) - V0 = np.array([c.V(a0) for c in constituents])[:, np.newaxis] - speed = np.array([c.speed(a0) for c in constituents])[:, np.newaxis] - u = [np.mod(np.array([c.u(a_i) for c in constituents])[:, np.newaxis], 360.0) - for a_i in a] - f = [np.mod(np.array([c.f(a_i) for c in constituents])[:, np.newaxis], 360.0) - for a_i in a] - - if radians: - speed = d2r*speed - V0 = d2r*V0 - u = [d2r*each for each in u] - return speed, u, f, V0 - - def at(self, t): - """ - Return the modelled tidal height at given times. - Arguments: - t -- array of times at which to evaluate the tidal height - """ - t0 = t[0] - hours = self._hours(t0, t) - partition = 240.0 - t = self._partition(hours, partition) - times = self._times(t0, [(i + 0.5)*partition for i in range(len(t))]) - speed, u, f, V0 = self.prepare(t0, times, radians = True) - H = self.model['amplitude'][:, np.newaxis] - p = d2r*self.model['phase'][:, np.newaxis] - - return np.concatenate([ - Tide._tidal_series(t_i, H, p, speed, u_i, f_i, V0) - for t_i, u_i, f_i in izip(t, u, f) - ]) - - def highs(self, *args): - """ - Generator yielding only the high tides. - Arguments: - see Tide.extrema() - """ - for t in ifilter(lambda e: e[2] == 'H', self.extrema(*args)): - yield t - - def lows(self, *args): - """ - Generator yielding only the low tides. - Arguments: - see Tide.extrema() - """ - for t in ifilter(lambda e: e[2] == 'L', self.extrema(*args)): - yield t - - def form_number(self): - """ - Returns the model's form number, a helpful heuristic for classifying tides. - """ - k1, o1, m2, s2 = ( - np.extract(self.model['constituent'] == c, self.model['amplitude']) - for c in [constituent._K1, constituent._O1, constituent._M2, constituent._S2] - ) - return (k1+o1)/(m2+s2) - - def classify(self): - """ - Classify the tide according to its form number - """ - form = self.form_number() - if 0 <= form <= 0.25: - return 'semidiurnal' - elif 0.25 < form <= 1.5: - return 'mixed (semidiurnal)' - elif 1.5 < form <= 3.0: - return 'mixed (diurnal)' - else: - return 'diurnal' - - def extrema(self, t0, t1 = None, partition = 2400.0): - """ - A generator for high and low tides. - Arguments: - t0 -- time after which extrema are sought - t1 -- optional time before which extrema are sought (if not given, the generator is infinite) - partition -- number of hours for which we consider the node factors to be constant (default: 2400.0) - """ - if t1: - #yield from in python 3.4 - for e in takewhile(lambda t: t[0] < t1, self.extrema(t0)): - yield e - else: - #We assume that extrema are separated by at least delta hours - delta = np.amin([ - 90.0 / c.speed(astro(t0)) for c in self.model['constituent'] - if not c.speed(astro(t0)) == 0 - ]) - #We search for stationary points from offset hours before t0 to - #ensure we find any which might occur very soon after t0. - offset = 24.0 - partitions = ( - Tide._times(t0, i*partition) for i in count()), (Tide._times(t0, i*partition) for i in count(1) - ) - - #We'll overestimate to be on the safe side; - #values outside (start,end) won't get yielded. - interval_count = int(np.ceil((partition + offset) / delta)) + 1 - amplitude = self.model['amplitude'][:, np.newaxis] - phase = d2r*self.model['phase'][:, np.newaxis] - - for start, end in izip(*partitions): - speed, [u], [f], V0 = self.prepare(start, Tide._times(start, 0.5*partition)) - #These derivatives don't include the time dependence of u or f, - #but these change slowly. - def d(t): - return np.sum(-speed*amplitude*f*np.sin(speed*t + (V0 + u) - phase), axis=0) - def d2(t): - return np.sum(-speed**2.0 * amplitude*f*np.cos(speed*t + (V0 + u) - phase), axis=0) - #We'll overestimate to be on the safe side; - #values outside (start,end) won't get yielded. - intervals = ( - delta * i -offset for i in range(interval_count)), (delta*(i+1) - offset for i in range(interval_count) - ) - for a, b in izip(*intervals): - if d(a)*d(b) < 0: - extrema = fsolve(d, (a + b) / 2.0, fprime = d2)[0] - time = Tide._times(start, extrema) - [height] = self.at([time]) - hilo = 'H' if d2(extrema) < 0 else 'L' - if start < time < end: - yield (time, height, hilo) - - @staticmethod - def _hours(t0, t): - """ - Return the hourly offset(s) of a (list of) time from a given time. - Arguments: - t0 -- time from which offsets are sought - t -- times to find hourly offsets from t0. - """ - if not isinstance(t, Iterable): - return Tide._hours(t0, [t])[0] - elif isinstance(t[0], datetime): - return np.array([(ti-t0).total_seconds() / 3600.0 for ti in t]) - else: - return t - - @staticmethod - def _partition(hours, partition = 3600.0): - """ - Partition a sorted list of numbers (or in this case hours). - Arguments: - hours -- sorted ndarray of hours. - partition -- maximum partition length (default: 3600.0) - """ - partition = float(partition) - relative = hours - hours[0] - total_partitions = np.ceil(relative[-1] / partition + 10*np.finfo(np.float).eps).astype('int') - return [hours[np.floor(np.divide(relative, partition)) == i] for i in range(total_partitions)] - - @staticmethod - def _times(t0, hours): - """ - Return a (list of) datetime(s) given an initial time and an (list of) hourly offset(s). - Arguments: - t0 -- initial time - hours -- hourly offsets from t0 - """ - if not isinstance(hours, Iterable): - return Tide._times(t0, [hours])[0] - elif not isinstance(hours[0], datetime): - return np.array([t0 + timedelta(hours=h) for h in hours]) - else: - return np.array(hours) - - @staticmethod - def _tidal_series(t, amplitude, phase, speed, u, f, V0): - return np.sum(amplitude*f*np.cos(speed*t + (V0 + u) - phase), axis=0) - - def normalize(self): - """ - Adapt self.model so that amplitudes are positive and phases are in [0,360) as per convention - """ - for i, (_, amplitude, phase) in enumerate(self.model): - if amplitude < 0: - self.model['amplitude'][i] = -amplitude - self.model['phase'][i] = phase + 180.0 - self.model['phase'][i] = np.mod(self.model['phase'][i], 360.0) - - @classmethod - def decompose( - cls, - heights, - t = None, - t0 = None, - interval = None, - constituents = constituent.noaa, - initial = None, - n_period = 2, - callback = None, - full_output = False - ): - """ - Return an instance of Tide which has been fitted to a series of tidal observations. - Arguments: - It is not necessary to provide t0 or interval if t is provided. - heights -- ndarray of tidal observation heights - t -- ndarray of tidal observation times - t0 -- datetime representing the time at which heights[0] was recorded - interval -- hourly interval between readings - constituents -- list of constituents to use in the fit (default: constituent.noaa) - initial -- optional Tide instance to use as first guess for least squares solver - n_period -- only include constituents which complete at least this many periods (default: 2) - callback -- optional function to be called at each iteration of the solver - full_output -- whether to return the output of scipy's leastsq solver (default: False) - """ - if t is not None: - if isinstance(t[0], datetime): - hours = Tide._hours(t[0], t) - t0 = t[0] - elif t0 is not None: - hours = t - else: - raise ValueError("t can be an array of datetimes, or an array " - "of hours since t0 in which case t0 must be " - "specified.") - elif None not in [t0, interval]: - hours = np.arange(len(heights)) * interval - else: - raise ValueError("Must provide t(datetimes), or t(hours) and " - "t0(datetime), or interval(hours) and t0(datetime) " - "so that each height can be identified with an " - "instant in time.") - - #Remove duplicate constituents (those which travel at exactly the same - #speed, irrespective of phase) - constituents = list(OrderedDict.fromkeys(constituents)) - - #No need for least squares to find the mean water level constituent z0, - #work relative to mean - constituents = [c for c in constituents if not c == constituent._Z0] - z0 = np.mean(heights) - heights = heights - z0 - - #Only analyse frequencies which complete at least n_period cycles over - #the data period. - constituents = [ - c for c in constituents - if 360.0 * n_period < hours[-1] * c.speed(astro(t0)) - ] - n = len(constituents) - - sort = np.argsort(hours) - hours = hours[sort] - heights = heights[sort] - - #We partition our time/height data into intervals over which we consider - #the values of u and f to assume a constant value (that is, their true - #value at the midpoint of the interval). Constituent - #speeds change much more slowly than the node factors, so we will - #consider these constant and equal to their speed at t0, regardless of - #the length of the time series. - - partition = 240.0 - - t = Tide._partition(hours, partition) - times = Tide._times(t0, [(i + 0.5)*partition for i in range(len(t))]) - - speed, u, f, V0 = Tide._prepare(constituents, t0, times, radians = True) - - #Residual to be minimised by variation of parameters (amplitudes, phases) - def residual(hp): - H, p = hp[:n, np.newaxis], hp[n:, np.newaxis] - s = np.concatenate([ - Tide._tidal_series(t_i, H, p, speed, u_i, f_i, V0) - for t_i, u_i, f_i in izip(t, u, f) - ]) - res = heights - s - if callback: - callback(res) - return res - - #Analytic Jacobian of the residual - this makes solving significantly - #faster than just using gradient approximation, especially with many - #measurements / constituents. - def D_residual(hp): - H, p = hp[:n, np.newaxis], hp[n:, np.newaxis] - ds_dH = np.concatenate([ - f_i*np.cos(speed*t_i+u_i+V0-p) - for t_i, u_i, f_i in izip(t, u, f)], - axis = 1) - - ds_dp = np.concatenate([ - H*f_i*np.sin(speed*t_i+u_i+V0-p) - for t_i, u_i, f_i in izip(t, u, f)], - axis = 1) - - return np.append(-ds_dH, -ds_dp, axis=0) - - #Initial guess for solver, haven't done any analysis on this since the - #solver seems to converge well regardless of the initial guess We do - #however scale the initial amplitude guess with some measure of the - #variation - amplitudes = np.ones(n) * (np.sqrt(np.dot(heights, heights)) / len(heights)) - phases = np.ones(n) - - if initial: - for (c0, amplitude, phase) in initial.model: - for i, c in enumerate(constituents): - if c0 == c: - amplitudes[i] = amplitude - phases[i] = d2r*phase - - initial = np.append(amplitudes, phases) - - lsq = leastsq(residual, initial, Dfun=D_residual, col_deriv=True, ftol=1e-7) - - model = np.zeros(1+n, dtype=cls.dtype) - model[0] = (constituent._Z0, z0, 0) - model[1:]['constituent'] = constituents[:] - model[1:]['amplitude'] = lsq[0][:n] - model[1:]['phase'] = lsq[0][n:] - - if full_output: - return cls(model = model, radians = True), lsq - return cls(model = model, radians = True) + dtype = np.dtype([ + ('constituent', object), + ('amplitude', float), + ('phase', float)]) + + def __init__( + self, + constituents=None, + amplitudes=None, + phases=None, + model=None, + radians=False): + """ + Initialise a tidal model. Provide constituents, amplitudes and phases + OR a model. + Arguments: + constituents -- list of constituents used in the model. + amplitudes -- list of amplitudes corresponding to constituents + phases -- list of phases corresponding to constituents + model -- an ndarray of type Tide.dtype representing the constituents, + amplitudes and phases. + radians -- boolean representing whether phases are in radians + (default False) + """ + if None not in [constituents, amplitudes, phases]: + if len(constituents) == len(amplitudes) == len(phases): + model = np.zeros(len(phases), dtype=Tide.dtype) + model['constituent'] = np.array(constituents) + model['amplitude'] = np.array(amplitudes) + model['phase'] = np.array(phases) + else: + raise ValueError("Constituents, amplitudes and phases should " + "all be arrays of equal length.") + elif model is not None: + if not model.dtype == Tide.dtype: + raise ValueError("Model must be a numpy array with " + "dtype == Tide.dtype") + else: + raise ValueError("Must be initialised with constituents, " + "amplitudes and phases; or a model.") + if radians: + model['phase'] = r2d * model['phase'] + self.model = model[:] + self.normalize() + + def prepare(self, *args, **kwargs): + return Tide._prepare(self.model['constituent'], *args, **kwargs) + + @staticmethod + def _prepare(constituents, t0, t=None, radians=True): + """ + Return constituent speed and equilibrium argument at a given time, and + constituent node factors at given times. + Arguments: + constituents -- list of constituents to prepare + t0 -- time at which to evaluate speed and equilibrium argument for each + constituent + t -- list of times at which to evaluate node factors for each + constituent (default: t0) + radians -- whether to return the angular arguments in radians or + degrees (default: True) + """ + # The equilibrium argument is constant and taken at the beginning of + # the time series (t0). The speed of the equilibrium argument changes + # very slowly, so again we take it to be constant over any length of + # data. The node factors change more rapidly. + if isinstance(t0, Iterable): + t0 = t0[0] + if t is None: + t = [t0] + if not isinstance(t, Iterable): + t = [t] + a0 = astro(t0) + a = [astro(t_i) for t_i in t] + + # For convenience give u, V0 (but not speed!) in [0, 360) + V0 = np.array([c.V(a0) for c in constituents])[:, np.newaxis] + speed = np.array([c.speed(a0) for c in constituents])[:, np.newaxis] + u = [np.mod(np.array([c.u(a_i) for c in constituents])[:, np.newaxis], + 360.0) for a_i in a] + f = [np.mod(np.array([c.f(a_i) for c in constituents])[:, np.newaxis], + 360.0) for a_i in a] + + if radians: + speed = d2r * speed + V0 = d2r * V0 + u = [d2r * each for each in u] + return speed, u, f, V0 + + def at(self, t): + """ + Return the modelled tidal height at given times. + Arguments: + t -- array of times at which to evaluate the tidal height + """ + t0 = t[0] + hours = self._hours(t0, t) + partition = 240.0 + t = self._partition(hours, partition) + times = self._times(t0, [(i + 0.5) * partition for i in range(len(t))]) + speed, u, f, V0 = self.prepare(t0, times, radians=True) + H = self.model['amplitude'][:, np.newaxis] + p = d2r * self.model['phase'][:, np.newaxis] + + return np.concatenate([ + Tide._tidal_series(t_i, H, p, speed, u_i, f_i, V0) + for t_i, u_i, f_i in izip(t, u, f) + ]) + + def highs(self, *args): + """ + Generator yielding only the high tides. + Arguments: + see Tide.extrema() + """ + for t in ifilter(lambda e: e[2] == 'H', self.extrema(*args)): + yield t + + def lows(self, *args): + """ + Generator yielding only the low tides. + Arguments: + see Tide.extrema() + """ + for t in ifilter(lambda e: e[2] == 'L', self.extrema(*args)): + yield t + + def form_number(self): + """ + Returns the model's form number, a helpful heuristic for classifying + tides. + """ + k1, o1, m2, s2 = ( + np.extract(self.model['constituent'] == c, self.model['amplitude']) + for c in [constituent._K1, constituent._O1, constituent._M2, + constituent._S2] + ) + return (k1 + o1) / (m2 + s2) + + def classify(self): + """ + Classify the tide according to its form number + """ + form = self.form_number() + if 0 <= form <= 0.25: + return 'semidiurnal' + elif 0.25 < form <= 1.5: + return 'mixed (semidiurnal)' + elif 1.5 < form <= 3.0: + return 'mixed (diurnal)' + else: + return 'diurnal' + + def extrema(self, t0, t1=None, partition=2400.0): + """ + A generator for high and low tides. + Arguments: + t0 -- time after which extrema are sought + t1 -- optional time before which extrema are sought (if not given, the + generator is infinite) + partition -- number of hours for which we consider the node factors to + be constant (default: 2400.0) + """ + if t1: + # yield from in python 3.4 + for e in takewhile(lambda t: t[0] < t1, self.extrema(t0)): + yield e + else: + # We assume that extrema are separated by at least delta hours + delta = np.amin([ + 90.0 / c.speed(astro(t0)) for c in self.model['constituent'] + if not c.speed(astro(t0)) == 0 + ]) + # We search for stationary points from offset hours before t0 to + # ensure we find any which might occur very soon after t0. + offset = 24.0 + partitions = ( + (Tide._times(t0, i * partition) for i in count()), + (Tide._times(t0, i * partition) for i in count(1)) + ) + + # We'll overestimate to be on the safe side; + # values outside (start,end) won't get yielded. + interval_count = int(np.ceil((partition + offset) / delta)) + 1 + amplitude = self.model['amplitude'][:, np.newaxis] + phase = d2r * self.model['phase'][:, np.newaxis] + + for start, end in izip(*partitions): + speed, [u], [f], V0 = self.prepare( + start, Tide._times(start, 0.5 * partition)) + # These derivatives don't include the time dependence of u or + # f, but these change slowly. + + def d(t): + return np.sum( + -speed * amplitude * f + * np.sin( + speed * t + (V0 + u) - phase), + axis=0) + + def d2(t): + return np.sum( + -speed ** 2.0 * amplitude * f + * np.cos( + speed * t + (V0 + u) - phase), + axis=0) + # We'll overestimate to be on the safe side; + # values outside (start,end) won't get yielded. + intervals = ( + (delta * i - offset for i in range(interval_count)), + (delta * (i + 1) - offset for i in range(interval_count)) + ) + for a, b in izip(*intervals): + if d(a) * d(b) < 0: + extrema = fsolve(d, (a + b) / 2.0, fprime=d2)[0] + time = Tide._times(start, extrema) + [height] = self.at([time]) + hilo = 'H' if d2(extrema) < 0 else 'L' + if start < time < end: + yield (time, height, hilo) + + @staticmethod + def _hours(t0, t): + """ + Return the hourly offset(s) of a (list of) time from a given time. + Arguments: + t0 -- time from which offsets are sought + t -- times to find hourly offsets from t0. + """ + if not isinstance(t, Iterable): + return Tide._hours(t0, [t])[0] + elif isinstance(t[0], datetime): + return np.array([(ti-t0).total_seconds() / 3600.0 for ti in t]) + else: + return t + + @staticmethod + def _partition(hours, partition=3600.0): + """ + Partition a sorted list of numbers (or in this case hours). + Arguments: + hours -- sorted ndarray of hours. + partition -- maximum partition length (default: 3600.0) + """ + partition = float(partition) + relative = hours - hours[0] + total_partitions = np.ceil( + relative[-1] / partition + 10 * np.finfo(np.float).eps + ).astype('int') + return [hours[np.floor(np.divide(relative, partition)) == i] + for i in range(total_partitions)] + + @staticmethod + def _times(t0, hours): + """ + Return a (list of) datetime(s) given an initial time and an (list of) + hourly offset(s). + Arguments: + t0 -- initial time + hours -- hourly offsets from t0 + """ + if not isinstance(hours, Iterable): + return Tide._times(t0, [hours])[0] + elif not isinstance(hours[0], datetime): + return np.array([t0 + timedelta(hours=h) for h in hours]) + else: + return np.array(hours) + + @staticmethod + def _tidal_series(t, amplitude, phase, speed, u, f, V0): + return np.sum(amplitude * f * np.cos(speed * t + (V0 + u) - phase), + axis=0) + + def normalize(self): + """ + Adapt self.model so that amplitudes are positive and phases are in + [0,360) as per convention + """ + for i, (_, amplitude, phase) in enumerate(self.model): + if amplitude < 0: + self.model['amplitude'][i] = -amplitude + self.model['phase'][i] = phase + 180.0 + self.model['phase'][i] = np.mod(self.model['phase'][i], 360.0) + + @classmethod + def decompose( + cls, + heights, + t=None, + t0=None, + interval=None, + constituents=constituent.noaa, + initial=None, + n_period=2, + callback=None, + full_output=False + ): + """ + Return an instance of Tide which has been fitted to a series of tidal + observations. + Arguments: + It is not necessary to provide t0 or interval if t is provided. + heights -- ndarray of tidal observation heights + t -- ndarray of tidal observation times + t0 -- datetime representing the time at which heights[0] was recorded + interval -- hourly interval between readings + constituents -- list of constituents to use in the fit + (default: constituent.noaa) + initial -- optional Tide instance to use as first guess for least + squares solver + n_period -- only include constituents which complete at least this many + periods (default: 2) + callback -- optional function to be called at each iteration of the + solver + full_output -- whether to return the output of scipy's leastsq solver + (default: False) + """ + if t is not None: + if isinstance(t[0], datetime): + hours = Tide._hours(t[0], t) + t0 = t[0] + elif t0 is not None: + hours = t + else: + raise ValueError("t can be an array of datetimes, or an array " + "of hours since t0 in which case t0 must be " + "specified.") + elif None not in [t0, interval]: + hours = np.arange(len(heights)) * interval + else: + raise ValueError("Must provide t(datetimes), or t(hours) and " + "t0(datetime), or interval(hours) and " + "t0(datetime) so that each height can be " + "identified with an instant in time.") + + # Remove duplicate constituents (those which travel at exactly the same + # speed, irrespective of phase) + constituents = list(OrderedDict.fromkeys(constituents)) + + # No need for least squares to find the mean water level constituent + # z0, work relative to mean + constituents = [c for c in constituents if not c == constituent._Z0] + z0 = np.mean(heights) + heights = heights - z0 + + # Only analyse frequencies which complete at least n_period cycles over + # the data period. + constituents = [ + c for c in constituents + if 360.0 * n_period < hours[-1] * c.speed(astro(t0)) + ] + n = len(constituents) + + sort = np.argsort(hours) + hours = hours[sort] + heights = heights[sort] + + # We partition our time/height data into intervals over which we + # consider the values of u and f to assume a constant value (that is, + # their true value at the midpoint of the interval). Constituent + # speeds change much more slowly than the node factors, so we will + # consider these constant and equal to their speed at t0, regardless of + # the length of the time series. + + partition = 240.0 + + t = Tide._partition(hours, partition) + times = Tide._times(t0, [(i + 0.5) * partition for i in range(len(t))]) + + speed, u, f, V0 = Tide._prepare(constituents, t0, times, radians=True) + + # Residual to be minimised by variation of parameters + # (amplitudes, phases) + def residual(hp): + H, p = hp[:n, np.newaxis], hp[n:, np.newaxis] + s = np.concatenate([ + Tide._tidal_series(t_i, H, p, speed, u_i, f_i, V0) + for t_i, u_i, f_i in izip(t, u, f) + ]) + res = heights - s + if callback: + callback(res) + return res + + # Analytic Jacobian of the residual - this makes solving significantly + # faster than just using gradient approximation, especially with many + # measurements / constituents. + def D_residual(hp): + H, p = hp[:n, np.newaxis], hp[n:, np.newaxis] + ds_dH = np.concatenate([ + f_i * np.cos(speed * t_i + u_i + V0 - p) + for t_i, u_i, f_i in izip(t, u, f)], + axis=1) + + ds_dp = np.concatenate([ + H * f_i * np.sin(speed * t_i + u_i + V0 - p) + for t_i, u_i, f_i in izip(t, u, f)], + axis=1) + + return np.append(-ds_dH, -ds_dp, axis=0) + + # Initial guess for solver, haven't done any analysis on this since the + # solver seems to converge well regardless of the initial guess We do + # however scale the initial amplitude guess with some measure of the + # variation + amplitudes = np.ones(n) * ( + np.sqrt(np.dot(heights, heights)) / len(heights)) + phases = np.ones(n) + + if initial: + for (c0, amplitude, phase) in initial.model: + for i, c in enumerate(constituents): + if c0 == c: + amplitudes[i] = amplitude + phases[i] = d2r * phase + + initial = np.append(amplitudes, phases) + + lsq = leastsq(residual, initial, Dfun=D_residual, col_deriv=True, + ftol=1e-7) + + model = np.zeros(1 + n, dtype=cls.dtype) + model[0] = (constituent._Z0, z0, 0) + model[1:]['constituent'] = constituents[:] + model[1:]['amplitude'] = lsq[0][:n] + model[1:]['phase'] = lsq[0][n:] + + if full_output: + return cls(model=model, radians=True), lsq + return cls(model=model, radians=True)