天文潮位を計算する (Rust)

参考: https://www1.kaiho.mlit.go.jp/kenkyu/report/rhr60/rhr60_r_01.pdf 潮汐表の計算について 海洋情報部研究報告 第60号 令和4年3月18日

↑の計算を Rust で実装してみました。資料だけみても読み解きにくい部分があるので、同様のことをされる方のために、Gistのかわりにここにコードを貼っておきます。以上です。

//! 参考: 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
// https://www.data.jma.go.jp/kaiyou/db/tide/suisan/harms60.php?stn=NG&year=2025&tyear=2025
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}");
    }
}