@@ -5,8 +5,6 @@ use crate::authoring::*;
5
5
use std:: f64:: consts:: FRAC_PI_2 ;
6
6
const EPS10 : f64 = 1e-10 ;
7
7
8
- // ----- C O M M O N -------------------------------------------------------------------
9
-
10
8
// ----- F O R W A R D -----------------------------------------------------------------
11
9
12
10
fn fwd ( op : & Op , _ctx : & dyn Context , operands : & mut dyn CoordinateSet ) -> usize {
@@ -40,49 +38,41 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
40
38
41
39
// The polar aspects are fairly simple
42
40
if north_polar || south_polar {
41
+ let sign = if north_polar { -1.0 } else { 1.0 } ;
43
42
for i in 0 ..n {
44
- let mut coord = operands. get_coord ( i) ;
45
- let sign = if north_polar { -1.0 } else { 1.0 } ;
46
-
47
- let lat = coord[ 1 ] ;
48
- let lon = coord[ 0 ] ;
43
+ let ( lon, lat) = operands. xy ( i) ;
49
44
let ( sin_lon, cos_lon) = ( lon - lon_0) . sin_cos ( ) ;
50
45
51
46
let q = ancillary:: qs ( lat. sin ( ) , e) ;
52
47
let rho = a * ( qp + sign * q) . sqrt ( ) ;
53
48
54
- coord [ 0 ] = x_0 + rho * sin_lon;
55
- coord [ 1 ] = y_0 + sign * rho * cos_lon;
56
- operands. set_coord ( i, & coord ) ;
49
+ let easting = x_0 + rho * sin_lon;
50
+ let northing = y_0 + sign * rho * cos_lon;
51
+ operands. set_xy ( i, easting , northing ) ;
57
52
successes += 1 ;
58
53
}
59
54
return successes;
60
55
}
61
56
57
+ // Either equatorial or oblique aspects
62
58
for i in 0 ..n {
63
- let mut coord = operands. get_coord ( i) ;
64
- let lon = coord[ 0 ] ;
65
- let lat = coord[ 1 ] ;
59
+ let ( lon, lat) = operands. xy ( i) ;
66
60
let ( sin_lon, cos_lon) = ( lon - lon_0) . sin_cos ( ) ;
67
61
68
62
// Authalic latitude, 𝜉
69
63
let xi = ( ancillary:: qs ( lat. sin ( ) , e) / qp) . asin ( ) ;
70
-
71
64
let ( sin_xi, cos_xi) = xi. sin_cos ( ) ;
65
+
72
66
let b = if oblique {
73
67
let factor = 1.0 + sin_xi_0 * sin_xi + ( cos_xi_0 * cos_xi * cos_lon) ;
74
68
rq * ( 2.0 / factor) . sqrt ( )
75
69
} else {
76
70
1.0
77
71
} ;
78
72
79
- // Easting
80
- coord[ 0 ] = x_0 + ( b * d) * ( cos_xi * sin_lon) ;
81
-
82
- // Northing
83
- coord[ 1 ] = y_0 + ( b / d) * ( cos_xi_0 * sin_xi - sin_xi_0 * cos_xi * cos_lon) ;
84
- operands. set_coord ( i, & coord) ;
85
-
73
+ let easting = x_0 + ( b * d) * ( cos_xi * sin_lon) ;
74
+ let northing = y_0 + ( b / d) * ( cos_xi_0 * sin_xi - sin_xi_0 * cos_xi * cos_lon) ;
75
+ operands. set_xy ( i, easting, northing) ;
86
76
successes += 1 ;
87
77
}
88
78
@@ -122,39 +112,33 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
122
112
let mut successes = 0_usize ;
123
113
let n = operands. len ( ) ;
124
114
125
- // The polar aspects are not as simple as in the forward case
115
+ // The polar aspects are not quite as simple as in the forward case
126
116
if north_polar || south_polar {
117
+ let sign = if north_polar { -1.0 } else { 1.0 } ;
127
118
for i in 0 ..n {
128
- let mut coord = operands. get_coord ( i) ;
129
- let sign = if north_polar { -1.0 } else { 1.0 } ;
130
-
131
- let x = coord[ 0 ] ;
132
- let y = coord[ 1 ] ;
119
+ let ( x, y) = operands. xy ( i) ;
133
120
let rho = ( x - x_0) . hypot ( y - y_0) ;
134
121
135
122
// The authalic latitude is a bit convoluted
136
123
let denom = a * a * ( 1.0 - ( ( 1.0 - es) / ( 2.0 * e) ) * ( ( 1.0 - e) / ( 1.0 + e) ) . ln ( ) ) ;
137
124
let xi = ( -sign) * ( 1.0 - rho * rho / denom) ;
138
125
139
- coord [ 0 ] = lon_0 + ( x - x_0) . atan2 ( sign * ( y - y_0) ) ;
140
- coord [ 1 ] = ellps. latitude_authalic_to_geographic ( xi, & authalic) ;
141
- operands. set_coord ( i, & coord ) ;
126
+ let lon = lon_0 + ( x - x_0) . atan2 ( sign * ( y - y_0) ) ;
127
+ let lat = ellps. latitude_authalic_to_geographic ( xi, & authalic) ;
128
+ operands. set_xy ( i, lon , lat ) ;
142
129
successes += 1 ;
143
130
}
144
131
return successes;
145
132
}
146
133
134
+ // Either equatorial or oblique aspects
147
135
for i in 0 ..n {
148
- let mut coord = operands. get_coord ( i) ;
149
- let x = coord[ 0 ] ;
150
- let y = coord[ 1 ] ;
136
+ let ( x, y) = operands. xy ( i) ;
151
137
let rho = ( ( x - x_0) / d) . hypot ( d * ( y - y_0) ) ;
152
138
153
139
// A bit of reality hardening ported from the PROJ implementation
154
140
if rho < EPS10 {
155
- coord[ 0 ] = lon_0;
156
- coord[ 1 ] = lat_0;
157
- operands. set_coord ( i, & coord) ;
141
+ operands. set_xy ( i, lon_0, lat_0) ;
158
142
successes += 1 ;
159
143
continue ;
160
144
}
@@ -163,23 +147,22 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
163
147
let asin_argument = 0.5 * rho / rq;
164
148
if asin_argument. abs ( ) > 1.0 {
165
149
debug ! ( "LAEA: ({x}, {y}) outside domain" ) ;
166
- coord[ 0 ] = f64:: NAN ;
167
- coord[ 1 ] = f64:: NAN ;
168
- operands. set_coord ( i, & coord) ;
150
+ operands. set_xy ( i, f64:: NAN , f64:: NAN ) ;
169
151
continue ;
170
152
}
171
153
172
154
let c = 2.0 * asin_argument. asin ( ) ;
173
155
let ( sin_c, cos_c) = c. sin_cos ( ) ;
156
+
174
157
// The authalic latitude, 𝜉
175
158
let xi = ( cos_c * sin_xi_0 + ( d * ( y - y_0) * sin_c * cos_xi_0) / rho) . asin ( ) ;
176
- coord[ 1 ] = ellps. latitude_authalic_to_geographic ( xi, & authalic) ;
159
+
160
+ let lat = ellps. latitude_authalic_to_geographic ( xi, & authalic) ;
177
161
178
162
let num = ( x - x_0) * sin_c;
179
163
let denom = d * rho * cos_xi_0 * cos_c - d * d * ( y - y_0) * sin_xi_0 * sin_c;
180
- coord[ 0 ] = num. atan2 ( denom) + lon_0;
181
- operands. set_coord ( i, & coord) ;
182
-
164
+ let lon = num. atan2 ( denom) + lon_0;
165
+ operands. set_xy ( i, lon, lat) ;
183
166
successes += 1 ;
184
167
}
185
168
@@ -219,12 +202,12 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error>
219
202
220
203
let polar = ( t - FRAC_PI_2 ) . abs ( ) < EPS10 ;
221
204
let north = polar && ( t > 0.0 ) ;
222
- let equatoreal = !polar && t < EPS10 ;
223
- let oblique = !polar && !equatoreal ;
224
- match ( polar, equatoreal , north) {
205
+ let equatorial = !polar && t < EPS10 ;
206
+ let oblique = !polar && !equatorial ;
207
+ match ( polar, equatorial , north) {
225
208
( true , _, true ) => params. boolean . insert ( "north_polar" ) ,
226
209
( true , _, false ) => params. boolean . insert ( "south_polar" ) ,
227
- ( _, true , _) => params. boolean . insert ( "equatoreal " ) ,
210
+ ( _, true , _) => params. boolean . insert ( "equatorial " ) ,
228
211
_ => params. boolean . insert ( "oblique" ) ,
229
212
} ;
230
213
@@ -247,7 +230,7 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error>
247
230
// D in the IOGP text
248
231
let d = if oblique {
249
232
a * ( cos_phi_0 / ( 1.0 - es * sin_phi_0 * sin_phi_0) . sqrt ( ) ) / ( rq * xi_0. cos ( ) )
250
- } else if equatoreal {
233
+ } else if equatorial {
251
234
rq. recip ( )
252
235
} else {
253
236
a
0 commit comments