diff --git a/resources/papers/1002.1417.pdf b/resources/papers/1002.1417.pdf new file mode 100644 index 000000000..46d486ab3 Binary files /dev/null and b/resources/papers/1002.1417.pdf differ diff --git a/src/GeographicPoint.php b/src/GeographicPoint.php index 0241c7a12..d35d8bd03 100644 --- a/src/GeographicPoint.php +++ b/src/GeographicPoint.php @@ -1797,12 +1797,16 @@ public function transverseMercator( $f = $ellipsoid->getFlattening(); $n = $f / (2 - $f); - $B = ($a / (1 + $n)) * (1 + $n ** 2 / 4 + $n ** 4 / 64); + $B = ($a / (1 + $n)) * (1 + $n ** 2 / 4 + $n ** 4 / 64 + $n ** 6 / 256 + (25 / 16384) * $n ** 8); - $h1 = $n / 2 - (2 / 3) * $n ** 2 + (5 / 16) * $n ** 3 + (41 / 180) * $n ** 4; - $h2 = (13 / 48) * $n ** 2 - (3 / 5) * $n ** 3 + (557 / 1440) * $n ** 4; - $h3 = (61 / 240) * $n ** 3 - (103 / 140) * $n ** 4; - $h4 = (49561 / 161280) * $n ** 4; + $h1 = $n / 2 - (2 / 3) * $n ** 2 + (5 / 16) * $n ** 3 + (41 / 180) * $n ** 4 - (127 / 288) * $n ** 5 + (7891 / 37800) * $n ** 6 + (72161 / 387072) * $n ** 7 - (18975107 / 50803200) * $n ** 8; + $h2 = (13 / 48) * $n ** 2 - (3 / 5) * $n ** 3 + (557 / 1440) * $n ** 4 + (281 / 630) * $n ** 5 - (1983433 / 1935360) * $n ** 6 + (13769 / 28800) * $n ** 7 + (148003883 / 174182400) * $n ** 8; + $h3 = (61 / 240) * $n ** 3 - (103 / 140) * $n ** 4 + (15061 / 26880) * $n ** 5 + (167603 / 181440) * $n ** 6 - (67102379 / 29030400) * $n ** 7 + (79682431 / 79833600) * $n ** 8; + $h4 = (49561 / 161280) * $n ** 4 - (179 / 168) * $n ** 5 + (6601661 / 7257600) * $n ** 6 + (97445 / 49896) * $n ** 7 - (40176129013 / 7664025600) * $n ** 8; + $h5 = (34729 / 80640) * $n ** 5 - (3418889 / 1995840) * $n ** 6 + (14644087 / 9123840) * $n ** 7 + (2605413599 / 622702080) * $n ** 8; + $h6 = (212378941 / 319334400) * $n ** 6 - (30705481 / 10378368) * $n ** 7 + (175214326799 / 58118860800) * $n ** 8; + $h7 = (1522256789 / 1383782400) * $n ** 7 - (16759934899 / 3113510400) * $n ** 8; + $h8 = (1424729850961 / 743921418240) * $n ** 8; if ($latitudeOrigin === 0.0) { $mO = 0; @@ -1818,7 +1822,11 @@ public function transverseMercator( $xiO2 = $h2 * sin(4 * $xiO0); $xiO3 = $h3 * sin(6 * $xiO0); $xiO4 = $h4 * sin(8 * $xiO0); - $xiO = $xiO0 + $xiO1 + $xiO2 + $xiO3 + $xiO4; + $xiO5 = $h5 * sin(10 * $xiO0); + $xiO6 = $h6 * sin(12 * $xiO0); + $xiO7 = $h7 * sin(14 * $xiO0); + $xiO8 = $h8 * sin(16 * $xiO0); + $xiO = $xiO0 + $xiO1 + $xiO2 + $xiO3 + $xiO4 + $xiO5 + $xiO6 + $xiO7 + $xiO8; $mO = $B * $xiO; } @@ -1834,8 +1842,16 @@ public function transverseMercator( $eta3 = $h3 * cos(6 * $xi0) * sinh(6 * $eta0); $xi4 = $h4 * sin(8 * $xi0) * cosh(8 * $eta0); $eta4 = $h4 * cos(8 * $xi0) * sinh(8 * $eta0); - $xi = $xi0 + $xi1 + $xi2 + $xi3 + $xi4; - $eta = $eta0 + $eta1 + $eta2 + $eta3 + $eta4; + $xi5 = $h5 * sin(10 * $xi0) * cosh(10 * $eta0); + $eta5 = $h5 * cos(10 * $xi0) * sinh(10 * $eta0); + $xi6 = $h6 * sin(12 * $xi0) * cosh(12 * $eta0); + $eta6 = $h6 * cos(12 * $xi0) * sinh(12 * $eta0); + $xi7 = $h7 * sin(14 * $xi0) * cosh(14 * $eta0); + $eta7 = $h7 * cos(14 * $xi0) * sinh(14 * $eta0); + $xi8 = $h8 * sin(16 * $xi0) * cosh(16 * $eta0); + $eta8 = $h8 * cos(16 * $xi0) * sinh(16 * $eta0); + $xi = $xi0 + $xi1 + $xi2 + $xi3 + $xi4 + $xi5 + $xi6 + $xi7 + $xi8; + $eta = $eta0 + $eta1 + $eta2 + $eta3 + $eta4 + $eta5 + $eta6 + $eta7 + $eta8; $easting = $falseEasting->asMetres()->getValue() + $kO * $B * $eta; $northing = $falseNorthing->asMetres()->getValue() + $kO * ($B * $xi - $mO); diff --git a/src/ProjectedPoint.php b/src/ProjectedPoint.php index 69b7584bd..12c880439 100644 --- a/src/ProjectedPoint.php +++ b/src/ProjectedPoint.php @@ -1893,12 +1893,16 @@ public function transverseMercator( $f = $ellipsoid->getFlattening(); $n = $f / (2 - $f); - $B = ($a / (1 + $n)) * (1 + $n ** 2 / 4 + $n ** 4 / 64); + $B = ($a / (1 + $n)) * (1 + $n ** 2 / 4 + $n ** 4 / 64 + $n ** 6 / 256 + (25 / 16384) * $n ** 8); - $h1 = $n / 2 - (2 / 3) * $n ** 2 + (37 / 96) * $n ** 3 - (1 / 360) * $n ** 4; - $h2 = (1 / 48) * $n ** 2 + (1 / 15) * $n ** 3 - (437 / 1440) * $n ** 4; - $h3 = (17 / 480) * $n ** 3 - (37 / 840) * $n ** 4; - $h4 = (4397 / 161280) * $n ** 4; + $h1 = $n / 2 - (2 / 3) * $n ** 2 + (5 / 16) * $n ** 3 + (41 / 180) * $n ** 4 - (127 / 288) * $n ** 5 + (7891 / 37800) * $n ** 6 + (72161 / 387072) * $n ** 7 - (18975107 / 50803200) * $n ** 8; + $h2 = (13 / 48) * $n ** 2 - (3 / 5) * $n ** 3 + (557 / 1440) * $n ** 4 + (281 / 630) * $n ** 5 - (1983433 / 1935360) * $n ** 6 + (13769 / 28800) * $n ** 7 + (148003883 / 174182400) * $n ** 8; + $h3 = (61 / 240) * $n ** 3 - (103 / 140) * $n ** 4 + (15061 / 26880) * $n ** 5 + (167603 / 181440) * $n ** 6 - (67102379 / 29030400) * $n ** 7 + (79682431 / 79833600) * $n ** 8; + $h4 = (49561 / 161280) * $n ** 4 - (179 / 168) * $n ** 5 + (6601661 / 7257600) * $n ** 6 + (97445 / 49896) * $n ** 7 - (40176129013 / 7664025600) * $n ** 8; + $h5 = (34729 / 80640) * $n ** 5 - (3418889 / 1995840) * $n ** 6 + (14644087 / 9123840) * $n ** 7 + (2605413599 / 622702080) * $n ** 8; + $h6 = (212378941 / 319334400) * $n ** 6 - (30705481 / 10378368) * $n ** 7 + (175214326799 / 58118860800) * $n ** 8; + $h7 = (1522256789 / 1383782400) * $n ** 7 - (16759934899 / 3113510400) * $n ** 8; + $h8 = (1424729850961 / 743921418240) * $n ** 8; if ($latitudeOrigin === 0.0) { $mO = 0; @@ -1914,10 +1918,23 @@ public function transverseMercator( $xiO2 = $h2 * sin(4 * $xiO0); $xiO3 = $h3 * sin(6 * $xiO0); $xiO4 = $h4 * sin(8 * $xiO0); - $xiO = $xiO0 + $xiO1 + $xiO2 + $xiO3 + $xiO4; + $xiO5 = $h5 * sin(10 * $xiO0); + $xiO6 = $h6 * sin(12 * $xiO0); + $xiO7 = $h7 * sin(14 * $xiO0); + $xiO8 = $h8 * sin(16 * $xiO0); + $xiO = $xiO0 + $xiO1 + $xiO2 + $xiO3 + $xiO4 + $xiO5 + $xiO6 + $xiO7 + $xiO8; $mO = $B * $xiO; } + $h1 = $n / 2 - (2 / 3) * $n ** 2 + (37 / 96) * $n ** 3 - (1 / 360) * $n ** 4 - (81 / 512) * $n ** 5 + (96199 / 604800) * $n ** 6 - (5406467 / 38707200) * $n ** 7 + (7944359 / 67737600) * $n ** 8; + $h2 = (1 / 48) * $n ** 2 + (1 / 15) * $n ** 3 - (437 / 1440) * $n ** 4 + (46 / 105) * $n ** 5 - (1118711 / 3870720) * $n ** 6 + (51841 / 1209600) * $n ** 7 + (24749483 / 348364800) * $n ** 8; + $h3 = (17 / 480) * $n ** 3 - (37 / 840) * $n ** 4 - (209 / 4480) * $n ** 5 + (5569 / 90720) * $n ** 6 + (9261899 / 58060800) * $n ** 7 - (6457463 / 17740800) * $n ** 8; + $h4 = (4397 / 161280) * $n ** 4 - (11 / 504) * $n ** 5 - (830251 / 7257600) * $n ** 6 + (466511 / 2494800) * $n ** 7 + (324154477 / 7664025600) * $n ** 8; + $h5 = (4583 / 161280) * $n ** 5 - (108847 / 3991680) * $n ** 6 - (8005831 / 63866880) * $n ** 7 + (22894433 / 124540416) * $n ** 8; + $h6 = (20648693 / 638668800) * $n ** 6 - (16363163 / 518918400) * $n ** 7 - (2204645983 / 12915302400) * $n ** 8; + $h7 = (219941297 / 5535129600) * $n ** 7 - (497323811 / 12454041600) * $n ** 8; + $h8 = (191773887257 / 3719607091200) * $n ** 8; + $eta = $easting / ($B * $kO); $xi = ($northing + $kO * $mO) / ($B * $kO); $xi1 = $h1 * sin(2 * $xi) * cosh(2 * $eta); @@ -1928,8 +1945,16 @@ public function transverseMercator( $eta3 = $h3 * cos(6 * $xi) * sinh(6 * $eta); $xi4 = $h4 * sin(8 * $xi) * cosh(8 * $eta); $eta4 = $h4 * cos(8 * $xi) * sinh(8 * $eta); - $xi0 = $xi - ($xi1 + $xi2 + $xi3 + $xi4); - $eta0 = $eta - ($eta1 + $eta2 + $eta3 + $eta4); + $xi5 = $h5 * sin(10 * $xi) * cosh(10 * $eta); + $eta5 = $h5 * cos(10 * $xi) * sinh(10 * $eta); + $xi6 = $h6 * sin(12 * $xi) * cosh(12 * $eta); + $eta6 = $h6 * cos(12 * $xi) * sinh(12 * $eta); + $xi7 = $h7 * sin(14 * $xi) * cosh(14 * $eta); + $eta7 = $h7 * cos(14 * $xi) * sinh(14 * $eta); + $xi8 = $h8 * sin(16 * $xi) * cosh(16 * $eta); + $eta8 = $h8 * cos(16 * $xi) * sinh(16 * $eta); + $xi0 = $xi - $xi1 - $xi2 - $xi3 - $xi4 - $xi5 - $xi6 - $xi7 - $xi8; + $eta0 = $eta - $eta1 - $eta2 - $eta3 - $eta4 - $eta5 - $eta6 - $eta7 - $eta8; $beta = self::asin(sin($xi0) / cosh($eta0)); diff --git a/tests/CoordinateOperation/AutoConversionTest.php b/tests/CoordinateOperation/AutoConversionTest.php index adcb59933..190a8defd 100644 --- a/tests/CoordinateOperation/AutoConversionTest.php +++ b/tests/CoordinateOperation/AutoConversionTest.php @@ -208,11 +208,11 @@ public function testBritishNationalGridToUTM(): void $to = $from->convert($toCRS); if (class_exists(OSTN15OSGM15Provider::class)) { - self::assertEqualsWithDelta(31326368.174145, $to->getEasting()->getValue(), 0.0001); - self::assertEqualsWithDelta(5708455.8991285, $to->getNorthing()->getValue(), 0.0001); + self::assertEqualsWithDelta(31326368.093447, $to->getEasting()->getValue(), 0.0001); + self::assertEqualsWithDelta(5708454.7196684, $to->getNorthing()->getValue(), 0.0001); } else { - self::assertEqualsWithDelta(31326366.159092, $to->getEasting()->getValue(), 0.0001); - self::assertEqualsWithDelta(5708455.7715515, $to->getNorthing()->getValue(), 0.0001); + self::assertEqualsWithDelta(31326366.078970, $to->getEasting()->getValue(), 0.0001); + self::assertEqualsWithDelta(5708454.600444621, $to->getNorthing()->getValue(), 0.0001); } } diff --git a/tests/GeographicPointTest.php b/tests/GeographicPointTest.php index f077e4ed6..b65303bf8 100644 --- a/tests/GeographicPointTest.php +++ b/tests/GeographicPointTest.php @@ -289,7 +289,7 @@ public function testDistanceDifferentCRSs(): void $to = GeographicPoint::create(Geographic2D::fromSRID(Geographic2D::EPSG_PZ_90), new Degree(51.507977), new Degree(-0.124588), null); if (class_exists(OSTN15OSGM15Provider::class)) { - self::assertEqualsWithDelta(3734.135, $from->calculateDistance($to)->getValue(), 0.001); + self::assertEqualsWithDelta(3735.308, $from->calculateDistance($to)->getValue(), 0.001); } else { self::assertEqualsWithDelta(3735.156, $from->calculateDistance($to)->getValue(), 0.001); } @@ -1117,8 +1117,8 @@ public function testOSTN15(): void $toCRS = Projected::fromSRID(Projected::EPSG_OSGB36_BRITISH_NATIONAL_GRID); $to = $from->OSTN15($toCRS, (new OSTN15OSGM15Provider())->provideGrid()); - self::assertEqualsWithDelta(651409.80373330, $to->getEasting()->asMetres()->getValue(), 0.00000001); - self::assertEqualsWithDelta(313177.44988696, $to->getNorthing()->asMetres()->getValue(), 0.00000001); + self::assertEqualsWithDelta(651409.804, $to->getEasting()->asMetres()->getValue(), 0.001); + self::assertEqualsWithDelta(313177.450, $to->getNorthing()->asMetres()->getValue(), 0.001); } public function testNADCON5Forward3DTransform3DPoint(): void