Skip to content

Instantly share code, notes, and snippets.

@ciscorn

ciscorn/tide_main.rs

Last active January 24, 2025 06:28
Show Gist options
  • Save ciscorn/5843c5d4b23432243d73efdb4492251e to your computer and use it in GitHub Desktop.
Save ciscorn/5843c5d4b23432243d73efdb4492251e to your computer and use it in GitHub Desktop.
Download ZIP
tide, tidal, 天文潮位
Raw
tide_main.rs
///
///
/// 参考: https://www1.kaiho.mlit.go.jp/kenkyu/report/rhr60/rhr60_r_01.pdf
/// 潮汐表の計算について 海洋情報部研究報告 第60号 令和4年3月18日
use chrono::{Datelike, NaiveDate};
/// 各分潮の角速度 (度/時)
const OMEGAS: [f64; 60] = [
0.0410686, 0.0821373, 0.5443747, 1.0158958, 1.0980331, 12.8542862, 12.9271398, 13.3986609,
13.4715145, 13.9430356, 14.0251729, 14.4920521, 14.5695476, 14.9178647, 14.9589314, 15.0000000,
15.0410686, 15.0821353, 15.1232059, 15.5125897, 15.5854433, 16.0569644, 16.1391017, 27.3416964,
27.4238337, 27.8953548, 27.9682084, 28.4397295, 28.5125831, 28.9019669, 28.9841042, 29.0662415,
29.4556253, 29.5284789, 29.9589333, 30.0000000, 30.0410667, 30.0821373, 30.5443747, 30.6265120,
31.0158958, 42.9271398, 43.4761563, 43.9430356, 44.0251729, 45.0410686, 57.4238337, 57.9682084,
58.4397295, 58.9841042, 59.0662415, 60.0000000, 60.0821373, 86.4079380, 86.9523127, 87.4238337,
87.9682084, 88.0503457, 88.9841042, 89.0662415,
];
// NAGOYA
const PARAMS: [[f64; 2]; 60] = [
[14.73, 145.84],
[1.47, 26.92],
[0.76, 81.97],
[0.37, 0.42],
[0.91, 115.08],
[0.48, 147.04],
[0.49, 156.30],
[3.64, 156.98],
[0.65, 161.63],
[17.78, 166.66],
[0.36, 205.04],
[0.88, 178.56],
[0.23, 174.23],
[0.45, 178.50],
[7.62, 182.74],
[0.87, 32.96],
[24.03, 186.49],
[0.07, 285.27],
[0.31, 167.51],
[0.28, 206.52],
[1.25, 208.62],
[0.35, 307.38],
[0.79, 239.49],
[0.10, 52.03],
[0.52, 192.88],
[1.50, 172.97],
[1.95, 197.10],
[11.08, 174.52],
[2.17, 173.21],
[0.16, 150.21],
[64.74, 177.70],
[0.36, 277.40],
[0.74, 166.05],
[2.24, 181.17],
[1.97, 204.37],
[30.35, 203.85],
[0.04, 69.45],
[8.51, 199.23],
[0.14, 349.01],
[0.55, 25.78],
[0.28, 10.85],
[0.67, 282.17],
[1.38, 197.55],
[0.77, 298.62],
[0.88, 283.96],
[0.76, 38.92],
[0.12, 19.60],
[0.05, 339.82],
[0.03, 345.27],
[0.10, 90.60],
[0.05, 57.23],
[0.35, 305.41],
[0.09, 252.79],
[0.08, 273.53],
[0.17, 270.95],
[0.06, 299.24],
[0.25, 300.12],
[0.06, 255.82],
[0.10, 331.95],
[0.07, 336.29],
];
/// ユリウス通日を求める
fn julian_day(date: NaiveDate) -> f32 {
let month = date.month();
let a = (14 - month) / 12;
let y = (date.year() as u32) + 4800 - a;
let m = month + 12 * a - 3;
let jd = date.day() + (153 * m + 2) / 5 + y * 365 + y / 4 - y / 100 + y / 400 - 32045;
jd as f32 - 0.5
}
pub struct TideEstimator {
v0: [f64; 60],
f: [f64; 60],
u: [f64; 60],
}
impl TideEstimator {
pub fn new(date_begin: NaiveDate, date_mid: NaiveDate) -> Self {
Self::setup_coefficients(date_begin, date_mid)
}
/// 天文潮位を推算する
pub fn estimate(&self, lng: f64, z0: f64, params: &[[f64; 2]; 60], t: f64) -> f64 {
z0 + (0..60)
.map(|i| {
let [hi, ki] = params[i];
let a0 = match i {
0..=4 => 0,
5..=22 => 1,
23..=40 => 2,
41..=45 => 3,
46..=52 => 4,
_ => 6,
} as f64;
self.f[i]
* hi
* (a0 * lng + OMEGAS[i] * (t - 9.) + self.v0[i] + self.u[i] - ki)
.to_radians()
.cos()
})
.sum::<f64>()
}
/// 天文潮位を推算するための係数を準備する
fn setup_coefficients(date_begin: NaiveDate, date_mid: NaiveDate) -> TideEstimator {
let t_md = (julian_day(date_begin) as f64 - 2378496.) / 36525.;
let t_sd = (julian_day(date_begin) as f64 - 2415021.) / 36525.;
let t_mm = (julian_day(date_mid) as f64 - 2378496.) / 36525.;
let s_d = 335.723436
+ 481267.887361 * t_md
+ 3.38888e-3 * t_md.powi(2)
+ 1.83333e-6 * t_md.powi(3);
let p_d =
225.397325 + 4069.053805 * t_md - 1.02869e-2 * t_md.powi(2) - 1.22222e-5 * t_md.powi(3);
let h_d = 280.6824 + 36000.769325 * t_sd + 7.2222e-4 * t_sd.powi(2);
let n_m =
33.272936 - 1934.144694 * t_mm + 2.08028e-3 * t_md.powi(2) + 2.08333e-6 * t_mm.powi(3);
let p_m =
225.397325 + 4069.053805 * t_mm - 1.02869e-2 * t_mm.powi(2) - 1.22222e-5 * t_mm.powi(3);
let cos_nm = n_m.to_radians().cos();
let cos_2nm = (2. * n_m).to_radians().cos();
let cos_3nm = (3. * n_m).to_radians().cos();
let sin_nm = n_m.to_radians().sin();
let sin_2nm = (2. * n_m).to_radians().sin();
let sin_3nm = (3. * n_m).to_radians().sin();
let f_mm = 1.0000 + (-0.1300 * cos_nm) + (0.0013 * cos_2nm) + (0.0000 * cos_3nm);
let u_mm = (0.00 * sin_nm) + (0.00 * sin_2nm) + (0.00 * sin_3nm);
let f_mf = 1.0429 + (0.4135 * cos_nm) + (-0.0040 * cos_2nm) + (0.0000 * cos_3nm);
let u_mf = (-23.74 * sin_nm) + (2.68 * sin_2nm) + (-0.38 * sin_3nm);
let f_o1 = 1.0089 + (0.1871 * cos_nm) + (-0.0147 * cos_2nm) + (0.0014 * cos_3nm);
let u_o1 = (10.80 * sin_nm) + (-1.34 * sin_2nm) + (0.19 * sin_3nm);
let f_k1 = 1.0060 + (0.1150 * cos_nm) + (-0.0088 * cos_2nm) + (0.0006 * cos_3nm);
let u_k1 = (-8.86 * sin_nm) + (0.68 * sin_2nm) + (-0.07 * sin_3nm);
let f_j1 = 1.0129 + (0.1676 * cos_nm) + (-0.0170 * cos_2nm) + (0.0016 * cos_3nm);
let u_j1 = (-12.94 * sin_nm) + (1.34 * sin_2nm) + (-0.19 * sin_3nm);
let f_oo1 = 1.1027 + (0.6504 * cos_nm) + (0.0317 * cos_2nm) + (-0.0014 * cos_3nm);
let u_oo1 = (-36.68 * sin_nm) + (4.02 * sin_2nm) + (-0.57 * sin_3nm);
let f_m2 = 1.0004 + (-0.0373 * cos_nm) + (0.0002 * cos_2nm) + (0.0000 * cos_3nm);
let u_m2 = (-2.14 * sin_nm) + (0.00 * sin_2nm) + (0.00 * sin_3nm);
let f_k2 = 1.0241 + (0.2863 * cos_nm) + (0.0083 * cos_2nm) + (-0.0015 * cos_3nm);
let u_k2 = (-17.74 * sin_nm) + (0.68 * sin_2nm) + (-0.04 * sin_3nm);
let (f_l2, u_l2) = {
let x = 1.
- 0.2505 * (2. * p_m).to_radians().cos()
- 0.1102 * (2. * p_m - n_m).to_radians().cos()
- 0.0156 * (2. * p_m - 2. * n_m).to_radians().cos()
- 0.0370 * (n_m).to_radians().cos();
let y = -0.2505 * (2. * p_m).to_radians().sin()
- 0.1102 * (2. * p_m - n_m).to_radians().sin()
- 0.0156 * (2. * p_m - 2. * n_m).to_radians().sin()
- 0.0370 * (n_m).to_radians().sin();
((x * x + y * y).sqrt(), f64::atan2(y, x).to_degrees())
};
let (f_m1, u_m1) = {
let x = 2. * p_m.to_radians().cos() + 0.4 * (p_m - n_m).to_radians().cos();
let y = p_m.to_radians().sin() + 0.2 * (p_m - n_m).to_radians().sin();
((x * x + y * y).sqrt(), f64::atan2(y, x).to_degrees())
};
let calc_v0 = |a1: i32, a2: i32, a3: i32, a4: i32| {
(a1 as f64 * s_d) + (a2 as f64 * h_d) + (a3 as f64 * p_d) + a4 as f64
};
let mut v0 = [0.0; 60];
let mut f = [1.0; 60];
let mut u = [0.0; 60];
// Sa 太陽年周潮
v0[0] = calc_v0(0, 1, 0, 0);
// Ssa 太陽半年周潮
v0[1] = calc_v0(0, 2, 0, 0);
// Mm 太陰月周潮
v0[2] = calc_v0(1, 0, -1, 0);
(f[2], u[2]) = (f_mm, u_mm);
// MSf 日月合成半月周潮
v0[3] = calc_v0(2, -2, 0, 0);
(f[3], u[3]) = (f_m2, -u_m2);
// Mf 太陰半月周潮
v0[4] = calc_v0(2, 0, 0, 0);
(f[4], u[4]) = (f_mf, u_mf);
// 2Q1 二次太陰楕率潮
v0[5] = calc_v0(-4, 1, 2, 270);
(f[5], u[5]) = (f_o1, u_o1);
// SIG1
v0[6] = calc_v0(-4, 3, 0, 270);
(f[6], u[6]) = (f_o1, u_o1);
// Q1 主太陰楕率潮
v0[7] = calc_v0(-3, 1, 1, 270);
(f[7], u[7]) = (f_o1, u_o1);
// RHO1 主太陰出差潮
v0[8] = calc_v0(-3, 3, -1, 270);
(f[8], u[8]) = (f_o1, u_o1);
// O1 主太陰日周潮
v0[9] = calc_v0(-2, 1, 0, 270);
(f[9], u[9]) = (f_o1, u_o1);
// MP1
v0[10] = calc_v0(-2, 3, 0, 90);
(f[10], u[10]) = (f_m2, u_m2);
// M1 副太陰楕率潮
v0[11] = calc_v0(-1, 1, 0, 90);
(f[11], u[11]) = (f_m1, u_m1);
// CHI1
v0[12] = calc_v0(-1, 3, -1, 90);
(f[12], u[12]) = (f_j1, u_j1);
// PI1 主太陽楕率潮
v0[13] = calc_v0(0, -2, 0, 193);
// P1 主太陽日周潮
v0[14] = calc_v0(0, -1, 0, 270);
// S1 気象日周潮
v0[15] = calc_v0(0, 0, 0, 180);
// K1 日月合成日周潮
v0[16] = calc_v0(0, 1, 0, 90);
(f[16], u[16]) = (f_k1, u_k1);
// PSI1 副太陽楕率潮
v0[17] = calc_v0(0, 2, 0, 167);
// PHI1 二次太陽日周潮
v0[18] = calc_v0(0, 3, 0, 90);
// THE1
v0[19] = calc_v0(1, -1, 1, 90);
(f[19], u[19]) = (f_j1, u_j1);
// J1 小太陰楕率潮
v0[20] = calc_v0(1, 1, -1, 90);
(f[20], u[20]) = (f_j1, u_j1);
// SO1
v0[21] = calc_v0(2, -1, 0, 90);
(f[21], u[21]) = (f_o1, -u_o1);
// OO1 二次太陰日周潮
v0[22] = calc_v0(2, 1, 0, 90);
(f[22], u[22]) = (f_oo1, u_oo1);
// OQ2
v0[23] = calc_v0(-5, 2, 1, 180);
(f[23], u[23]) = (f_o1 * f_o1, u_o1 * u_o1);
// MNS2
v0[24] = calc_v0(-5, 4, 1, 0);
(f[24], u[24]) = (f_m2 * f_m2, 2. * u_m2);
// 2N2 二次太陰楕率潮
v0[25] = calc_v0(-4, 2, 2, 0);
(f[25], u[25]) = (f_m2, u_m2);
// MU2 太陰二均差潮
v0[26] = calc_v0(-4, 4, 0, 0);
(f[26], u[26]) = (f_m2, u_m2);
// N2 主太陰楕率潮
v0[27] = calc_v0(-3, 2, 1, 0);
(f[27], u[27]) = (f_m2, u_m2);
// NU2 主太陰出差潮
v0[28] = calc_v0(-3, 4, -1, 0);
(f[28], u[28]) = (f_m2, u_m2);
// OP2
v0[29] = calc_v0(-2, 0, 0, 180);
(f[29], u[29]) = (f_o1, u_o1);
// M2 主太陰半日周潮
v0[30] = calc_v0(-2, 2, 0, 0);
(f[30], u[30]) = (f_m2, u_m2);
// MKS2
v0[31] = calc_v0(-2, 4, 0, 0);
(f[31], u[31]) = (f_m2 * f_k2, u_m2 + u_k2);
// LAM2 副太陰出差潮
v0[32] = calc_v0(-1, 0, 1, 180);
(f[32], u[32]) = (f_m2, u_m2);
// L2 副太陰楕率潮
v0[33] = calc_v0(-1, 2, -1, 180);
(f[33], u[33]) = (f_l2, u_l2);
// T2 主太陽楕率潮
v0[34] = calc_v0(0, -1, 0, 283);
// S2 主太陽半日周潮
v0[35] = calc_v0(0, 0, 0, 0);
// R2 副太陽楕率潮
v0[36] = calc_v0(0, 1, 0, 257);
// K2 日月合成半日周潮
v0[37] = calc_v0(0, 2, 0, 0);
(f[37], u[37]) = (f_k2, u_k2);
// MSN2
v0[38] = calc_v0(1, 0, -1, 0);
(f[38], u[38]) = (f_m2 * f_m2, 2. * u_m2);
// KJ2
v0[39] = calc_v0(1, 2, -1, 180);
(f[39], u[39]) = (f_k1 * f_j1, u_k1 + u_j1);
// 2SM2
v0[40] = calc_v0(2, -2, 0, 0);
(f[40], u[40]) = (f_m2, -u_m2);
// MO3
v0[41] = calc_v0(-4, 3, 0, 270);
(f[41], u[41]) = (f_m2 * f_o1, u_m2 + u_o1);
// M3
v0[42] = calc_v0(-3, 3, 0, 180);
(f[42], u[42]) = (f_m2.powf(1.5), 1.5 * u_m2);
// SO3
v0[43] = calc_v0(-2, 1, 0, 270);
(f[43], u[43]) = (f_o1, u_o1);
// MK3
v0[44] = calc_v0(-2, 3, 0, 90);
(f[44], u[44]) = (f_m2 * f_k1, u_m2 + u_k1);
// SK3
v0[45] = calc_v0(0, 1, 0, 90);
(f[45], u[45]) = (f_k1, u_k1);
// MN4
v0[46] = calc_v0(-5, 4, 1, 0);
(f[46], u[46]) = (f_m2 * f_m2, 2. * u_m2);
// M4 太陰1/4日周潮
v0[47] = calc_v0(-4, 4, 0, 0);
(f[47], u[47]) = (f_m2 * f_m2, 2. * u_m2);
// SN4
v0[48] = calc_v0(-3, 2, 1, 0);
(f[48], u[48]) = (f_m2, u_m2);
// MS4
v0[49] = calc_v0(-2, 2, 0, 0);
(f[49], u[49]) = (f_m2, u_m2);
// MK4
v0[50] = calc_v0(-2, 4, 0, 0);
(f[50], u[50]) = (f_m2 * f_k2, u_m2 + u_k2);
// S4 太陽1/4日周潮
v0[51] = calc_v0(0, 0, 0, 0);
// SK4
v0[52] = calc_v0(0, 2, 0, 0);
(f[52], u[52]) = (f_k2, u_k2);
// 2MN6
v0[53] = calc_v0(-7, 6, 1, 0);
(f[53], u[53]) = (f_m2.powi(3), 3. * u_m2);
// M6 太陰1/6日周潮
v0[54] = calc_v0(-6, 6, 0, 0);
(f[54], u[54]) = (f_m2.powi(3), 3. * u_m2);
// MSN6
v0[55] = calc_v0(-5, 4, 1, 0);
(f[55], u[55]) = (f_m2 * f_m2, 2. * u_m2);
// 2MS6
v0[56] = calc_v0(-4, 4, 0, 0);
(f[56], u[56]) = (f_m2 * f_m2, 2. * u_m2);
// 2MK6
v0[57] = calc_v0(-4, 6, 0, 0);
(f[57], u[57]) = (f_m2 * f_m2 * f_k2, 2. * u_m2 + u_k2);
// 2SM6
v0[58] = calc_v0(-2, 2, 0, 0);
(f[58], u[58]) = (f_m2, u_m2);
// MSK6
v0[59] = calc_v0(-2, 4, 0, 0);
(f[59], u[59]) = (f_m2 * f_k2, u_m2 + u_k2);
Self { v0, f, u }
}
}
fn main() {
let begin = NaiveDate::from_ymd_opt(2024, 1, 1).unwrap();
let mid = NaiveDate::from_ymd_opt(2024, 7, 1).unwrap();
let est = TideEstimator::new(begin, mid);
let day = NaiveDate::from_ymd_opt(2024, 1, 1).unwrap();
let lng = 136. + 52. / 60. + 51. / 3600.;
let z0 = 140.0;
for hour in 0..24 {
let h = est.estimate(
lng,
z0,
&PARAMS,
(day.signed_duration_since(begin).num_days() * 24 + hour) as f64,
);
println!("{hour}, {h}");
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment