<?php # -*- PHP -*-
if ($_REQUEST["SOURCE"] != "" || $_REQUEST["source"] != "") {
printf("\n<FONT SIZE=1>\n");
# highlight_file("kkj-funcs.php");
show_source("kkj-funcs.php");
printf("</FONT><P><HR><P>\n");
include("../include/sign-oh2mqk.inc");
printf("<P></TR></TABLE>\n");
include("../include/base.inc");
include("../include/foot.inc");
exit;
}
######### HTML things ###########
function intlinks() {
langfilt('
FI: <P>Lähdemuoto:
FI: <A HREF="#KKJ-PI">KKJ P/I</A>,
FI: <A HREF="#KKJ-LALO">KKJ Leveys/Pituus</A>,
FI: <A HREF="#KKJ-LALO">WGS84 Leveys/Pituus</A>,
FI: <A HREF="#Maidenhead">Maidenhead</A>,
FI: <A HREF="#explanations">Selityksiä</A>.
EN: <P>Source format:
EN: <A HREF="#KKJ-PI">KKJ P/I</A>,
EN: <A HREF="#KKJ-LALO">KKJ Lat/Long</A>,
EN: <A HREF="#WGS-LALO">WGS84 Lat/Long</A>,
EN: <A HREF="#Maidenhead">Maidenhead</A>,
EN: <A HREF="#explanations">Explanations</A>.
');
}
function form_mh_digits($MhDigits) {
printf("<INPUT TYPE=radio NAME=MHDIGITS VALUE=4");
if ($MhDigits == 4) printf(" CHECKED");
printf("> 4 \n");
printf("<INPUT TYPE=radio NAME=MHDIGITS VALUE=6");
if ($MhDigits == 6) printf(" CHECKED");
printf("> 6 \n");
printf("<INPUT TYPE=radio NAME=MHDIGITS VALUE=10");
if ($MhDigits == 10) printf(" CHECKED");
printf("> 20m \n");
printf("<INPUT TYPE=radio NAME=MHDIGITS VALUE=12");
if ($MhDigits == 12) printf(" CHECKED");
printf("> 2m \n");
printf("<INPUT TYPE=radio NAME=MHDIGITS VALUE=16");
if ($MhDigits == 16) printf(" CHECKED");
printf("> 10mm \n");
}
####### SUPPORT STUFF ########
#
# Input: $LALO[La]: Latitude, radians, north positive
# $LALO[Lo]: Longitude, radians, east positive
#
function WithinFinland($LALO,$print) {
$La = rad2deg($LALO[La]);
$Lo = rad2deg($LALO[Lo]);
$rc = ((18.5 < $Lo && $Lo < 32.0) &&
(59.0 < $La && $La < 70.5));
# printf(" WithinFinland: La = %f Lo = %f\n",$La,$Lo);
if ($print && !$rc) {
langfilt('FI:
FI: <FONT COLOR=RED>KKJ EI OLE MÄÄRITELTY SUOMEN ULKOPUOLELLA!</FONT>
EN: <FONT COLOR=RED>KKJ IS NOT DEFINED OUTSIDE FINLAND!</FONT>
');
}
return($rc);
}
# Parse the base60 notation value..
function parsedms($str) {
$d = 0.0;
$m = 0.0;
$s = 0.0;
sscanf($str,"%f %f %f", &$d, &$m, &$s);
# printf(" d/m/s = %f %f %f\n", $d, $m, $s);
$sign = 1;
if ($d < 0.0) {
$sign = -1;
$d = 0.0 - $d;
}
return $sign * ($d + $m/60.0 + $s / 3600.0);
}
#
# Input: Degrees: positive east (or north) IN RADIANS!
# latitude: non-zero -> N/S, zero -> E/W
# Output: String
#
function fmtdms($angle,$latitude) {
$res = "E ";
$deg = rad2deg($angle);
if ($latitude) { $res = "N"; }
if ($deg < 0.0) {
$deg = -$deg;
$res = "W ";
if ($latitude) { $res = "S"; }
}
$deg += 0.0000000001;
$d = floor($deg);
$deg = 60.0 * ($deg - $d);
$m = floor($deg);
$s = 60.0 * ($deg - $m);
return( sprintf("%s%3d°%02d'%02.4f\"",$res,$d,$m,$s) );
}
#
# Input: Degrees: positive east (or north) IN RADIANS!
# latitude: non-zero -> N/S, zero -> E/W
# Output: String
#
function fmtdmm($angle,$latitude) {
$res = "E ";
$deg = rad2deg($angle) + 0.0000000001;
if ($latitude) { $res = "N"; }
if ($deg < 0.0) {
$deg = -$deg;
$res = "W ";
if ($latitude) { $res = "S"; }
}
$d = floor($deg);
$m = 60.0 * ($deg - $d);
return( sprintf("%s%3d°%02.7f'",$res,$d,$m) );
}
#
# Input: Degrees: positive east (or north) IN RADIANS!
# latitude: non-zero -> N/S, zero -> E/W
# Output: String
#
function fmtddd($angle,$latitude) {
$res = "E ";
$deg = rad2deg($angle) + 0.0000000001;
if ($latitude) { $res = "N"; }
if ($deg < 0.0) {
$deg = -$deg;
$res = "W ";
if ($latitude) { $res = "S"; }
}
return( sprintf("%s%3.9f° ",$res,$deg) );
}
function dms2deg($deg,$min,$sec) {
$sign = 1;
if ($deg < 0) {
$sign = -1;
$deg = -$deg;
}
$deg = $deg + $min / 60.0 + $sec / 3600.0;
return ($deg * $sign);
}
# ----------- MH support stuff -----------
function printmhdimensions($dX, $dY) {
global $LANG;
if ($dY >= 10000.0 && $dX >= 10000.0) {
$dX /= 1000.0;
$dY /= 1000.0;
if ($LANG == "fi")
printf(" Kilometreinä : korkeus ± %.3f km\n".
" keskipisteestä: leveys ± %.3f km\n",
$dX, $dY);
else
printf(" In kilometers: height ± %.3f km\n".
" from centre : width ± %.3f km\n",
$dX, $dY);
} else {
if ($LANG == "fi")
printf(" Metreinä : korkeus ± %.3f m\n".
" keskipisteestä: leveys ± %.3f m\n",
$dX, $dY);
else
printf(" In meters : height ± %.3f m\n".
" from centre: width ± %.3f m\n",
$dX, $dY);
}
}
function highlight_mh($MH, $digits) {
$mh1 = substr($MH, 0, $digits);
$mh2 = substr($MH, $digits);
return "<B>" . $mh1 . "</B>" . $mh2;
}
function mh_edge_distance($WGS, $MH) {
global $LANG;
# GRS80 ellipsoid:
$EL[a] = 6378137.0;
$EL[f] = 1.0/298.25722210088;
# WGS84 ellipsoid:
$EL[a] = 6378137.0;
$EL[f] = 1.0/298.257223563;
$MH6 = strtoupper(substr($MH, 0, 6));
$MH6w = parsemaidenhead($MH6);
$MH4 = strtoupper(substr($MH, 0, 4));
$MH4w = parsemaidenhead($MH4);
if ($MH6 != $MH4) {
$M[La] = $WGS[La];
$M[Lo] = $MH6w[Lo] - $MH6w[dLo];
$R = ellipsoidal_great_circle($M, $WGS, $EL);
$distW = $R[dist];
$M[La] = $WGS[La];
$M[Lo] = $MH6w[Lo] + $MH6w[dLo];
$R = ellipsoidal_great_circle($M, $WGS, $EL);
$distE = $R[dist];
$M[La] = $MH6w[La] - $MH6w[dLa];
$M[Lo] = $WGS[Lo];
$R = ellipsoidal_great_circle($M, $WGS, $EL);
$distS = $R[dist];
$M[La] = $MH6w[La] + $MH6w[dLa];
$M[Lo] = $WGS[Lo];
$R = ellipsoidal_great_circle($M, $WGS, $EL);
$distN = $R[dist];
if ($LANG == "fi") {
printf("Etäisyys MH ruudun %s reunoihin:\n",$MH6);
printf(" Länteen: %5.2f m Itään: %5.2f m\n",
$distW, $distE);
printf(" Pohjoiseen: %5.2f m Etelään: %5.2f m\n",
$distN, $distS);
} else {
printf("Distance to the edges of the MH square %s:\n",$MH6);
printf(" West: %5.2f m East: %5.2f m\n",
$distW, $distE);
printf(" North: %5.2f m South: %5.2f m\n",
$distN, $distS);
}
}
$M[La] = $WGS[La];
$M[Lo] = $MH4w[Lo] - $MH4w[dLo];
$R = ellipsoidal_great_circle($M, $WGS, $EL);
$distW = $R[dist];
$M[La] = $WGS[La];
$M[Lo] = $MH4w[Lo] + $MH4w[dLo];
$R = ellipsoidal_great_circle($M, $WGS, $EL);
$distE = $R[dist];
$M[La] = $MH4w[La] - $MH4w[dLa];
$M[Lo] = $WGS[Lo];
$R = ellipsoidal_great_circle($M, $WGS, $EL);
$distS = $R[dist];
$M[La] = $MH4w[La] + $MH4w[dLa];
$M[Lo] = $WGS[Lo];
$R = ellipsoidal_great_circle($M, $WGS, $EL);
$distN = $R[dist];
if ($LANG == "fi") {
printf("Etäisyys MH suurruudun %s reunoihin:\n",$MH4);
printf(" Länteen: %3.3f km Itään: %3.3f km\n",
$distW/1000.0, $distE/1000.0);
printf(" Pohjoiseen: %3.3f km Etelään: %3.3f km\n",
$distN/1000.0, $distS/1000.0);
} else {
printf("Distance to the edges of the MH major square %s:\n",$MH4);
printf(" West: %3.3f km East: %3.3f km\n",
$distW/1000.0, $distE/1000.0);
printf(" North: %3.3f km South: %3.3f km\n",
$distN/1000.0, $distS/1000.0);
}
}
######## MATH THINGS #########
$Long0 = deg2rad(27.0);
function KKJ_to_WGSxy($KKJ) {
if ($KKJ[X] > 72000000.0) { # North
$a = 1.0000006841;
$b = 0.0000015179;
$dx = -132.881;
$dy = -185.993;
} else { # South
$a = 0.9999955294;
$b = -0.0000068348;
$dx = -125.056;
$dy = -107.850;
}
$ED50[X] = $dx + ($a * $KKJ[X] - $b * $KKJ[Y]);
$ED50[Y] = $dy + ($b * $KKJ[X] + $a * $KKJ[Y]);
return $ED50;
}
#
# function WGSxy_to_KKJ($ED50) {
#
# if ($ED50[X] > 72000000.0) { ## North
# $a = 0.9999993159;
# $b = -0.0000015179;
# $dx = 132.881;
# $dy = 185.993;
# } else { ## South
# $a = 1.0000044706;
# $b = 0.0000068348;
# $dx = 125.056;
# $dy = 107.850;
# }
#
# $KKJ[X] = $dx + ($a * $ED50[X] - $b * $ED50[Y]);
# $KKJ[Y] = $dy + ($b * $ED50[X] + $a * $ED50[Y]);
#
# return $KKJ;
# }
#
function WGSxy_to_KKJ($ED50) {
# Here the EUREF89 Y coordinate is WITH the zone number!
$a = 1.0000021251;
$b = 0.000004439;
$dx = 132.858;
$dy = 132.490;
$KKJ[X] = $dx + ($a * $ED50[X] - $b * $ED50[Y]);
$KKJ[Y] = $dy + ($b * $ED50[X] + $a * $ED50[Y]);
return $KKJ;
}
function ED50xy_to_KKJ($ED50) {
$a = 1.00000075150;
$b = -0.00000439333;
$dx = -61.5805;
$dy = 95.6691;
$ED50[Y] -= 3000000.0; # MAGIC KNOWLEDGE!
# Here the ED50 Y coordinate is WITHOUT zone number
$KKJ[X] = $dx + ($a * $ED50[X] - $b * $ED50[Y]);
$KKJ[Y] = $dy + ($b * $ED50[X] + $a * $ED50[Y]);
$KKJ[Y] += 3000000.0; # MAGIC KNOWLEDGE!
return $KKJ;
}
# Input: WGS[La] is Latitude in Radians, north is positive
# WGS[Lo] is Longitude in Radians, east is positive
# Output: KKJ[La] is Latitude in Radians, north is positive
# KKJ[Lo] is Longitude in Radians, east is positive
#
# Accuracy has some limitations in range of 0.2-1.0 meters
#
function WGSlalo_to_KKJlalo($WGS) {
$La = rad2deg( $WGS[La] );
$Lo = rad2deg( $WGS[Lo] );
$dLa = deg2rad( -0.124766E+01 +
0.269941E+00 * $La +
-0.191342E+00 * $Lo +
-0.356086E-02 * $La * $La +
0.122353E-02 * $La * $Lo +
0.335456E-03 * $Lo * $Lo ) / 3600.0;
$dLo = deg2rad( 0.286008E+02 +
-0.114139E+01 * $La +
0.581329E+00 * $Lo +
0.152376E-01 * $La * $La +
-0.118166E-01 * $La * $Lo +
-0.826201E-03 * $Lo * $Lo ) / 3600.0;
$KKJ[La] = $WGS[La] + $dLa;
$KKJ[Lo] = $WGS[Lo] + $dLo;
return ($KKJ);
}
# Input: KKJ[La] is Latitude in Radians, north is positive
# KKJ[Lo] is Longitude in Radians, east is positive
# Output: WGS[La] is Latitude in Radians, north is positive
# WGS[Lo] is Longitude in Radians, east is positive
#
# Accuracy has some limitations in range of 0.2-1.0 meters
#
function KKJlalo_to_WGSlalo($KKJ) {
$La = rad2deg( $KKJ[La] );
$Lo = rad2deg( $KKJ[Lo] );
$dLa = deg2rad( 0.124867E+01 +
-0.269982E+00 * $La +
0.191330E+00 * $Lo +
0.356119E-02 * $La * $La +
-0.122312E-02 * $La * $Lo +
-0.335514E-03 * $Lo * $Lo ) / 3600.0;
$dLo = deg2rad( -0.286111E+02 +
0.114183E+01 * $La +
-0.581428E+00 * $Lo +
-0.152421E-01 * $La * $La +
0.118177E-01 * $La * $Lo +
0.826646E-03 * $Lo * $Lo ) / 3600.0;
$WGS[La] = $KKJ[La] + $dLa;
$WGS[Lo] = $KKJ[Lo] + $dLo;
return ($WGS);
}
# Inputs: $ED50[Lo] = Longitude (Radians)
# $ED50[La] = Latitude (Radians)
# Outputs: $ED50[X], $ED50[Y] (meters)
function KKJLaLo_to_KKJxy($INP, $ZoneNumber, $Long0) {
# Following method is NOT excellent for more
# than about 2deg from midline.
# It might be usable within a meter for wider
# projections
$Lo = $INP[Lo] - $Long0;
$a = 6378388.0; # Hayford ellipsoid
$f = 1/297.0;
# $a = 6378137.0; # GRS80 -> WGS84
# $f = 1/298.257222101;
$b = (1.0 - $f) * $a;
$bb = $b * $b; # b²
$c = ($a / $b) * $a; # c = a²/b
$ee = ($a * $a - $bb) / $bb; # e'² = (a² - b²)/b²
$n = ($a - $b)/($a + $b); # n = (a-b)/(a+b)
$nn = $n * $n; # n²
$cosLa = cos($INP[La]);
$NN = $ee * $cosLa * $cosLa; # Eta² = e'² cos² Latitude
$LaF = atan(tan($INP[La]) /
cos($Lo * sqrt(1 + $NN)));
$cosLaF = cos($LaF);
$t = (tan($Lo) * $cosLaF) /
sqrt(1 + $ee * $cosLaF * $cosLaF);
$A = $a / ( 1 + $n );
$A1 = $A * (1 + $nn / 4 + $nn * $nn / 64);
$A2 = $A * 1.5 * $n * (1 - $nn / 8);
$A3 = $A * 0.9375 * $nn * (1 - $nn / 4);
$A4 = $A * 35/48.0 * $nn * $n;
$OUT[X] = $A1 * $LaF -
$A2 * sin(2*$LaF) +
$A3 * sin(4*$LaF) -
$A4 * sin(6*$LaF);
$OUT[Y] = $c * log($t + sqrt(1+$t*$t)) +
500000.0 + $ZoneNumber * 1000000.0;
#printf("X1 = %.1f X2 = %.1f X3 = %.1f\nX4 = %.1f".
# " Lo0 = %.1f LaF = %.6f\n",
# $A1 * $LaF,
# 0 - $A2 * sin(2*$LaF),
# 0 + $A3 * sin(4*$LaF),
# 0 - $A4 * sin(6*$LaF), rad2deg($Long0), $LaF);
#printf("n = %.9f nn = %.9f\n", $n, $nn);
#
#printf(" ED50: P = %.1f I = %.1f\n",$OUT[X],$OUT[Y]);
return ($OUT);
}
function WGSLaLo_to_WGSxy($INP, $ZoneNumber, $Long0) {
# Following method is NOT excellent for more
# than about 2deg from midline.
# It might be usable within a meter for wider
# projections
$Lo = $INP[Lo] - $Long0;
# $a = 6378388; # Hayford ellipsoid
# $f = 1/297.0;
$a = 6378137; # GRS80 -> WGS84
$f = 1/298.257223564;
$b = (1.0 - $f) * $a;
$bb = $b * $b; # b²
$c = ($a / $b) * $a; # c = a²/b
$ee = ($a * $a - $bb) / $bb; # e'² = (a² - b²)/b²
$n = ($a - $b)/($a + $b); # n = (a-b)/(a+b)
$nn = $n * $n; # n²
$cosLa = cos($INP[La]);
$NN = $ee * $cosLa * $cosLa; # Eta² = e'² cos² Latitude
$LaF = atan(tan($INP[La]) /
cos($Lo * sqrt(1 + $NN)));
$cosLaF = cos($LaF);
$t = (tan($Lo) * $cosLaF) /
sqrt(1 + $ee * $cosLaF * $cosLaF);
$A = $a / ( 1 + $n );
$A1 = $A * (1 + $nn / 4 + $nn * $nn / 64);
$A2 = $A * 1.5 * $n * (1 - $nn / 8);
$A3 = $A * 0.9375 * $nn * (1 - $nn / 4);
$A4 = $A * 35/48.0 * $nn * $n;
$OUT[X] = $A1 * $LaF -
$A2 * sin(2*$LaF) +
$A3 * sin(4*$LaF) -
$A4 * sin(6*$LaF);
$OUT[Y] = $c * log($t + sqrt(1+$t*$t)) +
500000.0 + $ZoneNumber * 1000000.0;
#printf("X1 = %.1f X2 = %.1f X3 = %.1f\nX4 = %.1f".
# " Lo0 = %.1f LaF = %.6f\n",
# $A1 * $LaF,
# 0 - $A2 * sin(2*$LaF),
# 0 + $A3 * sin(4*$LaF),
# 0 - $A4 * sin(6*$LaF), rad2deg($Long0), $LaF);
#printf("n = %.9f nn = %.9f\n", $n, $nn);
#
#printf(" WGS: P = %.1f I = %.1f\n",$OUT[X],$OUT[Y]);
return ($OUT);
}
#
# Scan iteratively the target area, until find matching
# KKJ coordinate value. Area is defined in WGS84.
#
### BUGGY RESULTS ###
function KKJxy_to_WGSlalo($KKJ, $INITIAL, $ZoneNumber, $Long0) {
$delta = deg2rad(0.1);
$MinLa = $INITIAL[La] - $delta;
$MaxLa = $INITIAL[La] + $delta;
$MinLo = $INITIAL[Lo] - $delta;
$MaxLo = $INITIAL[Lo] + $delta;
global $Debug;
for ($i = 1; $i < 35; ++$i) {
$DeltaLa = $MaxLa - $MinLa;
$DeltaLo = $MaxLo - $MinLo;
$WGS[La] = $MinLa + 0.5 * $DeltaLa;
$WGS[Lo] = $MinLo + 0.5 * $DeltaLo;
if ($Debug) {
printf("Test%2d: %s %s\n", $i,
fmtdms($WGS[La],1), fmtdms($WGS[Lo],0));
printf(" Max: %s %s\n",
fmtdms($MaxLa,1), fmtdms($MaxLo,0));
printf(" Min: %s %s\n",
fmtdms($MinLa,1), fmtdms($MinLo,0));
}
$ED5w = WGSLaLo_to_WGSxy($WGS, $ZoneNumber, $Long0);
# printf(" WGS: P = %.1f I = %.1f\n",$ED5w[X],$ED5w[Y]);
$KKJt = WGSxy_to_KKJ($ED5w);
if ($Debug) {
printf(" WGS: P = %7.3f I = %7.3f\n",
$KKJt[X], $KKJt[Y]);
printf(" Delta: %7.3f %7.3f\n",
$KKJt[X]-$KKJ[X], $KKJt[Y]-$KKJ[Y]);
}
if ($KKJt[X] < $KKJ[X]) {
$MinLa = $MinLa + 0.45 * $DeltaLa;
} else {
$MaxLa = $MinLa + 0.55 * $DeltaLa;
}
if ($KKJt[Y] < $KKJ[Y]) {
$MinLo = $MinLo + 0.45 * $DeltaLo;
} else {
$MaxLo = $MinLo + 0.55 * $DeltaLo;
}
}
return ($WGS);
}
#
# Scan iteratively the target area, until find matching
# KKJ coordinate value. Area is defined with Hayford Ellipsoid.
#
function KKJxy_to_KKJlalo($KKJ, $ZoneNumber, $Long0) {
$MinLa = deg2rad(59.0);
$MaxLa = deg2rad(70.5);
$MinLo = deg2rad(18.5);
$MaxLo = deg2rad(32.0);
global $Debug;
for ($i = 1; $i < 35; ++$i) {
$DeltaLa = $MaxLa - $MinLa;
$DeltaLo = $MaxLo - $MinLo;
$LALO[La] = $MinLa + 0.5 * $DeltaLa;
$LALO[Lo] = $MinLo + 0.5 * $DeltaLo;
if ($Debug) {
printf("Test%2d: %s %s\n", $i,
fmtdms($LALO[La],1), fmtdms($LALO[Lo],0));
printf(" Max: %s %s\n",
fmtdms($MaxLa,1), fmtdms($MaxLo,0));
printf(" Min: %s %s\n",
fmtdms($MinLa,1), fmtdms($MinLo,0));
}
$KKJt = KKJLaLo_to_KKJxy($LALO, $ZoneNumber, $Long0);
if ($Debug) {
printf(" KKJ: P = %7.3f I = %7.3f\n",
$KKJt[X], $KKJt[Y]);
printf(" Delta: %7.3f %7.3f\n",
$KKJt[X]-$KKJ[X], $KKJt[Y]-$KKJ[Y]);
}
if ($KKJt[X] < $KKJ[X]) {
$MinLa = $MinLa + 0.45 * $DeltaLa;
} else {
$MaxLa = $MinLa + 0.55 * $DeltaLa;
}
if ($KKJt[Y] < $KKJ[Y]) {
$MinLo = $MinLo + 0.45 * $DeltaLo;
} else {
$MaxLo = $MinLo + 0.55 * $DeltaLo;
}
}
return ($LALO);
}
?>
Matti Aarnio <matti.aarnio@zmailer.org>; OH2MQK
