source: subversion/sites/other/freemap-npe/latlong.php @ 5557

Last change on this file since 5557 was 1540, checked in by nick, 13 years ago

Freemap-NPE: initial version

File size: 13.6 KB
Line 
1<?php
2
3
4# Constants for the UK gridref system
5define ('N0',  -100000);
6define('E0', 400000);
7define( 'F0', 0.9996012717);
8define ('PHI0',    49.0);
9define ('LAMBDA0', -2.0);
10define('A',        6377563.396);
11define('B',        6356256.910);
12
13/*
14print_r ( wgs84_ll_to_gr(array ("lat"=>51.05,"long"=>-0.72)));
15echo "<br/>";
16print_r ( gr_to_wgs84_ll(array ("e"=>489600,"n"=>128500)));
17*/
18
19// NOTICE OF AUTHORSHIP: These are PHP versions of C functions from
20// the LGPL JEEPS library, of which author and licence info follows....
21
22/********************************************************************
23** @source JEEPS arithmetic/conversion functions
24**
25** @author Copyright (C) 1999 Alan Bleasby
26** @version 1.0
27** @modified Dec 28 1999 Alan Bleasby. First version
28** @@
29**
30** This library is free software; you can redistribute it and/or
31** modify it under the terms of the GNU Library General Public
32** License as published by the Free Software Foundation; either
33** version 2 of the License, or (at your option) any later version.
34**
35** This library is distributed in the hope that it will be useful,
36** but WITHOUT ANY WARRANTY; without even the implied warranty of
37** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
38** Library General Public License for more details.
39**
40** You should have received a copy of the GNU Library General Public
41** License along with this library; if not, write to the
42** Free Software Foundation, Inc., 59 Temple Place - Suite 330,
43** Boston, MA  02111-1307, USA.
44********************************************************************/
45
46
47// ll_to_gr()
48// Slightly modified version of:
49/* @func GPS_Math_Airy1830LatLonToNGEN **************************************
50**
51** Convert Airy 1830 datum latitude and longitude to UK Ordnance Survey
52** National Grid Eastings and Northings
53**
54** @param [r] phi [double] WGS84 latitude     (deg)
55** @param [r] lambda [double] WGS84 longitude (deg)
56** @param [w] E [double *] NG easting (metres)
57** @param [w] N [double *] NG northing (metres)
58**
59** @return [void]
60*/
61
62// Modifications:
63// Altered to take advantage of PHP associative arrays.
64
65function ll_to_gr ( $pos )
66{
67    $gridref = GPS_Math_LatLon_To_EN($pos['lat'],$pos['long'],N0,E0,PHI0,
68                                                                                LAMBDA0,F0, A,B);
69
70    return $gridref;
71}
72
73function wgs84_ll_to_gr ( $pos )
74{
75        $outlat = 0;
76        $outlon = 0;
77        $ht = 0;
78        GPS_Math_WGS84_To_Known_Datum_M($pos['lat'],$pos['long'],30,
79                                        $outlat,$outlon,$ht);
80        return ll_to_gr(array("lat"=>$outlat,"long"=>$outlon));
81}
82
83/* @func  GPS_Math_LatLon_To_EN **********************************
84**
85** Convert latitude and longitude to eastings and northings
86** Standard Gauss-Kruger Transverse Mercator
87**
88** @param [w] E [double *] easting (metres)
89** @param [w] N [double *] northing (metres)
90** @param [r] phi [double] latitude (deg)
91** @param [r] lambda [double] longitude (deg)
92** @param [r] N0 [double] true northing origin (metres)
93** @param [r] E0 [double] true easting  origin (metres)
94** @param [r] phi0 [double] true latitude origin (deg)
95** @param [r] lambda0 [double] true longitude origin (deg)
96** @param [r] F0 [double] scale factor on central meridian
97** @param [r] a [double] semi-major axis (metres)
98** @param [r] b [double] semi-minor axis (metres)
99**
100** @return [void]
101************************************************************************/
102function GPS_Math_LatLon_To_EN($phi, $lambda, $N0, $E0, $phi0,
103                                                                $lambda0, $F0, $a, $b)
104{
105        // Replaced these deg/rad conversions (originally done with a function)
106        // with simple multiplication
107    $phi0    *= (M_PI/180.0); 
108    $lambda0 *= (M_PI/180.0); 
109    $phi     *= (M_PI/180.0); 
110    $lambda  *= (M_PI/180.0); 
111   
112    $esq = (($a*$a)-($b*$b)) / ($a*$a);
113    $n   = ($a-$b) / ($a+$b);
114   
115    $tmp  = 1.0 - ($esq * sin($phi) * sin($phi));
116    $nu   = $a * $F0 * pow($tmp,-0.5);
117    $rho  = $a * $F0 * (1.0 - $esq) * pow($tmp,-1.5);
118    $etasq = ($nu / $rho) - 1.0;
119
120    $fdf   = 5.0 / 4.0;
121    $tmp   = 1.0 + $n + ($fdf * $n * $n) + ($fdf * $n * $n * $n);
122    $tmp  *= ($phi - $phi0);
123    $tmp2  = 3.0*$n + 3.0*$n*$n + (21.0/8.0)*$n*$n*$n;
124    $tmp2 *= (sin($phi-$phi0) * cos($phi+$phi0));
125    $tmp  -= $tmp2;
126
127    $fde   = (15.0 / 8.0);
128    $tmp2  = (($fde*$n*$n) + ($fde*$n*$n*$n)) * sin(2.0 * ($phi-$phi0));
129    $tmp2 *= cos(2.0 * ($phi+$phi0));
130    $tmp  += $tmp2;
131   
132    $tmp2  = (35.0/24.0) * $n * $n * $n;
133    $tmp2 *= sin(3.0 * ($phi-$phi0));
134    $tmp2 *= cos(3.0 * ($phi+$phi0));
135    $tmp  -= $tmp2;
136
137    $M     = $b * $F0 * $tmp;
138    $I     = $M + $N0;
139    $II    = ($nu / 2.0) * sin($phi) * cos($phi);
140    $III   = ($nu / 24.0) * sin($phi) * cos($phi) * cos($phi) * cos($phi);
141    $III  *= (5.0 - (tan($phi) * tan($phi)) + (9.0 * $etasq));
142    $IIIA  = ($nu / 720.0) * sin($phi) * pow(cos($phi),5.0);
143    $IIIA *= (61.0 - (58.0*tan($phi)*tan($phi)) +
144             pow(tan($phi),4.0));
145    $IV    = $nu * cos($phi);
146
147    $tmp   = pow(cos($phi),3.0);
148    $tmp  *= (($nu/$rho) - tan($phi) * tan($phi));
149    $V     = ($nu/6.0) * $tmp;
150
151    $tmp   = 5.0 - (18.0 * tan($phi) * tan($phi));
152    $tmp  += tan($phi)*tan($phi)*tan($phi)*tan($phi) + (14.0 * $etasq);
153    $tmp  -= (58.0 * tan($phi) * tan($phi) * $etasq);
154    $tmp2  = cos($phi)*cos($phi)*cos($phi)*cos($phi)*cos($phi) * $tmp;
155    $VI    = ($nu / 120.0) * $tmp2;
156   
157    $gridref['n'] = $I + $II*($lambda-$lambda0)*($lambda-$lambda0) +
158             $III*pow(($lambda-$lambda0),4.0) +
159             $IIIA*pow(($lambda-$lambda0),6.0);
160
161    $gridref['e']=$E0 + $IV*($lambda-$lambda0)+$V*pow(($lambda-$lambda0),3.0) +
162         $VI * pow(($lambda-$lambda0),5.0);
163
164    return $gridref;
165}
166
167// gr_to_ll()
168// Slightly modified version of:
169/* @func GPS_Math_NGENToAiry1830LatLon **************************************
170**
171** Convert  to UK Ordnance Survey National Grid Eastings and Northings to
172** Airy 1830 datum latitude and longitude
173**
174** @param [r] E [double] NG easting (metres)
175** @param [r] N [double] NG northing (metres)
176** @param [w] phi [double *] Airy latitude     (deg)
177** @param [w] lambda [double *] Airy longitude (deg)
178**
179** @return [void]
180************************************************************************/
181
182// Modifications:
183// Altered to take advantage of PHP associative arrays.
184function gr_to_ll($gridref)
185{
186    $point = GPS_Math_EN_To_LatLon($gridref['e'],$gridref['n'],
187                                                                        N0,E0,PHI0,LAMBDA0,F0,A,B);
188   
189    return $point;
190}
191
192function gr_to_wgs84_ll($gridref)
193{
194        $outlat = 0;
195        $outlon = 0;
196        $ht = 0;
197       
198        $latlon = gr_to_ll($gridref);
199        GPS_Math_Known_Datum_To_WGS84_M ($latlon['lat'],$latlon['long'],0,
200                                                        $outlat,$outlon,$ht);
201        return array ("lat"=>$outlat,"long"=>$outlon);
202}
203
204/* @func  GPS_Math_EN_To_LatLon **************************************
205**
206** Convert Eastings and Northings to latitude and longitude
207**
208** @param [w] E [double] NG easting (metres)
209** @param [w] N [double] NG northing (metres)
210** @param [r] phi [double *] Airy latitude     (deg)
211** @param [r] lambda [double *] Airy longitude (deg)
212** @param [r] N0 [double] true northing origin (metres)
213** @param [r] E0 [double] true easting  origin (metres)
214** @param [r] phi0 [double] true latitude origin (deg)
215** @param [r] lambda0 [double] true longitude origin (deg)
216** @param [r] F0 [double] scale factor on central meridian
217** @param [r] a [double] semi-major axis (metres)
218** @param [r] b [double] semi-minor axis (metres)
219**
220** @return [void]
221************************************************************************/
222function GPS_Math_EN_To_LatLon($E,  $N, 
223                           $N0, $E0,
224                           $phi0, $lambda0,
225                           $F0, $a, $b)
226{
227    $nphi=0.0;
228   
229        // Replaced these deg/rad conversions (originally done with a function)
230        // with simple multiplication
231    $phi0    *= (M_PI/180.0); 
232    $lambda0 *= (M_PI/180.0); 
233
234    $n     = ($a-$b) / ($a+$b);
235    $fdf   = 5.0 / 4.0;
236    $fde   = (15.0 / 8.0);
237
238    $esq = (($a*$a)-($b*$b)) / ($a*$a);
239
240
241    $phix = (($N-$N0)/($a*$F0)) + $phi0;
242   
243    $tmp  = 1.0 - ($esq * sin($phix) * sin($phix));
244    $nu   = $a * $F0 * pow($tmp,-0.5);
245    $rho  = $a * $F0 * (1.0 - $esq) * pow($tmp,-1.5);
246    $etasq = ($nu / $rho) - 1.0;
247
248    $M = -1e20;
249
250       
251        for($c=1; $c<=10; $c++)
252    //while($N-$N0-$M > 0.000001)
253    {
254        $nphi = $phix;
255       
256        $tmp   = 1.0 + $n + ($fdf * $n * $n) + ($fdf * $n * $n * $n);
257        $tmp  *= ($nphi - $phi0);
258        $tmp2  = 3.0*$n + 3.0*$n*$n +
259                (21.0/8.0)*$n*$n*$n;
260        $tmp2 *= (sin($nphi-$phi0) * cos($nphi+$phi0));
261        $tmp  -= $tmp2;
262
263
264        $tmp2  = (($fde*$n*$n)+($fde*$n*$n*$n))*sin(2.0*($nphi-$phi0));
265        $tmp2 *= cos(2.0 * ($nphi+$phi0));
266        $tmp  += $tmp2;
267   
268        $tmp2  = (35.0/24.0) * $n * $n * $n;
269        $tmp2 *= sin(3.0 * ($nphi-$phi0));
270        $tmp2 *= cos(3.0 * ($nphi+$phi0));
271        $tmp  -= $tmp2;
272
273        $M     = $b * $F0 * $tmp;
274
275        if($N-$N0-$M > 0.000001)
276            $phix = (($N-$N0-$M)/($a*$F0)) + $nphi;
277    }
278   
279
280    $VII  = tan($nphi) / (2.0 * $rho * $nu);
281
282    $tmp  = 5.0 + 3.0 * tan($nphi) * tan($nphi) + $etasq;
283    $tmp -= 9.0 * tan($nphi) * tan($nphi) * $etasq;
284    $VIII = (tan($nphi)*$tmp) / (24.0 * $rho * $nu*$nu*$nu);
285
286    $tmp  = 61.0 + 90.0 * tan($nphi) * tan($nphi);
287    $tmp += 45.0 * pow(tan($nphi),4.0);
288    $IX   = tan($nphi) / (720.0 * $rho * pow($nu,5.0)) * $tmp;
289
290    $X    = 1.0 / (cos($nphi) * $nu);
291
292    $tmp  = ($nu / $rho) + 2.0 * tan($nphi) * tan($nphi);
293    $XI   = (1.0 / (cos($nphi) * 6.0 * $nu*$nu*$nu)) * $tmp;
294
295    $tmp  = 5.0 + 28.0 * tan($nphi)*tan($nphi);
296    $tmp += 24.0 * pow(tan($nphi),4.0);
297    $XII  = (1.0 / (120.0 * pow($nu,5.0) * cos($nphi)))
298           * $tmp;
299
300    $tmp  = 61.0 + 662.0 * tan($nphi) * tan($nphi);
301    $tmp += 1320.0 * pow(tan($nphi),4.0);
302    $tmp += 720.0  * pow(tan($nphi),6.0);
303    $XIIA = (1.0 / (cos($nphi) * 5040.0 * pow($nu,7.0)))
304           * $tmp;
305
306    $point['lat'] = $nphi - $VII*pow(($E-$E0),2.0) + $VIII*pow(($E-$E0),4.0) -
307           $IX*pow(($E-$E0),6.0);
308   
309    $point['long'] = $lambda0 + $X*($E-$E0) - $XI*pow(($E-$E0),3.0) +
310              $XII*pow(($E-$E0),5.0) - $XIIA*pow(($E-$E0),7.0);
311
312        // Replaced these rad/deg conversions (originally done with a function)
313        // with simple multiplication
314    $point['lat']   *= (180.0/M_PI); 
315    $point['long'] *= (180.0/M_PI); 
316
317    return $point;
318}
319// jeeps funcs (c) Alan Bleasby Licence LGPL
320
321/* @func GPS_Math_Known_Datum_To_WGS84_M **********************************
322**
323** Transform datum to WGS84 using Molodensky
324**
325** @param [r] Sphi [double] source latitude (deg)
326** @param [r] Slam [double] source longitude (deg)
327** @param [r] SH   [double] source height  (metres)
328** @param [w] Dphi [$*] dest latitude (deg)
329** @param [w] Dlam [$*] dest longitude (deg)
330** @param [w] DH   [$*] dest height  (metres)
331** @param [r] n    [int32] datum number from GPS_Datum structure
332**
333** @return [void]
334************************************************************************/
335function GPS_Math_Known_Datum_To_WGS84_M($Sphi, $Slam, $SH,
336                                     &$Dphi, &$Dlam, &$DH)
337{
338    $Da  =  6378137.0;
339    $Dif =  298.257223563;
340   
341    $Sa   = 6378206.400; 
342    $Sif  = 294.9786982; 
343    $x    = -8; 
344    $y    = 160;
345    $z    = 176;
346
347    GPS_Math_Molodensky($Sphi,$Slam,$SH,$Sa,$Sif,$Dphi,$Dlam,$DH,$Da,$Dif,
348                                                        $x,$y,$z);
349}
350/* @func GPS_Math_WGS84_To_Known_Datum_M ********************************
351**
352** Transform WGS84 to other datum using Molodensky
353**
354** @param [r] Sphi [double] source latitude (deg)
355** @param [r] Slam [double] source longitude (deg)
356** @param [r] SH   [double] source height  (metres)
357** @param [w] Dphi [$*] dest latitude (deg)
358** @param [w] Dlam [$*] dest longitude (deg)
359** @param [w] DH   [$*] dest height  (metres)
360** @param [r] n    [int32] datum number from GPS_Datum structure
361**
362** @return [void]
363************************************************************************/
364function GPS_Math_WGS84_To_Known_Datum_M($Sphi, $Slam, $SH,
365                                     &$Dphi, &$Dlam, &$DH)
366{
367    $Sa  =  6378137.0;
368    $Sif =  298.257223563;
369   
370    $Da   = 6377563.396; 
371    $Dif  = 299.3249646;
372    $x    = -375;
373    $y    = 111;
374    $z    = -431;
375
376    GPS_Math_Molodensky($Sphi,$Slam,$SH,$Sa,$Sif,$Dphi,$Dlam,$DH,$Da,$Dif,
377                                                        $x,$y,$z);
378}
379/* @func GPS_Math_Molodensky *******************************************
380**
381** Transform one datum to another
382**
383** @param [r] Sphi [double] source latitude (deg)
384** @param [r] Slam [double] source longitude (deg)
385** @param [r] SH   [double] source height  (metres)
386** @param [r] Sa   [double] source semi-major axis (metres)
387** @param [r] Sif  [double] source inverse flattening
388** @param [w] Dphi [$*] dest latitude (deg)
389** @param [w] Dlam [$*] dest longitude (deg)
390** @param [w] DH   [$*] dest height  (metres)
391** @param [r] Da   [double]   dest semi-major axis (metres)
392** @param [r] Dif  [double]   dest inverse flattening
393** @param [r] dx  [double]   dx
394** @param [r] dy  [double]   dy
395** @param [r] dz  [double]   dz
396**
397** @return [void]
398************************************************************************/
399function GPS_Math_Molodensky($Sphi, $Slam, $SH, $Sa,
400                         $Sif, &$Dphi, &$Dlam,
401                         &$DH, $Da, $Dif, $dx,
402                         $dy, $dz)
403{
404    $Sf = 1.0 / $Sif;
405    $Df = 1.0 / $Dif;
406   
407    $esq = 2.0*$Sf - pow($Sf,2.0);
408    $bda = 1.0 - $Sf;
409    $Sphi = deg2rad($Sphi);
410    $Slam = deg2rad($Slam);
411       
412    $da = $Da - $Sa;
413    $df = $Df - $Sf;
414
415    $phis = sin($Sphi);
416    $phic = cos($Sphi);
417    $lams = sin($Slam);
418    $lamc = cos($Slam);
419   
420    $N = $Sa /  sqrt(1.0 - $esq*pow($phis,2.0));
421   
422    $tmp = (1.0-$esq) /pow((1.0-$esq*pow($phis,2.0)),1.5);
423    $M   = $Sa * $tmp;
424
425    $tmp  = $df * (($M/$bda)+$N*$bda) * $phis * $phic;
426    $tmp2 = $da * $N * $esq * $phis * $phic / $Sa;
427    $tmp2 += ((-$dx*$phis*$lamc-$dy*$phis*$lams) + $dz*$phic);
428    $dphi = ($tmp2 + $tmp) / ($M + $SH);
429   
430    $dlambda = (-$dx*$lams+$dy*$lamc) / (($N+$SH)*$phic);
431
432    $dheight = $dx*$phic*$lamc + $dy*$phic*$lams + $dz*$phis - $da*($Sa/$N) +
433        $df*$bda*$N*$phis*$phis;
434   
435    $Dphi = $Sphi + $dphi;
436    $Dlam = $Slam + $dlambda;
437    $DH   = $SH   + $dheight;
438   
439    $Dphi = rad2deg($Dphi);
440    $Dlam = rad2deg($Dlam);
441}
442?>
Note: See TracBrowser for help on using the repository browser.