高調波定数を使用して潮汐を予測する方法

Question

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

私はブリッジポートから（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 

Answer 1

あなたがなりC 1のアプローチを使用している場合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; 
} 

Answer 2

National Tidal Datum EpochがUnix epochには関係ありません。その期間の平均潮位の大きさ参照です。 Bridgeportデータには、各コンポーネントの大きさと位相が含まれており、最後の列のコンポーネントの頻度が示されます。フェーズはGMT（GMTは0度）を基準にしているため、時間をGMTで指定するか、タイムゾーンに従ってオフセットします。最後の列は、場所ではなくコンポーネントの関数なので、列は任意の場所で同じです。だから問題は、フーリエ解析のような、異なる周波数の正弦波を合計することです。結果の関数は、コンポーネントがすべて（24時間以外の）異なる期間を持つため、ではなくに24時間の期間があります。私は参考としてthisを使用しました。

以下のコードでは、いくつかのオフセット引数を追加しました（あなたの3.16999のように、どこから来たかわかりません）。私は読解と解析のコードを計算コードから分離しました。 get_tidesから返される関数には、測定データがclosureに含まれています。あなたはそうする必要はなく、代わりにクラスを使うことができます。

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 

EDIT：今、特定の日付、または時刻の場合 ：

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))) 

Answer 3

「フェーズは度であり、GMTを基準あなたのリンク上では、ブリッジポートの潮は、グリニッジ子午線の理論的な「平衡」の潮汐から同期していない「段階的な」角度であることを意味します。

また、時間の関数として平衡潮汐が必要です。 1つの方法は、http://coastalengineeringmanual.tpub.com/Part-II-Chap5/Part-II-Chap50024.htmの2010列の定数をコピーし、2010-01-01T00：00：00を時間データ/エポック/ゼロとして使用することです。

さらに、月の節の歳差運動に応じて多数の成分を調節する節点因子が存在する。

テーブルからコピーするのではなく、tide_facを編集して実行できます。fプログラムをhttp://adcirc.org/home/related-software/adcirc-utility-programs/から選択して、自分の選択した時間データ/エポック/ゼロで節点因子と平衡潮汐成分のセットを生成します。

はここ2013-01-01T00のために結節と平衡潮相の引数を生成するtide_fac.f修正プログラムからの出力です：00：00：

TIDAL FACTORS STARTING: HR- 0.00, DAY- 1, MONTH- 1 YEAR- 2013 

FOR A RUN LASTING 365.00 DAYS 


CONST NODE  EQ ARG (ref GM) 
NAME FACTOR (DEG) 


M2 1.02718  270.21 
S2 1.00000  0.00 
N2 1.02718  16.12 
K1 0.92336  17.70 
M4 1.05510  180.42 
O1 0.87509  248.94 
M6 1.08377  90.63 
MK3 0.94846  287.91 
S4 1.00000  0.00 
MN4 1.05510  286.33 
NU2 1.02718  73.02 
S6 1.00000  0.00 
MU2 1.02718  178.93 
2N2 1.02718  122.03 
OO1 0.63334  333.60 
LAMB 1.02718  287.40 
S1 1.00000  180.00 
M1 0.00000  159.37 
J1 0.89282  275.36 
MM 1.09369  254.09 
SSA 1.00000  201.62 
SA 1.00000  280.81 
MSF 1.02718  91.28 
MF 0.74476  312.33 
RHO1 0.87509  51.75 
Q1 0.87509  354.85 
T2 1.00000  2.35 
R2 1.00000  177.65 
2Q1 0.87509  100.76 
P1 1.00000  349.19 
2SM2 1.02718  89.79 
M3 1.04114  225.31 
L2 0.00000  348.15 
2MK3 0.97424  162.72 
K2 0.81662  214.58 
M8 1.11323  0.84 
MS4 1.02718  270.21 

Answer 4

私はちょうど潮の分析と予測のためPytidesと呼ばれる小さなPythonライブラリを書いています。それはまだ初期段階ですが、あなたはそれを使いやすくすることができます。私はこれを行う方法のan exampleを書きました。

Daveが述べたように、あなたの方法がうまくいかない理由は、NOAA段階が構成要素の平衡議論に関して公開されているからです。したがって、これらの平衡議論を行うための何らかの方法が必要です（私はPytides君は！）。