gistools/proj/project/
airy.rs

1use crate::proj::{
2    CoordinateStep, LAT_B, NO_CUT, Proj, ProjMode, ProjValue, ProjectCoordinates,
3    TransformCoordinates,
4};
5use alloc::rc::Rc;
6use core::{cell::RefCell, f64::consts::FRAC_PI_2};
7use libm::{cos, fabs, log, sin, tan};
8
9/// Airy variables
10#[derive(Debug, Default, Clone, Copy, PartialEq)]
11pub struct Airy {
12    p_halfpi: f64,
13    sinph0: f64,
14    cosph0: f64,
15    cb: f64,
16    mode: ProjMode,
17    // do not cut at hemisphere limit
18    no_cut: bool,
19}
20
21const EPS: f64 = 1e-10;
22
23/// Airy Projection
24#[derive(Debug, Clone, PartialEq)]
25pub struct AiryProjection {
26    proj: Rc<RefCell<Proj>>,
27    store: RefCell<Airy>,
28}
29impl ProjectCoordinates for AiryProjection {
30    fn code(&self) -> i64 {
31        -1
32    }
33    fn name(&self) -> &'static str {
34        "Airy"
35    }
36    fn names() -> &'static [&'static str] {
37        &["Airy", "airy"]
38    }
39}
40impl CoordinateStep for AiryProjection {
41    fn new(proj: Rc<RefCell<Proj>>) -> Self {
42        let mut store = Airy {
43            no_cut: proj.borrow().params.get(&NO_CUT).unwrap_or(&ProjValue::default()).bool(),
44            ..Default::default()
45        };
46        {
47            let proj = &mut proj.borrow_mut();
48            let beta =
49                0.5 * (FRAC_PI_2 - proj.params.get(&LAT_B).unwrap_or(&ProjValue::default()).f64());
50            if fabs(beta) < EPS {
51                store.cb = -0.5;
52            } else {
53                store.cb = 1. / tan(beta);
54                store.cb = store.cb * store.cb * log(cos(beta));
55            }
56
57            if fabs(fabs(proj.phi0) - FRAC_PI_2) < EPS {
58                if proj.phi0 < 0. {
59                    store.p_halfpi = -FRAC_PI_2;
60                    store.mode = ProjMode::SPole;
61                } else {
62                    store.p_halfpi = FRAC_PI_2;
63                    store.mode = ProjMode::NPole;
64                }
65            } else if fabs(proj.phi0) < EPS {
66                store.mode = ProjMode::Equit;
67            } else {
68                store.mode = ProjMode::Obliq;
69                store.sinph0 = sin(proj.phi0);
70                store.cosph0 = cos(proj.phi0);
71            }
72            proj.es = 0.;
73        }
74        AiryProjection { proj, store: store.into() }
75    }
76    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
77        airy_s_forward(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
78    }
79    fn inverse<P: TransformCoordinates>(&self, _p: &mut P) {
80        // There is no inverse
81    }
82}
83
84/// Ellipsoid/spheroid, forward
85pub fn airy_s_forward<P: TransformCoordinates>(airy: &mut Airy, _proj: &Proj, p: &mut P) {
86    let sinlam = sin(p.lam());
87    let coslam = cos(p.lam());
88    match airy.mode {
89        ProjMode::Equit | ProjMode::Obliq => {
90            let sinphi = sin(p.phi());
91            let cosphi = cos(p.phi());
92            let mut cosz = cosphi * coslam;
93            if airy.mode == ProjMode::Obliq {
94                cosz = airy.sinph0 * sinphi + airy.cosph0 * cosz;
95            }
96            if !airy.no_cut && cosz < -EPS {
97                panic!("Coordinates are outside the projection domain");
98            }
99            let s = 1. - cosz;
100            let k_rho = if fabs(s) > EPS {
101                let t = 0.5 * (1. + cosz);
102                if t == 0. {
103                    panic!("Coordinates are outside the projection domain");
104                }
105                -log(t) / s - airy.cb / t
106            } else {
107                0.5 - airy.cb
108            };
109            p.set_x(k_rho * cosphi * sinlam);
110            if airy.mode == ProjMode::Obliq {
111                p.set_y(k_rho * (airy.cosph0 * sinphi - airy.sinph0 * cosphi * coslam));
112            } else {
113                p.set_y(k_rho * sinphi);
114            }
115        }
116        ProjMode::SPole | ProjMode::NPole => {
117            p.set_phi(fabs(airy.p_halfpi - p.phi()));
118            if !airy.no_cut && (p.phi() - EPS) > FRAC_PI_2 {
119                panic!("Coordinates are outside the projection domain");
120            }
121            p.set_phi(p.phi() * 0.5);
122            if p.phi() > EPS {
123                let t = tan(p.phi());
124                let k_rho = -2. * (log(cos(p.phi())) / t + t * airy.cb);
125                p.set_x(k_rho * sinlam);
126                p.set_y(k_rho * coslam);
127                if airy.mode == ProjMode::NPole {
128                    p.set_y(-p.y());
129                }
130            } else {
131                p.set_x(0.);
132                p.set_y(0.);
133            }
134        }
135    }
136}
gistools/proj/project/airy.rs

gistools/proj/project/
airy.rs
