A small program to calculate the area of a triangle on a spherical surface using spherical geometry.

spherical_coordinate_system.rs

use nalgebra::Point3;
use std::f64::consts::PI;

/// Reference: http://sshmathgeom.private.coocan.jp/geometryonsphere.pdf

pub struct SphericalCoordinate {
    pub theta: f64, // latitude: -half pi <= theta <= half pi
    pub phi: f64,   // longitude: -2pi <= phi <= 2pi
}

impl SphericalCoordinate {
    /// Returns the angle ∠BAC
    /// define the angle between ∠BAC by tangent line in arc AB and tangent line in arc AC
    /// cos(θ) = (OA x OB) . (OA x OC) / |OA x OB||OB x OC|
    pub fn angle(&self, b: &Self, c: &Self) -> f64 {
        let a = self.point(1.).coords;
        let b = b.point(1.).coords;
        let c = c.point(1.).coords;

        // normal by oab plane
        let ab = a.cross(&b);

        // normal by oac plane
        let ac = a.cross(&c);

        let cos = ab.dot(&ac) / (ab.magnitude() * ac.magnitude());
        cos.acos()
    }

    pub fn point(&self, r: f64) -> Point3<f64> {
        let t_cos = self.theta.cos();
        let p_cos = self.phi.cos();
        let t_sin = self.theta.sin();
        let p_sin = self.phi.sin();
        let x = r * t_cos * p_cos;
        let y = r * t_cos * p_sin;
        let z = r * t_sin;
        Point3::new(x, y, z)
    }
}

pub struct SphericalTriangle {
    a: SphericalCoordinate,
    b: SphericalCoordinate,
    c: SphericalCoordinate,
}

impl SphericalTriangle {
    /// Returns the area of the triangle ABC
    /// (∠BAC + ∠ABC + ∠ACB - π) * r^2
    pub fn area(&self, r: f64) -> f64 {
        let a = self.a.angle(&self.b, &self.c);
        let b = self.b.angle(&self.a, &self.c);
        let c = self.c.angle(&self.a, &self.b);
        (a + b + c - PI) * r * r
    }
}

#[cfg(test)]
mod tests {
    use std::f64::consts::FRAC_PI_2;

    use nalgebra::Point3;

    use crate::{SphericalCoordinate, SphericalTriangle};

    #[test]
    fn test() {
        // let r = 6378.137;
        let r = 1.;
        let tokyo = SphericalCoordinate {
            theta: (35_f64 + 40. / 60.).to_radians(),
            phi: (139_f64 + 45. / 60.).to_radians(),
        };
        let diff = tokyo.point(r) - Point3::new(-0.6200675211, 0.5249259043, 0.5830686615);
        assert!(diff.magnitude() < 1e-8);
        let nyc = SphericalCoordinate {
            theta: (40_f64 + 47. / 60.).to_radians(),
            phi: (-73_f64 - 58. / 60.).to_radians(),
        };
        dbg!(nyc.point(r));

        let southpole = SphericalCoordinate {
            theta: -FRAC_PI_2,
            phi: 0.,
        };
        dbg!(southpole.point(r));

        let angle = tokyo.angle(&nyc, &southpole);
        dbg!(angle.to_degrees());
        let diff = angle.to_degrees() - 154.91600061;
        assert!(diff.abs() < 1e-8);

        let triangle = SphericalTriangle {
            a: tokyo,
            b: nyc,
            c: southpole,
        };
        let area = triangle.area(r);
        dbg!(area);
        let diff = area - 4.7846893065;
        assert!(diff.abs() < 1e-8);
    }
}