python - 如何使用谐波常数预测潮汐

Question

我正在尝试了解如何使用调和常数创建潮汐预测表。

我使用了 Bridgeport 的谐波常数 (http://tidesandcurrents.noaa.gov/data_menu.shtml?stn=8467150%20Bridgeport,%20CT&type=Harmonic%20Constituents)

并使用这个 python 脚本总结潮汐成分 -

import math
import time

tidalepoch = 0
epoch = time.mktime(time.gmtime()) - tidalepoch
f = open('bridgeport.txt', 'r')
M_PI = 3.14159
lines = f.readlines()
t = epoch - 24 * 3600
i = -24
while t < epoch:
  height = 0
  for line in lines:
    x = line.split()
    A = float(x[2]) # amplitude
    B = float(x[3]) # phase
    B *=  M_PI / 180.0
    C = float(x[4]) # speed
    C *= M_PI / 648000

    # h = R cost (wt - phi)
    height += A * math.cos(C * t - B)

  print str(i) + " " + str(height + 3.61999)
  i += 1
  t += 3600

“今天”每小时打印一个高度。生成的高度大约在我预期的范围内，即 -0.5 到 7.5 英尺，但与日期不符。

我走在正确的轨道上吗？如何确定潮汐纪元？在维基百科 (http://en.wikipedia.org/wiki/Arthur_Thomas_Doodson) 上的 Doodsen 示例中，他们使用 0 来获得 1991 年 9 月 1 日的结果。但是它们的谐波值与当前发布的谐波值不同，而且那个日期似乎对我不起作用。

这是我的 bridgeport.txt 文件的内容 -

1 M2                       3.251 109.6  28.9841042 
2 S2                       0.515 135.9  30.0000000 
3 N2                       0.656  87.6  28.4397295 
4 K1                       0.318 191.6  15.0410686 
5 M4                       0.039 127.4  57.9682084 
6 O1                       0.210 219.5  13.9430356 
7 M6                       0.044 353.9  86.9523127 
8 MK3                      0.023 198.8  44.0251729 
9 S4                       0.000   0.0  60.0000000 
10 MN4                      0.024  97.2  57.4238337 
11 NU2                      0.148  89.8  28.5125831 
12 S6                       0.000   0.0  90.0000000 
13 MU2                      0.000   0.0  27.9682084 
14 2N2                      0.077  65.6  27.8953548 
15 OO1                      0.017 228.7  16.1391017 
16 LAM2                     0.068 131.1  29.4556253 
17 S1                       0.031 175.5  15.0000000 
18 M1                       0.024 264.4  14.4966939 
19 J1                       0.021 237.0  15.5854433 
20 MM                       0.000   0.0   0.5443747 
21 SSA                      0.072  61.2   0.0821373 
22 SA                       0.207 132.0   0.0410686 
23 MSF                      0.000   0.0   1.0158958 
24 MF                       0.000   0.0   1.0980331 
25 RHO                      0.015 258.1  13.4715145 
26 Q1                       0.059 205.7  13.3986609 
27 T2                       0.054 106.4  29.9589333 
28 R2                       0.004 136.9  30.0410667 
29 2Q1                      0.014 238.8  12.8542862 
30 P1                       0.098 204.1  14.9589314 
31 2SM2                     0.000   0.0  31.0158958 
32 M3                       0.012 200.1  43.4761563 
33 L2                       0.162 134.1  29.5284789 
34 2MK3                     0.015 203.7  42.9271398 
35 K2                       0.150 134.7  30.0821373 
36 M8                       0.000   0.0 115.9364166 
37 MS4                      0.000   0.0  58.9841042

Answer 1

如果您使用 C，一种方法是使用 libTCD http://www.flaterco.com/xtide/libtcd.html存储组成数据，这为访问数据提供了一个很好的框架。

要进行计算，您需要将纪元调整为当年的开始。此代码使用 libTCD 中的函数，并基于 xTide 中使用的算法。

TIDE_RECORD record;
int readStatus = read_tide_record(self.stationID, &record);
int epochStartSeconds = staringSecondForYear(year);

for (unsigned constituentNumber=0; constituentNumber<constituentCount; constituentNumber++) {
    float constituentAmplitude = record.amplitude[constituentNumber];
    float nodeFactor = get_node_factor(constituentNumber, year);
    amplitudes[constituentNumber] =  constituentAmplitude * nodeFactor;

    float constituentEpoch =  deg2rad(-record.epoch[constituentNumber]);
    float equilibrium = deg2rad(get_equilibrium(constituentNumber, year));
    phases[constituentNumber] = constituentEpoch + equilibrium;
}

然后计算自年初以来偏移量的潮汐:

- (float)getHeightForTimeSince1970:(long)date {
  //calculate the time to use for this calculation
  int secondsSinceEpoch = date - epochStartSeconds;

  float height = 0;
  for(int i = 0; i < constituentCount; i++) {
    float degreesPerHour = get_speed(i);
    float radiansPerSecond = degreesPerHour * M_PI / 648000.0;
    height += amplitudes[i] * cos(radiansPerSecond * secondsSinceEpoch + phases[i]);
  }

  height += record.datum_offset;
  return height;
}