source: subversion/applications/editors/osm-editor/qt3/llgr.cpp @ 16590

Last change on this file since 16590 was 1158, checked in by nick, 14 years ago

moved current version to qt3 directory in preparation for qt4

File size: 10.8 KB
Line 
1
2// NOTICE OF AUTHORSHIP: These functions are taken from the LGPL JEEPS library,
3// author and licence info follows....
4
5/********************************************************************
6** @source JEEPS arithmetic/conversion functions
7**
8** @author Copyright (C) 1999 Alan Bleasby
9** @version 1.0
10** @modified Dec 28 1999 Alan Bleasby. First version
11** @@
12**
13** This library is free software; you can redistribute it and/or
14** modify it under the terms of the GNU Library General Public
15** License as published by the Free Software Foundation; either
16** version 2 of the License, or (at your option) any later version.
17**
18** This library is distributed in the hope that it will be useful,
19** but WITHOUT ANY WARRANTY; without even the implied warranty of
20** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
21** Library General Public License for more details.
22**
23** You should have received a copy of the GNU Library General Public
24** License along with this library; if not, write to the
25** Free Software Foundation, Inc., 59 Temple Place - Suite 330,
26** Boston, MA  02111-1307, USA.
27********************************************************************/
28
29#include "llgr.h"
30#include <cmath>
31#include <sstream>
32#include <cctype>
33
34namespace OpenStreetMap
35{
36
37// ll_to_gr()
38// Slightly modified version of:
39/* @func modGPS_Math_Airy1830LatLonToNGEN **************************************
40**
41** Convert Airy 1830 datum latitude and longitude to UK Ordnance Survey
42** National Grid Eastings and Northings
43**
44** @param [r] phi [double] WGS84 latitude     (deg)
45** @param [r] lambda [double] WGS84 longitude (deg)
46** @param [w] E [double *] NG easting (metres)
47** @param [w] N [double *] NG northing (metres)
48**
49** @return [void]
50*/
51
52// Modifications:
53// - Returns a Location structure.
54// - Parameters renamed from "phi" and "lambda" to "lat" and "lng".
55
56EarthPoint ll_to_gr ( const EarthPoint& ll )
57{
58        return ll_to_gr(ll.y,ll.x);
59}
60
61EarthPoint ll_to_gr ( double lat, double lon ) 
62{
63
64    double N0      = -100000;
65    double E0      =  400000;
66    double F0      = 0.9996012717;
67    double phi0    = 49.;
68    double lambda0 = -2.;
69    double a       = 6377563.396;
70    double b       = 6356256.910;
71        double e, n;
72
73    modGPS_Math_LatLon_To_EN(&e,&n,lat,lon,N0,E0,phi0,lambda0,F0,
74                                                        a,b);
75
76    return EarthPoint(e,n);
77}
78
79/* @func  modGPS_Math_LatLon_To_EN **********************************
80**
81** Convert latitude and longitude to eastings and northings
82** Standard Gauss-Kruger Transverse Mercator
83**
84** @param [w] E [double *] easting (metres)
85** @param [w] N [double *] northing (metres)
86** @param [r] phi [double] latitude (deg)
87** @param [r] lambda [double] longitude (deg)
88** @param [r] N0 [double] true northing origin (metres)
89** @param [r] E0 [double] true easting  origin (metres)
90** @param [r] phi0 [double] true latitude origin (deg)
91** @param [r] lambda0 [double] true longitude origin (deg)
92** @param [r] F0 [double] scale factor on central meridian
93** @param [r] a [double] semi-major axis (metres)
94** @param [r] b [double] semi-minor axis (metres)
95**
96** @return [void]
97************************************************************************/
98void modGPS_Math_LatLon_To_EN(double *E, double *N, double phi,
99                           double lambda, double N0, double E0,
100                           double phi0, double lambda0,
101                           double F0, double a, double b)
102{
103    double esq;
104    double n;
105    double etasq;
106    double nu;
107    double rho;
108    double M;
109    double I;
110    double II;
111    double III;
112    double IIIA;
113    double IV;
114    double V;
115    double VI;
116   
117    double tmp;
118    double tmp2;
119    double fdf;
120    double fde;
121   
122        // Replaced these deg/rad conversions (originally done with a function)
123        // with simple multiplication
124    phi0    *= (M_PI/180.0); 
125    lambda0 *= (M_PI/180.0); 
126    phi     *= (M_PI/180.0); 
127    lambda  *= (M_PI/180.0); 
128   
129    esq = ((a*a)-(b*b)) / (a*a);
130    n   = (a-b) / (a+b);
131   
132    tmp  = (double)1.0 - (esq * sin(phi) * sin(phi));
133    nu   = a * F0 * pow(tmp,(double)-0.5);
134    rho  = a * F0 * ((double)1.0 - esq) * pow(tmp,(double)-1.5);
135    etasq = (nu / rho) - (double)1.0;
136
137    fdf   = (double)5.0 / (double)4.0;
138    tmp   = (double)1.0 + n + (fdf * n * n) + (fdf * n * n * n);
139    tmp  *= (phi - phi0);
140    tmp2  = (double)3.0*n + (double)3.0*n*n + ((double)21./(double)8.)*n*n*n;
141    tmp2 *= (sin(phi-phi0) * cos(phi+phi0));
142    tmp  -= tmp2;
143
144    fde   = ((double)15.0 / (double)8.0);
145    tmp2  = ((fde*n*n) + (fde*n*n*n)) * sin((double)2.0 * (phi-phi0));
146    tmp2 *= cos((double)2.0 * (phi+phi0));
147    tmp  += tmp2;
148   
149    tmp2  = ((double)35.0/(double)24.0) * n * n * n;
150    tmp2 *= sin((double)3.0 * (phi-phi0));
151    tmp2 *= cos((double)3.0 * (phi+phi0));
152    tmp  -= tmp2;
153
154    M     = b * F0 * tmp;
155    I     = M + N0;
156    II    = (nu / (double)2.0) * sin(phi) * cos(phi);
157    III   = (nu / (double)24.0) * sin(phi) * cos(phi) * cos(phi) * cos(phi);
158    III  *= ((double)5.0 - (tan(phi) * tan(phi)) + ((double)9.0 * etasq));
159    IIIA  = (nu / (double)720.0) * sin(phi) * pow(cos(phi),(double)5.0);
160    IIIA *= ((double)61.0 - ((double)58.0*tan(phi)*tan(phi)) +
161             pow(tan(phi),(double)4.0));
162    IV    = nu * cos(phi);
163
164    tmp   = pow(cos(phi),(double)3.0);
165    tmp  *= ((nu/rho) - tan(phi) * tan(phi));
166    V     = (nu/(double)6.0) * tmp;
167
168    tmp   = (double)5.0 - ((double)18.0 * tan(phi) * tan(phi));
169    tmp  += tan(phi)*tan(phi)*tan(phi)*tan(phi) + ((double)14.0 * etasq);
170    tmp  -= ((double)58.0 * tan(phi) * tan(phi) * etasq);
171    tmp2  = cos(phi)*cos(phi)*cos(phi)*cos(phi)*cos(phi) * tmp;
172    VI    = (nu / (double)120.0) * tmp2;
173   
174    *N = I + II*(lambda-lambda0)*(lambda-lambda0) +
175             III*pow((lambda-lambda0),(double)4.0) +
176             IIIA*pow((lambda-lambda0),(double)6.0);
177
178    *E = E0 + IV*(lambda-lambda0) + V*pow((lambda-lambda0),(double)3.0) +
179         VI * pow((lambda-lambda0),(double)5.0);
180
181    return;
182}
183
184// gr_to_ll()
185// Slightly modified version of:
186/* @func modGPS_Math_NGENToAiry1830LatLon **************************************
187**
188** Convert  to UK Ordnance Survey National Grid Eastings and Northings to
189** Airy 1830 datum latitude and longitude
190**
191** @param [r] E [double] NG easting (metres)
192** @param [r] N [double] NG northing (metres)
193** @param [w] phi [double *] Airy latitude     (deg)
194** @param [w] lambda [double *] Airy longitude (deg)
195**
196** @return [void]
197************************************************************************/
198
199// Modifications:
200// - Grid reference passed as a struct of type Location.
201// - Lat/long parameters renamed from "phi" and "lambda" to "lat" and "lng".
202
203EarthPoint gr_to_ll(const EarthPoint& gridref)
204{
205    double N0      = -100000;
206    double E0      =  400000;
207    double F0      = 0.9996012717;
208    double phi0    = 49.;
209    double lambda0 = -2.;
210    double a       = 6377563.396;
211    double b       = 6356256.910;
212
213        EarthPoint retval;
214
215    modGPS_Math_EN_To_LatLon(gridref.x,gridref.y,&retval.y,&retval.x,N0,E0,
216                                        phi0,lambda0,F0, a,b);
217   
218    return retval;
219}
220
221/* @func  modGPS_Math_EN_To_LatLon **************************************
222**
223** Convert Eastings and Northings to latitude and longitude
224**
225** @param [w] E [double] NG easting (metres)
226** @param [w] N [double] NG northing (metres)
227** @param [r] phi [double *] Airy latitude     (deg)
228** @param [r] lambda [double *] Airy longitude (deg)
229** @param [r] N0 [double] true northing origin (metres)
230** @param [r] E0 [double] true easting  origin (metres)
231** @param [r] phi0 [double] true latitude origin (deg)
232** @param [r] lambda0 [double] true longitude origin (deg)
233** @param [r] F0 [double] scale factor on central meridian
234** @param [r] a [double] semi-major axis (metres)
235** @param [r] b [double] semi-minor axis (metres)
236**
237** @return [void]
238************************************************************************/
239void modGPS_Math_EN_To_LatLon(double E, double N, double *phi,
240                           double *lambda, double N0, double E0,
241                           double phi0, double lambda0,
242                           double F0, double a, double b)
243{
244    double esq;
245    double n;
246    double etasq;
247    double nu;
248    double rho;
249    double M;
250    double VII;
251    double VIII;
252    double IX;
253    double X;
254    double XI;
255    double XII;
256    double XIIA;
257    double phix;
258    double nphi=0.0;
259   
260    double tmp;
261    double tmp2;
262    double fdf;
263    double fde;
264
265        // Replaced these deg/rad conversions (originally done with a function)
266        // with simple multiplication
267    phi0    *= (M_PI/180.0); 
268    lambda0 *= (M_PI/180.0); 
269
270    n     = (a-b) / (a+b);
271    fdf   = (double)5.0 / (double)4.0;
272    fde   = ((double)15.0 / (double)8.0);
273
274    esq = ((a*a)-(b*b)) / (a*a);
275
276
277    phix = ((N-N0)/(a*F0)) + phi0;
278   
279    tmp  = (double)1.0 - (esq * sin(phix) * sin(phix));
280    nu   = a * F0 * pow(tmp,(double)-0.5);
281    rho  = a * F0 * ((double)1.0 - esq) * pow(tmp,(double)-1.5);
282    etasq = (nu / rho) - (double)1.0;
283
284    M = (double)-1e20;
285
286    while(N-N0-M > (double)0.000001)
287    {
288        nphi = phix;
289       
290        tmp   = (double)1.0 + n + (fdf * n * n) + (fdf * n * n * n);
291        tmp  *= (nphi - phi0);
292        tmp2  = (double)3.0*n + (double)3.0*n*n +
293                ((double)21./(double)8.)*n*n*n;
294        tmp2 *= (sin(nphi-phi0) * cos(nphi+phi0));
295        tmp  -= tmp2;
296
297
298        tmp2  = ((fde*n*n) + (fde*n*n*n)) * sin((double)2.0 * (nphi-phi0));
299        tmp2 *= cos((double)2.0 * (nphi+phi0));
300        tmp  += tmp2;
301   
302        tmp2  = ((double)35.0/(double)24.0) * n * n * n;
303        tmp2 *= sin((double)3.0 * (nphi-phi0));
304        tmp2 *= cos((double)3.0 * (nphi+phi0));
305        tmp  -= tmp2;
306
307        M     = b * F0 * tmp;
308
309        if(N-N0-M > (double)0.000001)
310            phix = ((N-N0-M)/(a*F0)) + nphi;
311    }
312   
313
314    VII  = tan(nphi) / ((double)2.0 * rho * nu);
315
316    tmp  = (double)5.0 + (double)3.0 * tan(nphi) * tan(nphi) + etasq;
317    tmp -= (double)9.0 * tan(nphi) * tan(nphi) * etasq;
318    VIII = (tan(nphi)*tmp) / ((double)24.0 * rho * nu*nu*nu);
319
320    tmp  = (double)61.0 + (double)90.0 * tan(nphi) * tan(nphi);
321    tmp += (double)45.0 * pow(tan(nphi),(double)4.0);
322    IX   = tan(nphi) / ((double)720.0 * rho * pow(nu,(double)5.0)) * tmp;
323
324    X    = (double)1.0 / (cos(nphi) * nu);
325
326    tmp  = (nu / rho) + (double)2.0 * tan(nphi) * tan(nphi);
327    XI   = ((double)1.0 / (cos(nphi) * (double)6.0 * nu*nu*nu)) * tmp;
328
329    tmp  = (double)5.0 + (double)28.0 * tan(nphi)*tan(nphi);
330    tmp += (double)24.0 * pow(tan(nphi),(double)4.0);
331    XII  = ((double)1.0 / ((double)120.0 * pow(nu,(double)5.0) * cos(nphi)))
332           * tmp;
333
334    tmp  = (double)61.0 + (double)662.0 * tan(nphi) * tan(nphi);
335    tmp += (double)1320.0 * pow(tan(nphi),(double)4.0);
336    tmp += (double)720.0  * pow(tan(nphi),(double)6.0);
337    XIIA = ((double)1.0 / (cos(nphi) * (double)5040.0 * pow(nu,(double)7.0)))
338           * tmp;
339
340    *phi = nphi - VII*pow((E-E0),(double)2.0) + VIII*pow((E-E0),(double)4.0) -
341           IX*pow((E-E0),(double)6.0);
342   
343    *lambda = lambda0 + X*(E-E0) - XI*pow((E-E0),(double)3.0) +
344              XII*pow((E-E0),(double)5.0) - XIIA*pow((E-E0),(double)7.0);
345
346        // Replaced these rad/deg conversions (originally done with a function)
347        // with simple multiplication
348    *phi   *= (180.0/M_PI); 
349    *lambda *= (180.0/M_PI); 
350
351    return;
352}
353
354}
Note: See TracBrowser for help on using the repository browser.