gistools/proj/project/
labrd.rs

1use crate::proj::{
2    AZIMUTH_PROJECTION_CENTRE, CoordinateStep, LABORDE, Proj, ProjValue, ProjectCoordinates,
3    TransformCoordinates,
4};
5use alloc::rc::Rc;
6use core::cell::RefCell;
7use libm::{atan, cos, exp, fabs, log, sin, sqrt, tan};
8
9const M_FORTPI: f64 = 4. / core::f64::consts::PI;
10const EPS: f64 = 1.0e-10;
11
12/// Laborde variables
13#[derive(Debug, Default, Clone, Copy, PartialEq)]
14pub struct LabordeData {
15    k_rg: f64,
16    p0s: f64,
17    a: f64,
18    c: f64,
19    ca: f64,
20    cb: f64,
21    cc: f64,
22    cd: f64,
23}
24
25/// Laborde Projection
26#[derive(Debug, Clone, PartialEq)]
27pub struct LabordeProjection {
28    proj: Rc<RefCell<Proj>>,
29    store: RefCell<LabordeData>,
30}
31impl ProjectCoordinates for LabordeProjection {
32    fn code(&self) -> i64 {
33        LABORDE
34    }
35    fn name(&self) -> &'static str {
36        "Laborde"
37    }
38    fn names() -> &'static [&'static str] {
39        &["Laborde", "Laborde Oblique Mercator", "labrd"]
40    }
41}
42impl CoordinateStep for LabordeProjection {
43    fn new(proj: Rc<RefCell<Proj>>) -> Self {
44        let mut store = LabordeData::default();
45        {
46            let proj = &mut proj.borrow_mut();
47            if proj.phi0 == 0. {
48                panic!("Invalid value for lat_0: lat_0 should be different from 0");
49            }
50
51            let az = proj
52                .params
53                .get(&AZIMUTH_PROJECTION_CENTRE) // (lat_ts)
54                .unwrap_or(&ProjValue::default())
55                .f64();
56            let sinp = sin(proj.phi0);
57            let mut t = 1. - proj.es * sinp * sinp;
58            let _n = 1. / sqrt(t);
59            let _r = proj.one_es * _n / t;
60            store.k_rg = proj.k0 * sqrt(_n * _r);
61            store.p0s = atan(sqrt(_r / _n) * tan(proj.phi0));
62            store.a = sinp / sin(store.p0s);
63            t = proj.e * sinp;
64            store.c = 0.5 * proj.e * store.a * log((1. + t) / (1. - t))
65                + -store.a * log(tan(M_FORTPI + 0.5 * proj.phi0))
66                + log(tan(M_FORTPI + 0.5 * store.p0s));
67            t = az + az;
68            store.cb = 1. / (12. * store.k_rg * store.k_rg);
69            store.ca = (1. - cos(t)) * store.cb;
70            store.cb *= sin(t);
71            store.cc = 3. * (store.ca * store.ca - store.cb * store.cb);
72            store.cd = 6. * store.ca * store.cb;
73        }
74        LabordeProjection { proj, store: store.into() }
75    }
76    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
77        labrd_e_forward(&self.store.borrow(), &self.proj.borrow(), p);
78    }
79    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
80        labrd_e_inverse(&self.store.borrow(), &self.proj.borrow(), p);
81    }
82}
83
84/// Laborde Ellipsoidal forward project
85pub fn labrd_e_forward<P: TransformCoordinates>(laborde: &LabordeData, proj: &Proj, p: &mut P) {
86    let phi = p.phi();
87    let lam = p.lam();
88
89    let mut v1 = laborde.a * log(tan(M_FORTPI + 0.5 * phi));
90    let mut t = proj.e * sin(phi);
91    let mut v2 = 0.5 * proj.e * laborde.a * log((1. + t) / (1. - t));
92    let ps = 2. * (atan(exp(v1 - v2 + laborde.c)) - M_FORTPI);
93    let i1 = ps - laborde.p0s;
94    let cosps = cos(ps);
95    let cosps2 = cosps * cosps;
96    let sinps = sin(ps);
97    let sinps2 = sinps * sinps;
98    let i4 = laborde.a * cosps;
99    let i2 = 0.5 * laborde.a * i4 * sinps;
100    let i3 = i2 * laborde.a * laborde.a * (5. * cosps2 - sinps2) / 12.;
101    let mut i6 = i4 * laborde.a * laborde.a;
102    let i5 = i6 * (cosps2 - sinps2) / 6.;
103    i6 *= laborde.a * laborde.a * (5. * cosps2 * cosps2 + sinps2 * (sinps2 - 18. * cosps2)) / 120.;
104    t = lam * lam;
105    let mut x = laborde.k_rg * lam * (i4 + t * (i5 + t * i6));
106    let mut y = laborde.k_rg * (i1 + t * (i2 + t * i3));
107    let x2 = x * x;
108    let y2 = y * y;
109    v1 = 3. * x * y2 - x * x2;
110    v2 = y * y2 - 3. * x2 * y;
111    x += laborde.ca * v1 + laborde.cb * v2;
112    y += laborde.ca * v2 - laborde.cb * v1;
113
114    p.set_x(x);
115    p.set_y(y);
116}
117
118/// Laborde Ellipsoidal inverse project
119pub fn labrd_e_inverse<P: TransformCoordinates>(laborde: &LabordeData, proj: &Proj, p: &mut P) {
120    // t = 0.0 optimization is to avoid a false positive cppcheck warning
121    let mut t;
122    // (cppcheck git beaf29c15867984aa3c2a15cf15bd7576ccde2b3). Might no
123    // longer be necessary with later versions.
124    let mut x = p.x();
125    let mut y = p.y();
126
127    let mut x2 = x * x;
128    let y2 = y * y;
129    let mut v1 = 3. * x * y2 - x * x2;
130    let mut v2 = y * y2 - 3. * x2 * y;
131    let v3 = x * (5. * y2 * y2 + x2 * (-10. * y2 + x2));
132    let v4 = y * (5. * x2 * x2 + y2 * (-10. * x2 + y2));
133    x += -laborde.ca * v1 - laborde.cb * v2 + laborde.cc * v3 + laborde.cd * v4;
134    y += laborde.cb * v1 - laborde.ca * v2 - laborde.cd * v3 + laborde.cc * v4;
135    let ps = laborde.p0s + y / laborde.k_rg;
136    let mut pe = ps + proj.phi0 - laborde.p0s;
137
138    for _ in 0..20 {
139        v1 = laborde.a * log(tan(M_FORTPI + 0.5 * pe));
140        let tpe = proj.e * sin(pe);
141        v2 = 0.5 * proj.e * laborde.a * log((1. + tpe) / (1. - tpe));
142        t = ps - 2. * (atan(exp(v1 - v2 + laborde.c)) - M_FORTPI);
143        pe += t;
144        if fabs(t) < EPS {
145            break;
146        }
147    }
148
149    t = proj.e * sin(pe);
150    t = 1. - t * t;
151    let re = proj.one_es / (t * sqrt(t));
152    t = tan(ps);
153    let t2 = t * t;
154    let s = laborde.k_rg * laborde.k_rg;
155    let mut d = re * proj.k0 * laborde.k_rg;
156    let i7 = t / (2. * d);
157    let i8 = t * (5. + 3. * t2) / (24. * d * s);
158    d = cos(ps) * laborde.k_rg * laborde.a;
159    let i9 = 1. / d;
160    d *= s;
161    let i10 = (1. + 2. * t2) / (6. * d);
162    let i11 = (5. + t2 * (28. + 24. * t2)) / (120. * d * s);
163    x2 = x * x;
164
165    p.set_phi(pe + x2 * (-i7 + i8 * x2));
166    p.set_lam(x * (i9 + x2 * (-i10 + x2 * i11)));
167}
gistools/proj/project/labrd.rs

gistools/proj/project/
labrd.rs
