Automatically exported from code.google.com/p/planningalerts
Non puoi selezionare più di 25 argomenti Gli argomenti devono iniziare con una lettera o un numero, possono includere trattini ('-') e possono essere lunghi fino a 35 caratteri.
 
 
 
 
 
 

776 righe
24 KiB

  1. <?php
  2. //--------------------------------------------------------------------------
  3. // PHPcoord
  4. // phpcoord.php
  5. //
  6. // (c) 2005 Jonathan Stott
  7. //
  8. // Created on 11-Aug-2005
  9. //
  10. // 2.2 - 11 Feb 2006
  11. // - Used different algorithm for calculating distance between latitudes
  12. // and longitudes - fixes a number of problems with distance calculations
  13. // 2.1 - 22 Dec 2005
  14. // - Added getOSRefFromSixFigureReference function
  15. // 2.0 - 21 Dec 2005
  16. // - Completely different object design - conversion functions now through
  17. // objects rather than static functions
  18. // - Updated comments and documentation
  19. // 1.1 - 11 Sep 2005
  20. // - Added OSGB36/WGS84 data conversions
  21. // 1.0 - 11 Aug 2005
  22. // - Initial version
  23. //--------------------------------------------------------------------------
  24. // ================================================================== LatLng
  25. class LatLng {
  26. var $lat;
  27. var $lng;
  28. /**
  29. * Create a new LatLng object from the given latitude and longitude
  30. *
  31. * @param lat latitude
  32. * @param lng longitude
  33. */
  34. function LatLng($lat, $lng) {
  35. $this->lat = $lat;
  36. $this->lng = $lng;
  37. }
  38. /**
  39. * Return a string representation of this LatLng object
  40. *
  41. * @return a string representation of this LatLng object
  42. */
  43. function toString() {
  44. return "(" . $this->lat . ", " . $this->lng . ")";
  45. }
  46. /**
  47. * Calculate the surface distance between this LatLng object and the one
  48. * passed in as a parameter.
  49. *
  50. * @param to a LatLng object to measure the surface distance to
  51. * @return the surface distance
  52. */
  53. function distance($to) {
  54. $er = 6366.707;
  55. $latFrom = deg2rad($this->lat);
  56. $latTo = deg2rad($to->lat);
  57. $lngFrom = deg2rad($this->lng);
  58. $lngTo = deg2rad($to->lng);
  59. $x1 = $er * cos($lngFrom) * sin($latFrom);
  60. $y1 = $er * sin($lngFrom) * sin($latFrom);
  61. $z1 = $er * cos($latFrom);
  62. $x2 = $er * cos($lngTo) * sin($latTo);
  63. $y2 = $er * sin($lngTo) * sin($latTo);
  64. $z2 = $er * cos($latTo);
  65. $d = acos(sin($latFrom)*sin($latTo) + cos($latFrom)*cos($latTo)*cos($lngTo-$lngFrom)) * $er;
  66. return $d;
  67. }
  68. /**
  69. * Convert this LatLng object from OSGB36 datum to WGS84 datum.
  70. */
  71. function OSGB36ToWGS84() {
  72. $airy1830 = new RefEll(6377563.396, 6356256.909);
  73. $a = $airy1830->maj;
  74. $b = $airy1830->min;
  75. $eSquared = $airy1830->ecc;
  76. $phi = deg2rad($this->lat);
  77. $lambda = deg2rad($this->lng);
  78. $v = $a / (sqrt(1 - $eSquared * sinSquared($phi)));
  79. $H = 0; // height
  80. $x = ($v + $H) * cos($phi) * cos($lambda);
  81. $y = ($v + $H) * cos($phi) * sin($lambda);
  82. $z = ((1 - $eSquared) * $v + $H) * sin($phi);
  83. $tx = 446.448;
  84. $ty = -124.157;
  85. $tz = 542.060;
  86. $s = -0.0000204894;
  87. $rx = deg2rad( 0.00004172222);
  88. $ry = deg2rad( 0.00006861111);
  89. $rz = deg2rad( 0.00023391666);
  90. $xB = $tx + ($x * (1 + $s)) + (-$rx * $y) + ($ry * $z);
  91. $yB = $ty + ($rz * $x) + ($y * (1 + $s)) + (-$rx * $z);
  92. $zB = $tz + (-$ry * $x) + ($rx * $y) + ($z * (1 + $s));
  93. $wgs84 = new RefEll(6378137.000, 6356752.3141);
  94. $a = $wgs84->maj;
  95. $b = $wgs84->min;
  96. $eSquared = $wgs84->ecc;
  97. $lambdaB = rad2deg(atan($yB / $xB));
  98. $p = sqrt(($xB * $xB) + ($yB * $yB));
  99. $phiN = atan($zB / ($p * (1 - $eSquared)));
  100. for ($i = 1; $i < 10; $i++) {
  101. $v = $a / (sqrt(1 - $eSquared * sinSquared($phiN)));
  102. $phiN1 = atan(($zB + ($eSquared * $v * sin($phiN))) / $p);
  103. $phiN = $phiN1;
  104. }
  105. $phiB = rad2deg($phiN);
  106. $this->lat = $phiB;
  107. $this->lng = $lambdaB;
  108. }
  109. /**
  110. * Convert this LatLng object from WGS84 datum to OSGB36 datum.
  111. */
  112. function WGS84ToOSGB36() {
  113. $wgs84 = new RefEll(6378137.000, 6356752.3141);
  114. $a = $wgs84->maj;
  115. $b = $wgs84->min;
  116. $eSquared = $wgs84->ecc;
  117. $phi = deg2rad($this->lat);
  118. $lambda = deg2rad($this->lng);
  119. $v = $a / (sqrt(1 - $eSquared * sinSquared($phi)));
  120. $H = 0; // height
  121. $x = ($v + $H) * cos($phi) * cos($lambda);
  122. $y = ($v + $H) * cos($phi) * sin($lambda);
  123. $z = ((1 - $eSquared) * $v + $H) * sin($phi);
  124. $tx = -446.448;
  125. $ty = 124.157;
  126. $tz = -542.060;
  127. $s = 0.0000204894;
  128. $rx = deg2rad(-0.00004172222);
  129. $ry = deg2rad(-0.00006861111);
  130. $rz = deg2rad(-0.00023391666);
  131. $xB = $tx + ($x * (1 + $s)) + (-$rx * $y) + ($ry * $z);
  132. $yB = $ty + ($rz * $x) + ($y * (1 + $s)) + (-$rx * $z);
  133. $zB = $tz + (-$ry * $x) + ($rx * $y) + ($z * (1 + $s));
  134. $airy1830 = new RefEll(6377563.396, 6356256.909);
  135. $a = $airy1830->maj;
  136. $b = $airy1830->min;
  137. $eSquared = $airy1830->ecc;
  138. $lambdaB = rad2deg(atan($yB / $xB));
  139. $p = sqrt(($xB * $xB) + ($yB * $yB));
  140. $phiN = atan($zB / ($p * (1 - $eSquared)));
  141. for ($i = 1; $i < 10; $i++) {
  142. $v = $a / (sqrt(1 - $eSquared * sinSquared($phiN)));
  143. $phiN1 = atan(($zB + ($eSquared * $v * sin($phiN))) / $p);
  144. $phiN = $phiN1;
  145. }
  146. $phiB = rad2deg($phiN);
  147. $this->lat = $phiB;
  148. $this->lng = $lambdaB;
  149. }
  150. /**
  151. * Convert this LatLng object into an OSGB grid reference. Note that this
  152. * function does not take into account the bounds of the OSGB grid -
  153. * beyond the bounds of the OSGB grid, the resulting OSRef object has no
  154. * meaning
  155. *
  156. * @return the converted OSGB grid reference
  157. */
  158. function toOSRef() {
  159. $airy1830 = new RefEll(6377563.396, 6356256.909);
  160. $OSGB_F0 = 0.9996012717;
  161. $N0 = -100000.0;
  162. $E0 = 400000.0;
  163. $phi0 = deg2rad(49.0);
  164. $lambda0 = deg2rad(-2.0);
  165. $a = $airy1830->maj;
  166. $b = $airy1830->min;
  167. $eSquared = $airy1830->ecc;
  168. $phi = deg2rad($this->lat);
  169. $lambda = deg2rad($this->lng);
  170. $E = 0.0;
  171. $N = 0.0;
  172. $n = ($a - $b) / ($a + $b);
  173. $v = $a * $OSGB_F0 * pow(1.0 - $eSquared * sinSquared($phi), -0.5);
  174. $rho =
  175. $a * $OSGB_F0 * (1.0 - $eSquared) * pow(1.0 - $eSquared * sinSquared($phi), -1.5);
  176. $etaSquared = ($v / $rho) - 1.0;
  177. $M =
  178. ($b * $OSGB_F0)
  179. * (((1 + $n + ((5.0 / 4.0) * $n * $n) + ((5.0 / 4.0) * $n * $n * $n))
  180. * ($phi - $phi0))
  181. - (((3 * $n) + (3 * $n * $n) + ((21.0 / 8.0) * $n * $n * $n))
  182. * sin($phi - $phi0)
  183. * cos($phi + $phi0))
  184. + ((((15.0 / 8.0) * $n * $n) + ((15.0 / 8.0) * $n * $n * $n))
  185. * sin(2.0 * ($phi - $phi0))
  186. * cos(2.0 * ($phi + $phi0)))
  187. - (((35.0 / 24.0) * $n * $n * $n)
  188. * sin(3.0 * ($phi - $phi0))
  189. * cos(3.0 * ($phi + $phi0))));
  190. $I = $M + $N0;
  191. $II = ($v / 2.0) * sin($phi) * cos($phi);
  192. $III =
  193. ($v / 24.0)
  194. * sin($phi)
  195. * pow(cos($phi), 3.0)
  196. * (5.0 - tanSquared($phi) + (9.0 * $etaSquared));
  197. $IIIA =
  198. ($v / 720.0)
  199. * sin($phi)
  200. * pow(cos($phi), 5.0)
  201. * (61.0 - (58.0 * tanSquared($phi)) + pow(tan($phi), 4.0));
  202. $IV = $v * cos($phi);
  203. $V = ($v / 6.0) * pow(cos($phi), 3.0) * (($v / $rho) - tanSquared($phi));
  204. $VI =
  205. ($v / 120.0)
  206. * pow(cos($phi), 5.0)
  207. * (5.0
  208. - (18.0 * tanSquared($phi))
  209. + (pow(tan($phi), 4.0))
  210. + (14 * $etaSquared)
  211. - (58 * tanSquared($phi) * $etaSquared));
  212. $N =
  213. $I
  214. + ($II * pow($lambda - $lambda0, 2.0))
  215. + ($III * pow($lambda - $lambda0, 4.0))
  216. + ($IIIA * pow($lambda - $lambda0, 6.0));
  217. $E =
  218. $E0
  219. + ($IV * ($lambda - $lambda0))
  220. + ($V * pow($lambda - $lambda0, 3.0))
  221. + ($VI * pow($lambda - $lambda0, 5.0));
  222. return new OSRef($E, $N);
  223. }
  224. /**
  225. * Convert a latitude and longitude to an UTM reference
  226. *
  227. * @return the converted UTM reference
  228. */
  229. function toUTMRef() {
  230. $wgs84 = new RefEll(6378137, 6356752.314);
  231. $UTM_F0 = 0.9996;
  232. $a = $wgs84->maj;
  233. $eSquared = $wgs84->ecc;
  234. $longitude = $this->lng;
  235. $latitude = $this->lat;
  236. $latitudeRad = $latitude * (pi() / 180.0);
  237. $longitudeRad = $longitude * (pi() / 180.0);
  238. $longitudeZone = (int) (($longitude + 180.0) / 6.0) + 1;
  239. // Special zone for Norway
  240. if ($latitude >= 56.0
  241. && $latitude < 64.0
  242. && $longitude >= 3.0
  243. && $longitude < 12.0) {
  244. $longitudeZone = 32;
  245. }
  246. // Special zones for Svalbard
  247. if ($latitude >= 72.0 && $latitude < 84.0) {
  248. if ($longitude >= 0.0 && $longitude < 9.0) {
  249. $longitudeZone = 31;
  250. } else if ($longitude >= 9.0 && $longitude < 21.0) {
  251. $longitudeZone = 33;
  252. } else if ($longitude >= 21.0 && $longitude < 33.0) {
  253. $longitudeZone = 35;
  254. } else if ($longitude >= 33.0 && $longitude < 42.0) {
  255. $longitudeZone = 37;
  256. }
  257. }
  258. $longitudeOrigin = ($longitudeZone - 1) * 6 - 180 + 3;
  259. $longitudeOriginRad = $longitudeOrigin * (pi() / 180.0);
  260. $UTMZone = getUTMLatitudeZoneLetter($latitude);
  261. $ePrimeSquared = ($eSquared) / (1 - $eSquared);
  262. $n = $a / sqrt(1 - $eSquared * sin($latitudeRad) * sin($latitudeRad));
  263. $t = tan($latitudeRad) * tan($latitudeRad);
  264. $c = $ePrimeSquared * cos($latitudeRad) * cos($latitudeRad);
  265. $A = cos($latitudeRad) * ($longitudeRad - $longitudeOriginRad);
  266. $M =
  267. $a
  268. * ((1
  269. - $eSquared / 4
  270. - 3 * $eSquared * $eSquared / 64
  271. - 5 * $eSquared * $eSquared * $eSquared / 256)
  272. * $latitudeRad
  273. - (3 * $eSquared / 8
  274. + 3 * $eSquared * $eSquared / 32
  275. + 45 * $eSquared * $eSquared * $eSquared / 1024)
  276. * sin(2 * $latitudeRad)
  277. + (15 * $eSquared * $eSquared / 256
  278. + 45 * $eSquared * $eSquared * $eSquared / 1024)
  279. * sin(4 * $latitudeRad)
  280. - (35 * $eSquared * $eSquared * $eSquared / 3072)
  281. * sin(6 * $latitudeRad));
  282. $UTMEasting =
  283. (double) ($UTM_F0
  284. * $n
  285. * ($A
  286. + (1 - $t + $c) * pow($A, 3.0) / 6
  287. + (5 - 18 * $t + $t * $t + 72 * $c - 58 * $ePrimeSquared)
  288. * pow($A, 5.0)
  289. / 120)
  290. + 500000.0);
  291. $UTMNorthing =
  292. (double) ($UTM_F0
  293. * ($M
  294. + $n
  295. * tan($latitudeRad)
  296. * ($A * $A / 2
  297. + (5 - $t + (9 * $c) + (4 * $c * $c)) * pow($A, 4.0) / 24
  298. + (61 - (58 * $t) + ($t * $t) + (600 * $c) - (330 * $ePrimeSquared))
  299. * pow($A, 6.0)
  300. / 720)));
  301. // Adjust for the southern hemisphere
  302. if ($latitude < 0) {
  303. $UTMNorthing += 10000000.0;
  304. }
  305. return new UTMRef($UTMEasting, $UTMNorthing, $UTMZone, $longitudeZone);
  306. }
  307. }
  308. // =================================================================== OSRef
  309. // References given with OSRef are accurate to 1m.
  310. class OSRef {
  311. var $easting;
  312. var $northing;
  313. /**
  314. * Create a new OSRef object representing an OSGB grid reference. Note
  315. * that the parameters for this constructor require eastings and
  316. * northings with 1m accuracy and need to be absolute with respect to
  317. * the whole of the British Grid. For example, to create an OSRef
  318. * object from the six-figure grid reference TG514131, the easting would
  319. * be 651400 and the northing would be 313100.
  320. *
  321. * Grid references with accuracy greater than 1m can be represented
  322. * using floating point values for the easting and northing. For example,
  323. * a value representing an easting or northing accurate to 1mm would be
  324. * given as 651400.0001.
  325. *
  326. * @param easting the easting of the reference (with 1m accuracy)
  327. * @param northing the northing of the reference (with 1m accuracy)
  328. */
  329. function OSRef($easting, $northing) {
  330. $this->easting = $easting;
  331. $this->northing = $northing;
  332. }
  333. /**
  334. * Convert this grid reference into a string showing the exact values
  335. * of the easting and northing.
  336. *
  337. * @return
  338. */
  339. function toString() {
  340. return "(" . $this->easting . ", " . $this->northing . ")";
  341. }
  342. /**
  343. * Convert this grid reference into a string using a standard six-figure
  344. * grid reference including the two-character designation for the 100km
  345. * square. e.g. TG514131.
  346. *
  347. * @return
  348. */
  349. function toSixFigureString() {
  350. $hundredkmE = floor($this->easting / 100000);
  351. $hundredkmN = floor($this->northing / 100000);
  352. $firstLetter = "";
  353. if ($hundredkmN < 5) {
  354. if ($hundredkmE < 5) {
  355. $firstLetter = "S";
  356. } else {
  357. $firstLetter = "T";
  358. }
  359. } else if ($hundredkmN < 10) {
  360. if ($hundredkmE < 5) {
  361. $firstLetter = "N";
  362. } else {
  363. $firstLetter = "O";
  364. }
  365. } else {
  366. $firstLetter = "H";
  367. }
  368. $secondLetter = "";
  369. $index = 65 + ((4 - ($hundredkmN % 5)) * 5) + ($hundredkmE % 5);
  370. $ti = $index;
  371. if ($index >= 73) $index++;
  372. $secondLetter = chr($index);
  373. $e = floor(($this->easting - (100000 * $hundredkmE)) / 100);
  374. $n = floor(($this->northing - (100000 * $hundredkmN)) / 100);
  375. $es = $e;
  376. if ($e < 100) $es = "0$es";
  377. if ($e < 10) $es = "0$es";
  378. $ns = $n;
  379. if ($n < 100) $ns = "0$ns";
  380. if ($n < 10) $ns = "0$ns";
  381. return $firstLetter . $secondLetter . $es . $ns;
  382. }
  383. /**
  384. * Convert this grid reference into a latitude and longitude
  385. *
  386. * @return
  387. */
  388. function toLatLng() {
  389. $airy1830 = new RefEll(6377563.396, 6356256.909);
  390. $OSGB_F0 = 0.9996012717;
  391. $N0 = -100000.0;
  392. $E0 = 400000.0;
  393. $phi0 = deg2rad(49.0);
  394. $lambda0 = deg2rad(-2.0);
  395. $a = $airy1830->maj;
  396. $b = $airy1830->min;
  397. $eSquared = $airy1830->ecc;
  398. $phi = 0.0;
  399. $lambda = 0.0;
  400. $E = $this->easting;
  401. $N = $this->northing;
  402. $n = ($a - $b) / ($a + $b);
  403. $M = 0.0;
  404. $phiPrime = (($N - $N0) / ($a * $OSGB_F0)) + $phi0;
  405. do {
  406. $M =
  407. ($b * $OSGB_F0)
  408. * (((1 + $n + ((5.0 / 4.0) * $n * $n) + ((5.0 / 4.0) * $n * $n * $n))
  409. * ($phiPrime - $phi0))
  410. - (((3 * $n) + (3 * $n * $n) + ((21.0 / 8.0) * $n * $n * $n))
  411. * sin($phiPrime - $phi0)
  412. * cos($phiPrime + $phi0))
  413. + ((((15.0 / 8.0) * $n * $n) + ((15.0 / 8.0) * $n * $n * $n))
  414. * sin(2.0 * ($phiPrime - $phi0))
  415. * cos(2.0 * ($phiPrime + $phi0)))
  416. - (((35.0 / 24.0) * $n * $n * $n)
  417. * sin(3.0 * ($phiPrime - $phi0))
  418. * cos(3.0 * ($phiPrime + $phi0))));
  419. $phiPrime += ($N - $N0 - $M) / ($a * $OSGB_F0);
  420. } while (($N - $N0 - $M) >= 0.001);
  421. $v = $a * $OSGB_F0 * pow(1.0 - $eSquared * sinSquared($phiPrime), -0.5);
  422. $rho =
  423. $a
  424. * $OSGB_F0
  425. * (1.0 - $eSquared)
  426. * pow(1.0 - $eSquared * sinSquared($phiPrime), -1.5);
  427. $etaSquared = ($v / $rho) - 1.0;
  428. $VII = tan($phiPrime) / (2 * $rho * $v);
  429. $VIII =
  430. (tan($phiPrime) / (24.0 * $rho * pow($v, 3.0)))
  431. * (5.0
  432. + (3.0 * tanSquared($phiPrime))
  433. + $etaSquared
  434. - (9.0 * tanSquared($phiPrime) * $etaSquared));
  435. $IX =
  436. (tan($phiPrime) / (720.0 * $rho * pow($v, 5.0)))
  437. * (61.0
  438. + (90.0 * tanSquared($phiPrime))
  439. + (45.0 * tanSquared($phiPrime) * tanSquared($phiPrime)));
  440. $X = sec($phiPrime) / $v;
  441. $XI =
  442. (sec($phiPrime) / (6.0 * $v * $v * $v))
  443. * (($v / $rho) + (2 * tanSquared($phiPrime)));
  444. $XII =
  445. (sec($phiPrime) / (120.0 * pow($v, 5.0)))
  446. * (5.0
  447. + (28.0 * tanSquared($phiPrime))
  448. + (24.0 * tanSquared($phiPrime) * tanSquared($phiPrime)));
  449. $XIIA =
  450. (sec($phiPrime) / (5040.0 * pow($v, 7.0)))
  451. * (61.0
  452. + (662.0 * tanSquared($phiPrime))
  453. + (1320.0 * tanSquared($phiPrime) * tanSquared($phiPrime))
  454. + (720.0
  455. * tanSquared($phiPrime)
  456. * tanSquared($phiPrime)
  457. * tanSquared($phiPrime)));
  458. $phi =
  459. $phiPrime
  460. - ($VII * pow($E - $E0, 2.0))
  461. + ($VIII * pow($E - $E0, 4.0))
  462. - ($IX * pow($E - $E0, 6.0));
  463. $lambda =
  464. $lambda0
  465. + ($X * ($E - $E0))
  466. - ($XI * pow($E - $E0, 3.0))
  467. + ($XII * pow($E - $E0, 5.0))
  468. - ($XIIA * pow($E - $E0, 7.0));
  469. return new LatLng(rad2deg($phi), rad2deg($lambda));
  470. }
  471. }
  472. // ================================================================== UTMRef
  473. class UTMRef {
  474. var $easting;
  475. var $northing;
  476. var $latZone;
  477. var $lngZone;
  478. /**
  479. * Create a new object representing a UTM reference.
  480. *
  481. * @param easting
  482. * @param northing
  483. * @param latZone
  484. * @param lngZone
  485. */
  486. function UTMRef($easting, $northing, $latZone, $lngZone) {
  487. $this->easting = $easting;
  488. $this->northing = $northing;
  489. $this->latZone = $latZone;
  490. $this->lngZone = $lngZone;
  491. }
  492. /**
  493. * Return a string representation of this UTM reference
  494. *
  495. * @return
  496. */
  497. function toString() {
  498. return $this->lngZone . $this->latZone . " " .
  499. $this->easting . " " . $this->northing;
  500. }
  501. /**
  502. * Convert this UTM reference to a latitude and longitude
  503. *
  504. * @return the converted latitude and longitude
  505. */
  506. function toLatLng() {
  507. $wgs84 = new RefEll(6378137, 6356752.314);
  508. $UTM_F0 = 0.9996;
  509. $a = $wgs84->maj;
  510. $eSquared = $wgs84->ecc;
  511. $ePrimeSquared = $eSquared / (1.0 - $eSquared);
  512. $e1 = (1 - sqrt(1 - $eSquared)) / (1 + sqrt(1 - $eSquared));
  513. $x = $this->easting - 500000.0;;
  514. $y = $this->northing;
  515. $zoneNumber = $this->lngZone;
  516. $zoneLetter = $this->latZone;
  517. $longitudeOrigin = ($zoneNumber - 1.0) * 6.0 - 180.0 + 3.0;
  518. // Correct y for southern hemisphere
  519. if ((ord($zoneLetter) - ord("N")) < 0) {
  520. $y -= 10000000.0;
  521. }
  522. $m = $y / $UTM_F0;
  523. $mu =
  524. $m
  525. / ($a
  526. * (1.0
  527. - $eSquared / 4.0
  528. - 3.0 * $eSquared * $eSquared / 64.0
  529. - 5.0
  530. * pow($eSquared, 3.0)
  531. / 256.0));
  532. $phi1Rad =
  533. $mu
  534. + (3.0 * $e1 / 2.0 - 27.0 * pow($e1, 3.0) / 32.0) * sin(2.0 * $mu)
  535. + (21.0 * $e1 * $e1 / 16.0 - 55.0 * pow($e1, 4.0) / 32.0)
  536. * sin(4.0 * $mu)
  537. + (151.0 * pow($e1, 3.0) / 96.0) * sin(6.0 * $mu);
  538. $n =
  539. $a
  540. / sqrt(1.0 - $eSquared * sin($phi1Rad) * sin($phi1Rad));
  541. $t = tan($phi1Rad) * tan($phi1Rad);
  542. $c = $ePrimeSquared * cos($phi1Rad) * cos($phi1Rad);
  543. $r =
  544. $a
  545. * (1.0 - $eSquared)
  546. / pow(
  547. 1.0 - $eSquared * sin($phi1Rad) * sin($phi1Rad),
  548. 1.5);
  549. $d = $x / ($n * $UTM_F0);
  550. $latitude = (
  551. $phi1Rad
  552. - ($n * tan($phi1Rad) / $r)
  553. * ($d * $d / 2.0
  554. - (5.0
  555. + (3.0 * $t)
  556. + (10.0 * $c)
  557. - (4.0 * $c * $c)
  558. - (9.0 * $ePrimeSquared))
  559. * pow($d, 4.0)
  560. / 24.0
  561. + (61.0
  562. + (90.0 * $t)
  563. + (298.0 * $c)
  564. + (45.0 * $t * $t)
  565. - (252.0 * $ePrimeSquared)
  566. - (3.0 * $c * $c))
  567. * pow($d, 6.0)
  568. / 720.0)) * (180.0 / pi());
  569. $longitude = $longitudeOrigin + (
  570. ($d
  571. - (1.0 + 2.0 * $t + $c) * pow($d, 3.0) / 6.0
  572. + (5.0
  573. - (2.0 * $c)
  574. + (28.0 * $t)
  575. - (3.0 * $c * $c)
  576. + (8.0 * $ePrimeSquared)
  577. + (24.0 * $t * $t))
  578. * pow($d, 5.0)
  579. / 120.0)
  580. / cos($phi1Rad)) * (180.0 / pi());
  581. return new LatLng($latitude, $longitude);
  582. }
  583. }
  584. // ================================================================== RefEll
  585. class RefEll {
  586. var $maj;
  587. var $min;
  588. var $ecc;
  589. /**
  590. * Create a new RefEll object to represent a reference ellipsoid
  591. *
  592. * @param maj the major axis
  593. * @param min the minor axis
  594. */
  595. function RefEll($maj, $min) {
  596. $this->maj = $maj;
  597. $this->min = $min;
  598. $this->ecc = (($maj * $maj) - ($min * $min)) / ($maj * $maj);
  599. }
  600. }
  601. // ================================================== Mathematical Functions
  602. function sinSquared($x) {
  603. return sin($x) * sin($x);
  604. }
  605. function cosSquared($x) {
  606. return cos($x) * cos($x);
  607. }
  608. function tanSquared($x) {
  609. return tan($x) * tan($x);
  610. }
  611. function sec($x) {
  612. return 1.0 / cos($x);
  613. }
  614. /**
  615. * Take a string formatted as a six-figure OS grid reference (e.g.
  616. * "TG514131") and return a reference to an OSRef object that represents
  617. * that grid reference. The first character must be H, N, S, O or T.
  618. * The second character can be any uppercase character from A through Z
  619. * excluding I.
  620. *
  621. * @param ref
  622. * @return
  623. * @since 2.1
  624. */
  625. function getOSRefFromSixFigureReference($ref) {
  626. $char1 = substr($ref, 0, 1);
  627. $char2 = substr($ref, 1, 1);
  628. $east = substr($ref, 2, 3) * 100;
  629. $north = substr($ref, 5, 3) * 100;
  630. if ($char1 == 'H') {
  631. $north += 1000000;
  632. } else if ($char1 == 'N') {
  633. $north += 500000;
  634. } else if ($char1 == 'O') {
  635. $north += 500000;
  636. $east += 500000;
  637. } else if ($char1 == 'T') {
  638. $east += 500000;
  639. }
  640. $char2ord = ord($char2);
  641. if ($char2ord > 73) $char2ord--; // Adjust for no I
  642. $nx = (($char2ord - 65) % 5) * 100000;
  643. $ny = (4 - floor(($char2ord - 65) / 5)) * 100000;
  644. return new OSRef($east + $nx, $north + $ny);
  645. }
  646. /**
  647. * Work out the UTM latitude zone from the latitude
  648. *
  649. * @param latitude
  650. * @return
  651. */
  652. function getUTMLatitudeZoneLetter($latitude) {
  653. if ((84 >= $latitude) && ($latitude >= 72)) return "X";
  654. else if (( 72 > $latitude) && ($latitude >= 64)) return "W";
  655. else if (( 64 > $latitude) && ($latitude >= 56)) return "V";
  656. else if (( 56 > $latitude) && ($latitude >= 48)) return "U";
  657. else if (( 48 > $latitude) && ($latitude >= 40)) return "T";
  658. else if (( 40 > $latitude) && ($latitude >= 32)) return "S";
  659. else if (( 32 > $latitude) && ($latitude >= 24)) return "R";
  660. else if (( 24 > $latitude) && ($latitude >= 16)) return "Q";
  661. else if (( 16 > $latitude) && ($latitude >= 8)) return "P";
  662. else if (( 8 > $latitude) && ($latitude >= 0)) return "N";
  663. else if (( 0 > $latitude) && ($latitude >= -8)) return "M";
  664. else if (( -8 > $latitude) && ($latitude >= -16)) return "L";
  665. else if ((-16 > $latitude) && ($latitude >= -24)) return "K";
  666. else if ((-24 > $latitude) && ($latitude >= -32)) return "J";
  667. else if ((-32 > $latitude) && ($latitude >= -40)) return "H";
  668. else if ((-40 > $latitude) && ($latitude >= -48)) return "G";
  669. else if ((-48 > $latitude) && ($latitude >= -56)) return "F";
  670. else if ((-56 > $latitude) && ($latitude >= -64)) return "E";
  671. else if ((-64 > $latitude) && ($latitude >= -72)) return "D";
  672. else if ((-72 > $latitude) && ($latitude >= -80)) return "C";
  673. else return 'Z';
  674. }
  675. ?>