Перевод координат из Sk42 в Wsg84: отладка

[whoim]

Добрый день всем! Полный новичок, и при этом далеко не математик. Возникла необходимость на php переводить координаты из SK42 в WSG84, нашел где то на просторах этого форума скрипт на VB, делающий это. Тщательно перевел на PHP. И вот: ошибки.

Скрипт:

Код: Выделить всё

<?
#преобразование координат из SK42 в WSG84 и наоборот
# по методу http://gps-club.ru/gps_think/detail.php?ID=4886

#константы
$Pi = 3.14159265358979; #Число Пи 
$ro = 206264.8062; #Число угловых секунд в радиане

#Эллипсоид Красовского 
$aP = 6378245; #Большая полуось
$alP = 1 / 298.3; #Сжатие 
$e2P =  2 * $alP - $alP ^ 2; #Квадрат эксцентриситета;

#Эллипсоид WGS84 (GRS80, эти два эллипсоида сходны по большинству параметров)
$aW  = 6378137; #Большая полуось
$alW  = 1 / 298.257223563; #Сжатие
$e2W  = 2 * $alW - $alW ^ 2; #Квадрат эксцентриситета

#Вспомогательные значения для преобразования эллипсоидов
$a  = ($aP + $aW) / 2;
$e2  = ($e2P + $e2W) / 2;
$da  = $aW - $aP;
$de2  = $e2W - $e2P;

#Линейные элементы трансформирования, в метрах
$dx  = 23.92;
$dy  = -141.27;
$dz  = -80.9;
#Угловые элементы трансформирования, в секундах
$wx  = 0;
$wy  = 0;
$wz  = 0;
#Дифференциальное различие масштабов
$ms  = 0;

function WGS84_SK42_Lat($Bd, $Ld, $H) {
  return ($Bd - dB_($Bd, $Ld, $H) / 3600);
}

function SK42_WGS84_Lat($Bd, $Ld, $H) {
  return ($Bd + dB_($Bd, $Ld, $H) / 3600);
}

function WGS84_SK42_Long($Bd, $Ld, $H) {
  return ($Ld - dL_($Bd, $Ld, $H) / 3600);
}

function SK42_WGS84_Long($Bd, $Ld, $H) {
  return ($Ld + dL_($Bd, $Ld, $H) / 3600);
}

function dB_($Bd, $Ld, $H) {
  $B = $Bd * $Pi / 180;
  $L = $Ld * $Pi / 180;
  $M = $a * (1 - $e2) / (1 - $e2 * Sin($B) ^ 2) ^ 1.5;
  $N = $a * (1 - $e2 * Sin($B) ^ 2) ^ -0.5;
  return $ro / ($M + $H) * ($N / $a * $e2 * Sin($B) * Cos($B) * $da + ($N ^ 2 / $a ^ 2 + 1) * $N * Sin($B) * Cos($B) * $de2 / 2 - ($dx * Cos($L) + $dy * Sin($L)) * Sin($B) + $dz * Cos($B)) - $wx * Sin($L) * (1 + $e2 * Cos(2 * $B)) + $wy * Cos($L) * (1 + $e2 * Cos(2 * $B)) - $ro * $ms * $e2 * Sin($B) * Cos($B) ;
}

function dL_($Bd, $Ld, $H) { 
  $B = $Bd * $Pi / 180;
  $L = $Ld * $Pi / 180;
  $N = $a * (1 - $e2 * Sin($B) ^ 2) ^ -0.5;
  return $ro / (($N + $H) * Cos($B)) * (($dx * -1) * Sin($L) + $dy * Cos($L)) + Tan($B) * (1 - $e2) * ($wx * Cos($L) + $wy * Sin($L)) - $wz ;
}

function WGS84Alt($Bd, $Ld, $H) { 
  $B = $Bd * $Pi / 180;
  $L = $Ld * $Pi / 180;
  $N = $a * (1 - $e2 * Sin($B) ^ 2) ^ -0.5;
  $dH = ($a * -1) / $N * $da + $N * Sin($B) ^ 2 * $de2 / 2 + ($dx * Cos($L) + $dy * Sin($L)) * Cos($B) + $dz * Sin($B) - $N * $e2 * Sin($B) * Cos($B) * ($wx / $ro * Sin($L) - $wy / $ro * Cos($L)) + ($a ^ 2 / $N + $H) * $ms;
  return ($H + $dH);
}

$h_ = 44.830063;
$v_ = 38.831393;
$alt_ = 82.3;

echo 'sk42 H:'.$h_.' V: '.$v_.' ALT: '.$alt_;
echo '<hr>';
echo 'wsg84 H:'.SK42_WGS84_Lat($h_,$v_,$alt_).' V: '.SK42_WGS84_Long($h_,$v_,$alt_).' ALT: '.WGS84Alt($h_,$v_,$alt_);
?>

Ошибки:

Код: Выделить всё

Warning: Division by zero in Z:\home\fototrack.ru\www\inc\sk42_to_wsg84.php on line 57

Warning: Division by zero in Z:\home\fototrack.ru\www\inc\sk42_to_wsg84.php on line 57

Warning: Division by zero in Z:\home\fototrack.ru\www\inc\sk42_to_wsg84.php on line 71

Warning: Division by zero in Z:\home\fototrack.ru\www\inc\sk42_to_wsg84.php on line 71

Warning: Division by zero in Z:\home\fototrack.ru\www\inc\sk42_to_wsg84.php on line 71

Warning: Division by zero in Z:\home\fototrack.ru\www\inc\sk42_to_wsg84.php on line 71

Происходит деление на ноль в результирующих формулах функций dL_ и dB_. Прошу подсказать - отчего

[Максим Дубинин]

вы форумы не перепутали? у вас источник - gps-club, хотя они свою статью у нас скопировали, с тех пор были какие-то исправления, точно не помню какие
http://gis-lab.info/qa/wgs84-sk42-wgs84-formula.html

[whoim]

Максим Дубинин, здравствуйте! )
Нет, не перепутал - в копирайте стоял именно Ваш форум, и статью я находил на нем но потом потерял. Нашел только у них. Думаю, лучше к первоисточнику.

Спасибо за ссылку! Пошел курить.

[Александр Мурый]

Может, вам поможет proj4php?

[whoim]

Покурил. Никаких различий к сожалению. Продолжаю искать какую либо синтаксическую или организационную ошибку, но таких пока нет.

В некоторых местах кода стречается нижнее подчеркивание: Cos(L)) _ + Tan(B)

Я прочитал что это символ переноса строки, решил что он нужен чтобы связать формулы на разных строчках (удобночитаемость) и исключил их. Верно сделал?

"Может, вам поможет proj4php?"
Не хотелось ставить пушку чтобы стрелять по воробьям Посмотрю, что это такое.

[Voltron]

В некоторых местах кода стречается нижнее подчеркивание: Cos(L)) _ + Tan(B)

В Basic'е нижнее подчеркиванием используется для разрыва длинных строк. Т.е. это исключительно для удобочитаемости, так что вы все правильно сделали.

[whoim]

спасибо.
Может быть, входные значения должны быть другого формата?
$h_ = 44.830063;
$v_ = 38.831393;
$alt_ = 82.3;

это взято из трека Ozi Explorer

[Максим Дубинин]

PHP-скрипты для пересчета системы формата координат

[whoim]

Максим Дубинин, огромное спасибо. Чего ж поиском то я не смог найти..

Но только все равно есть проблема Готового скрипта для моих нужд там нет, а разобраться КАК надо делать.. Если бы я разобрался в деталях, то, думаю, сюда бы не написал. Крыша едет от длинных формул, да и теории навигационных дел в голове - нет..

В общем, несколько путей вырисовывалось, буду думать. Спасибо!

[Максим Дубинин]

не сдавайтесь - все получится, думаю комбинация вашего и по ссылке должна проблему решить и без "въезжания" в геодезию

если все получится и тесты совпадут - давайте все-таки опубликуем результат на странице формул

[whoim]

Если все получится - публиковать можно будет кому угодно и где угодно
Дело в том, что хочется сохранить идентичность уже отработанному Вами алгоритму. Просто представлять его на другом языке. Сдаваться, в общем, не собираюсь

[whoim]

=WGS84_SK42_Lat(50;50;0) и получить результат 49.99980414

пока не выходит. Ошибка деления на ноль..

[whoim]

"ТАВО РОТ" как говорит мой друг. Переменные, обозначенные так, внутри функций не видны..
Счас допилю )

[whoim]

Полностью рабочий вариант:

Код: Выделить всё

<?
#преобразование координат из SK42 в wgs84 и наоборот
# по методу http://gis-lab.info/qa/wgs84-sk42-wgs84-formula.html

#константы
define(Pi, 3.14159265358979); #Число Пи 
define(ro, 206264.806); #Число угловых секунд в радиане 206264.8062

#Эллипсоид Красовского 
define(aP, 6378245); #Большая полуось
define(alP, 1 / 298.3); #Сжатие 
define(e2P,  2 * alP - alP ^ 2); #Квадрат эксцентриситета;

#Эллипсоид WGS84 (GRS80, эти два эллипсоида сходны по большинству параметров)
define(aW, 6378137); #Большая полуось
define(alW, 1 / 298.257223563); #Сжатие
define(e2W, 2 * alW - alW ^ 2); #Квадрат эксцентриситета

#Вспомогательные значения для преобразования эллипсоидов
define(a, (aP + aW) / 2);
define(e2, (e2P + e2W) / 2);
define(da, aW - aP);
define(de2, e2W - e2P);

#Линейные элементы трансформирования, в метрах
define(dx, 23.92);
define(dy, -141.27);
define(dz, -80.9);
#Угловые элементы трансформирования, в секундах
define(wx, 0);
define(wy, 0);
define(wz, 0);
#Дифференциальное различие масштабов
define(ms, 0);

function WGS84_SK42_Lat($Bd, $Ld, $H) {
  return ($Bd - dB_($Bd, $Ld, $H) / 3600);
}

function SK42_WGS84_Lat($Bd, $Ld, $H) {
  return ($Bd + dB_($Bd, $Ld, $H) / 3600);
}

function WGS84_SK42_Long($Bd, $Ld, $H) {
  return ($Ld - dL_($Bd, $Ld, $H) / 3600);
}

function SK42_WGS84_Long($Bd, $Ld, $H) {
  return ($Ld + dL_($Bd, $Ld, $H) / 3600);
}

function dB_($Bd, $Ld, $H) {
  $B = $Bd * Pi / 180;
  $L = $Ld * Pi / 180;
  $M = a * (1 - e2) / (1 - e2 * Sin($B) ^ 2) ^ 1.5;
  $N = a * (1 - e2 * Sin($B) ^ 2) ^ -0.5;
  return ( ro / ($M + $H) * ($N / a * e2 * Sin($B) * Cos($B) * da + ($N ^ 2 / a ^ 2 + 1) * $N * Sin($B) * Cos($B) * de2 / 2 - (dx * Cos($L) + dy * Sin($L)) * Sin($B) + dz * Cos($B)) - wx * Sin($L) * (1 + e2 * Cos(2 * $B)) + wy * Cos($L) * (1 + e2 * Cos(2 * $B)) - ro * ms * e2 * Sin($B) * Cos($B) );
}

function dL_($Bd, $Ld, $H) { 
  $B = $Bd * Pi / 180;
  $L = $Ld * Pi / 180;
  $N = a * (1 - e2 * Sin($B) ^ 2) ^ -0.5;
  return ro / (($N + $H) * Cos($B)) * (-dx * Sin($L) + dy * Cos($L)) + Tan($B) * (1 - e2) * (wx * Cos($L) + wy * Sin($L)) - wz ;
}

function WGS84Alt($Bd, $Ld, $H) { 
  $B = $Bd * Pi / 180;
  $L = $Ld * Pi / 180;
  $N = a * (1 - e2 * Sin($B) ^ 2) ^ -0.5;
  $dH = -a / $N * da + $N * Sin($B) ^ 2 * de2 / 2 + (dx * Cos($L) + dy * Sin($L)) * Cos($B) + dz * Sin($B) - $N * e2 * Sin($B) * Cos($B) * (wx / ro * Sin($L) - wy / ro * Cos($L)) + (a ^ 2 / $N + $H) * ms;
  return ($H + $dH);
}

$h_ = 50;
$v_ = 50;
$alt_ = 0;

echo 'sk42 H:'.$h_;
echo '<hr>';
echo 'wgs84 H:'.WGS84_SK42_Lat($h_,$v_,$alt_);
?>

[whoim]

Если вдруг будет синтаксическая ошибка define (на старых php) - заключите имя константы - первый параметр define - в кавычки

Код: Выделить всё

define("Pi", 3.14159265358979); #Число Пи