sjtsk2wgs84

  protected function _sjtsk2wgs84($xx, $yy, $hh)
  {
    # vstup X a Y souřadnice v S-JTKS a nadmořská výška
    # vystupem jsou souradnice WGS84, tedy longitude a latitude a upravena nadmořská výška

    # nadmorska vyska
    $hh += 45;

    # Vypocet zemepisnych souradnic z rovinnych souradnic
    $a = 6377397.15508;
    $e = 0.081696831215303;
    $n = 0.97992470462083;
    $konst_u_ro = 12310230.12797036;
    $sinUQ = 0.863499969506341;
    $cosUQ = 0.504348889819882;
    $sinVQ = 0.420215144586493;
    $cosVQ = 0.907424504992097;
    $alfa = 1.000597498371542;
    $k = 1.003419163966575;
    $ro = sqrt($xx*$xx + $yy*$yy);
    $epsilon = 2*atan($yy/($ro+$xx));
    $D = $epsilon/$n;
    $S = 2*atan(exp(1/$n*log($konst_u_ro/$ro)))-1*(M_PI)/2;
    $sinS = sin($S);
    $cosS = cos($S);
    $sinU = $sinUQ*$sinS-$cosUQ*$cosS*cos($D);
    $cosU = sqrt(1-$sinU*$sinU);
    $sinDV = sin($D)*$cosS/$cosU;
    $cosDV = sqrt(1-$sinDV*$sinDV);
    $sinV = $sinVQ*$cosDV-$cosVQ*$sinDV;
    $cosV = $cosVQ*$cosDV+$sinVQ*$sinDV;
    $Ljtsk = 2*atan($sinV/(1+$cosV))/$alfa;
    $t = exp(2/$alfa*log((1+$sinU)/$cosU/$k));
    $pom = ($t-1)/($t+1);

    $sinB = '';

    while (abs($pom-$sinB)>1E-15)
    {
      $sinB = $pom;
      $pom = $t*exp($e*log((1+$e*$sinB)/(1-$e*$sinB)));
      $pom = ($pom-1)/($pom+1);
    }
    $Bjtsk = atan($pom/sqrt(1-$pom*$pom));
    # Pravoúhlé souřadnice ve S-JTSK
    $a = 6377397.15508;
    $f_1 = 299.152812853;
    $e2 = 1-(1-1/$f_1)*(1-1/$f_1);
    $ro = $a/sqrt(1-$e2*sin($Bjtsk)*sin($Bjtsk));
    $x = ($ro+$hh)*cos($Bjtsk)*cos($Ljtsk);
    $y = ($ro+$hh)*cos($Bjtsk)*sin($Ljtsk);
    $z = ((1-$e2)*$ro+$hh)*sin($Bjtsk);

    # Pravoúhlé souřadnice v WGS-84
    $dx = 570.69;
    $dy = 85.69;
    $dz = 462.84;
    $wz = -5.2611/3600*M_PI/180;
    $wy = -1.58676/3600*M_PI/180;
    $wx = -4.99821/3600*M_PI/180;
    $m = 3.543e-6;
    $xn = $dx+(1+$m)*($x+$wz*$y-$wy*$z);
    $yn = $dy+(1+$m)*(-$wz*$x+$y+$wx*$z);
    $zn = $dz+(1+$m)*($wy*$x-$wx*$y+$z);

    # Geodetické souřadnice v systému WGS-84
    $a = 6378137.0;
    $f_1 = 298.257223563;
    $a_b = $f_1/($f_1-1);
    $p = sqrt($xn*$xn+$yn*$yn);
    $e2 = 1-(1-1/$f_1)*(1-1/$f_1);
    $theta = atan($zn*$a_b/$p);
    $st = sin($theta);
    $ct = cos($theta);
    $t = ($zn+$e2*$a_b*$a*$st*$st*$st)/($p-$e2*$a*$ct*$ct*$ct);
    $B = atan($t);
    $L = 2*atan($yn/($p+$xn));
    $H = sqrt(1+$t*$t)*($p-$a/sqrt(1+(1-$e2)*$t*$t));
    # Formát výstupních hodnot
    $B = $B/M_PI*180;
    $sirka = "N";
    if ($B<0) { $B = -$B; $sirka = "S"; };
    $stsirky = floor($B);
    $B = ($B-$stsirky)*60;
    $minsirky = floor($B);
    $B = ($B-$minsirky)*60;
    $vtsirky = round($B*1000)/1000;
    $sirka = $sirka . $stsirky . "°" . $minsirky . "'" . $vtsirky;
    $L = $L/M_PI*180;
    $delka = "E";
    if ($L<0) { $L = -$L; $delka = "W"; };
    $stdelky = floor($L);
    $L = ($L-$stdelky)*60;
    $mindelky = floor($L);
    $L = ($L-$mindelky)*60;
    $vtdelky = round($L*1000)/1000;
    $delka = $delka . $stdelky . "°" . $mindelky . "'" . $vtdelky;
    $vyska = round(($H-45)*100)/100;
    $vyska .= " m";

    # vratime vysledek, nejlepe asi jako pole tri stringu

    return array($sirka,$delka,$vyska);
  }