地図投影座標変換 (XYCONV)

１．機　能
  ・xyconv／nxycnv ルーチンは，座標原点の緯度・経度と計算点の緯度・経度を与えて
    国土調査法に基づく平面直角座標（JTM座標と呼ぶ）の座標値を計算する． xyconv 
    は 2002年3月まで使用された旧東京測地系，nxycnv は 2002年4月以降の「日本測地
    系2000」に基づいた計算を行う
  ・ikconv／nikcnv ルーチンは，xyconv／nxycnv ルーチンの逆変換を行う．
  ・utm ルーチンは，緯度・経度を与えて指定のゾーン番号における UTM座標の値を計
    算する．旧東京測地系の場合 1〜60 の値で，日本測地系2000 の場合 200 を加算し
    た 201〜260 の値で，そのゾーン番号を指定する．
  ・ikutm ルーチンは，utm ルーチンの逆変換を行う．
 (注) これらはガウスクルィーゲル図法（横メルカトール図法）に基づく座標計算を行
    うものであり，地球楕円体モデルとして，旧東京測地系ではベッセル原子，日本測
    地系2000 では GRS80楕円体 を用いており，中央経線上での縮尺係数は，UTM座標で
    は 0.9996，JTM座標では 0.9999 である．JTM座標では，その原点で X＝Y＝0. であ
    るが，UTM座標では，赤道で X＝0.，中央経線で Y＝500.(km) である．また，元来
    の UTM座標では，南半球の X座標を 10000.(km) 加算した正の値で表すこととされ
    ているが，本ルーチンの計算では加算を行わず，負の値で表現している．

２．呼出し手順と引数

      call xyconv(gi, gk, fi, fk, x, y)
      call nxycnv(gi, gk, fi, fk, x, y)
      call utm(nc, fi, fk, x, y)
  gi, gk  [float] 座標原点の緯度・経度を分単位で与える．
  nc      [int]   UTM座標のゾーン番号（1〜60）または それに 200を加算したもの．
                  1〜60 であれば，東経 (nc*6-183)度を中央経線とし，旧東京測地系
                  で計算し，201〜260 であれば，日本測地系2000 で計算する．
  fi, fk  [float] 計算点の緯度・経度を分単位で与える．
  x,  y   [float] 計算結果．X座標(北向き)，Y座標(東向き)の値（km単位）を得る．

      call ikconv(gi, gk, x, y, fi, fk)
      call nikcnv(gi, gk, x, y, fi, fk)
      call utmik(nc, x, y, fi, fk)
  gi, gk  [float] 座標原点の緯度・経度を分単位で与える．
  nc      [int]   UTM座標のゾーン番号を示す．（1〜60／201〜260）
  x,  y   [float] 計算点の X座標値(北向き)，Y座標値(東向き)を km単位で与える．
  fi, fk  [float] 計算結果．緯度・経度の値（分単位）を得る．

１．機　能
    各種の地図投影法に基づく緯度・経度から直角座標への変換，およびその逆変換の
  計算を行う．cvinit／cvdinit ルーチンで，図法番号とそのパラメータを指定したの
  ち，cviken／cvdiken ルーチンで緯度・経度から直角座標への変換，もしくは cvenik
  ／cvdenik ルーチンで直角座標から緯度・経度への変換の計算を行う．
    地球モデルとしては，ベッセル原子の楕円体（長半径が 6377.397155 km で，偏平
  率が 1／299.152813；いずれも定義値）と GRS80楕円体（長半径[定義値]が 6378.137
  km で，偏平率[誘導定数]が 1／298.257222）を用いる．
    現在，用意されている図法／座標系とその番号は，ベッセル楕円体に対して，次の
  通りであり，GRS楕円体に対しては 200 を加算したものである．

  （図法番号）   （図法／座標系）
         0    国土調査法に基づく平面直角座標[JTM座標]（ガウスクルィーゲル図法）
      1-60    UTM座標（ガウスクルィーゲル図法）
        61    北極UPS座標（Universal Polar Stereographic Projection）
        62    南極UPS座標（   ＝ 国際極心平射図法                  ）
        65    UTM座標と同様であるが非標準の中央経線を指定できる
        70    メルカトール図法
        71    ランベルト等角円錐図法（1基本緯線）
        72    ランベルト等角円錐図法（2基本緯線）
       100    ランベルト等積方位図法（楕円体と表面積の等しい球からの投影）
       109    ランベルト等積方位図法（楕円体の長半径を半径とする球から投影）
       199    緯経度の分単位の値をそのまま km 単位の北・東向き座標値とみなす

  (注1) 球からの等積投影では，半径として，楕円体と表面積が等しい球の半径を採
      用すべきであり，図法番号 100 は，その半径を用いているが，図法番号 109 
      では，"Magnetic Anomaly Map of East Asia" の投影法に合わせて，約 0.11％
      大きい長半径を半径とした計算に対応させてある．
  (注2) 直角座標はすべて縮尺 1：1 の場合の km単位を採用しており，一般の図法で
      は，座標展開の原点で X＝Y＝0.であるが，UTM座標においては，X座標(北向き)
      が赤道で 0.，Y座標(東向き)が中央経線上で 500.(km) と決められており，UPS
      座標においては，極（原点）で X＝Y＝2000.(km) である．

２．呼出し手順と引数

      call cvinit(nc, gi, gk, sp1, sp2)
      call cvdinit(nc, dgi, dgk, dp1, dp2)
  nc        [int]     図法番号（前記）を与える．
  gi,  gk   [float]   座標原点（X＝Y＝0.）の緯度・経度を分単位で与える．
                       （UTM座標／UPS座標 の場合は参照されない．）
  dgi, dgk  [double]   [同上]
  sp1, sp2  [float]   2基本緯線タイプのランベルト等角円錐図法の場合にその2基本
                      緯線の緯度を分単位で与える．
                       （他の図法の場合には参照されない．）
  dp1, dp2  [double]   [同上]

      call cviken(fi, fk, ye, xn)
      call cvdiken(dfi, dfk, dye, dxn)
  fi,  fk   [float]   計算点の緯度・経度を分単位で与える．
  dfi, dfk  [double]   [同上]
  ye,  xn   [float]   計算結果．東向き座標値，北向き座標値（km単位）を得る．
  dye, dxn  [double]   [同上]

      call cvenik(ye, xn, fi, fk)
      call cvdenik(dye, dxn, dfi, dfk)
  ye,  xn   [float]   計算点の東向き座標値，北向き座標値をkm単位で与える．
  dye, dxn  [double]   [同上]
  fi,  fk   [float]   計算結果．緯度・経度の値（分単位）を得る．
  dfi, dfk  [double]   [同上]

    なお，cviken／cvenik／cvdiken／cvdenik ルーチンにおいては，直角座標値引数の
  配置順序が，utm／utmik／xyconv／ikconv／nxycnv／nikcnv ルーチンの場合と異なっ
  ているので，注意されたい．