師匠の散歩

きままにPerlでも

Hubeny standard eq.

標準ヒュベニィの式で測地線距離を求める/ subDistanceHubeny.cgi

Hubenyの式で測地線距離を求めます。

当初は、カシミール3Dに掲載されていた日本測地系で使用される式を用いていましたが、VBAで勉強した世界測地系でのFubeny Standard EquationをPerlに移植しました。

つくば 36.10056    140.09111;
東京  35.655  139.74472;
測地線長:58502.34248 m
最初の地点の経緯度を入力してください(初期値:富士山−筑波山)
地点1: 地点2: 項数:

スクリプト

以下は参考資料です。正式版は、Maplib.pmで、Sub distanceHubenyStandard をご覧ください。

#  ------------------------------------------------------------------
#  ヒュベニの距離計算式(旧日本測地系)
#  http:# www.kashmir3d.com/kash/manual/std_siki.htm
#  緯度・経度はラジアンで与え、長さmを戻す
#  ------------------------------------------------------------------
sub subSokuchisen{
  my ($lat1,$lng1,$lat2,$lng2) = ;
  my $P     = ($lat1+$lat2)/2;  #  P : 2点の平均緯度
  my $dP    =  $lat2-$lat1;     #  dP: 2点の緯度差
  my $dR    =  $lng2-$lng1;     #  dR: 2点の経度差
  my $M     = 6334834/sqrt( (1-0.006674*sin($P)**2)**3 ) ;
  my $N     = 6377397/sqrt(1-0.006674*sin($Pt)**2);
  my $D     = sqrt(($M * $dP)*($M * $dP)+($N*cos($P)*$dR)**2);
return $D;
}
sub distanceHubenyStandard{
  my ($x1,$y1,$x2,$y2,$kou) = @_;
  $x1 = doRad($x1);
  $y1 = doRad($y1);
  $x2 = doRad($x2);
  $y2 = doRad($y2);
  $kou=10 if (! $kou);
    
  #f = 1 / 298.257222101;                                      # GRS80 扁平率 
  #a = 6378137 GRS80                                           # 長半期
  my $e2 = $Henpei * (2 - $Henpei );                        # 第一離心率の2乗
  my $eD2 = $e2 / (1 - $e2);                                # 第二離心率の2乗
  
  my $T = ($x1 + $x2) / 2;                                  # 2点の平均緯度
  my $dT = $x2 - $x1;                                       # 2点の緯度差
  my $dG = $y2 - $y1;
  my $W = sqrt(1 - $e2 * sin($T) ** 2);
  my $M = $R_long * (1 - $e2) / $W ** 3;                    # M: 子午線曲率半径
  my $N = $R_long / $W;                                     # N: 卯酉線曲率半径
  
  my $i2 = $eD2 * cos($T) ** 2;
  my $t2 = tan($T) ** 2;
  
  my $cost = cos($T);

# 項数セットで、目的項数までのHubenyの式を取り込む。#kouが無記入の場合は=10が入るので、全ての項数で計算する。
# $sSinAの8項と9項は合わせて8項と置いた。
# Switch文が使えるとシンプルになる

  my $sSinA = $N * $cost * $dG;
  my $sCosA = $N / (1 + $i2) * $dT;
  if($kou<1) { $kou=1;}
  if($kou>10){$kou=10;}

  if($kou<2){
    $sSinA += $N * $cost / 24 * (1 - $i2 + $i2 ** 2 - $i2 ** 3 - 9 * $t2 * $i2 + 18 * $t2 * $i2 ** 2 - 27 * $t2 * $i2 ** 3) * $dT ** 2 * $dG ;
    $sCosA += $N / 24 * (3 * $i2 - 6 * $i2 ** 2 + 9 * $i2 ** 3 - 3 * $t2 * $i2 + 21 * $t2 * $i2 ** 2 - 54 * $t2 * $i2 ** 3) * $dT ** 3 ;
  }
  if($kou<3){
    $sSinA += $N * $cost ** 3 / 24 * (-1 * $t2) * $dG ** 3;
    $sCosA += $N * $cost ** 2 / 24 * (-2 - 3 * $t2 + 3 * $t2 * $i2 - 3 * $t2 * $i2 ** 2 + 3 * $t2 * $i2 ** 3) * $dT * $dG ** 2 ;
  }
  if($kou<4){
    $sSinA += $N * $cost / 5760 * (7 + 10 * $i2 - 27 * $i2 ** 2 - 54 * $t2 * $i2 - 642 * $t2 * $i2 ** 2 + 675 * $t2 * $i2 ** 3) * $dT ** 4 * $dG;
    $sCosA += $N / 5760 * (-36 * $i2 + 207 * $i2 ** 2 + 36 * $t2 * $i2 - 1062 * $t2 * $i2 ** 2 + 135 * $t2 ** 2 * $i2 ** 2) * $dT ** 5 ;
  }
  if($kou<5){
    $sSinA += $N * $cost ** 3 / 5760 * (-16 - 70 * $t2 - 158 * $t2 * $i2 + 158 * $t2 * $i2 ** 2 + 90 * $t2 ** 2 * $i2 - 180 * $t2 ** 2 * $i2 ** 2) * $dT ** 2 * $dG ** 3;
    $sCosA += $N * $cost ** 2 / 5760 * (-16 - 60 * $t2 + 4 * $i2 - 4 * $i2 ** 2 + 102 * $t2 * $i2 + 48 * $t2 * $i2 ** 2 + 90 * $t2 ** 2 * $i2 - 630 * $t2 ** 2 * $i2 ** 2) * $dT ** 3 * $dG ** 2 ;
  }
  if($kou<6){
    $sSinA += $N * $cost ** 5 / 5760 * (-24 * $t2 + 3 * $t2 ** 2 - 27 * $t2 * $i2) * $dG ** 5 ;
    $sCosA += $N * $cost ** 4 / 5760 * (-8 - 20 * $t2 + 15 * $t2 ** 2 - 8 * $i2 + 96 * $t2 * $i2 - 15 * $t2 ** 2 * $i2 + 15 * $t2 ** 2 * $i2 ** 2) * $dT * $dG ** 4 ;
  }
  if($kou<7){
    $sSinA += $N * $cost / 1935360 * 62 * $dT ** 6 * $dG ;
    $sCosA += $N * $cost ** 2 / 1935360 * (-192 - 2016 * $t2) * $dT ** 5 * $dG ** 2 ;
  }
  if($kou<8){
    $sSinA += $N * $cost ** 3 / 1935360 * (-416 - 2954 * $t2) * $dT ** 4 * $dG ** 3 ;
    $sCosA += $N * $cost ** 4 / 1935360 * (256 + 784 * $t2 + 4200 * $t2 ** 2) * $dT ** 3 * $dG ** 4 ;
  }
  if($kou<9){
    $sSinA += $N * $cost ** 5 / 1935360 * (-192 - 1680 * $t2 + 2652 * $t2 ** 2) * $dT ** 2 * $dG ** 5 ;
    $sSinA += $N * $cost ** 7 / 1935360 * (-816 * $t2 + 528 * $t2 ** 2 - 6 * $t2 ** 3) * $dG ** 7;
    $sCosA += $N * $cost ** 6 / 1935360 * (-64 - 224 * $t2 + 1148 * $t2 ** 2 - 42 * $t2 ** 3) * $dT * $dG ** 6;;
  }
  return sqrt( $sSinA ** 2 + $sCosA ** 2);
}

subDistanceHubeny.cgi // Topに戻る // indexに戻る
Copyright(C) 2009-2021 Grandmaster Last up : 2020/09/18