Skip to content

Commit

Permalink
refactor: use mean lunar time as independent variable
Browse files Browse the repository at this point in the history 
refactor: modify table to use original Doodson coefficients
fix: flake8 checks
  • Loading branch information
@tsutterley
tsutterley authored and root committed Jan 11, 2024
1 parent 9e1ad0d commit a4a6c8d
Show file tree
Hide file tree
Showing 3 changed files with 98 additions and 87 deletions.
157 changes: 83 additions & 74 deletions pyTMD/arguments.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,10 +138,8 @@ def arguments(
nt = len(np.atleast_1d(MJD))
# initial time conversions
hour = 24.0*np.mod(MJD, 1)
# convert from hours solar time into degrees
# note: Doodson uses mean lunar time as the independent variable
# this conversion occurs with the coefficients in the table (tau - s + h)
tau = 15.0*hour
# convert from hours solar time into mean lunar time in degrees
tau = 15.0*hour - s + h
# variable for multiples of 90 degrees (Ray technical note 2017)
k = 90.0 + np.zeros((nt))
# degrees to radians
Expand Down Expand Up @@ -556,9 +554,6 @@ def doodson_number(
j, = [j for j,val in enumerate(cindex) if (val == c)]
# remove phase dependence "k"
coef = coefficients[:-1,j]
# revert changes for mean lunar time (tau = t - s + h)
coef[1] += coef[0] # add "s" values removed in table
coef[2] -= coef[0] # remove "h" values added in table
# add 5 to values following Doodson convention (prevent negatives)
coef[1:] += 5
# convert to single number and round off floating point errors
Expand Down Expand Up @@ -596,70 +591,77 @@ def _arguments_table(**kwargs):

# modified Doodson coefficients for constituents
# using 7 index variables: tau, s, h, p, n, pp, k
# tau: mean lunar time
# s: mean longitude of moon
# h: mean longitude of sun
# p: mean longitude of lunar perigee
# n: mean longitude of ascending lunar node
# pp: mean longitude of solar perigee
# k: 90-degree phase
coef = np.zeros((7, 60))
coef[:,0] = [0.0, 0.0, 1.0, 0.0, 0.0, -1.0, 0.0] # Sa
coef[:,1] = [0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 0.0] # Ssa
coef[:,2] = [0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0] # Mm
coef[:,3] = [0.0, 2.0, -2.0, 0.0, 0.0, 0.0, 0.0] # MSf
coef[:,4] = [0.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0] # Mf
coef[:,5] = [0.0, 3.0, 0.0, -1.0, 0.0, 0.0, 0.0] # Mt
coef[:,6] = [1.0, -5.0, 3.0, 1.0, 0.0, 0.0, -1.0] # alpha1
coef[:,7] = [1.0, -4.0, 1.0, 2.0, 0.0, 0.0, -1.0] # 2q1
coef[:,8] = [1.0, -4.0, 3.0, 0.0, 0.0, 0.0, -1.0] # sigma1
coef[:,9] = [1.0, -3.0, 1.0, 1.0, 0.0, 0.0, -1.0]# q1
coef[:,10] = [1.0, -3.0, 3.0, -1.0, 0.0, 0.0, -1.0] # rho1
coef[:,11] = [1.0, -2.0, 1.0, 0.0, 0.0, 0.0, -1.0] # o1
coef[:,12] = [1.0, -2.0, 3.0, 0.0, 0.0, 0.0, 1.0] # tau1
coef[:,13] = [1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0] # m1
coef[:,14] = [1.0, -1.0, 3.0, -1.0, 0.0, 0.0, 1.0] # chi1
coef[:,15] = [1.0, 0.0, -2.0, 0.0, 0.0, 1.0, -1.0] # pi1
coef[:,16] = [1.0, 0.0, -1.0, 0.0, 0.0, 0.0, -1.0] # p1
coef[:,6] = [1.0, -4.0, 2.0, 1.0, 0.0, 0.0, -1.0] # alpha1
coef[:,7] = [1.0, -3.0, 0.0, 2.0, 0.0, 0.0, -1.0] # 2q1
coef[:,8] = [1.0, -3.0, 2.0, 0.0, 0.0, 0.0, -1.0] # sigma1
coef[:,9] = [1.0, -2.0, 0.0, 1.0, 0.0, 0.0, -1.0]# q1
coef[:,10] = [1.0, -2.0, 2.0, -1.0, 0.0, 0.0, -1.0] # rho1
coef[:,11] = [1.0, -1.0, 0.0, 0.0, 0.0, 0.0, -1.0] # o1
coef[:,12] = [1.0, -1.0, 2.0, 0.0, 0.0, 0.0, 1.0] # tau1
coef[:,13] = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0] # m1
coef[:,14] = [1.0, 0.0, 2.0, -1.0, 0.0, 0.0, 1.0] # chi1
coef[:,15] = [1.0, 1.0, -3.0, 0.0, 0.0, 1.0, -1.0] # pi1
coef[:,16] = [1.0, 1.0, -2.0, 0.0, 0.0, 0.0, -1.0] # p1
if kwargs['corrections'] in ('OTIS', 'ATLAS', 'TMD3', 'netcdf'):
coef[:,17] = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0] # s1
coef[:,17] = [1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 1.0] # s1
elif kwargs['corrections'] in ('GOT', 'FES'):
coef[:,17] = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0] # s1 (Doodson's phase)
coef[:,18] = [1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0] # k1
coef[:,19] = [1.0, 0.0, 2.0, 0.0, 0.0, -1.0, 1.0] # psi1
coef[:,20] = [1.0, 0.0, 3.0, 0.0, 0.0, 0.0, 1.0] # phi1
coef[:,21] = [1.0, 1.0, -1.0, 1.0, 0.0, 0.0, 1.0] # theta1
coef[:,22] = [1.0, 1.0, 1.0, -1.0, 0.0, 0.0, 1.0] # j1
coef[:,23] = [1.0, 2.0, 1.0, 0.0, 0.0, 0.0, 1.0] # oo1
coef[:,24] = [2.0, -4.0, 2.0, 2.0, 0.0, 0.0, 0.0] # 2n2
coef[:,25] = [2.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0] # mu2
coef[:,26] = [2.0, -3.0, 2.0, 1.0, 0.0, 0.0, 0.0] # n2
coef[:,27] = [2.0, -3.0, 4.0, -1.0, 0.0, 0.0, 0.0] # nu2
coef[:,28] = [2.0, -2.0, 1.0, 0.0, 0.0, 1.0, 0.0] # m2a
coef[:,29] = [2.0, -2.0, 2.0, 0.0, 0.0, 0.0, 0.0] # m2
coef[:,30] = [2.0, -2.0, 3.0, 0.0, 0.0, -1.0, 0.0] # m2b
coef[:,31] = [2.0, -1.0, 0.0, 1.0, 0.0, 0.0, 2.0] # lambda2
coef[:,32] = [2.0, -1.0, 2.0, -1.0, 0.0, 0.0, 2.0] # l2
coef[:,33] = [2.0, 0.0, -1.0, 0.0, 0.0, 1.0, 0.0] # t2
coef[:,34] = [2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # s2
coef[:,35] = [2.0, 0.0, 1.0, 0.0, 0.0, -1.0, 2.0] # r2
coef[:,36] = [2.0, 0.0, 2.0, 0.0, 0.0, 0.0, 0.0] # k2
coef[:,37] = [2.0, 1.0, 2.0, 0.0, 0.0, -1.0, 0.0] # eta2
coef[:,38] = [2.0, -5.0, 4.0, 1.0, 0.0, 0.0, 0.0] # mns2
coef[:,39] = [2.0, 2.0, -2.0, 0.0, 0.0, 0.0, 0.0] # 2sm2
coef[:,40] = [3.0, -3.0, 3.0, 0.0, 0.0, 0.0, 0.0] # m3
coef[:,17] = [1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 2.0] # s1 (Doodson's phase)
coef[:,18] = [1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0] # k1
coef[:,19] = [1.0, 1.0, 1.0, 0.0, 0.0, -1.0, 1.0] # psi1
coef[:,20] = [1.0, 1.0, 2.0, 0.0, 0.0, 0.0, 1.0] # phi1
coef[:,21] = [1.0, 2.0, -2.0, 1.0, 0.0, 0.0, 1.0] # theta1
coef[:,22] = [1.0, 2.0, 0.0, -1.0, 0.0, 0.0, 1.0] # j1
coef[:,23] = [1.0, 3.0, 0.0, 0.0, 0.0, 0.0, 1.0] # oo1
coef[:,24] = [2.0, -2.0, 0.0, 2.0, 0.0, 0.0, 0.0] # 2n2
coef[:,25] = [2.0, -2.0, 2.0, 0.0, 0.0, 0.0, 0.0] # mu2
coef[:,26] = [2.0, -1.0, 0.0, 1.0, 0.0, 0.0, 0.0] # n2
coef[:,27] = [2.0, -1.0, 2.0, -1.0, 0.0, 0.0, 0.0] # nu2
coef[:,28] = [2.0, 0.0, -1.0, 0.0, 0.0, 1.0, 0.0] # m2a
coef[:,29] = [2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # m2
coef[:,30] = [2.0, 0.0, 1.0, 0.0, 0.0, -1.0, 0.0] # m2b
coef[:,31] = [2.0, 1.0, -2.0, 1.0, 0.0, 0.0, 2.0] # lambda2
coef[:,32] = [2.0, 1.0, 0.0, -1.0, 0.0, 0.0, 2.0] # l2
coef[:,33] = [2.0, 2.0, -3.0, 0.0, 0.0, 1.0, 0.0] # t2
coef[:,34] = [2.0, 2.0, -2.0, 0.0, 0.0, 0.0, 0.0] # s2
coef[:,35] = [2.0, 2.0, -1.0, 0.0, 0.0, -1.0, 2.0] # r2
coef[:,36] = [2.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0] # k2
coef[:,37] = [2.0, 3.0, 0.0, 0.0, 0.0, -1.0, 0.0] # eta2
coef[:,38] = [2.0, -3.0, 2.0, 1.0, 0.0, 0.0, 0.0] # mns2
coef[:,39] = [2.0, 4.0, -4.0, 0.0, 0.0, 0.0, 0.0] # 2sm2
coef[:,40] = [3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # m3
coef[:,41] = coef[:,18] + coef[:,29] # mk3
coef[:,42] = [3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # s3
coef[:,42] = [3.0, 3.0, -3.0, 0.0, 0.0, 0.0, 0.0] # s3
coef[:,43] = coef[:,26] + coef[:,29] # mn4
coef[:,44] = [4.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0] # m4
coef[:,44] = [4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # m4
coef[:,45] = coef[:,29] + coef[:,34] # ms4
coef[:,46] = coef[:,29] + coef[:,36] # mk4
coef[:,47] = [4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # s4
coef[:,48] = [5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # s5
coef[:,49] = [6.0, -6.0, 6.0, 0.0, 0.0, 0.0, 0.0] # m6
coef[:,50] = [6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # s6
coef[:,51] = [7.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # s7
coef[:,52] = [8.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # s8
coef[:,47] = [4.0, 4.0, -4.0, 0.0, 0.0, 0.0, 0.0] # s4
coef[:,48] = [5.0, 5.0, -5.0, 0.0, 0.0, 0.0, 0.0] # s5
coef[:,49] = [6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # m6
coef[:,50] = [6.0, 6.0, -6.0, 0.0, 0.0, 0.0, 0.0] # s6
coef[:,51] = [7.0, 7.0, -7.0, 0.0, 0.0, 0.0, 0.0] # s7
coef[:,52] = [8.0, 8.0, -8.0, 0.0, 0.0, 0.0, 0.0] # s8
# shallow water constituents
coef[:,53] = [8.0, -8.0, 8.0, 0.0, 0.0, 0.0, 0.0] # m8
coef[:,53] = [8.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # m8
coef[:,54] = coef[:,29] + coef[:,36] - coef[:,34] # mks2
coef[:,55] = [0.0, 4.0, -2.0, 0.0, 0.0, 0.0, 0.0] # msqm
coef[:,56] = [0.0, 3.0, 0.0, -1.0, 0.0, 0.0, 0.0] # mtm
coef[:,57] = [4.0, -6.0, 4.0, 2.0, 0.0, 0.0, 0.0] # n4
coef[:,58] = [2.0, -5.0, 4.0, 1.0, 0.0, 0.0, 0.0] # eps2
coef[:,57] = [4.0, -2.0, 0.0, 2.0, 0.0, 0.0, 0.0] # n4
coef[:,58] = [2.0, -3.0, 2.0, 1.0, 0.0, 0.0, 0.0] # eps2
# mean sea level
coef[:,59] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # z0
# return the coefficient table
Expand All @@ -686,26 +688,33 @@ def _minor_table(**kwargs):
"""
# modified Doodson coefficients for constituents
# using 7 index variables: tau, s, h, p, n, pp, k
# tau: mean lunar time
# s: mean longitude of moon
# h: mean longitude of sun
# p: mean longitude of lunar perigee
# n: mean longitude of ascending lunar node
# pp: mean longitude of solar perigee
# k: 90-degree phase
coef = np.zeros((7, 20))
coef[:,0] = [1.0, -4.0, 1.0, 2.0, 0.0, 0.0, -1.0] # 2q1
coef[:,1] = [1.0, -4.0, 3.0, 0.0, 0.0, 0.0, -1.0] # sigma1
coef[:,2] = [1.0, -3.0, 3.0, -1.0, 0.0, 0.0, -1.0] # rho1
coef[:,3] = [1.0, -1.0, 1.0, -1.0, 0.0, 0.0, 1.0] # m1
coef[:,4] = [1.0, -1.0, 1.0, 1.0, 0.0, 0.0, 1.0] # m1
coef[:,5] = [1.0, -1.0, 3.0, -1.0, 0.0, 0.0, 1.0] # chi1
coef[:,6] = [1.0, 0.0, -2.0, 0.0, 0.0, 1.0, -1.0] # pi1
coef[:,7] = [1.0, 0.0, 3.0, 0.0, 0.0, 0.0, 1.0] # phi1
coef[:,8] = [1.0, 1.0, -1.0, 1.0, 0.0, 0.0, 1.0] # theta1
coef[:,9] = [1.0, 1.0, 1.0, -1.0, 0.0, 0.0, 1.0] # j1
coef[:,10] = [1.0, 2.0, 1.0, 0.0, 0.0, 0.0, 1.0] # oo1
coef[:,11] = [2.0, -4.0, 2.0, 2.0, 0.0, 0.0, 0.0] # 2n2
coef[:,12] = [2.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0] # mu2
coef[:,13] = [2.0, -3.0, 4.0, -1.0, 0.0, 0.0, 0.0] # nu2
coef[:,14] = [2.0, -1.0, 0.0, 1.0, 0.0, 0.0, 2.0]# lambda2
coef[:,15] = [2.0, -1.0, 2.0, -1.0, 0.0, 0.0, 2.0] # l2
coef[:,16] = [2.0, -1.0, 2.0, 1.0, 0.0, 0.0, 0.0] # l2
coef[:,17] = [2.0, 0.0, -1.0, 0.0, 0.0, 1.0, 0.0] # t2
coef[:,18] = [2.0, -5.0, 4.0, 1.0, 0.0, 0.0, 0.0] # eps2
coef[:,19] = [2.0, 1.0, 2.0, 0.0, 0.0, -1.0, 0.0] # eta2
coef[:,0] = [1.0, -3.0, 0.0, 2.0, 0.0, 0.0, -1.0] # 2q1
coef[:,1] = [1.0, -3.0, 2.0, 0.0, 0.0, 0.0, -1.0] # sigma1
coef[:,2] = [1.0, -2.0, 2.0, -1.0, 0.0, 0.0, -1.0] # rho1
coef[:,3] = [1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 1.0] # m1
coef[:,4] = [1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0] # m1
coef[:,5] = [1.0, 0.0, 2.0, -1.0, 0.0, 0.0, 1.0] # chi1
coef[:,6] = [1.0, 1.0, -3.0, 0.0, 0.0, 1.0, -1.0] # pi1
coef[:,7] = [1.0, 1.0, 2.0, 0.0, 0.0, 0.0, 1.0] # phi1
coef[:,8] = [1.0, 2.0, -2.0, 1.0, 0.0, 0.0, 1.0] # theta1
coef[:,9] = [1.0, 2.0, 0.0, -1.0, 0.0, 0.0, 1.0] # j1
coef[:,10] = [1.0, 3.0, 0.0, 0.0, 0.0, 0.0, 1.0] # oo1
coef[:,11] = [2.0, -2.0, 0.0, 2.0, 0.0, 0.0, 0.0] # 2n2
coef[:,12] = [2.0, -2.0, 2.0, 0.0, 0.0, 0.0, 0.0] # mu2
coef[:,13] = [2.0, -1.0, 2.0, -1.0, 0.0, 0.0, 0.0] # nu2
coef[:,14] = [2.0, 1.0, -2.0, 1.0, 0.0, 0.0, 2.0]# lambda2
coef[:,15] = [2.0, 1.0, 0.0, -1.0, 0.0, 0.0, 2.0] # l2
coef[:,16] = [2.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0] # l2
coef[:,17] = [2.0, 2.0, -3.0, 0.0, 0.0, 1.0, 0.0] # t2
coef[:,18] = [2.0, -3.0, 2.0, 1.0, 0.0, 0.0, 0.0] # eps2
coef[:,19] = [2.0, 3.0, 0.0, 0.0, 0.0, -1.0, 0.0] # eta2
# return the coefficient table
return coef
20 changes: 9 additions & 11 deletions pyTMD/predict.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -355,19 +355,17 @@ def infer_minor(
zmin[:,18] = 0.53285*z[:,8] - 0.03304*z[:,4]# eps2
zmin[:,19] = -0.0034925*z[:,5] + 0.0831707*z[:,7]# eta2

# initial time conversions
hour = 24.0*np.mod(MJD, 1)
# convert from hours solar time into degrees
# note: Doodson uses mean lunar time as the independent variable
# this conversion occurs with the coefficients in the table (tau - s + h)
tau = 15.0*hour
# variable for multiples of 90 degrees (Ray technical note 2017)
k = 90.0 + np.zeros((n))
# set function for astronomical longitudes
ASTRO5 = True if kwargs['corrections'] in ('GOT', 'FES') else False
# convert from Modified Julian Dates into Ephemeris Time
s, h, p, n, pp = pyTMD.astro.mean_longitudes(MJD + kwargs['deltat'],
ASTRO5=ASTRO5)
# initial time conversions
hour = 24.0*np.mod(MJD, 1)
# convert from hours solar time into mean lunar time in degrees
tau = 15.0*hour - s + h
# variable for multiples of 90 degrees (Ray technical note 2017)
k = 90.0 + np.zeros((n))

# determine equilibrium arguments
fargs = np.c_[tau, s, h, p, n, pp, k]
Expand Down Expand Up @@ -424,9 +422,9 @@ def infer_minor(
xi[xi > np.pi] -= 2.0*np.pi
nu = at1 - at2
I2 = np.tan(II/2.0)
Ra1 = np.sqrt(1.0 - 12.0*(I2**2)*np.cos(2.0*(P - xi)) + 36.0*(I2**4))
P2 = np.sin(2.0*(P - xi))
Q2 = 1.0/(6.0*(I2**2)) - np.cos(2.0*(P - xi))
Ra1 = np.sqrt(1.0 - 12.0*(I2**2)*np.cos(2.0*(p - xi)) + 36.0*(I2**4))
P2 = np.sin(2.0*(p - xi))
Q2 = 1.0/(6.0*(I2**2)) - np.cos(2.0*(p - xi))
R = np.arctan(P2/Q2)

f[:,0] = np.sin(II)*(np.cos(II/2.0)**2)/0.38 # 2Q1
Expand Down
8 changes: 6 additions & 2 deletions test/test_arguments.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,8 @@ def test_arguments(corrections):
# convert from hours into degrees
t1 = 15.0*hour
t2 = 30.0*hour
# convert from hours solar time into mean lunar time in degrees
tau = 15.0*hour - s + h
# variable for multiples of 90 degrees (Ray technical note 2017)
k = 90.0 + np.zeros((nt))

Expand Down Expand Up @@ -104,7 +106,7 @@ def test_arguments(corrections):
arg[:,59] = 0.0 # Z0

# determine equilibrium arguments using table
fargs = np.c_[t1, s, h, p, omega, pp, k]
fargs = np.c_[tau, s, h, p, omega, pp, k]
test = np.dot(fargs, _arguments_table(corrections=corrections))

# validate arguments between methods
Expand All @@ -125,6 +127,8 @@ def test_minor():
# convert from hours into degrees
t1 = 15.0*hour
t2 = 30.0*hour
# convert from hours solar time into mean lunar time in degrees
tau = 15.0*hour - s + h
# variable for multiples of 90 degrees (Ray technical note 2017)
k = 90.0 + np.zeros((nt))

Expand Down Expand Up @@ -152,7 +156,7 @@ def test_minor():
arg[:,19] = t2 + s + 2.0*h - pp # eta2

# determine equilibrium arguments
fargs = np.c_[t1, s, h, p, n, pp, k]
fargs = np.c_[tau, s, h, p, n, pp, k]
test = np.dot(fargs, _minor_table())

# validate arguments between methods
Expand Down

0 comments on commit a4a6c8d

Please sign in to comment.