feat: 添加 lagran 和 yint 模块

- lagran: Lagrange 三点插值函数
- yint: 二次插值函数（数组版本）

这两个纯数学函数用于数值插值。

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

src/math/lagran.rs

@@ -0,0 +1,81 @@
+//! Lagrange 三点插值。
+//!
+//! 重构自 TLUSTY `lagran.f`
+
+/// Lagrange 三点插值。
+///
+/// 给定三个点 (x0,y0), (x1,y1), (x2,y2)，计算 x 处的插值 y。
+///
+/// # 参数
+/// * `x0`, `y0` - 第一个点
+/// * `x1`, `y1` - 第二个点
+/// * `x2`, `y2` - 第三个点
+/// * `x` - 需要插值的位置
+///
+/// # 返回值
+/// 返回插值结果 y
+///
+/// # 注意
+/// 调用者需确保 x0, x1, x2 互不相等
+pub fn lagran(x0: f64, x1: f64, x2: f64, y0: f64, y1: f64, y2: f64, x: f64) -> f64 {
+    // Lagrange 基函数
+    let xl0 = (x - x1) * (x - x2) / ((x0 - x1) * (x0 - x2));
+    let xl1 = (x - x0) * (x - x2) / ((x1 - x0) * (x1 - x2));
+    let xl2 = (x - x0) * (x - x1) / ((x2 - x0) * (x2 - x1));
+
+    // 插值结果
+    y0 * xl0 + y1 * xl1 + y2 * xl2
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use approx::assert_relative_eq;
+
+    #[test]
+    fn test_lagran_at_points() {
+        // 测试在插值节点上的值
+        let y = lagran(0.0, 1.0, 2.0, 1.0, 2.0, 5.0, 0.0);
+        assert_relative_eq!(y, 1.0, epsilon = 1e-15);
+
+        let y = lagran(0.0, 1.0, 2.0, 1.0, 2.0, 5.0, 1.0);
+        assert_relative_eq!(y, 2.0, epsilon = 1e-15);
+
+        let y = lagran(0.0, 1.0, 2.0, 1.0, 2.0, 5.0, 2.0);
+        assert_relative_eq!(y, 5.0, epsilon = 1e-15);
+    }
+
+    #[test]
+    fn test_lagran_interpolation() {
+        // 测试线性插值
+        let y = lagran(0.0, 1.0, 2.0, 0.0, 1.0, 2.0, 0.5);
+        assert_relative_eq!(y, 0.5, epsilon = 1e-15);
+
+        // 测试二次插值 (y = x^2)
+        let y = lagran(0.0, 1.0, 2.0, 0.0, 1.0, 4.0, 0.5);
+        assert_relative_eq!(y, 0.25, epsilon = 1e-15);
+
+        let y = lagran(0.0, 1.0, 2.0, 0.0, 1.0, 4.0, 1.5);
+        assert_relative_eq!(y, 2.25, epsilon = 1e-15);
+    }
+
+    #[test]
+    fn test_lagran_extrapolation() {
+        // 测试外推
+        let y = lagran(0.0, 1.0, 2.0, 0.0, 1.0, 4.0, -1.0);
+        assert_relative_eq!(y, 1.0, epsilon = 1e-15); // (-1)^2 = 1
+
+        let y = lagran(0.0, 1.0, 2.0, 0.0, 1.0, 4.0, 3.0);
+        assert_relative_eq!(y, 9.0, epsilon = 1e-15); // 3^2 = 9
+    }
+
+    #[test]
+    fn test_lagran_negative_x() {
+        // 测试负 x 值
+        let y = lagran(-2.0, 0.0, 2.0, 4.0, 0.0, 4.0, 0.0);
+        assert_relative_eq!(y, 0.0, epsilon = 1e-15);
+
+        let y = lagran(-2.0, 0.0, 2.0, 4.0, 0.0, 4.0, -1.0);
+        assert_relative_eq!(y, 1.0, epsilon = 1e-15);
+    }
+}

src/math/mod.rs

@@ -112,6 +112,7 @@ mod intlem;
 mod intxen;
 mod irc;
 mod interpolate;
+mod lagran;
 mod laguer;
 mod lemini;
 mod levsol;
@@ -267,6 +268,7 @@ mod voigte;
 mod wn;
 mod wnstor;
 mod xk2dop;
+mod yint;
 mod ylintp;
 mod zmrho;
 

src/math/yint.rs

@@ -0,0 +1,80 @@
+//! 二次插值函数。
+//!
+//! 重构自 TLUSTY `yint.f`
+
+/// 二次插值。
+///
+/// 给定三个点 (xl[0],yl[0]), (xl[1],yl[1]), (xl[2],yl[2])，
+/// 计算 xl0 处的插值结果。
+///
+/// # 参数
+/// * `xl` - x 坐标数组（长度为 3）
+/// * `yl` - y 坐标数组（长度为 3）
+/// * `xl0` - 需要插值的位置
+///
+/// # 返回值
+/// 返回插值结果
+///
+/// # 注意
+/// 调用者需确保 xl[0], xl[1], xl[2] 互不相等
+pub fn yint(xl: &[f64], yl: &[f64], xl0: f64) -> f64 {
+    // Fortran 索引转换：XL(1) -> xl[0], XL(2) -> xl[1], XL(3) -> xl[2]
+    let a0 = (xl[1] - xl[0]) * (xl[2] - xl[1]) * (xl[2] - xl[0]);
+    let a1 = (xl0 - xl[1]) * (xl0 - xl[2]) * (xl[2] - xl[1]);
+    let a2 = (xl0 - xl[0]) * (xl[2] - xl0) * (xl[2] - xl[0]);
+    let a3 = (xl0 - xl[0]) * (xl0 - xl[1]) * (xl[1] - xl[0]);
+
+    (yl[0] * a1 + yl[1] * a2 + yl[2] * a3) / a0
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use approx::assert_relative_eq;
+
+    #[test]
+    fn test_yint_at_points() {
+        let xl = [0.0, 1.0, 2.0];
+        let yl = [1.0, 2.0, 5.0];
+
+        // 测试在插值节点上的值
+        assert_relative_eq!(yint(&xl, &yl, 0.0), 1.0, epsilon = 1e-15);
+        assert_relative_eq!(yint(&xl, &yl, 1.0), 2.0, epsilon = 1e-15);
+        assert_relative_eq!(yint(&xl, &yl, 2.0), 5.0, epsilon = 1e-15);
+    }
+
+    #[test]
+    fn test_yint_quadratic() {
+        // y = x^2
+        let xl = [0.0, 1.0, 2.0];
+        let yl = [0.0, 1.0, 4.0];
+
+        assert_relative_eq!(yint(&xl, &yl, 0.5), 0.25, epsilon = 1e-15);
+        assert_relative_eq!(yint(&xl, &yl, 1.5), 2.25, epsilon = 1e-15);
+    }
+
+    #[test]
+    fn test_yint_extrapolation() {
+        let xl = [0.0, 1.0, 2.0];
+        let yl = [0.0, 1.0, 4.0];
+
+        // 外推: y = x^2
+        assert_relative_eq!(yint(&xl, &yl, -1.0), 1.0, epsilon = 1e-15);
+        assert_relative_eq!(yint(&xl, &yl, 3.0), 9.0, epsilon = 1e-15);
+    }
+
+    #[test]
+    fn test_yint_vs_lagran() {
+        // yint 和 lagran 应该给出相同的结果
+        use super::super::lagran::lagran;
+
+        let xl = [0.0, 1.0, 2.0];
+        let yl = [1.0, 3.0, 2.0];
+
+        for x in [0.0, 0.5, 1.0, 1.5, 2.0, 3.0].iter() {
+            let y1 = yint(&xl, &yl, *x);
+            let y2 = lagran(xl[0], xl[1], xl[2], yl[0], yl[1], yl[2], *x);
+            assert_relative_eq!(y1, y2, epsilon = 1e-15);
+        }
+    }
+}