xybl.inc.php

下記ファイルの内容を表示しています。 ダウンロードを行うにはファイル名をクリックしてください。

<?php

//============================================================================
// 緯度経度⇔平面座標 変換
//============================================================================

// 原点座標位置を戻す
function gentenset($prj, &$b0,&$L0) {

	switch ($prj) {
		case  1:$b0=33.0; $L0=129.5;             break;
		case  2:$b0=33.0; $L0=131.0;             break;
		case  3:$b0=36.0; $L0=132.1666666666666; break;
		case  4:$b0=33.0; $L0=133.5;             break;
		case  5:$b0=36.0; $L0=134.3333333333333; break;
		case  6:$b0=36.0; $L0=136.0;             break;
		case  7:$b0=36.0; $L0=137.1666666666666; break;
		case  8:$b0=36.0; $L0=138.5;             break;
		case  9:$b0=36.0; $L0=139.833333333333;  break;
		case 10:$b0=40.0; $L0=140.833333333333;  break;
		case 11:$b0=44.0; $L0=140.25;            break;
		case 12:$b0=44.0; $L0=142.25;            break;
		case 13:$b0=44.0; $L0=144.25;            break;
		case 14:$b0=26.0; $L0=142.0;             break;
		case 15:$b0=26.0; $L0=127.5;             break;
		case 16:$b0=26.0; $L0=124.0;             break;
		case 17:$b0=26.0; $L0=131.0;             break;
		case 18:$b0=20.0; $L0=136.0;             break;
		case 19:$b0=26.0; $L0=154.0;             break;
	}

}

// BLtoXY 緯度経度を平面座標(測地座標)に変換
//
//      K :19座標系の系番号
// 	Lp:経度の度点
// 	bp:緯度の度点
//  var rX:19座標系のX座標
//  var rY:19座標系のY座標
// 	戻り値=子午線収差(度点表記)
//
function BLtoXY($K, $Lp,$bp, &$rX,&$rY) {

	gentenset($K,$b0p,$L0p);

	$ra=6377397.155;
	$rb=6356078.965;
	$k1=1.00503730604855;
	$k2=5.0478492403e-3;
	$k3=1.0563786831e-5;
	$k4=2.0633322e-8;
	$m0=0.9999;

	$b=$bp; $L=$Lp; $b0=$b0p; $L0=$L0p;

	$r0=180/pi();
	$b0=$b0/$r0;
	$L0=$L0/$r0;
	$b=$b/$r0;
	$L=$L/$r0;

	$e =sqrt(1-pow(($rb/$ra),2));
	$e1=sqrt(pow(($ra/$rb),2)-1);
	$w=sqrt(1-pow(($e*sin($b)),2));
	$n=$ra/$w;
	$a1=$k1*($b-$b0);
	$a2=$k2*(sin(2*$b)-sin(2*$b0))/2;
	$a3=$k3*(sin(4*$b)-sin(4*$b0))/4;
	$a4=$k4*(sin(6*$b)-sin(6*$b0))/6;
	$db=$ra*(1-pow($e,2))*($a1-$a2+$a3-$a4);
	$t=tan($b);
	$h=$e1*cos($b);
	$dl=$L-$L0;

	$x=$db+pow($dl,2)*$n/2*sin($b)*cos($b);
	$x=$x+pow($dl,4)*$n/24*sin($b)*pow(cos($b),3)*(5-pow($t,2)+9*pow($h,2)+4*pow($h,4));
	$x=$x+pow($dl,6)*$n/720*sin($b)*pow(cos($b),5)*(61-58*pow($t,2)+pow($t,4)+270*pow($h,2)-330*pow(($h*$t),2));
	$y=$dl*$n*cos($b)+pow($dl*cos($b),3)*$n/6*(1-pow($t,2)+pow($h,2));
	$y=$y+pow(($dl*cos($b)),5)*$n/120*(5-18*pow($t,2)+pow($t,4)+14*pow($h,2)-58*pow(($h*$t),2));

	$g=$dl*sin($b)+pow($dl,3)/3*sin($b)*pow(cos($b),2)*(1+3*pow($h,2)+2*pow($h,4))+pow($dl,5)/15*sin($b)*pow(cos($b),4)*(2-pow($t,2));

	$x=$x*$m0;
	$y=$y*$m0;
	$g=$g*$r0;
	$rX=$y;
	$rY=$x;

	return($g);

}

// XYtoBL 測地座標を緯度経度に変換する
//
//      K :19座標系の系番号
// 	yp:測地X座標
// 	xp:測地Y座標
//  var	oB:経度
//  var	oL:緯度
// 	戻り値=子午線収差
//
function XYtoBL($K, $yp,$xp, &$oB,&$oL) {

	gentenset($K,$b0p,$L0p);

	$ra=6377397.155;
	$rb=6356078.963;
	$k1=1.005037306048555;
	$k2=5.0478492403e-3;
	$k3=1.0563786831e-5;
	$k4=2.0633322e-8;
	$m0=0.9999;

	$x=$xp; $y=$yp; $b0=$b0p; $L0=$L0p;

	$r0=180/pi();
	$b0=$b0/$r0;
	$L0=$L0/$r0;

	$e =sqrt(1-pow(($rb/$ra),2));
	$e1=sqrt(pow(($ra/$rb),2)-1);
	$a=$ra*(1-pow($e,2));
	$d=$a*($k1*$b0 - $k2/2*sin(2*$b0) + $k3/4*sin(4*$b0) - $k4/6*sin(6*$b0)) + $x/$m0;
	for ($i=1; $i<=6; $i++) {
		if ($i==1) {
			$p[$i]=$b0;
		} else {
			$p[$i]=$p[$i-1]-($s[$i-1]-$d)/$m[$i-1];
		}
		$s[$i]=$a*($k1*$p[$i] - $k2/2*sin(2*$p[$i]) + $k3/4*sin(4*$p[$i]) - $k4/6*sin(6*$p[$i]));
		$m[$i]=$a/(pow((1-pow(($e*sin($p[$i])),2)),1.5));
	}
	$pp=$p[5];
	$y1=$y/$m0;
	$t =tan($pp);
	$h =$e1*cos($pp);
	$n =$ra/(pow((1-pow(($e*sin($pp)),2)),0.5));
	$mm=$m[5];

	$LB=$pp-(pow($y1,2)*$t/(2*$mm*$n)-pow($y1,4)*$t*(5+3*pow($t,2)
	   +pow($h,2)-9*pow(($t*$h),2))/(24*$mm*pow($n,3)));
	$LB=$LB-(pow($y1,6)*$t*(61+90*pow($t,2)+45*pow($t,4))/(720*$mm*pow($n,5)));

	$L=($y1/($n*cos($pp))-pow($y1,3)*(1+2*pow($t,2)+pow($h,2))/(6*pow($n,3)*cos($pp)));
	$L=$L+ (pow($y1,5)*(5+28*pow($t,2)+24*pow($t,4))/(120*pow($n,5)*cos($pp)));
	$L=$L+$L0;

	$G=($y1*$t/$n-pow($y1,3)*$t*(1+pow($t,2)-pow($h,2))/(3*pow($n,3))
	  +pow(($y1/$n),5)*$t*(2+5*pow($t,2)+3*pow($t,4))/15);

	$L=$L*$r0;
	$LB=$LB*$r0;
	$G=$G*$r0;
	$oL=$LB;
	$oB=$L;

	return($G);

}

?>