15

調和定数を使用して潮汐予測表を作成する方法を理解しようとしています。

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

これは、「今日」の 1 時間ごとに 1 つの高さを出力します。結果の高さは、-0.5 から 7.5 フィートの範囲であると予想されますが、日付には正しくありません。

私は正しい軌道に乗っていますか？潮汐の時代をどのように決定するのですか？ウィキペディア ( http://en.wikipedia.org/wiki/Arthur_Thomas_Doodson )の Doodsen の例では、1991 年 9 月 1 日の結果を得るために 0 を使用しました。私にはうまくいかないようです。

これが私の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
4

4 に答える 4

4

National Tidal Datum EpochはUnixのエポックとは関係ありません。これは、その期間の平均潮位のマグニチュード リファレンスです。Bridgeport データには、各成分の振幅と位相が含まれており、最後の列に成分の周波数が示されています。位相は GMT (GMT は 0 度) に対して相対的であるため、時間を GMT で指定するか、タイムゾーンに従ってオフセットします。最後の列は場所ではなくコンポーネントの関数であるため、その列はどの場所でも同じです。問題は、フーリエ解析のように、異なる周波数の正弦波を合計することです。コンポーネントはすべて異なる (24 時間ではない) 期間を持つため、結果として得られる関数には24 時間の期間が含まれない可能性が高くなります。私はこれを使いました参考までに。

以下のコードでは、いくつかのオフセット修飾子を追加するためにいくつかのオプションの引数を含めました (あなたの 3.16999 のように、それがどこからのものかはわかりません)。読み取りと解析のコードを計算コードから分離しました。から返される関数にget_tidesは、測定データが含まれています。そのようにする必要はなく、代わりにクラスを使用できます。

from __future__ import division
import math
import time

def get_tides(fn):
    """Read tide data from file, return a function that
    gives tide level for specified hour."""
    with open(fn,'r') as fh:
        lines = fh.readlines()

    measured = []
    for line in lines:
        index,cons,ampl,phase,speed = line.split()
        measured.append(tuple(float(x) for x in (ampl,phase,speed)))

    def find_level(hour,total_offset=0,time_offset=0):
        total = total_offset
        for ampl,phase,speed in measured:
            # you could break out the constituents here
            # to see the individual contributions
            cosarg = speed * (hour + time_offset) + phase
            total += ampl * math.cos(cosarg * math.pi/180) # do rad conversion here
        return total

    return find_level


bridgeport = get_tides('bridgeport.txt')

for hour in range(-24,25,1):
    print("%3sh %7.4f" % (hour,bridgeport(hour,total_offset=3.61999)))

そして出力：

-24h -0.5488
-23h -1.8043
-22h -1.7085
-21h -0.3378
-20h  1.8647
-19h  4.4101
-18h  6.8374
-17h  8.5997
-16h  9.1818
-15h  8.4168
-14h  6.5658
-13h  4.1003
-12h  1.5669
-11h -0.3936
-10h -1.2009
 -9h -0.6705
 -8h  0.9032
 -7h  3.0316
 -6h  5.2485
 -5h  7.0432
 -4h  7.8633
 -3h  7.4000
 -2h  5.8028
 -1h  3.5317
  0h  1.1223
  1h -0.8425
  2h -1.7748
  3h -1.3620
  4h  0.2196
  5h  2.4871
  6h  4.9635
  7h  7.2015
  8h  8.6781
  9h  8.9524
 10h  7.9593
 11h  6.0188
 12h  3.6032
 13h  1.2440
 14h -0.4444
 15h -0.9554
 16h -0.2106
 17h  1.4306
 18h  3.4834
 19h  5.5091
 20h  7.0246
 21h  7.5394
 22h  6.8482
 23h  5.1795
 24h  2.9995

編集: 特定の日付、または現在の時刻:

import datetime

epoch = datetime.datetime(1991,9,1,0,0,0)
now = datetime.datetime.utcnow()
hourdiff = (now-epoch).days*24

for hour in range(0,25,1):
    print("%3sh %7.4f" % (hour,bridgeport(hour+hourdiff)))
于 2013-02-11T19:46:34.777 に答える