位置計算装置、位置計算方法、位置計算プログラム、及び位置計算プログラムの記憶媒体
【課題】 平面直角座標上の座標ｘ，ｙから、地球楕円体上の緯度、経度を精度よく算出する位置計算装置を提供する。
【解決手段】 演算部６が、赤道から当該座標値ｘまでの子午線弧長Ｗを算出し、子午線曲率半径を以って緯度φj にて赤道から垂線緯度φj まで積分することにより子午線弧長Ｓ（φj ）を算出する。演算部６は、緯度φj を与えることによって
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）を計算して緯度φj+1 を算出する。緯度φj+1 と緯度φj との差の絶対値が基準値以内の場合、
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 式３を計算し、求めた緯度φk+1 を緯度φn とし、演算部６に当該緯度φn における縮尺係数ｍt+1 を算出させ、当該縮尺係数ｍt+1 を用いて求めた、子午線弧長Ｗにて、上記式１の計算を行わせる。
【解決手段】 演算部６が、赤道から当該座標値ｘまでの子午線弧長Ｗを算出し、子午線曲率半径を以って緯度φj にて赤道から垂線緯度φj まで積分することにより子午線弧長Ｓ（φj ）を算出する。演算部６は、緯度φj を与えることによって
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）を計算して緯度φj+1 を算出する。緯度φj+1 と緯度φj との差の絶対値が基準値以内の場合、
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 式３を計算し、求めた緯度φk+1 を緯度φn とし、演算部６に当該緯度φn における縮尺係数ｍt+1 を算出させ、当該縮尺係数ｍt+1 を用いて求めた、子午線弧長Ｗにて、上記式１の計算を行わせる。
【発明の詳細な説明】
【技術分野】
【０００１】
本願発明は、位置計算装置、位置計算方法、位置計算プログラム、及び位置計算プログラムの記憶媒体に関するものである。
【背景技術】
【０００２】
地球は、楕円体（地球楕円体）であり、地理学的経緯度を以って、地球上の位置を特定することができる。しかし、取り扱いの便から、一般的な測量等において、平面直角座標が利用され、座標によって、位置の特定がされる（測量法第11条）。この座標平面については、回転楕円体（地球楕円体）から、直接平面に投影するガウス・クリューゲルの等角投影（横メルカトール投影）を採用して、原点における縮尺係数を0.99990 としている。
尚、地球楕円体を測量の基準にするためには、楕円体の中心を実際の地球上のどの位置に、またその楕円体の座標軸が実際の地球のどこを通るかということを決める必要があり、この位置と方向が決められた地球楕円体が準拠楕円体とされる。
【０００３】
地球上の点の水平位置は、厳密には準拠楕円体上の地理学的経緯度によって表されるべきであるが、位置・方向・距離等を、上記の通り平面上に投影して測量計算を行うことは曲面上に比べ非常に簡単になり便利である。また、公共測量のように測量範囲が狭い場合には、十分正確に表すことができるとされている。
日本で用いられている平面直角座標は、ガウス・クリューゲルの等角投影法によるもので、座標系原点を通る子午線は等長に、図形は等角の相似形に投影される。しかし、距離については、基準子午線から東西に離れるに従って平面距離が増大していくため、投影距離の誤差を相対的に1/10,000以内に収めるよう座標原点に上記の縮尺係数（0.9999）を与え、かつ、基準子午線より東西130km 以内を適用範囲とした座標系（平面直角座標系は日本全国を１９の座標系に区分している。）を設けている（非特許文献１）。
なお、座標系のＸ軸は、原点において子午線に一致する軸とし、原点から真北に向かう値を正としている。Ｙ軸は原点においてＸ軸に直交する軸とし、真東に向かう値を正としている。この点、三次元直交座標（Ｘ，Ｙ, Ｚ）とは定義が異なる。又、各座標系原点の値はＸ＝Ｙ＝0.000 メートルとしている。
【０００４】
測量結果の成果表において、距離については楕円体面上の距離が記載され、座標については緯度・経度のほか、平面直角座標（ｘ，ｙ）が記載される。楕円体面と平面の関係については、準拠楕円体に被せられた楕円体筒（仮想筒）を広げて平面座標（ｘ，ｙ）を得る。この場合、基準子午線における平面距離は、楕円体上の距離より１万分の１だけ短くなり、縮尺係数は、0.9999である。基準子午線から東西に約90ｋｍ離れた場所では、平面と楕円体面の距離が等しくなり、縮尺係数は1.0000である。基準子午線から東西に130 ｋｍ離れると、平面距離は、楕円体面距離より、１万分の１だけ長くなり、縮尺係数は、1.0001である。
【０００５】
わが国土について、平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙから、その緯度・経度を算出する手段として、国土地理院のホームページにおいて、「緯度・経度への換算」（非特許文献３）のプログラムが利用できる。
このプログラムを利用することにより、ｗｅｂ上において、上記座標値ｘ，ｙの入力により、簡便に、その緯度・経度を算出することができる。
ところが、このプログラムにおいて、座標値ｘ，ｙが属する座標系の基準子午線から、東西に離れるに従って、その誤差が大きくなる。
この誤差は、国土地理院ホームページの当該「緯度・経度への換算」プログラムで得た緯度・経度を、同じく国土地理院ホームページ「平面直角座標への換算」（非特許文献４）プログラムにて、簡単に確認することができる。
【０００６】
例えば、国土地理院ホームページ「緯度・経度への換算」（非特許文献３）において、測地系を世界測地系とし、平面直角座標系を３として、上記「平面直角座標から緯度・経度への換算」プログラムへ、平面直角座標上の任意の５点の座標値Ｘ，Ｙを与えて、ａ−１）〜ａ−５）に示す結果（緯度・経度）を得た。
ａ−１）は当該座標系３の座標原点（０，０）を通る子午線上の点であり、ａ−２）〜ａ−５）において、入力する座標値Ｙの値を基準子午線から西方へ徐々に遠ざかるものとしている。
そして上記「平面直角座標への換算」プログラムにおいて、上記ａ−１）〜ａ−５）で入力した座標値と、換算結果が一致するように、試行錯誤で北緯と東経を入力した（北緯と東経の値を種々変えて、出力値が、「平面直角座標への換算」プログラム）における入力値と一致するものを探った）。そして、ｂ−１）〜ｂ−５）に示す結果（平面直角座標上の座標値）を得た。
【０００７】
ａ−１）
入力値 Ｘ座標 -207504.143 ｍ Ｙ座標 0 ｍ
換算結果 緯度 34°07'45.93167" 経度 132°10'00.00000"
真北方向角＋ 0°00'00.00" 縮尺係数 0.99990000
ａ−２）
入力値 Ｘ座標 -207494.938 ｍ Ｙ座標 -59997.338 ｍ
換算結果 緯度 34°07'40.02853" 経度 131°30'58.32911"
真北方向角＋ 0°21'53.81" 縮尺係数 0.99994436
ａ−３）
入力値 Ｘ座標 -207483.393 ｍ Ｙ座標 -90082.603 ｍ
換算結果 緯度 34°07'32.62476" 経度 131°11'24.23439"
真北方向角＋ 0°32'52.52" 縮尺係数 1.00000000
ａ−４）
入力値 Ｘ座標 -207467.727 ｍ Ｙ座標 -119340.074 ｍ
換算結果 緯度 34°07'22.57910" 経度 130°52'22.58340"
真北方向角＋ 0°43'32.98" 縮尺係数 1.00007551
ａ−５）
入力値 Ｘ座標 -207462.646 ｍ Ｙ座標 -127394.933 ｍ
換算結果 緯度 34°07'19.32100" 経度 130°47'08.30605"
真北方向角＋ 0°46'29.28" 縮尺係数 1.00010000
【０００８】
ｂ−１）
入力値 北緯 34°07'45.93167" 東経 132°10'00.00000"
換算結果 Ｘ座標 -207504.143 ｍ Ｙ座標 0.0000 ｍ
真北方向角＋ 0°00'00.00" 縮尺係数 0.99990000
ｂ−２）
入力値 北緯 34°07'40.02854" 東経 131°30'58.32910"
換算結果 Ｘ座標 -207494.938 ｍ Ｙ座標 -59997.338ｍ
真北方向角＋ 0°21'53.81" 縮尺係数 0.99994436
ｂ−３）
入力値 北緯 34°07'32.62478" 東経 131°11'24.23440"
換算結果 Ｘ座標 -207483.393 ｍ Ｙ座標 -90082.603ｍ
真北方向角＋ 0°32'52.52" 縮尺係数 1.00000000
ｂ−４）
入力値 北緯 34°07'22.57909" 東経 130°52'22.58341"
換算結果 Ｘ座標 -207467.727 ｍ Ｙ座標 -119340.074 ｍ
真北方向角＋ 0°43'32.98" 縮尺係数 1.00007551
ｂ−５）
入力値 北緯 34°07'19.32099" 東経 130°47'08.30603"
換算結果 Ｘ座標 -207462.646 ｍ Ｙ座標 -127394.933 ｍ
真北方向角＋ 0°46'29.28" 縮尺係数 1.00010000
【０００９】
上記の結果を見ると、ａ−１）の入力値と一致する、出力値が得られるようにｂ−１）において探って得た入力値は、ａ−１）における換算結果と一致している（上記ａ−１）における換算結果の緯度34°07'45.93167"と経度 132°10'00.00000"は、ｂ−１）における入力値の北緯・東経と一致している）。即ち、当該原点を通る基準子午線上において、「平面直角座標への換算」プログラムは、少なくとも計算可能な桁数内において、誤差を生じていないといえる。
ところが、ａ−２）について、上記と同様の作業を行ったｂ−２）を見ると、ａ−２）における換算結果と、ｂ−２）における入力値とは一致していない。即ち、前者ａ−２）が緯度を34°07'40.02853"とし経度を 131°30'58.32911"とするのに対して、後者ｂ−２）が北緯を34°07'40.02854"とし東経を 131°30'58.32910"とし、緯度及び経度について、夫々小数点以下第５位（最下位）において食い違いが生じている。
このような誤差の発生について、ａ−３〜ａ−５）と対応するｂ−３）〜ｂ−５）と見れば、基準子午線から西方へ離れるに従って、特に経度についての当該誤差が大きくなっていることが分かる。
例えば、基準子午線から西方へ約１２７ｋｍ離れた地点であるａ−５）における換算結果（緯度34°07'19.32100"、経度 130°47'08.30605"）と、対応するｂ−５）における入力値（北緯34°07'19.32099"、東経 130°47'08.30603"）とでは、特に、緯度について小数点以下第３位にて既に食い違いが生じている。
【００１０】
【非特許文献１】国土地理院ホームページ http://vldb.gsi.go.jp/sokuchi/datum/tokyodatum.html
【非特許文献２】国土地理院ホームページ「赤道からの子午線弧長を与えて緯度を求める計算」http://vldb.gsi.go.jp/sokuchi/surveycalc/algorithm/s2b/s2b.htm
【非特許文献３】国土地理院ホームページ「緯度・経度への換算」http://vldb.gsi.go.jp/sokuchi/surveycalc/xy2blf.html
【非特許文献４】国土地理院ホームページ「平面直角座標への換算」http://vldb.gsi.go.jp/sokuchi/surveycalc/bl2xyf.html
【００１１】
従って、このプログラムでは、入力座標値Ｘ，Ｙが属する座標系原点の基準子午線付近においては、問題がないものの、基準子午線から東西へ離れるに従って、換算結果である緯度及び経度の精度が著しく低下し、厳密なデータとして取り扱うことができなくなる。
上記において、そのような誤差が生じる原因は、プログラムの内容が公開されていないことから、本願発明者において正確なことは分からない。
また、発明者にて可能な範囲で、先行技術文献として、特許公報を調べたが、同種の発明が記載された公報を発見するに至っていない。
そこで、発明者は、この点を熟考し、推測するに、上記国土地理院のプログラムにおいて、緯度について例えば、国土地理院ホームページ「赤道からの子午線弧長を与えて緯度を求める計算」（非特許文献２）などに示される公式を用いて、緯度についての収束計算を行っているのであろうが、基準子午線から離れた場合の前述の縮尺係数の変化が十分に反映するものとなっておらず、その結果、収束（反復計算）により安定状態となった緯度の精度が十分なものでない（ずれた位置に収束してしまっている）と考える。
具体的には、国土地理院の上記プログラムは、非特許文献２に示す国土地理院ホームページの「６．赤道からの子午線弧長を与えて緯度を求める計算」で示されているφn+1 を求める式を用いて反復計算を行うものであると考えられる（この式を用いた結果が上記と一致した）。この式を用いて反復計算を行うと上述の通り、大きな誤差を持った値で数値が安定（収束）する。その理由を検討すると、この式の計算において、縮尺係数は、原点における縮尺係数ｍ0 （＝0.9999）に固定されており、基準子午線以外の、求めようとする座標位置における縮尺係数は、原点の縮尺係数ｍ0 （＝0.9999）と異なるものであるという事実は、無視されている（尚、国土地理院の式で用いているＭは、赤道からの座標系原点までの子午線弧長を示し、この明細書で用いる子午線曲率半径を示すＭとは異なる）。
現状では、推測の域を出るものではないが、上記プログラムを用いた結果は、歴然としている。
従って、当該プログラムは、著しく信頼性の欠くものとなっており、厳密性が要求される状況において、その利用が困難となっている。
【発明の開示】
【発明が解決しようとする課題】
【００１２】
本願発明は、等角投影法により地球楕円体を投影する平面直角座標Ｘ，Ｙ上の座標値から地球楕円体上の少なくとも緯度を求めることが可能な位置計算装置について、縮尺係数の取り扱いを柔軟に行うことにて、誤差の低い新規な装置を提供することにより、上記課題の解決を図った。
【課題を解決するための手段】
【００１３】
本願第１の発明は、等角投影法により地球楕円体を投影する平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙを入力することにより、少なくとも当該座標値に対応する地球楕円体上の緯度φn を出力することが可能な位置計算装置として、次の構成を備えたものを提供する。
即ち、この位置計算装置は、座標入力部と、出力部と、演算部と、判定部と、回帰制御部とを備える。座標入力部は、上記座標値ｘ，ｙの入力を受け付ける。演算部は、直角座標の座標値ｘが与えられることにより、地球楕円体の赤道から当該座標値ｘ，ｙが属する座標系原点までの子午線弧長に対し、当該座標値ｘを仮の縮尺係数ｍt で除した値を、加算して、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとする。
そして、この子午線弧長Ｗに対応する基準子午線上の緯度、即ち、上記入力座標から基準子午線上に下ろした垂線の緯度を算出するため、演算部は、地球楕円体の長半径と第１離心率と緯度とから定まる子午線曲率半径について、地球楕円体の長半径ａと第１離心率ｅと共に任意の緯度φj を初期値として与え、別途に子午線曲率半径を以って赤道から上記緯度φj まで積分することにより、赤道からの基準子午線上の子午線弧長Ｓ（φj ）を算出する。更に、演算部は、基準子午線上の上記緯度φj を与えることによって、式１を計算して基準子午線上の緯度φj+1 を漸近値として算出する。
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）…式１
判定部は、上記式１にて算出した基準子午線上の上記漸近緯度φj+1 と基準子午線上の緯度φj との差の絶対値と、基準値との比較を行うものである。判定部にて、基準子午線上の、上記漸近緯度φj+1 と緯度φj との差の絶対値が基準値を超えると判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするには十分収束していないと判定した場合、回帰制御部は、この漸近緯度φj+1 を上記基準子午線上の緯度φj とし、演算部に基準子午線上の赤道からの子午線弧長Ｓ（φj ）を算出させると共に当該子午線弧長Ｓ（φj ）から上記式１の計算にて再度漸近緯度φj+1 を算出させ、判定部に再度上記比較を行わせる。また、判定部が、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするのに十分収束したと判定した場合、回帰制御部は、当該漸近緯度φj+1 を、入力座標から基準子午線上に下ろした垂線の緯度φk として、仮の縮尺係数ｍt と、当該垂線緯度φk と、当該垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式３を計算して緯度φk+1 を演算部に算出させる。
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …式３
回帰制御部は、上記緯度φk+1 を入力座標に対応する緯度φn として、演算部に、縮尺係数ｍt+1 を算出させ、前記仮の縮尺係数ｍt に代え当該縮尺係数ｍt+1 を用いて求めた、基準子午線上の子午線弧長Ｗにて、上記式１の計算を行わせる。
ここでいう地球楕円体には、各種の準拠楕円体を含む。
基準子午線とは、当該座標が属する座標系の原点を通る子午線をいう。
また、垂線緯度とは、入力座標から基準子午線に下ろした垂線と、当該基準子午線との交点の緯度をいう。
また、積分には、その近似計算を含む。
上記において、回帰制御部は、演算部に対して直接、計算の指令を下すものに限定するものではなく、間接的に、演算部へ計算を行わせるものも含む。
【００１４】
本願第２の発明は、上記本願第１の発明にあって、次の特徴を備えた位置計算装置を提供する。即ち、上記の仮の縮尺係数ｍt とは、演算部に、直角座標の座標値ｙを与えることにより、原点における縮尺係数ｍ0 又はその代替値と、当該座標が属する座標系の原点における平均曲率半径Ｒ0 とから、式２の計算により、演算部が算出した縮尺係数ｍt である。
ｍt ＝｛１＋（ｙ2 ／２ｍ02Ｒ02）＋（ｙ4 ／24ｍ04Ｒ04）｝×0.9999…式２
上記の原点における平均曲率半径Ｒ0 は、当該座標が属する座標系の原点の緯度φ0 又はその代替値と共に、地球楕円体の長半径ａ及び短半径ｂ、或いは当該半径ａ，ｂから導かれる第１離心率ｅ、第２離心率ｅ’を用いて定めたものである。
本願第３の発明は、上記本願第２の発明にあって、次の特徴を備えた位置計算装置を提供する。即ち、上記縮尺係数ｍt+1 は、原点縮尺係数ｍ0 に代え上記縮尺係数ｍt を用い、且つ、原点における平均曲率半径Ｒ0 に代え演算部に式３にて算出させた緯度φn における平均曲率半径Ｒφn を用いて、演算部が、上記式２により算出するものである。
【００１５】
本願第４の発明は、上記本願第１乃至３の何れかの発明にあって、次の特徴を備えた位置計算装置を提供する。即ち、上記にあって、判定部が、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合に、回帰制御部は、上記緯度φn に対応する縮尺係数ｍt+1 と上記仮の縮尺係数ｍt との差の絶対値と、所定値との比較を、更に判定部にさせ、判定部が当該縮尺係数ｍt+1 と上記仮の縮尺係数ｍt との差の絶対値が所定値を超えると判定した場合にのみ、上記の通り、前記の仮の縮尺係数ｍt に代え、当該縮尺係数ｍt+1 を用いて求めた基準子午線上の子午線弧長Ｗにて、上記式１の演算部に計算を行わせるものである。
【００１６】
本願第５の発明は、上記本願第１乃至４の何れかの発明にあって、次の特徴を備えた位置計算装置を提供する。即ち、上記にあって、判定部が、基準子午線上の、緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、回帰制御部は、垂線緯度φk として上記漸近緯度φj+1 を演算部に与えるものであり、演算部は、上記仮の縮尺係数ｍt と、上記垂線緯度φk における卯酉線曲率半径Ｎk と、当該垂線緯度φk と第２離心率とによりｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式４にて、経度差分Δλを算出することが可能であり、
Δλ＝（ｙ／ｍt ）／Ｎk cos φk
−（ｙ／ｍt ）3 （１＋２ｔk 2 ＋ηk 2 ）／６Ｎk 3cosφk
＋（ｙ／ｍt ）5 （５＋２８ｔk 2 ＋２４ｔk 4 ）／１２０Ｎk 5cosφk …式４
演算部は、当該座標系が属する原点における経度λ0 に対し当該経度差分Δλを加算することにより、上記座標値ｘ，ｙにおける経度λを算出するものである。
【００１７】
本願第６の発明は、上記本願第１乃至５の何れかの発明にあって、次の特徴を備えた位置計算装置を提供する。即ち、上記にあって、判定部が、基準子午線上の、緯度φj+1 と緯度φj との差の絶対値を基準値以内と判定した場合、回帰制御部は、垂線緯度φk として緯度φj+1 を演算部に与えるものであり、演算部は、上記の縮尺係数ｍt と、当該緯度φk における、子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、当該垂線緯度φk と第２離心率とによりｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式５にて、座標ｘ，ｙにおける子午線収差γを算出するものである。
γ＝（ｙ／ｍt ）ｔk ／Ｎk
−（ｙ／ｍt ）3 ｔk （１＋ｔk 2 −ηk 2 ）／３Ｎk 3
＋（ｙ／ｍt ）5 ｔk （２＋５ｔk 2 ＋３ｔk 4 ）／１５Ｎk 5
…式５
【００１８】
本願第７の発明は、等角投影法により地球楕円体を投影する平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙを入力することにより、少なくとも当該座標値に対応する地球楕円体上の緯度φn を出力することが可能な位置計算プログラムについて、次の特徴を備えたものを提供する。
即ち、この位置計算プログラムは、座標入力部と、出力部と、演算部と、判定部と、回帰制御部とを備えたコンピュータに導入される。このプログラムが導入されたコンピュータは、座標入力部にて、上記座標値ｘ，ｙの入力を受け付け、演算部に、直角座標の座標値ｘを与えて、地球楕円体の赤道から当該座標値ｘ，ｙが属する座標系原点までの子午線弧長に対し、当該座標値ｘを仮の縮尺係数ｍt で除した値を加算して、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとする。
そして、この子午線弧長Ｗに対応する基準子午線上の緯度、即ち、上記入力座標値から基準子午線上に下ろした垂線の緯度を算出するため、演算部にて、地球楕円体の長半径と第１離心率と緯度とから定まる子午線曲率半径について、地球楕円体の長半径ａと第１離心率ｅと共に任意の緯度φj を初期値として与え、別途に子午線曲率半径を以って赤道から上記緯度φj まで積分することにより、赤道からの基準子午線上の子午線弧長Ｓ（φj ）を算出し、更に、演算部に、基準子午線上の緯度φj を与えることによって、式１を計算して基準子午線上の漸近緯度φj+1 を漸近値として算出させる。
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）…式１
そして、判定部に、上記式１にて算出した基準子午線上の漸近緯度φj+1 と基準子午線上の緯度φj との差の絶対値と、基準値との比較を行わせる。判定部にて、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値を超えると判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするには十分収束していないと判定した場合、回帰制御部により、漸近緯度φj+1 を上記基準子午線上の緯度φj とし、演算部に基準子午線上の赤道からの子午線弧長Ｓ（φj ）を算出させると共に当該子午線弧長Ｓ（φj ）から上記の式１の計算にて再度漸近緯度φj+1 を算出させ、判定部に再度上記比較を行わせるものである。また、判定部にて、基準子午線上の漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線とするのに十分収束したと判定した場合、当該漸近緯度φj+1 を入力座標からおろした垂線の緯度φk として、回帰制御部により、仮の縮尺係数ｍt と、当該垂線緯度φk と、当該垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式３を計算して緯度φk+1 を演算部に算出させる。
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …式３
更に、回帰制御部により、上記緯度φk+1 を入力座標に対応する緯度φn として、演算部に、縮尺係数ｍt+1 を算出させ、前記仮の縮尺係数ｍt に代え当該縮尺係数ｍt+1 を用いて求めた、基準子午線上の子午線弧長Ｗにて、上記式１の計算を行わせるものである。
【００１９】
本願第８の発明は、上記第７の発明に係る位置計算プログラムの記録媒体を提供する。
【００２０】
本願第９の発明は、等角投影法により地球楕円体を投影する平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙを入力することにより、少なくとも当該座標値に対応する地球楕円体上の緯度φn を出力することが可能な位置計算方法について、次の特徴を備えた方法を提供する。
即ち、この方法は、座標入力部と、出力部と、演算部と、判定部と、回帰制御部とを備えたコンピュータを用いるものである。そして、座標入力部にて、上記座標値ｘ，ｙの入力を受け付ける。演算部に、直角座標の座標値ｘを与えて、地球楕円体の赤道から当該座標値ｘ，ｙが属する座標系原点までの子午線弧長に対し、当該座標値ｘを仮の縮尺係数ｍt で除した値を加算して、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとし、演算部にて、地球楕円体の長半径と第１離心率と緯度とから定まる子午線曲率半径について、地球楕円体の長半径ａと第１離心率ｅと共に任意の緯度φj を初期値として与え、別途に子午線曲率半径を以って赤道から上記緯度φj まで積分することにより、赤道からの基準子午線上の子午線弧長Ｓ（φj ）を算出し、演算部に、基準子午線上の緯度φj を与えることによって、式１を計算して基準子午線上の漸近緯度φj+1 漸近値として算出させる。
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）…式１
判定部にて、上記式１にて算出した基準子午線上の漸近緯度φj+1 と基準子午線上の緯度φj との差の絶対値と、基準値との比較を行わせるものである。そして、判定部にて、上記基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値を超えると判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線に下ろした垂線の緯度とするには十分収束していないと判定した場合、回帰制御部により、漸近緯度φj+1 を上記基準子午線上の緯度φj とし、演算部に基準子午線上の赤道からの子午線弧長Ｓ（φj ）を算出させると共に当該子午線弧長Ｓ（φj ）から上記の式１の計算にて再度漸近緯度φj+1 を算出させ、判定部に再度上記比較を行わせる。判定部にて、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線に下ろした垂線の緯度とするのに十分収束したと判定した場合、回帰制御部は、当該漸近緯度φj+1 を、入力座標から基準子午線上に下ろした垂線の緯度φk として、仮の縮尺係数ｍt と、上記垂線緯度φk と、当該垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式３を計算して緯度φk+1 を演算部に算出させる。
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …式３
そして、回帰制御部により、上記緯度φk+1 を入力座標に対応する緯度φn として、演算部に、縮尺係数ｍt+1 を算出させ、前記仮の縮尺係数ｍt に代え当該縮尺係数ｍt+1 を用いて求めた、基準子午線上の子午線弧長Ｗにて、上記式１の計算を行わせるものである。
【発明の効果】
【００２１】
上記本願の各発明は、直角座標上の座標点の入力により、自動的に、少なくとも当該緯度を算出するものであって、当該緯度計算中、必要とする垂線緯度の反復計算（収束）において、緯度に対応して縮尺係数の収束も反映させることが可能な装置・方法・プログラムやその記録媒体を提供し得た。これによって、極めて精度の高い、緯度計算を行うことができる。
即ち、この発明では、入力座標から基準子午線上に下ろした垂線の緯度（垂線緯度）を算出するに当たり、緯度の収束計算に適した式１をイテレーション式として用いて緯度の計算（反復計算）を行い、得た垂線緯度から式３を用いて当該点（入力座標値に対応する）緯度を算出し直すことにより、より精度の高い緯度を得ることができる。更に、このようにして求めた緯度から縮尺係数を算出して当該縮尺係数の収束を加味した後処理を行うことにより、極めて精度の高い当該点の緯度の算出を可能とした。
また、特に、本願第４の発明では、上記の計算を可能としつつも、基準子午線に近い場合など十分に誤差の小さい場合にまで、後処理として、式３を経て得た上記の緯度で算出した縮尺係数による、式１の計算を強要するものではなく、縮尺係数の収束状況で、後処理を行うか否かを決定し、無駄な反復計算を抑制し、迅速に計算を完了する手段を備えた装置を提供し得た。
更に本願第５の発明では、上記の当該点の緯度と共に、経度λを、また本願第６の発明では、子午線収差を、夫々算出し得る位置計算装置を提供した。
【発明を実施するための最良の形態】
【００２２】
以下、本願発明の好ましい実施の形態について、説明する。
【００２３】
１．本願発明に係る計算方法の概要
本願発明は、平面直角座標上の座標ｘ，ｙから、対応する地球楕円体上の緯度、更には、経度、子午線収差を精度よく算出する位置計算方法を提供するものである。即ち、ガウス・クリューゲルの等角投影法により地球楕円体を投影する平面直角座標の座標点の座標値から、当該地球楕円上の対応する点の緯度、経度及び子午線収差を、高精度に算出する方法を提供するものである。
この方法は、垂線緯度の算出のための収束計算において、縮尺係数の収束を加味して、計算結果の精度を向上させたものである。
【００２４】
即ち、この方法は、主工程において、下記の式１にて垂線緯度の反復計算を行い、反復計算後の値から下記の式３を計算して緯度を算出する。そして、式３を経て得た値を用いて下記の式２にて縮尺係数を算出し、算出した縮尺係数から導いた基準子午線上の子午線弧長を用いて再度式１にて垂線緯度の反復計算を、更に、再度の式１による反復計算後の値にて式３の計算を、後処理工程として行うものである。
また、上記後処理工程に先立ち、算出した縮尺係数が収束したかを調べ、その結果当該縮尺係数が収束していた場合、後処理工程を行わずに、計算を終了することもできる。
従って、主工程の処理に続いて後処理工程を行い上記にて算出した縮尺係数を加味した結果を得るものとしてもよく、また、主工程の結果に対して、算出した縮尺係数の収束の判定を行い、当該縮尺係数が収束していない場合にのみ、上記の後処理工程を行うものとしてもよい。
このように、主工程後必ず後処理工程を行う場合も、主工程後の縮尺係数の収束の判定により後処理工程を行う場合も、何れも、垂線緯度の収束計算において、縮尺係数の収束を加味したものとなり、緯度の計算結果について必要な精度を確保することができる。尚、縮尺係数の収束を調べて後処理工程を行うか否かの判断を行う場合、無駄な反復計算（後処理工程）を行う必要がない。例えば、基準子午線付近では、主工程のみにて緯度及び経度の十分な精度を確保している場合がある。このような場合にまで後処理工程を行う必要がないからである。
但し、近年のコンピュータの性能の向上により、主工程後、必ず後処理工程を行うものとしても、実用上、計算に要する時間の増加や、ハードウエアの負荷による問題はない。
また、後処理工程は、１回に限らない。更に、上記と同様、後処理工程を行う毎に縮尺係数の収束の判定を行い、その結果により更に後処理工程を行うか否か判断するものとしても実施可能である。
【００２５】
以下、主工程後の縮尺係数の収束の判定により後処理工程を行うものについてその概要を説明する。
この位置計算方法は、上記座標値ｘ，ｙの入力を受け付ける。そして当該座標値ｘが与えられることにより、地球楕円体の赤道から当該座標値ｘ，ｙが属する座標系原点までの子午線弧長に対し、当該座標値ｘを０以外の仮の縮尺係数ｍt で除した値を、加算して、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとする。また、地球楕円体の長半径と第１離心率と緯度とから定まる子午線曲率半径について、地球楕円体の長半径ａと第１離心率ｅと共に任意の数値（緯度φj ）を初期値として与え、別途に子午線曲率半径を以って赤道から上記緯度φj まで積分（近似計算）することにより、赤道からの基準子午線上の子午線弧長Ｓ（φj ）を算出する。そして、基準子午線上の上記漸近緯度φj から、式１を計算して基準子午線上の緯度φj+1 を漸近値（漸近緯度φj+1 ）として算出する。
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）…式１
上記式１にて算出した基準子午線上の緯度φj+1 と基準子午線上の緯度φj との差の絶対値と、基準値との比較を行い、上記基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値を超えると判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするには十分収束していないと判定した場合、この漸近緯度φj+1 を上記基準子午線上の緯度φj とし、基準子午線上の赤道からの子午線弧長Ｓ（φj ）を算出すると共に当該子午線弧長Ｓ（φj ）から上記式１の計算にて再度漸近緯度φj+1 を算出し、再度上記比較を行う。その結果、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、即ち、この漸近緯度φj+1 （以下必要に応じて最近似垂線緯度と呼ぶ。）を、上記入力座標から基準子午線上に下ろした垂線とするのに十分収束したと判定した場合、当該漸近緯度φj+1 を、入力座標から基準子午線上に下ろした垂線の緯度φk として、仮の縮尺係数ｍt と、当該垂線緯度φk と、当該垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式３を計算して緯度φk+1 を算出する。
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …式３
上記にて求めた緯度φk+1 を、当該（入力座標値に対応する）緯度φn として縮尺係数ｍt+1 を算出し、縮尺係数ｍt+1 と上記仮の縮尺係数ｍt との差の絶対値と、所定値との比較を行い（以上主工程）、当該縮尺係数ｍt+1 と上記仮の縮尺係数ｍt との差の絶対値が所定値を超えると判定した場合に、当該縮尺係数ｍt+1 を用いて求めた、基準子午線上の子午線弧長Ｗにて、上記式１の計算をやり直す（後処理工程）。
その結果得られた緯度を計算結果として出力する。また、緯度の計算の過程にて、経度及び子午線収差を計算する。
上記の式１をイテレーション式（反復式）として行う反復計算は、図３に示す平面座標において、基準子午線（Ｘ軸）上の漸近緯度φj+1 を、入力座標点（ｘ，ｙ）から基準子午線（Ｘ軸）に下ろした垂線ｈと、基準子午線（Ｘ軸）との交点Ｇに近似させる（近づける）計算である。基準子午線上において縮尺係数は0.9999である。式１による反復計算により、漸近緯度φj+1 を上記交点Ｇと取り扱うに十分近接した場合、当該漸近緯度φj+1 を最近似垂線緯度とし、式３の計算により、座標値（ｘ，ｙ）に対応する緯度φk+1 を算出するのである。
また、上記の仮の縮尺係数とは、上記直角座標の座標値ｙと、原点における縮尺係数ｍ0 又はその代替値と、当該座標が属する座標系の原点における平均曲率半径Ｒ0 とから、式２の計算により、算出した縮尺係数ｍt である。
ｍt ＝｛１＋（ｙ2 ／２ｍ02Ｒ02）＋（ｙ4 ／24ｍ04Ｒ04）｝×0.9999…式２
上記の原点における平均曲率半径Ｒ0 は、当該座標が属する座標系の原点の緯度φ0 又はその代替値と共に、地球楕円体の長半径ａ及び短半径ｂ、或いは当該半径ａ，ｂから導かれる第１離心率ｅ、第２離心率ｅ’を用いて定めたものである。
上記の仮の縮尺係数には、式２にて算出するものに限定するものではなく、０以外の任意の数値を採用することが可能である。但し、無駄な後処理工程を減ずる観点から、上記式２を用いて算出したものを採用するのが好ましい。
また上記において、式２を用いる場合も、原点における縮尺係数ｍ0 の代替値として０以外の任意の数値を用いることが可能であるが、上記と同様の観点から、このような代替値ではなく、当該原点における縮尺係数ｍ0 を用いて計算を行うのが好ましい。
また、上記の上記縮尺係数ｍt+1 は、原点縮尺係数ｍ0 に代え上記縮尺係数ｍt を用い、且つ、原点における平均曲率半径Ｒ0 に代え式３にて算出した緯度φn における平均曲率半径Ｒφn （後述のＲk2）を用いて、上記式２により算出するものである（言い換えると、式２において、ｍt をｍ0 に代入し、Ｒφn をＲ0 に代入して計算する）。
【００２６】
上記の通り、主工程は、式１にて垂線緯度を反復収束計算させ、その後式３、式４（式１９）、式５、及び式２を経ることにより、緯度、経度、子午線収差角及び縮尺係数を算出する工程である。垂線緯度の反復収束計算は、上記の式１にて、他の値は変更せずに、漸近緯度φj+1 を緯度φj とし直して新たに漸近緯度φj+1 を求め、当該計算を｜φj+1 −φj ｜が収束するまで行う。）。そして、主工程で｜φj+1 −φj ｜が収束した際、即ち当該漸近緯度φj+1 が最近似垂線緯度となった際、漸近緯度φj+1 を垂線緯度φk として用い、式３を計算する。式３の計算にて得られた値を緯度φk +1をφn として、上記の通り当該緯度φn における平均曲率半径から縮尺係数ｍt+1 を求める。
上記にて算出した縮尺係数ｍt+1 と、それまで各値を算出するのに用いた縮尺係数ｍt とについて、｜ｍt+1 −ｍt ｜が収束しているか判定する。その結果、収束していれば、計算を終了し、式３、式４（式１９）、式５及び式２を経て得た上記の緯度φn 、経度λ、子午線収差角γ、及び縮尺係数ｍt+1 を、目的の数値として出力する。
上記において｜ｍt+1 −ｍt ｜が収束していないと判定した場合、後処理工程を行う。後処理工程において、上記主工程で得た縮尺係数ｍt+1 を用い、他の値は、（主工程で用いたものと）変更することなく、主工程と同様の処理を行う。当該後処理工程において｜ｍt+1 −ｍt ｜が収束していなかった場合、当該後処理工程で得た縮尺係数ｍt+1 を用いて、更なる後処理工程を行う。
要するに、式２において、上記縮尺係数ｍt を原点縮尺係数ｍ0 の代わりに用いて、式３にて計算した緯度φn における平均曲率半径Ｒφn （後述のＲk2）を原点における平均曲率半径Ｒ0 の代わりに用いることにより、当該点の縮尺係数を算出する。
これを反復させると、当該点の縮尺係数は収束する性質を持ち、完全に収束した場合の縮尺係数は、当該点（入力座標）における縮尺係数である。
【００２７】
２．本願発明に係る装置の概要とその構成
図１は、この装置の説明図である。図２は、その処理のフローを示す説明図である。
２−１）位置計算装置の概要
この装置は、オペレータから平面直角座標上の座標値ｘ，ｙの入力を受付け、自動的に（出力まで人手を要することなく専ら装置のみにて）計算を行い、対応する地球楕円体（準拠楕円体を指す。ここではＷＧＳ８４楕円体を採用する。但し、ベッセル、ＧＲＳ８０等の他の楕円体を採用しても実施可能である。）上の緯度及び経度を出力するものである。
【００２８】
２−２）ハードウエア構成
この位置計算装置は、位置計算プログラムをコンピュータへ導入することにより、構築される。このコンピュータは、入力装置と、演算装置と、制御装置と、記憶装置（内部記憶装置及び外部記憶装置）とを備えた一般的なＰＣ（パーソナルコンピュータ）やＷＳ（ワークステーション）を採用して実施することができる。即ち、キーボード、マウスに代表される入力装置、ディスプレイやプリンタに代表される出力装置、ＣＰＵ（中央上方処理装置）に代表されるされる演算装置や制御装置（演算制御装置）、ＲＡＭ（ランダムアクセスメモリ）に代表される内部記憶装置、ハードディスクやフレキシブルディスクその他の磁気ディスク或いは光学ディスクに代表される外部記憶装置を備えた一般的なコンピュータを採用して実施することができる。
【００２９】
２−３）ソフトウエア
上記の位置計算プログラムは、コンピュータの資源を全て占有する専用のソフトウエアであっても実施可能であるが、ＭＳ−ＤＯＳ（登録商標）、Ｗｉｎｄｏｗｓ（登録商標）やＵＮＩＸ（登録商標）に代表される周知のＯＳ（オペレーティングシステム）上で動作するものとしても実施可能である。
また、位置計算プログラムは、インタプリンタ型のものであっても、アセンブラやコンパイラを使用して生成されたものであっても何れでもよい。
この実施の形態において、位置計算プログラムは、ＯＳが導入されたコンピュータ上に、導入されて位置計算装置を実現する。以下、ＯＳと、位置計算プログラムとを位置計算ソフトウエアと呼ぶ（必要に応じて単にソフトと呼ぶ）。
【００３０】
２−４）計算装置の構成
この実施の形態に示す位置計算装置は、等角投影法により地球楕円体（準拠楕円体）を投影する平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙを入力することにより、少なくとも当該座標値に対応する地球楕円体上の緯度φn 、経度λ、子午線収差γを出力するものである。
図１に示す通り、この装置は、データ置部１０を有する座標入力部１と、出力部２と、定数設定部３１と、初期値設定部３２と、条件設定部３３と、第１仮緯度収容部４１と、第２仮緯度収容部４２と、第１仮縮尺係数収容部５１と、第２仮縮尺係数収容部５２と、第１〜第７の７つの演算部６１，６２，６３，６４，６５，６６，６７を有する演算部６と、演算制御部７と、主判定部８１及び副判定部８２を有する判定部８と、回帰制御部９とを備える。
この実施の形態において、座標入力部１は、入力されたデータを一時的に保持するデータ置部１０を備える。また、定数設定部３１は、算出値置部３０を備える。
座標入力部１は、上記ソフトと、コンピュータが有する、上記の入力装置と、演算制御装置とにて構築される。出力部２は、上記ソフトと、コンピュータが有する上記出力装置と、演算制御装置とにて構築される。
定数設定部３１、初期値設定部３２及び条件設定部３３の夫々は、上記ソフトと、コンピュータが有する、上記の演算制御装置、内部記憶装置及び外部記憶装置とにて構築される。
第１仮緯度収容部４１、第２仮緯度収容部４２、第１仮縮尺係数収容部５１及び第２仮縮尺係数収容部５２は、夫々、上記ソフトと、コンピュータが有する、上記の演算制御装置と、内部記憶装置とにて構築される。
第１〜第７の７つの演算部６１，６２，６３，６４，６５，６６，６７と、演算制御部７と、主判定部８１と、副判定部８２と、回帰制御部９とは、夫々、上記ソフトと、コンピュータが有する、上記の演算制御装置と、内部記憶装置とにて構築される。
【００３１】
２−５）計算装置の構成の詳細
上記の座標入力部１は、オペレータからの上記座標値ｘ，ｙの入力を受け付け。データ置部１０に収容する。
定数設定部３１は、上記地球楕円体の、長半径ａ、短半径ｂ、第１離心率ｅ、第２離心率ｅ’の夫々、又はこれらに準じた値と、上記の座標ｘ，ｙが属する座標系の原点の、緯度φ0 又はその代替値、経度λ0 の夫々と、当該原点における縮尺係数ｍ0 又はその代替値と、後述する定数Ａ〜Ｆとが、既知の数値として設定されている。
ここで設定とは、この装置を用いた位置計算の度にオペレータがこれらの値（定数）を上記入力装置を通じて入力するもの、オペレータが予め外部記憶装置にこれらの値（定数）を収容しておくもの、位置計算プログラムの一部としてプログラムが備えプログラムの導入によって自動的に外部記憶装置に収容されるものの何れも含む。
ここでは、必要な定数は、プログラムに含まれるものとして外部記憶装置に導入され、導入後オペレータによってコンピュータが備える上記入力装置を通じて、書き換え可能とする（書き換えられるのは、主として、上記の、原点の緯度φ0 又はその代替値と、経度λ0 と、原点の縮尺係数ｍ0 又はその代替値である）。
定数設定部３１は、算出値置部３０を備える。算出値置部３０は、演算部６が算出した値を収容する。
【００３２】
演算部６について説明する。尚、ここでは（この２−５）の説明では）、演算部６を構成する第１〜第７の演算部６１〜６７夫々について、他の演算部との関係を無視し（演算部６１〜６７間のデータの流れについては無視し）、個々の演算部の内容を中心に説明する。従って、この２−５）において、演算部夫々の説明に用いる記号は、特に断りがない限り、他の演算部で用いるものと独立したものとする（他の演算部が、同じ記号で表された定数やパラメータを使うものであっても、特に断りがない限り、当該記号が示す数値は別のものである）。
夫々の演算にて、処理される加減乗除及び、正弦、余弦、正接は、周知の方法によって、コンピュータを用いたソフトウエアによる情報処理（数値計算）にて行われる（近似計算も含む）。
【００３３】
第１演算部６１は、少なくとも次の４つの演算１）〜４）を行うことが可能である。
演算１）
ｃ＝ａ2 ／ｂ＝ａ／（１−ｅ2 ）1/2 ＝ａ（１＋ｅ'2）1/2 …式６
を、地球楕円体の長半径ａ、地球楕円体の短半径ｂ、地球楕円体の第１離心率ｅ、地球楕円体の第２離心率ｅ’の各値から計算し、地球楕円体の極の曲率半径ｃを求めること。尚、上記第１及び第２の離心率は、地球楕円体の長半径ａと短半径ｂとから算定することができる。また、この地球楕円体の極の曲率半径ｃは、定数として算出によらずに、定数設定部３１が保持するものとしてもよい。
演算２）
Ｖ＝（１＋ｅ'2cos2φ）1/2 …式７
上記式７を、緯度φと第２離心率ｅ’の各値から計算し、Ｖを求めること。
演算３）
Ｒ＝ｃ／Ｖ2 …式８
上記式８を、地球楕円体の極の曲率半径ｃと、上記Ｖの値とから計算し、緯度φにおける平均曲率半径を求めること。
演算４）
ｍn+1 ＝（１＋ｙ2 ／２ｍn 2 Ｒ2 ＋ｙ4 ／２４ｍn 4 Ｒ4 ）×0.9999
…式９
上記式９を、平面直角座標上の座標値ｘ，ｙのｙの値と、仮の縮尺係数ｍn とから計算し、縮尺係数ｍn+1 を求めること。
この演算４）は、特許請求の範囲に示す式２を処理する機能である（式９の縮尺係数ｍn を原点縮尺係数ｍ0 とし、平均曲率半径Ｒを原点における平均曲率半径Ｒ0 としたものが、請求項２に示す式２である）。
【００３４】
第２演算部６２は、次の演算５）を行うことが可能である。
演算５）
Ｗ＝ａ（１−ｅ2 ）・｛Ａ・φ−（Ｂ／２）・（sin ２φ）
＋（Ｃ／４）・（sin ４φ）−（Ｄ／６）・（sin ６φ）
＋（Ｅ／８）・（sin ８φ）−（Ｆ／１０）・（sin 10φ）｝
＋ｘ／ｍ …式１０
上記式１０を、第２演算部６２は、地球楕円体の長半径ａの値と、定数Ａ〜Ｆの各値と、緯度φの値と、平面直角座標上の座標値ｘ，ｙのｘの値と、０以外の任意の実数ｍとから計算し、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗを求めること。尚、各定数Ａ〜Ｆは、第１離心率ｅによって定まるものである（後述のステップ１０７参照）。
【００３５】
第３演算部６３は、次の演算６）を行うことが可能である。
演算６）
Ｓ (φ )＝ａ（１−ｅ2 ）・｛Ａ・φ−（Ｂ／２）・（sin ２φ）
＋（Ｃ／４）・（sin ４φ）−（Ｄ／６）・（sin ６φ）
＋（Ｅ／８）・（sin ８φ）−（Ｆ／１０）・（sin10 φ）｝
…式１１
上記の式１１を、地球楕円体の長半径ａと、第１離心率ｅと、上記の各定数Ａ〜Ｆと、緯度φとから計算し、赤道から緯度φまでの子午線弧長Ｓ (φ )を求めること。
【００３６】
第４演算部６４は、次の演算７）を行うことが可能である。
演算７）
φn+1 ＝φn
＋（Ｗ−Ｓ (φn ) ）（１−ｅ2 sin2φn ）3/2 ／ａ（１−ｅ2 ）
…式１２
上記式１２を、仮の垂線緯度φn と、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗと、仮の垂線緯度φn における赤道からの子午線弧長Ｓ (φn ) と、地球楕円体の長半径ａと、第１離心率ｅの、各値から計算し垂線緯度φn+1 を求めること。
この式１２は、特許請求の範囲に示す式１と対応する。
【００３７】
第５演算部６５は、次の５つの演算８）〜１２）を行うことが可能である。
演算８）
ｔ＝tan φ …式１３
上記式１３を、与えられた緯度φ値から、計算し、ｔを求めること。
演算９）
η＝ｅ’cos φ …式１４
上記式１４を、地球楕円体の第２離心率ｅ’と、緯度φの、夫々の値から、計算し、ηを求めること。
演算１０）
Ｍ＝ｃ／Ｖ 3 …式１５
上記式１５を、地球楕円体の極の曲率半径ｃの値と、緯度φにて定まるＶの値とから計算し、当該緯度φについての子午線曲率半径Ｍを求めること。
演算１１）
Ｎ＝ｃ／Ｖ …式１６
上記式１６を、地球楕円体の極の曲率半径ｃの値と、緯度φにて定まるＶの値とから計算し、当該緯度φについての卯酉線曲率半径Ｎを求めること。
演算１２）
φn+1 ＝φn −ｔn （ｙ／ｍt ）2 ／２Ｍn Ｎn ＋ｔn （ｙ／ｍt ）4 （５＋３ｔn 2 ＋ηn 2 −９ｔn 2 ηn 2 −４ηn 4 ）／24Ｍn Ｎn 3 −ｔn （ｙ／ｍt ）6 （61＋90ｔn 2 ＋45ｔn 4 ）／720 Ｍn Ｎn 5 …式１７
上記式１７を、平面直角座標上の座標値ｘ，ｙのｙの値、最近似垂線緯度φn 、縮尺係数ｍt 、最近似垂線緯度φn についてのｔn 、最近似垂線緯度φn についてのηn 、最近似垂線緯度φn についての子午線曲率半径Ｍn 、及び、最近似垂線緯度φn についての卯酉線曲率半径Ｎn の、夫々の値から計算し、当該緯度φn+1 を求めること。
この式１７は、特許請求の範囲に示す式３に対応する。
【００３８】
第６演算部６６は、次の２つの演算１３）及び１４）を行うことが可能である。
演算１３）
Δλ＝（ｙ／ｍt ）／Ｎn cos φn
−（ｙ／ｍt ）3 （１＋２ｔn 2 ＋ηn 2 ）／６Ｎn 3cosφn
＋（ｙ／ｍt ）5 （５＋２８ｔn 2 ＋２４ｔn 4 ）／１２０Ｎn 5cosφn
…式１８
上記式１８を、平面直角座標上の座標値ｘ，ｙのｙの値、最近似垂線緯度φn 、縮尺係数ｍt 、最近似垂線緯度φn についてのｔn 、最近似垂線緯度φn についてのηn 、及び、最近似垂線緯度φn についての卯酉線曲率半径Ｎn の、夫々の値から計算し、経度差分Δλを求めること。
この式１８は、特許請求の範囲（請求項４）の式４に対応する。
演算１４）
λ＝λ0 ＋Δλ …式１９
上記式１９を、平面直角座標上の座標値ｘ，ｙの属する原点の経度λ0 と、経度差分の夫々の値とから、計算し、経度λを求めること。
【００３９】
第７演算部６７は、次の演算１５）を行うことが可能である。
演算１５）
γ＝（ｙ／ｍt ）／Ｎn
−（ｙ／ｍt ）3 ｔn （１＋ｔn 2 −ηn 2 ）／３Ｎn 3
＋（ｙ／ｍt ）5 ｔn （２＋５ｔn 2 ＋３ｔn 4 ）／１５Ｎn 5
…式２０
上記式２０を、平面直角座標上の座標値ｘ，ｙのｙの値、縮尺係数ｍt 、最近似垂線緯度φn についてのｔn 、最近似垂線緯度φn についてのηn 、及び、最近似垂線緯度φn についての卯酉線曲率半径Ｎn の、夫々の値から計算し、子午線収差γを求めること。
この式２０は、特許請求の範囲（請求項５）の式５に対応する。
【００４０】
演算制御部７は、演算部６の各演算部６１〜６７へ演算の指令を行う。各演算部６１〜６７から、演算及び演算に必要な数値の所在（置き場）の指令を受けて、演算制御部７の指令の順に各演算部６１〜６７は、指令された演算を行う。
また、演算制御部７は、演算部６１〜６７の演算の結果得た緯度及び縮尺係数の収束状況について、判定部８に判定を行わせる。
判定部８は、垂線緯度の収束状態の判定を行う主判定部８１と、縮尺係数の判定状態の判定を行う副判定部８２とを備える。
回帰制御部９は、判定部８の判定を受けて、演算制御部７に、判定に対応した指示を行う。当該指示に従って、演算制御部７は、各演算部６１〜６７に指令を出す。即ち、回帰演算部９は、判定部８による判定の結果に従って、演算制御部７を通じて、各演算部６１〜６７を制御する。
【００４１】
初期値設定部３２は、（演算部６が行う）垂線緯度計算の初期値を決定するための数値（以下、必要に応じて初期値決定用助数と呼ぶ。）が設定され、当該数値を保持する。
条件設定部３３は、判定部８の判定（比較）に必要な基準値が設定され、設定された基準値を保持する。
【００４２】
第１仮緯度収容部４１は、演算部６が算出した緯度の値を収容する。第１仮緯度収容部４１は、新たに演算部６が算出した緯度の値を受け取ることにより、それまで収容していた緯度の値を、第２仮緯度収容部４２へ渡し、それまで収容した緯度の値を新たに受け取った緯度の値に書き換えて、新たな緯度の値を収容する。
第２仮緯度収容部４２は、第１仮緯度収容部４１から、第１仮緯度収容部４１が収容していた緯度の値を受け取り、当該緯度の値を収容する。第２仮緯度収容部４２は、更に第１仮緯度収容部４１から緯度の値を受け取ることにより、それまで収容していた緯度の値を受け取った後の緯度の値に書き換えて、当該後の緯度の値を収容する。
【００４３】
第１仮縮尺係数収容部５１は、演算部６が算出した縮尺係数の値を収容する。第１仮縮尺係数収容部５１は、新たに演算部６が算出した縮尺係数の値を受け取ることにより、それまで収容していた縮尺係数の値を、第２仮縮尺係数収容部５２へ渡し、それまで収容した縮尺係数の値を新たに受け取った縮尺係数の値に書き換えて、新たな縮尺係数の値を収容する。
第２仮縮尺係数収容部５２は、第１仮縮尺係数収容部５１から、第１仮縮尺係数収容部５１が収容していた縮尺係数の値を受け取り、当該縮尺係数の値を収容する。第２仮緯度収容部４２は、更に第１仮縮尺係数収容部５１から縮尺係数の値を受け取ることにより、それまで収容していた縮尺係数の値を受け取った後の縮尺係数の値に書き換えて、当該後の縮尺係数の値を収容する。
出力部２は、演算制御部７と回帰演算部９の指令を受けて、演算部６の演算結果を出力する。
【００４４】
３．計算の手順
以下、この装置を利用した緯度・経度・子午線収差の算出の流れについて、図２を用いて、説明する。
この計算方法は、主工程１００を遂行し、更に縮尺係数の収束状況の判定結果により後処理工程を遂行するものである。
また、上記の主工程１００の遂行に先立ち、事前に、前処理工程を行っておく。
【００４５】
３−０）前処理工程
前処理工程は、オペレータにより、位置計算前に行われる種々の設定工程である。
詳しくは、前処理工程において、オペレータにより、初期値設定部３２にて緯度計算の初期値を決定するための数値（初期値決定用助数）が設定される。このような助数として、１”（１秒）をrad （ラジアン）に換算した値である 1／206265が、条件設定部３３において設定される。但しこのような数値は変更可能である。
又、前処理工程において、オペレータにより、条件設定部３３に、垂線緯度の収束条件ＩＰＳ（しきい値）が設定される。
この収束条件ＩＰＳは、計算の精度を確保する上で重要である。即ち、演算部における各式にて的確な計算がなされても、収束条件ＩＰＳが広いと、十分な精度を確保するのが困難になるからである。
ここでは、収束条件ＩＰＳは、２”× 1０-5／206265rad と設定する。
更に、前処理工程において、オペレータにより、初期値設定部３２にて、縮尺係数の収束条件ＩＰＳ２が設定される。ここでは、ＩＰＳ２は、座標値の有効数字（最大９桁）を考慮して、４× 1０-12 と設定してある。
但し、各収束条件ＩＰＳ、ＩＰＳ２の各数値は、上記の値に限定するものではなく、変更可能である。
主工程１００以降の工程は、オペレータから座標値ｘ，ｙを受け付けたこの位置計算装置が、位置計算を実際に行う工程である。
主工程は、次のステップ１０１〜１２５を順次遂行する。
【００４６】
３−１）ステップ１０１（座標値受付工程）
オペレータは、座標入力部１により、直角座標の座標値ｘ，ｙ、及び当該座標値が属する座標系Ｚを入力する。座標系Ｚの入力については、日本に設定されている１９の座標系のうち何れに属するかを特定するものである。
上記の入力作業は、オペレータが行う（これ以降の他の作業は、位置計算装置が行う。尚、この実施の形態において、後述する初期値設定部３２の設定、及び定数設定部３１内の各代替値の設定は、位置計算前に、オペレータにより行われる）。
座標入力部１は、入力された座標値ｘ，ｙを、データ置部１０に収容し、演算制御部７に座標値の入力があったことを通知する。
演算制御部７は、上記通知を受けて、第１演算部６１に、定数設定部３１を参照させる。即ち、定数設定部３１は、各座標系の原点の緯度及び経度の数値を保持しており、第１演算部６１は、入力された（座標値ｘ，ｙが属する）座標系Ｚの、原点の緯度φ0 及び経度λ0 を特定する。これにて、以下、各演算部が定数設定部３１にて参照する原点の緯度φ0 及び経度λ0 は、確定する。
【００４７】
３−２）ステップ１０２（極の曲率半径計算工程）
このステップ１０２において、演算制御部７は、第１演算部６１に、定数設定部３１を参照させて、定数設定部３１が保持している、地球楕円体の長半径ａと短半径ｂから、或いは、地球楕円体の長半径ａと第１離心率ｅから、或いは、地球楕円体の長半径ａと第２離心率ｅ’から、下記の計算１を行わせ、地球楕円体の極の曲率半径ｃを算出させる。第１演算部６１は、算出した極の曲率半径ｃを、算出値置部３０に収容する。
ｃ＝ａ2 ／ｂ＝ａ／（１−ｅ2 ）1/2 ＝ａ（１＋ｅ'2）1/2 …計算１
計算１は、第１演算部６１の演算１）の式６の計算機能にて実行される。
【００４８】
３−３）ステップ１０３（Ｖ0 計算工程）
ステップ１０３において、演算制御部７は、第１演算部６１に、定数設定部３１を参照させて、第２離心率ｅ’と当該座標が属する原点の緯度φ0 （以下単に原点緯度φ0 と呼ぶ。）とから、計算２を行わせ、Ｖ0 を算出させる。第１演算部６１は、算出したＶφ0 を、算出値置部３０に収容する。
Ｖ0 ＝（１＋ｅ'2cos2φ0 ）1/2 …計算２
計算２は、第１演算部６１の演算２）の式７の計算機能にて実行される。
【００４９】
３−４）ステップ１０４（原点の平均曲率半径計算工程）
ステップ１０４において、演算制御部７は、第１演算部６１に、定数設定部３１（の算出値置部３０）を参照させて、極の曲率半径ｃとＶ0 の値とから、計算３を行わせ、当該座標系原点の平均曲率半径Ｒ0 を算出させる。第１演算部６１は、算出した平均曲率半径Ｒ0 を算出値置部３０に収容する。
Ｒ0 ＝ｃ／Ｖ02…計算３
計算３は、第１演算部６１の演算３）の式８の計算機能にて実行される。
【００５０】
３−５）ステップ１０５（縮尺係数計算工程）
ステップ１０５において、演算制御部７は、定数設定部３１の原点における縮尺係数ｍ0 （以下原点縮尺係数ｍ0 と呼ぶ。）を、第１仮縮尺係数収容部５１に収容させる（定数設定部３１の上記原点縮尺係数ｍ0 を第１仮縮尺係数収容部５１へコピーする）。この原点縮尺係数ｍ0 は、0.9999という数値を採る。
そして、演算制御部７は、第１演算部６１に、座標入力部１のデータ置部１０及び定数設定部３１（又は第１仮縮尺係数収容部５１）を参照させて、入力された座標値ｙの値と、原点縮尺係数ｍ0 と、上記平均曲率半径Ｒ0 とから、計算４を行わせ、縮尺係数ｍt を算出させる。
ｍt ＝（１＋ｙ2 ／２ｍ02Ｒ02＋ｙ4 ／２４ｍ04Ｒ04）×0.9999…計算４
（右辺の0.9999は、数値として原点縮尺係数ｍ0 と同じであるが、これは基準子午線上の縮尺係数を0.9999とするための一定の縮尺係数である。）
計算４は、第１演算部６１の演算４）の式９の計算機能にて実行される。
第１演算部６１は、算出した縮尺係数ｍt を第１仮縮尺係数収容部５１に収容する。第１仮縮尺係数収容部５１は、収容していた原点縮尺係数ｍ0 を第２仮縮尺係数収容部５２へ渡し、原点縮尺係数ｍ0 に代え上記の通り上記縮尺係数ｍt を収容する。
【００５１】
３−６）ステップ１０６（垂線緯度の初期値設定工程）
ステップ１０６において、演算制御部７は、第１演算部６１に、定数設定部３１と、初期値設定部３２とを参照させて、原点緯度φ0 から助数を減じた値、即ち、原点緯度φ0 − 1／206265を計算させる。第１演算部６１は、算出した値を緯度φj として、第１仮緯度収容部４１に収容する。
【００５２】
３−７）ステップ１０７（赤道から当該座標値ｘまでの基準子午線上の子午線弧長計算工程）
このステップ１０７において、演算制御部７は、第２演算部６２に、定数設定部３１と、座標入力部１のデータ置部１０と、第１仮縮尺係数収容部５１とを参照させて、長半径ａと、第１離心率ｅと、原点緯度φ0 と、座標値ｘと、縮尺係数ｍt とから、計算５をさせ、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとして算出させる。
Ｗ＝ａ（１−ｅ2 ）・｛Ａ・φ0 −（Ｂ／２）・（sin ２φ0 ）
＋（Ｃ／４）・（sin ４φ0 ）−（Ｄ／６）・（sin ６φ0 ）
＋（Ｅ／８）・（sin ８φ0 ）−（Ｆ／10）・（sin10 φ0 ）｝
＋ｘ／ｍt …計算５
計算５の、Ａ〜Ｆは、次の通りである。
Ａ＝１＋（３／４）ｅ2 ＋（45／64）ｅ4 ＋（175 ／256 ）ｅ6
＋（11025 ／16384 ）ｅ8 ＋（43659 ／65536 ）ｅ10
Ｂ＝ （３／４）ｅ2 ＋（15／16）ｅ4 ＋（525 ／512 ）ｅ6
＋（2205／2048）ｅ8 ＋（72765 ／65536 ）ｅ10
Ｃ＝ （15／64）ｅ4 ＋（105 ／256 ）ｅ6
＋（2205／4096）ｅ8 ＋（10395 ／16384 ）ｅ10
Ｄ＝ （35／512 ）ｅ6
＋（315 ／2048）ｅ8 ＋（31185 ／131072）ｅ10
Ｅ＝ （315 ／16384 ）ｅ8 ＋（3465／65536 ）ｅ10
Ｆ＝ （693 ／131072）ｅ10
計算５は、第２演算部６２の演算５）の式１０の計算機能にて実行される。
前述の通り、この定数Ａ〜Ｆは、定数設定部３１に予め収容されており、第２演算部６２や他の演算部が参照して利用する。
ここで、計算５の右辺のｘ／ｍt 以外の部分が、地球楕円体の赤道から座標値ｘ，ｙが属する座標系原点までの子午線弧長を示している。また、この縮尺係数ｍt は、特許請求の範囲でいう、仮の縮尺係数である。
第２演算部６２は、算出した、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗを定数設定部３１の算出値置部３０に収容する。
【００５３】
３−８）ステップ１０８（緯度に対する赤道からの子午線弧長算出工程）
このステップ１０８において、演算制御部７は、第３演算部６３に、定数設定部３１と、第１仮緯度収容部４１とを参照させ、長半径ａ、第１離心率ｅ、定数Ａ〜Ｆと、緯度φj とから、計算６を行わせ、緯度φj に対する赤道からの子午線弧長Ｓ (φj ) を算出させる。
Ｓ (φj ) ＝ａ（１−ｅ2 ）・｛Ａ・φj −（Ｂ／２）・（sin ２φj ）
＋（Ｃ／４）・（sin ４φj ）−（Ｄ／６）・（sin ６φj ）
＋（Ｅ／８）・（sin ８φj ）−（Ｆ／１０）・（sin10 φj ）｝
…計算６
即ち、このステップ１０８は、子午線曲率半径Ｍを以って赤道から垂線緯度φj まで積分するステップであり、上記の通り計算６は、第３演算部６３にて当該積分の近似計算である（ここでいう積分とは、当該近似計算を指している）。
計算６は、第３演算部６３の演算６）の式１１の計算機能にて実行される。
第３演算部６３は、上記にて算出した子午線弧長Ｓ (φj ) を、上記算出値置部３０に収容する。
【００５４】
３−９）ステップ１０９（初期の垂線緯度の収束計算工程）
ステップ１０９において、演算制御部７は、第４演算部６４に、定数設定部３１と、第１仮緯度収容部４１とを参照させ、緯度φj 、長半径ａ、第１離心率ｅ、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗ、子午線弧長Ｓ (φj ) とから、計算７を行わせ、垂線緯度φj+1 を算出させる。
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ） …計算７
計算７は、第４演算部６４の演算７）の式１２の計算機能にて実行される。
第４演算部６４は、算出した緯度φj+1 を、第１仮緯度収容部４１へ収容する。第４演算部６４から上記の緯度φj+1 を受けた第１仮緯度収容部４１は、収容していた緯度φj を第２仮緯度収容部４２へ渡し、緯度φj に代え上記緯度φj+1 の値を収容する。
【００５５】
３−１０）ステップ１１０（垂線緯度の収束判定工程）
ステップ１１０において、判定部８の主判定部８１は、第１仮緯度収容部４１と第２仮緯度収容部４２と条件設定部３３とを参照し、緯度φj+1 と緯度φj の差の絶対値即ち、｜φj+1 −φj ｜と、緯度の収束条件ＩＰＳ（しきい値）とを比較する。そして、｜φj+1 −φj ｜＞ＩＰＳと判定した場合、判定部８（主判定部８１）が、その結果を回帰制御部９に通知する。
この通知を受けた回帰制御部９は、第１仮緯度収容部４１の緯度φj+1 の値を緯度φj の値として（第１仮緯度収容部４１の緯度φj+1 を緯度φj に書き換えて）、装置の処理を再びステップ１０８に移す。そして、回帰制御部９は、ステップ１０８からステップ１１０を、判定部８（主判定部８１）が｜φj+1 −φj ｜≦ＩＰＳと判定するまで繰り返させる（垂線緯度の反復計算をさせる）。この反復計算において、垂線緯度φj+1 を垂線緯度φj に置き換える（書き換える）のみであり、演算部６（の各演算部）が参照する他の数値を変更しない（当該他の数値については、当初の値を用いて最初から計算される）。
以上により、ステップ１１０において、判定部８（の主判定部８１）は、｜φj+1 −φj ｜≦ＩＰＳと判定した場合、その結果を回帰制御部９に通知する。
回帰制御部９は、当該収束の通知を受けて、次のステップ１１１へ処理を移行させる。即ち、回帰制御部９は、演算制御部７に次のステップ１１１の処理を行わせる。
【００５６】
３−１１）ステップ１１１（Ｖk 計算工程）
ステップ１１１において、演算制御部７は、第１演算部６１に、定数設定部３１と、第１仮緯度収容部４１とを参照させて、第２離心率ｅ’と垂線緯度φj+1 とから、計算８を行わせ、Ｖk を算出させる。尚、説明の便宜上、この垂線緯度φj+1 を緯度φk として説明する。
Ｖk ＝（１＋ｅ'2cos2φk ）1/2 …計算８
計算８は、第１演算部６１の演算２）の式７の計算機能にて実行される。
第１演算部６１は、算出したＶk を、算出値置部３０に収容する。
【００５７】
３−１２）ステップ１１２（第２の値計算工程）
ステップ１１２において、演算制御部７は、第５演算部６５に、第１仮緯度収容部４１を参照させ、緯度φk の値を用いて（即ち、緯度φj+1 の値を緯度φk の値として用いて）、計算９を行わせ、ｔk （説明の便宜上第２の値ｔk と呼ぶ。）を算出させる。
ｔk ＝tan φk …計算９
計算９は、第５演算部６５の演算８）の式１３の計算機能にて実行される。
第５演算部６５は、算出した第２の値ｔk を、算出値置部３０へ収容する。
【００５８】
３−１３）ステップ１１３（第１の値計算工程）
ステップ１１３において、演算制御部７は、第５演算部６５に、定数設定部３１と第１仮緯度収容部４１を参照させ、第２離心率ｅ’と共に、緯度φk の値を用いて、計算１０を行わせ、ηk （説明の便宜上第１の値ηk と呼ぶ。）を算出させる。
ηk ＝ｅ’cos φk …計算１０
計算１０は第５演算部６５の演算９）の式１４の計算機能にて実行される。
第５演算部６５は、算出した第１の値ηk を、算出値置部３０へ収容する。
【００５９】
３−１４）ステップ１１４（子午線曲率半径の計算工程）
ステップ１１４において、演算制御部７は、第５演算部６５に、定数設定部３１を参照させ、極の曲率半径ｃとＶk とから、計算１１にて、緯度φk における子午線曲率半径Ｍk を算出させる。
Ｍk ＝ｃ／Ｖk 3 …計算１１
計算１１は、第５演算部６５の演算１０）の式１５の計算機能にて実行される。
第５演算部６５は、算出した子午線曲率半径Ｍk の値を、算出値置部３０へ収容する。
【００６０】
３−１５）ステップ１１５（卯酉線曲率半径の計算工程）
ステップ１１５において、演算制御部７は、第５演算部６５に、定数設定部３１を参照させ、極の曲率半径ｃとＶk とから、計算１２にて、緯度φk における卯酉線曲率半径Ｎk を算出させる。
Ｎk ＝ｃ／Ｖk …計算１２
計算１２は、第５演算部６５の演算１１）の式１６の計算機能にて実行される。
第５演算部６５は、算出した卯酉線曲率半径Ｎk の値を、算出値置部３０へ収容する。
【００６１】
３−１６）ステップ１１６（入力座標に対応する緯度の計算工程）
ステップ１１６において、回帰制御部９は、演算部制御部７を通じて、第５演算部６５に、緯度φk を与え、計算１３を行わせる。具体的には、演算制御部７は、回帰制御部９からの指令を受けて、第５演算部６５に、座標入力部１のデータ置部１０、第１仮緯度収容部４１、第１仮縮尺収容部５１及び定数設定部３１を参照させて、上記にて、第４演算部６４に算出させた最近似垂線緯度φk と、当該点の縮尺係数ｍt と、最近似垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、計算１３の計算にて、緯度φk+1 を計算させる。
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …計算１３
計算１３は、第５演算部６５の演算１２）の式１７の計算機能にて実行される。
第５演算部６５は、算出した緯度φk+1 を、当該点（入力座標点）の緯度φn として、第１仮緯度収容部４１へ、それまで第１仮緯度収容部４１に収容されていた値と置き換えて収容する。このとき、第１仮緯度収容部４１は、それまで保持していた緯度φk の値を第２仮緯度収容部４２へ渡す。
【００６２】
３−１７）ステップ１１７（経度差分計算工程）
このステップ１１７は、上記のステップ１１６と平行として行われる。即ち、ステップ１１７は、ステップ１１６で得られた結果（緯度φk+1 ）を用いるのではなく、ステップ１１６の処理に用いたデータを用いて、実行される。
ステップ１１７において、演算制御部７は、第６演算部６６に、座標入力部１のデータ置部１０、第１仮縮尺係数収容部５１及び算出値置部３０を参照させて、最近似垂線緯度φk （ステップ１１１に至るまでの一連の処理で得られた緯度、即ちステップ１１１〜１１６の処理のために用いた緯度φk ）と、当該点の縮尺係数ｍt と、最近似垂線緯度φk における子午線曲率半径Ｍk （ステップ１１４の結果）及び卯酉線曲率半径Ｎk （ステップ１１５の結果）と、ｅ’cos φk にて定まる第１の値ηk （ステップ１１３の結果）と、tan φk にて定まる第２の値ｔk （ステップ１１２の結果）と、上記座標値ｙとから、計算１４を行わせ、経度差分Δλを算出させる。
Δλ＝（ｙ／ｍt ）／Ｎk cos φk
−（ｙ／ｍt ）3 （１＋２ｔk 2 ＋ηk 2 ）／６Ｎk 3cosφk
＋（ｙ／ｍt ）5 （５＋２８ｔk 2 ＋２４ｔk 4 ）／１２０Ｎk 5cosφk
…計算１４
第６演算部６６は、算出した経度差分Δλを、定数設定部３１の算出値置部３０に、収容する。
【００６３】
３−１８）ステップ１１８（経度計算工程）
このステップ１１８と当該ステップ１１８に続くステップ１１９は、ステップ１１７と同様、上記のステップ１１６と平行として行われる。
ステップ１１８において、演算制御部７は、第６演算部６６に、定数設定部３１を参照させて、原点経度λ0 と経度差分Δλの和の計算を即ち、計算１５をさせて、経度λを算出させる。
第６演算部６６は、算出した経度λの値を算出値置部３０に収容する。
λ＝λ0 ＋Δλ …計算１５
計算１５は、第６演算部６６の演算１４）の式１９の計算機能にて実行される。
【００６４】
３−１９）ステップ１１９（子午線収差の計算工程）
このステップ１１９は、上記の通り、ステップ１１８に続く工程であり、ステップ１１６と平行して行われる。即ち、ステップ１１６で得られた結果（緯度φk+1 ）を用いるのではなく、ステップ１１６の処理に用いたデータを用いて、実行される。
ステップ１１９において、演算制御部７は、第７演算部６７に、座標入力部１のデータ置部１０、第１仮縮尺係数収容部５１及び算出値置部３０を参照させて、最近似垂線緯度φk （ステップ１１１に至るまでの一連の処理で得られた緯度、即ちステップ１１１〜１１６の処理のために用いた緯度φk ）と、当該点の縮尺係数ｍt と、最近似垂線緯度φk における卯酉線曲率半径Ｎk （ステップ１１５の結果）と、ｅ’cos φk にて定まる第１の値ηk （ステップ１１３の結果）と、tan φk にて定まる第２の値ｔk （ステップ１１２の結果）と、上記座標値ｙとから、計算１６を行わせ、子午線収差γを算出させる。
γ＝（ｙ／ｍt ）／Ｎk
−（ｙ／ｍt ）3 ｔk （１＋ｔk 2 −ηk 2 ）／３Ｎk 3
＋（ｙ／ｍt ）5 ｔk （２＋５ｔk 2 ＋３ｔk 4 ）／１５Ｎk 5
…計算１６
計算１６は、第７演算部６７の演算１５）の式２０の計算機能にて実行される。
第７演算部６７は、算出した子午線収差γを算出値置部３０に収容する。
【００６５】
３−２０）ステップ１２０（Ｖk2計算工程）
ステップ１２０は、ステップ１１６に続く工程である。
ステップ１２０において、演算制御部７は、第１演算部６１に、定数設定部３１と、第１仮緯度収容部４１とを参照させて、第２離心率ｅ’とステップ１１６で得た当該点の緯度φk+1 とから、計算１７を行わせ、Ｖk2を算出させる。尚、説明の便宜上、この緯度φk+1 を緯度φk2として説明する。
Ｖk2＝（１＋ｅ'2cos2φk2）1/2 …計算１７
計算１７は、第１演算部６１の演算２）の式７の計算機能にて実行される。
第１演算部６１は、算出したＶk2を、算出値置部３０に収容する。
【００６６】
３−２１）ステップ１２１（子午線曲率半径の計算工程）
ステップ１２１において、演算制御部７は、第５演算部６５に、定数設定部３１を参照させ、極の曲率半径ｃとＶk2とから、計算１８にて、緯度φk2における子午線曲率半径Ｍk2を算出させる。
Ｍk2＝ｃ／Ｖk23 …計算１８
計算１８は、第５演算部６５の演算１０）の式１５の計算機能にて実行される。
第５演算部６５は、算出した子午線曲率半径Ｍk2の値を、算出値置部３０へ収容する。
【００６７】
３−２２）ステップ１２２（卯酉線曲率半径の計算工程）
ステップ１２２において、演算制御部７は、第５演算部６５に、定数設定部３１を参照させ、極の曲率半径ｃとＶk2とから、計算１９にて、緯度φk2における卯酉線曲率半径Ｎk2を算出させる。
Ｎk2＝ｃ／Ｖk2 …計算１９
計算１９は、第５演算部６５の演算１１）の式１６の計算機能にて実行される。
第５演算部６５は、算出した卯酉線曲率半径Ｎk2の値を、算出値置部３０へ収容する。
【００６８】
３−２３）ステップ１２３（平均曲率半径の計算工程）
ステップ１２３において、演算制御部７は、第１演算部６１に、定数設定部３１を参照させて、極の曲率半径ｃとＶk2の値とから、計算２０を行わせ、当該座標（入力座標）の平均曲率半径Ｒk2を算出させる。第１演算部６１は、算出した平均曲率半径Ｒk2を算出値置部３０に収容する。
Ｒk2＝ｃ／Ｖk22 …計算２０
計算２０は、第１演算部６１の演算３）の式８の計算機能にて実行される。
【００６９】
３−２４）ステップ１２４（縮尺係数計算工程）
ステップ１２４において、演算制御部７は、第１演算部６１に、座標入力部１のデータ置部１０及び第１仮縮尺係数収容部５１を参照させて、入力された座標値ｙの値と、縮尺係数ｍt1（第１仮縮尺係数収容部５１に収容されている前記縮尺係数ｍt を縮尺係数ｍt1として説明する。以下同じ。）と、上記平均曲率半径Ｒk2とから、計算２１を行わせ、縮尺係数ｍt2を算出させる。
ｍt2＝（１＋ｙ2 ／２ｍt12 Ｒk22 ＋ｙ4 ／２４ｍt14 Ｒk24 ）×0.9999
…計算２１
計算２１は、第１演算部６１の演算４）の式９の計算機能にて実行される。
第１演算部６１は、算出した縮尺係数ｍt2を第１仮縮尺係数収容部５１に収容する。第１仮縮尺係数収容部５１は、収容していた原点縮尺係数ｍt を第２仮縮尺係数収容部５２へ渡し、上記の通り縮尺係数ｍt に代え上記縮尺係数ｍt2を収容する。
【００７０】
３−２５）ステップ１２５（縮尺係数の収束判定工程）
ステップ１２５において、判定部８の副判定部８２は、第１仮縮尺係数収容部５１と、第２仮縮尺係数収容部５２と、条件設定部３３とを参照し、縮尺係数ｍt1と縮尺係数ｍt2の差の絶対値即ち｜ｍt2−ｍt1｜と、縮尺係数の収束条件ＩＰＳ２とを比較する。そして、縮尺係数ｍt2（と縮尺係数ｍt1の差の絶対値）が収束したと判定された場合、即ち、副判定部８２が｜ｍt2−ｍt1｜≦ＩＰＳ２と判定した場合、判定部８は、その結果を回帰制御部９に通知する。
【００７１】
３−２６）ステップ１２６（出力工程）
上記ステップ１２５にて、縮尺係数の収束の通知を受けた回帰制御部９は、出力部２により、少なくとも、第１仮緯度収容部４１が収容している当該点の緯度φk+1 の値（φn ) と、算出値置部３０が収容している経度λの値及び子午線収差γの値とを、出力させる。このとき、必要に応じて、座標入力部１のデータ置部１０の入力座標値ｘ，ｙや、定数設定部３１が保持しその算出値置部３０が収容している、上記入力座標が属する原点の緯度や経度の数値、地球楕円体の長半径や短半径などの定数、或いは計算工程中に算出した中間の値を、併せて、出力部２に出力させるものとしても、上記の結果がどのような座標値に対応するか、また、どのような経緯を経て算出されたか、オペレータに確認させることができ、便利である。
【００７２】
一方、ステップ１２５において、判定部８（副判定部８２）が｜ｍt2−ｍt1｜＞ＩＰＳ２と判定した場合、ステップ１２６に移行せずに、前述の後処理工程を遂行する。
即ち、上記にて、｜ｍt2−ｍt1｜＞ＩＰＳ２と判定された場合、判定部８は、その結果を回帰制御部９に通知する。この通知を受けた回帰制御部９は、演算制御部７を通じて、上記主工程１００の各ステップ（ステップ１０２〜１２５）を繰り返させる。この後処理工程において、主工程１００と異なるのは、ステップ１０５の処理において、前述の原点縮尺係数ｍ0 に代え、第１仮縮尺係数収容部５１に収容されている縮尺係数ｍt2を第１演算部６１に参照させて、計算４の計算を行わせるということである（原点縮尺係数ｍ0 として、上記主工程１００を経て得た縮尺係数ｍt2を用いる。即ち、後処理工程で原点縮尺係数ｍ0 として用いる値は、上記主工程１００を経て得た縮尺係数ｍt2である）。他の値については、上記主工程１００と同様のものを用いる。即ち、後処理工程は、縮尺係数ｍt2を用いて、主工程１００をやり直す工程である。従って、後処理工程において、主工程１００の結果として加味されるのは、主工程１００で算出された縮尺係数のみであり、主工程１００で算出した緯度等の値は用いず、再度主工程１００で用いたのと同様の値を用いて計算する（縮尺係数を使わない計算や、縮尺係数を用いる計算であっても、縮尺係数以外の数値については、主工程１００で用いたものと同様の数値を用いて計算を行う）。
この後処理工程でのステップ１２５において、｜ｍt2−ｍt1｜＞ＩＰＳ２と判定された場合、判定部８は、その結果を回帰制御部９に通知する。そして、当該後処理工程を経て得た縮尺係数ｍt2を原点縮尺係数ｍ0 として、再度後処理工程（主工程１００のステップ１０２からのやり直し）を行う。このように、１回の後処理工程を経ても、縮尺係数ｍt2（と縮尺係数ｍt1の差の絶対値）が収束したと判定されない場合、収束するまで、即ち、副判定部８２が｜ｍt2−ｍt1｜≦ＩＰＳ２となるまで、上記の後処理工程を繰り返す。
後処理工程におけるステップ１２５にて、縮尺係数ｍt2が収束した即ち｜ｍt2−ｍt1｜≦ＩＰＳ２と副判定部８２に判定された場合、上記の通り、回帰制御部９は、演算制御部７に、ステップ１２６を遂行させる。
【００７３】
上記の通り、後処理工程（主工程１００のやり直し工程）が繰り返されると、ステップ１０７における計算５のｘ／ｍt が変化し、座標値ｘまでの基準子午線上の子午線弧長Ｗの値が変化し、その結果、計算７の計算結果がより精度の高いものとなるのである。
上記の実施の形態において、後処理工程は、判定部８（副判定部８２）による縮尺係数の収束不十分、即ち、｜ｍt2−ｍt1｜＞ＩＰＳ２の判定により、行うものとした。
この他、前述の通り、このような縮尺係数の収束の判定を行わずに、設定した回数必ず後処理工程を実行するものとしても実施可能である。例えば、後処理工程を必ず１回（或いは２回以上）行うものとして実施することが可能である。
また、上記の実施の形態において、ステップ１０７で用いる仮の数値として、ステップ１０５で算出した縮尺係数ｍt を用いたが、初回のステップ１０７（後処理工程でのステップ１０７は除く。）については、０以外の任意の実数を当該縮尺係数ｍt の代替値として用いて計算を行うものとしても実施可能である。
但し、ステップ１０７で用いる仮の縮尺係数として、ステップ１０５で算出した縮尺係数ｍt を用いることにより、収束が効率良く行われるので、例えば、ステップ１１０における判定の結果ステップ１０８への反復（小ループ）を３回程度に抑えることができ、また、ステップ１２５における後処理工程の反復（大ループ）を２回以内に抑えることができる。従って、この場合、後処理工程をステップ１２５の判定無しに２回に固定しても、事実上良好な結果が得られる。
【００７４】
ステップ１２５における縮尺係数の判定を行わずに後処理工程を必ず行うものとした場合、通常、主工程１００のステップ１０５において、仮の数値として、原点縮尺係数ｍ0 を用いて計算４をした縮尺係数ｍt を採用する場合、入力した座標値ｘ，ｙがその座標系のどの位置にあっても、主工程１００において、ステップ１１０の判定による反復計算（ステップ１０８〜１１０）を繰り返し、後処理工程を１回（当該後処理工程中においてもテップ１１０の判定による反復計算（ステップ１０８〜１１０）を繰り返す。）行えば、十分な精度の結果が得られる。
また、ステップ１０６において、原点緯度φ0 に対して、１秒（−1/206265ラジアン）ずらすものとした。このずらす値については、原点緯度φ0 に対して±５秒を超えない範囲で行うのが好ましい。この範囲であれば、上記の通り、主工程１００のステップ１１０の判定による反復計算（ステップ１０８〜１１０）を最大３回とし後処理工程を１回として十分な精度が得られることが発明者によって確認されている。従って、このステップ１０８〜１１０の反復計算を３回を限度として実施するのが好ましい。但し、反復計算は、収束条件ＩＰＳに依存しているため、このような回数に限定するものではない。
【００７５】
平面上の点を楕円体上の点とを対応させるのには、近似値を求める計算（反復計算即ち収束計算）が必要であり、精度良く対応させるには本願発明のような工夫が必要となるが、逆に、楕円体上の点を平面上の点に正確に対応させるのは、対応する点を簡単に算出することができる（反復計算は不要である）。即ち、周知の計算方法や装置で、緯度・経度から平面直角座標上の座標値を逆算して正確な値を算出することができる。
従って、このような逆算によって、本願発明に係る方法で得られた計算結果の精度を検証することができる。
【００７６】
検証法の例について簡単に説明する。
本願発明に係る上記の方法にて得られた結果（出力値）を緯度φ（緯度φn ）、経度λとして、下記の計算を行う。
計算２２にて、経度差分Δλを逆算する。
Δλ＝λ−λ0 …計算２２
上記結果の緯度φから、
計算２３にて、第１の値ηを求める。
η＝ｅ’cos φ…計算２３
また、上記の結果の緯度φから、計算２４にて、第２の値ｔを求める。
ｔ＝tan φ…計算２４
更に、上記の結果の緯度φから、計算２５にて、Ｖの値を求める。
Ｖ＝（１＋ｅ'2cos2φ）1/2 …計算２５
求めた当該Ｖの値から、計算２６にて子午線曲率半径Ｍを、計算２７にて、卯酉線曲率半径Ｎを、更に計算２８から平均曲率半径Ｒを、夫々求める。
Ｍ＝ｃ／Ｖ3 …計算２６
Ｎ＝ｃ／Ｖ…計算２７
Ｒ＝ｃ／Ｖ2 …計算２８
上記結果の緯度φ及び経度λ、更に、上記にて求めた第１の値η、第２の値ｔとから、計算２９にて、当該緯度φ及び経度λにおける縮尺係数ｍを求める。
ｍ＝｛（１＋Δλ2cos2 φ（１＋η2 ）／２＋Δλ4cos4 φ（５−４ｔ2 ）／24）｝×0.9999…計算２９
そして、上記の結果の緯度φと、当該φについての上記子午線曲率半径Ｍ、卯酉線曲率半径Ｎ、平均曲率半径Ｒ、上記の経度差分Δλ、第１の値η、第２の値ｔ、縮尺係数ｍとから、計算３０にて、当該緯度φ及び経度λに対応する平面直角座標上の座標値ｙを求めることができる（計算２９の右辺の0.9999は、数値として原点縮尺係数ｍ0 と同じであるが、前述の通り、これは基準子午線上の縮尺係数を0.9999とするための一定の縮尺係数である）。
ｙ＝｛ΔλＮcos φ＋Δλ3 Ｎcos 3 φ（１−ｔ2 ＋η2 ）／６
＋Δλ5 Ｎcos 5 φ（５−１８ｔ2 ＋ｔ4 ＋１４η2 −５８η2 ｔ2 ）／120 ｝×ｍ…計算３０
【００７７】
そして、地球楕円体の長半径ａと、第１離心率ｅと、上記結果の緯度φと、座標値ｘと、前述の定数Ａ〜Ｆとから、計算３１にて、上記結果の緯度φに対する赤道からの子午線弧長Ｓ (φ )を求める。
Ｓ (φ )＝ａ（１−ｅ2 ）・｛Ａ・φ−（Ｂ／２）・（sin ２φ）
＋（Ｃ／４）・（sin ４φ）−（Ｄ／６）・（sin ６φ）
＋（Ｅ／８）・（sin ８φ）−（Ｆ／１０）・（sin10 φ）｝
…計算３１
上記で求めた緯度φに対する赤道からの子午線弧長Ｓ (φ )と、原点緯度φ0 に対する赤道からの子午線弧長Ｓ (φ0 ) とにより、計算３２にて、その差Ｂを求める。
Ｂ＝Ｓ (φ )−Ｓ (φ0 ) …計算３２
更に、この差Ｂと、上記の結果の緯度φと、当該φについての上記子午線曲率半径Ｍ、卯酉線曲率半径Ｎ、平均曲率半径Ｒ、上記の経度差分Δλ、第１の値η、第２の値ｔ、縮尺係数ｍとから、計算３３にて、にて、当該緯度φ及び経度λに対応する平面直角座標上の座標値ｘを求めることができる。
ｘ＝｛Ｂ＋Δλ2sinφcos φ／２＋Δλ4 Ｎsin φcos3φ（５−ｔ2 ＋９η2 ＋４η4 ）／24＋Δλ6 Ｎsin φcos 5 φ（61−58ｔ2 ＋ｔ4 ＋ 270η2 −330 η2 ｔ2 ）／720 ｝×ｍ…計算３３
また、上記の結果の緯度φと、上記の経度差分Δλ、第１の値η、第２の値ｔとから、計算３４にて、当該緯度φ及び経度λにおける子午線収差γを求めることができる。
γ＝Δλsin φ＋Δλ3 sin φcos2φ（１＋３η2 ＋２η4 ）／３
＋Δλ5 sin φcos4φ（２−ｔ2 ）／15…計算３４
【００７８】
本明細書の背景技術の欄で、国土地理院のプログラムによる計算に使用した値について、本願発明に係る計算装置に入力して計算し、その結果の値を上記の検算法で検算したところ、国土地理院のプログラムが設定してある有効数字の範囲内（算出値の桁数内）において、本願発明に係る計算装置に入力した値（ｘ，ｙ）と、検算で得た値（ｘ，ｙ）とは、完全に一致することが分った。
更に、具体的な実施例について、以下に掲げる。
【実施例】
【００７９】
本願発明に係る装置を用いて、座標系３に属する５箇所の座標値ｘ，ｙ（入力）から、夫々の緯度φ・経度λを算出（出力）した（結果データ１〜５）。ここでは、座標値ｘを固定し、結果データ１から結果データ５へかけて、座標値ｙが原点から遠ざかる点を調べた。結果データ５は、縮尺係数が、1.0001となる点である。
結果データ１〜５で得た緯度・経度から、上記の検算法（検算プログラム）にて、座標値についての検証データ１〜５を得た（検証データの番号は、検証した結果データの番号と対応している）。
各データの計算において、ＷＧＳ８４地球楕円体の長半径ａ、地球楕円体の極の曲率半径ｃ（＝ａ2 ／ｂ）、離心率（第１離心率）ｅの２乗、第２離心率ｅ’の２乗、座標系Ｚ、当該座標系Ｚの原点緯度φ0 、当該座標系Ｚの原点経度λ0 、原点縮尺係数ｍ0 について、次の数値を定数として用いた。
ａ＝6378137.000
ｃ＝6399593.62576
ｅ2 ＝ 0.006694379990141
ｅ’2 ＝0.006739496742276
Ｚ＝3
φ0 ＝36- 0- 0.000
λ0 ＝ 132-10- 0.000
ｍ0 ＝0.9999
また、以下において、
ZAHYOCHI X( 1)は、入力座標値ｘを示し、単位はメートルである。
ZAHYOCHI Y( 1)は、入力座標値ｙを示し、単位はメートルである。
PHAIは、最近似垂線緯度φj+1 を示す。
LATITUDEは、当該点の緯度φn （北緯）を示す。例えば 34- 7-45.93167 は、 34 度 7分45.93167秒（34°7'45".93167 ）を示している。
LONGITUDE は、当該点の経度λ（東経）を示す。
GAMMA ANGLEは、当該点の子午線収差（真北方向角）γを示す。
MERIDIAN RADIUSは、当該点の子午線曲率半径Ｍを示す。
ORTHOG. MERIDIAN RADIUS は、当該点の卯酉線曲率半径Ｎを示す。
MEAN RADIUS は、当該点の平均曲率半径Ｒを示す。
SCALE FACTOR は、当該点の縮尺係数を示す。
【００８０】
（結果データ１）
１ ZAHYOCHI X( 1) = -207504.143
ZAHYOCHI Y( 1) = 0.000
PHAI = 34- 7-45.93167
LATITUDE = 34- 7-45.93167
LONGITUDE =132-10- 0.00000
GAMMA ANGLE = 0- 0- 0.00000
MERIDIAN RADIUS = 6355518.63673
ORTHOG. MERIDIAN RADIUS =6384868.10519
MEAN RADIUS = 6370176.46816
SCALE FACTOR = 0.9999000000
【００８１】
（結果データ２）
2 ZAHYOCHI X( 2) = -207504.143
ZAHYOCHI Y( 2) = -60000.000
PHAI = 34- 7-46.23045
LATITUDE = 34- 7-40.02854
LONGITUDE = 131-30-58.32910
GAMMA ANGLE = - 0-21-53.81142
MERIDIAN RADIUS = 6355516.93662
ORTHOG. MERIDIAN RADIUS =6384867.53587
MEAN RADIUS = 6370175.33215
SCALE FACTOR = 0.9999443586
【００８２】
（結果データ３）
3 ZAHYOCHI X( 3) = -207504.143
ZAHYOCHI Y( 3) = -90091.612
PHAI = 34- 7-46.60518
LATITUDE = 34-7-32.62478
LONGITUDE = 131-11-24.23440
GAMMA ANGLE = -0-32-52.51513
MERIDIAN RADIUS = 6355514.80440
ORTHODG. MERIDIAN RADIUS = 6384866.82185
MEAN RADIUS = 6370173.90739
SCALE FACTOR = 1.0000000000
【００８３】
（結果データ４）
4 ZAHYOCHI X( 4) = -207504.143
ZAHYOCHI Y( 4) = -119361.021
PHAI = 34-7-47.111365
LATITUDE = 34-7-22.57909
LONGITUDE = 130-52-22.58341
GAMMA ANGLE = -0-43-32.98171
MERIDIAN RADIUS = 6355511.91142
ORTHOG. MERIDIAN RADIUS = 6384865.85306
MEAN RADIUS = 6370171.97429
SCALE FACTOR = 1.0000755078
【００８４】
（結果データ５）
5 ZAHYOCHI X( 5) = -207504.143
ZAHYOCHI Y( 5) = -127420.415
PHAI = 34-7-47.27856
LATITUDE = 34-7-19.32099
LONGITUDE = 130-47-8.30603
GAMMA ANGLE = -0-46-29.28380
MERIDIAN RADIUS = 6355510.97317
ORTHOG. MERIDIAN RADIUS = 6384865.53887
MEAN RADIUS = 6370171.34734
SCALE FACTOR = 1.0001000000
【００８５】
（検証データ１）
１ ZAHYOCHI X( 1) = -207504.143
ZAHYOCHI Y( 1) = 0.000
LATITUDE = 34-7-45.93167
LONGITUDE = 132-10-0.00000
GAMMA ANGLE = 0-0-0.00000
MERIDIAN RADIUS = 6355518.63672
ORTHODG. MERIDIAN RADIUS =6384868.10519
MEAN RADIUS = 6370176.46816
SCALE FACTOR = 0.9999000000
【００８６】
（検証データ２）
2 ZAHYOCHI X( 2) = -207504.143
ZAHYOCHI Y( 2) = -60000.000
LATITUDE = 34-7-40.02854
LONGITUDE = 131-30-58.32910
GAMMA ANGLE = -0-21-53.81142
MERIDIAN RADIUS = 6355516.93662
ORTHODG. MERIDIAN RADIUS =6384867.53587
MEAN RADIUS = 6370175.33215
SCALE FACTOR = 0.9999443586
【００８７】
（検証データ３）
3 ZAHYOCHI X( 3) = -207504.143
ZAHYOCHI Y( 3) = -90091.612
LATITUDE = 34-7-32.62478
LONGITUDE = 131-11-24.23440
GAMMA ANGLE = -0-32-52.51514
MERIDIAN RADIUS = 6355514.80440
ORTHODG. MERIDIAN RADIUS = 6384866.82185
MEAN RADIUS = 6370173.90739
SCALE FACTOR = 0.9999999999
【００８８】
（検証データ４）
4 ZAHYOCHI X( 4) = -207504.143
ZAHYOCHI Y( 4) = -119361.021
LATITUDE = 34-7-22.57909
LONGITUDE = 130-52-22.58341
GAMMA ANGLE = -0-43-32.98172
MERIDIAN RADIUS = 6355511.91142
ORTHOG. MERIDIAN RADIUS = 6384865.85306
MEAN RADIUS = 6370171.97429
SCALE FACTOR = 1.0000755076
【００８９】
（検証データ５）
5 ZAHYOCHI X( 5) = -207504.143
ZAHYOCHI Y( 5) = -127420.415
LATITUDE = 34-7-19.32099
LONGITUDE = 130-47-8.30603
GAMMA ANGLE = -0-46-29.28381
MERIDIAN RADIUS = 6355510.97317
ORTHOG. MERIDIAN RADIUS= 6384865.53887
MEAN RADIUS = 6370171.34734
SCALE FACTOR = 1.0000999998
【００９０】
上記の通り、何れも、結果データの入力値と、検証データ出力値は、完全に一致している。ここで求めた緯度・経度は、小数点以下第５位であり、この範囲で誤差が生じている国土地理院のプログラムとでは、その精度の差は歴然としている。
また、結果データ５を見れば、縮尺係数を示す SCALE FACTORも、1.0001に、精度よく収束しているのが分かる。
【図面の簡単な説明】
【００９１】
【図１】本願発明に係る位置計算装置の一実施の形態の説明図である。
【図２】本願発明に係る位置計算方法のフローを示す説明図である。
【図３】本願発明に係る位置計算方法の、垂線緯度の計算についての説明図である。
【符号の説明】
【００９２】
１ 座標入力部
２ 主力部
６ 演算部
７ 演算制御部
８ 判定部
９ 回帰制御部
【技術分野】
【０００１】
本願発明は、位置計算装置、位置計算方法、位置計算プログラム、及び位置計算プログラムの記憶媒体に関するものである。
【背景技術】
【０００２】
地球は、楕円体（地球楕円体）であり、地理学的経緯度を以って、地球上の位置を特定することができる。しかし、取り扱いの便から、一般的な測量等において、平面直角座標が利用され、座標によって、位置の特定がされる（測量法第11条）。この座標平面については、回転楕円体（地球楕円体）から、直接平面に投影するガウス・クリューゲルの等角投影（横メルカトール投影）を採用して、原点における縮尺係数を0.99990 としている。
尚、地球楕円体を測量の基準にするためには、楕円体の中心を実際の地球上のどの位置に、またその楕円体の座標軸が実際の地球のどこを通るかということを決める必要があり、この位置と方向が決められた地球楕円体が準拠楕円体とされる。
【０００３】
地球上の点の水平位置は、厳密には準拠楕円体上の地理学的経緯度によって表されるべきであるが、位置・方向・距離等を、上記の通り平面上に投影して測量計算を行うことは曲面上に比べ非常に簡単になり便利である。また、公共測量のように測量範囲が狭い場合には、十分正確に表すことができるとされている。
日本で用いられている平面直角座標は、ガウス・クリューゲルの等角投影法によるもので、座標系原点を通る子午線は等長に、図形は等角の相似形に投影される。しかし、距離については、基準子午線から東西に離れるに従って平面距離が増大していくため、投影距離の誤差を相対的に1/10,000以内に収めるよう座標原点に上記の縮尺係数（0.9999）を与え、かつ、基準子午線より東西130km 以内を適用範囲とした座標系（平面直角座標系は日本全国を１９の座標系に区分している。）を設けている（非特許文献１）。
なお、座標系のＸ軸は、原点において子午線に一致する軸とし、原点から真北に向かう値を正としている。Ｙ軸は原点においてＸ軸に直交する軸とし、真東に向かう値を正としている。この点、三次元直交座標（Ｘ，Ｙ, Ｚ）とは定義が異なる。又、各座標系原点の値はＸ＝Ｙ＝0.000 メートルとしている。
【０００４】
測量結果の成果表において、距離については楕円体面上の距離が記載され、座標については緯度・経度のほか、平面直角座標（ｘ，ｙ）が記載される。楕円体面と平面の関係については、準拠楕円体に被せられた楕円体筒（仮想筒）を広げて平面座標（ｘ，ｙ）を得る。この場合、基準子午線における平面距離は、楕円体上の距離より１万分の１だけ短くなり、縮尺係数は、0.9999である。基準子午線から東西に約90ｋｍ離れた場所では、平面と楕円体面の距離が等しくなり、縮尺係数は1.0000である。基準子午線から東西に130 ｋｍ離れると、平面距離は、楕円体面距離より、１万分の１だけ長くなり、縮尺係数は、1.0001である。
【０００５】
わが国土について、平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙから、その緯度・経度を算出する手段として、国土地理院のホームページにおいて、「緯度・経度への換算」（非特許文献３）のプログラムが利用できる。
このプログラムを利用することにより、ｗｅｂ上において、上記座標値ｘ，ｙの入力により、簡便に、その緯度・経度を算出することができる。
ところが、このプログラムにおいて、座標値ｘ，ｙが属する座標系の基準子午線から、東西に離れるに従って、その誤差が大きくなる。
この誤差は、国土地理院ホームページの当該「緯度・経度への換算」プログラムで得た緯度・経度を、同じく国土地理院ホームページ「平面直角座標への換算」（非特許文献４）プログラムにて、簡単に確認することができる。
【０００６】
例えば、国土地理院ホームページ「緯度・経度への換算」（非特許文献３）において、測地系を世界測地系とし、平面直角座標系を３として、上記「平面直角座標から緯度・経度への換算」プログラムへ、平面直角座標上の任意の５点の座標値Ｘ，Ｙを与えて、ａ−１）〜ａ−５）に示す結果（緯度・経度）を得た。
ａ−１）は当該座標系３の座標原点（０，０）を通る子午線上の点であり、ａ−２）〜ａ−５）において、入力する座標値Ｙの値を基準子午線から西方へ徐々に遠ざかるものとしている。
そして上記「平面直角座標への換算」プログラムにおいて、上記ａ−１）〜ａ−５）で入力した座標値と、換算結果が一致するように、試行錯誤で北緯と東経を入力した（北緯と東経の値を種々変えて、出力値が、「平面直角座標への換算」プログラム）における入力値と一致するものを探った）。そして、ｂ−１）〜ｂ−５）に示す結果（平面直角座標上の座標値）を得た。
【０００７】
ａ−１）
入力値 Ｘ座標 -207504.143 ｍ Ｙ座標 0 ｍ
換算結果 緯度 34°07'45.93167" 経度 132°10'00.00000"
真北方向角＋ 0°00'00.00" 縮尺係数 0.99990000
ａ−２）
入力値 Ｘ座標 -207494.938 ｍ Ｙ座標 -59997.338 ｍ
換算結果 緯度 34°07'40.02853" 経度 131°30'58.32911"
真北方向角＋ 0°21'53.81" 縮尺係数 0.99994436
ａ−３）
入力値 Ｘ座標 -207483.393 ｍ Ｙ座標 -90082.603 ｍ
換算結果 緯度 34°07'32.62476" 経度 131°11'24.23439"
真北方向角＋ 0°32'52.52" 縮尺係数 1.00000000
ａ−４）
入力値 Ｘ座標 -207467.727 ｍ Ｙ座標 -119340.074 ｍ
換算結果 緯度 34°07'22.57910" 経度 130°52'22.58340"
真北方向角＋ 0°43'32.98" 縮尺係数 1.00007551
ａ−５）
入力値 Ｘ座標 -207462.646 ｍ Ｙ座標 -127394.933 ｍ
換算結果 緯度 34°07'19.32100" 経度 130°47'08.30605"
真北方向角＋ 0°46'29.28" 縮尺係数 1.00010000
【０００８】
ｂ−１）
入力値 北緯 34°07'45.93167" 東経 132°10'00.00000"
換算結果 Ｘ座標 -207504.143 ｍ Ｙ座標 0.0000 ｍ
真北方向角＋ 0°00'00.00" 縮尺係数 0.99990000
ｂ−２）
入力値 北緯 34°07'40.02854" 東経 131°30'58.32910"
換算結果 Ｘ座標 -207494.938 ｍ Ｙ座標 -59997.338ｍ
真北方向角＋ 0°21'53.81" 縮尺係数 0.99994436
ｂ−３）
入力値 北緯 34°07'32.62478" 東経 131°11'24.23440"
換算結果 Ｘ座標 -207483.393 ｍ Ｙ座標 -90082.603ｍ
真北方向角＋ 0°32'52.52" 縮尺係数 1.00000000
ｂ−４）
入力値 北緯 34°07'22.57909" 東経 130°52'22.58341"
換算結果 Ｘ座標 -207467.727 ｍ Ｙ座標 -119340.074 ｍ
真北方向角＋ 0°43'32.98" 縮尺係数 1.00007551
ｂ−５）
入力値 北緯 34°07'19.32099" 東経 130°47'08.30603"
換算結果 Ｘ座標 -207462.646 ｍ Ｙ座標 -127394.933 ｍ
真北方向角＋ 0°46'29.28" 縮尺係数 1.00010000
【０００９】
上記の結果を見ると、ａ−１）の入力値と一致する、出力値が得られるようにｂ−１）において探って得た入力値は、ａ−１）における換算結果と一致している（上記ａ−１）における換算結果の緯度34°07'45.93167"と経度 132°10'00.00000"は、ｂ−１）における入力値の北緯・東経と一致している）。即ち、当該原点を通る基準子午線上において、「平面直角座標への換算」プログラムは、少なくとも計算可能な桁数内において、誤差を生じていないといえる。
ところが、ａ−２）について、上記と同様の作業を行ったｂ−２）を見ると、ａ−２）における換算結果と、ｂ−２）における入力値とは一致していない。即ち、前者ａ−２）が緯度を34°07'40.02853"とし経度を 131°30'58.32911"とするのに対して、後者ｂ−２）が北緯を34°07'40.02854"とし東経を 131°30'58.32910"とし、緯度及び経度について、夫々小数点以下第５位（最下位）において食い違いが生じている。
このような誤差の発生について、ａ−３〜ａ−５）と対応するｂ−３）〜ｂ−５）と見れば、基準子午線から西方へ離れるに従って、特に経度についての当該誤差が大きくなっていることが分かる。
例えば、基準子午線から西方へ約１２７ｋｍ離れた地点であるａ−５）における換算結果（緯度34°07'19.32100"、経度 130°47'08.30605"）と、対応するｂ−５）における入力値（北緯34°07'19.32099"、東経 130°47'08.30603"）とでは、特に、緯度について小数点以下第３位にて既に食い違いが生じている。
【００１０】
【非特許文献１】国土地理院ホームページ http://vldb.gsi.go.jp/sokuchi/datum/tokyodatum.html
【非特許文献２】国土地理院ホームページ「赤道からの子午線弧長を与えて緯度を求める計算」http://vldb.gsi.go.jp/sokuchi/surveycalc/algorithm/s2b/s2b.htm
【非特許文献３】国土地理院ホームページ「緯度・経度への換算」http://vldb.gsi.go.jp/sokuchi/surveycalc/xy2blf.html
【非特許文献４】国土地理院ホームページ「平面直角座標への換算」http://vldb.gsi.go.jp/sokuchi/surveycalc/bl2xyf.html
【００１１】
従って、このプログラムでは、入力座標値Ｘ，Ｙが属する座標系原点の基準子午線付近においては、問題がないものの、基準子午線から東西へ離れるに従って、換算結果である緯度及び経度の精度が著しく低下し、厳密なデータとして取り扱うことができなくなる。
上記において、そのような誤差が生じる原因は、プログラムの内容が公開されていないことから、本願発明者において正確なことは分からない。
また、発明者にて可能な範囲で、先行技術文献として、特許公報を調べたが、同種の発明が記載された公報を発見するに至っていない。
そこで、発明者は、この点を熟考し、推測するに、上記国土地理院のプログラムにおいて、緯度について例えば、国土地理院ホームページ「赤道からの子午線弧長を与えて緯度を求める計算」（非特許文献２）などに示される公式を用いて、緯度についての収束計算を行っているのであろうが、基準子午線から離れた場合の前述の縮尺係数の変化が十分に反映するものとなっておらず、その結果、収束（反復計算）により安定状態となった緯度の精度が十分なものでない（ずれた位置に収束してしまっている）と考える。
具体的には、国土地理院の上記プログラムは、非特許文献２に示す国土地理院ホームページの「６．赤道からの子午線弧長を与えて緯度を求める計算」で示されているφn+1 を求める式を用いて反復計算を行うものであると考えられる（この式を用いた結果が上記と一致した）。この式を用いて反復計算を行うと上述の通り、大きな誤差を持った値で数値が安定（収束）する。その理由を検討すると、この式の計算において、縮尺係数は、原点における縮尺係数ｍ0 （＝0.9999）に固定されており、基準子午線以外の、求めようとする座標位置における縮尺係数は、原点の縮尺係数ｍ0 （＝0.9999）と異なるものであるという事実は、無視されている（尚、国土地理院の式で用いているＭは、赤道からの座標系原点までの子午線弧長を示し、この明細書で用いる子午線曲率半径を示すＭとは異なる）。
現状では、推測の域を出るものではないが、上記プログラムを用いた結果は、歴然としている。
従って、当該プログラムは、著しく信頼性の欠くものとなっており、厳密性が要求される状況において、その利用が困難となっている。
【発明の開示】
【発明が解決しようとする課題】
【００１２】
本願発明は、等角投影法により地球楕円体を投影する平面直角座標Ｘ，Ｙ上の座標値から地球楕円体上の少なくとも緯度を求めることが可能な位置計算装置について、縮尺係数の取り扱いを柔軟に行うことにて、誤差の低い新規な装置を提供することにより、上記課題の解決を図った。
【課題を解決するための手段】
【００１３】
本願第１の発明は、等角投影法により地球楕円体を投影する平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙを入力することにより、少なくとも当該座標値に対応する地球楕円体上の緯度φn を出力することが可能な位置計算装置として、次の構成を備えたものを提供する。
即ち、この位置計算装置は、座標入力部と、出力部と、演算部と、判定部と、回帰制御部とを備える。座標入力部は、上記座標値ｘ，ｙの入力を受け付ける。演算部は、直角座標の座標値ｘが与えられることにより、地球楕円体の赤道から当該座標値ｘ，ｙが属する座標系原点までの子午線弧長に対し、当該座標値ｘを仮の縮尺係数ｍt で除した値を、加算して、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとする。
そして、この子午線弧長Ｗに対応する基準子午線上の緯度、即ち、上記入力座標から基準子午線上に下ろした垂線の緯度を算出するため、演算部は、地球楕円体の長半径と第１離心率と緯度とから定まる子午線曲率半径について、地球楕円体の長半径ａと第１離心率ｅと共に任意の緯度φj を初期値として与え、別途に子午線曲率半径を以って赤道から上記緯度φj まで積分することにより、赤道からの基準子午線上の子午線弧長Ｓ（φj ）を算出する。更に、演算部は、基準子午線上の上記緯度φj を与えることによって、式１を計算して基準子午線上の緯度φj+1 を漸近値として算出する。
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）…式１
判定部は、上記式１にて算出した基準子午線上の上記漸近緯度φj+1 と基準子午線上の緯度φj との差の絶対値と、基準値との比較を行うものである。判定部にて、基準子午線上の、上記漸近緯度φj+1 と緯度φj との差の絶対値が基準値を超えると判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするには十分収束していないと判定した場合、回帰制御部は、この漸近緯度φj+1 を上記基準子午線上の緯度φj とし、演算部に基準子午線上の赤道からの子午線弧長Ｓ（φj ）を算出させると共に当該子午線弧長Ｓ（φj ）から上記式１の計算にて再度漸近緯度φj+1 を算出させ、判定部に再度上記比較を行わせる。また、判定部が、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするのに十分収束したと判定した場合、回帰制御部は、当該漸近緯度φj+1 を、入力座標から基準子午線上に下ろした垂線の緯度φk として、仮の縮尺係数ｍt と、当該垂線緯度φk と、当該垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式３を計算して緯度φk+1 を演算部に算出させる。
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …式３
回帰制御部は、上記緯度φk+1 を入力座標に対応する緯度φn として、演算部に、縮尺係数ｍt+1 を算出させ、前記仮の縮尺係数ｍt に代え当該縮尺係数ｍt+1 を用いて求めた、基準子午線上の子午線弧長Ｗにて、上記式１の計算を行わせる。
ここでいう地球楕円体には、各種の準拠楕円体を含む。
基準子午線とは、当該座標が属する座標系の原点を通る子午線をいう。
また、垂線緯度とは、入力座標から基準子午線に下ろした垂線と、当該基準子午線との交点の緯度をいう。
また、積分には、その近似計算を含む。
上記において、回帰制御部は、演算部に対して直接、計算の指令を下すものに限定するものではなく、間接的に、演算部へ計算を行わせるものも含む。
【００１４】
本願第２の発明は、上記本願第１の発明にあって、次の特徴を備えた位置計算装置を提供する。即ち、上記の仮の縮尺係数ｍt とは、演算部に、直角座標の座標値ｙを与えることにより、原点における縮尺係数ｍ0 又はその代替値と、当該座標が属する座標系の原点における平均曲率半径Ｒ0 とから、式２の計算により、演算部が算出した縮尺係数ｍt である。
ｍt ＝｛１＋（ｙ2 ／２ｍ02Ｒ02）＋（ｙ4 ／24ｍ04Ｒ04）｝×0.9999…式２
上記の原点における平均曲率半径Ｒ0 は、当該座標が属する座標系の原点の緯度φ0 又はその代替値と共に、地球楕円体の長半径ａ及び短半径ｂ、或いは当該半径ａ，ｂから導かれる第１離心率ｅ、第２離心率ｅ’を用いて定めたものである。
本願第３の発明は、上記本願第２の発明にあって、次の特徴を備えた位置計算装置を提供する。即ち、上記縮尺係数ｍt+1 は、原点縮尺係数ｍ0 に代え上記縮尺係数ｍt を用い、且つ、原点における平均曲率半径Ｒ0 に代え演算部に式３にて算出させた緯度φn における平均曲率半径Ｒφn を用いて、演算部が、上記式２により算出するものである。
【００１５】
本願第４の発明は、上記本願第１乃至３の何れかの発明にあって、次の特徴を備えた位置計算装置を提供する。即ち、上記にあって、判定部が、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合に、回帰制御部は、上記緯度φn に対応する縮尺係数ｍt+1 と上記仮の縮尺係数ｍt との差の絶対値と、所定値との比較を、更に判定部にさせ、判定部が当該縮尺係数ｍt+1 と上記仮の縮尺係数ｍt との差の絶対値が所定値を超えると判定した場合にのみ、上記の通り、前記の仮の縮尺係数ｍt に代え、当該縮尺係数ｍt+1 を用いて求めた基準子午線上の子午線弧長Ｗにて、上記式１の演算部に計算を行わせるものである。
【００１６】
本願第５の発明は、上記本願第１乃至４の何れかの発明にあって、次の特徴を備えた位置計算装置を提供する。即ち、上記にあって、判定部が、基準子午線上の、緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、回帰制御部は、垂線緯度φk として上記漸近緯度φj+1 を演算部に与えるものであり、演算部は、上記仮の縮尺係数ｍt と、上記垂線緯度φk における卯酉線曲率半径Ｎk と、当該垂線緯度φk と第２離心率とによりｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式４にて、経度差分Δλを算出することが可能であり、
Δλ＝（ｙ／ｍt ）／Ｎk cos φk
−（ｙ／ｍt ）3 （１＋２ｔk 2 ＋ηk 2 ）／６Ｎk 3cosφk
＋（ｙ／ｍt ）5 （５＋２８ｔk 2 ＋２４ｔk 4 ）／１２０Ｎk 5cosφk …式４
演算部は、当該座標系が属する原点における経度λ0 に対し当該経度差分Δλを加算することにより、上記座標値ｘ，ｙにおける経度λを算出するものである。
【００１７】
本願第６の発明は、上記本願第１乃至５の何れかの発明にあって、次の特徴を備えた位置計算装置を提供する。即ち、上記にあって、判定部が、基準子午線上の、緯度φj+1 と緯度φj との差の絶対値を基準値以内と判定した場合、回帰制御部は、垂線緯度φk として緯度φj+1 を演算部に与えるものであり、演算部は、上記の縮尺係数ｍt と、当該緯度φk における、子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、当該垂線緯度φk と第２離心率とによりｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式５にて、座標ｘ，ｙにおける子午線収差γを算出するものである。
γ＝（ｙ／ｍt ）ｔk ／Ｎk
−（ｙ／ｍt ）3 ｔk （１＋ｔk 2 −ηk 2 ）／３Ｎk 3
＋（ｙ／ｍt ）5 ｔk （２＋５ｔk 2 ＋３ｔk 4 ）／１５Ｎk 5
…式５
【００１８】
本願第７の発明は、等角投影法により地球楕円体を投影する平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙを入力することにより、少なくとも当該座標値に対応する地球楕円体上の緯度φn を出力することが可能な位置計算プログラムについて、次の特徴を備えたものを提供する。
即ち、この位置計算プログラムは、座標入力部と、出力部と、演算部と、判定部と、回帰制御部とを備えたコンピュータに導入される。このプログラムが導入されたコンピュータは、座標入力部にて、上記座標値ｘ，ｙの入力を受け付け、演算部に、直角座標の座標値ｘを与えて、地球楕円体の赤道から当該座標値ｘ，ｙが属する座標系原点までの子午線弧長に対し、当該座標値ｘを仮の縮尺係数ｍt で除した値を加算して、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとする。
そして、この子午線弧長Ｗに対応する基準子午線上の緯度、即ち、上記入力座標値から基準子午線上に下ろした垂線の緯度を算出するため、演算部にて、地球楕円体の長半径と第１離心率と緯度とから定まる子午線曲率半径について、地球楕円体の長半径ａと第１離心率ｅと共に任意の緯度φj を初期値として与え、別途に子午線曲率半径を以って赤道から上記緯度φj まで積分することにより、赤道からの基準子午線上の子午線弧長Ｓ（φj ）を算出し、更に、演算部に、基準子午線上の緯度φj を与えることによって、式１を計算して基準子午線上の漸近緯度φj+1 を漸近値として算出させる。
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）…式１
そして、判定部に、上記式１にて算出した基準子午線上の漸近緯度φj+1 と基準子午線上の緯度φj との差の絶対値と、基準値との比較を行わせる。判定部にて、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値を超えると判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするには十分収束していないと判定した場合、回帰制御部により、漸近緯度φj+1 を上記基準子午線上の緯度φj とし、演算部に基準子午線上の赤道からの子午線弧長Ｓ（φj ）を算出させると共に当該子午線弧長Ｓ（φj ）から上記の式１の計算にて再度漸近緯度φj+1 を算出させ、判定部に再度上記比較を行わせるものである。また、判定部にて、基準子午線上の漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線とするのに十分収束したと判定した場合、当該漸近緯度φj+1 を入力座標からおろした垂線の緯度φk として、回帰制御部により、仮の縮尺係数ｍt と、当該垂線緯度φk と、当該垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式３を計算して緯度φk+1 を演算部に算出させる。
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …式３
更に、回帰制御部により、上記緯度φk+1 を入力座標に対応する緯度φn として、演算部に、縮尺係数ｍt+1 を算出させ、前記仮の縮尺係数ｍt に代え当該縮尺係数ｍt+1 を用いて求めた、基準子午線上の子午線弧長Ｗにて、上記式１の計算を行わせるものである。
【００１９】
本願第８の発明は、上記第７の発明に係る位置計算プログラムの記録媒体を提供する。
【００２０】
本願第９の発明は、等角投影法により地球楕円体を投影する平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙを入力することにより、少なくとも当該座標値に対応する地球楕円体上の緯度φn を出力することが可能な位置計算方法について、次の特徴を備えた方法を提供する。
即ち、この方法は、座標入力部と、出力部と、演算部と、判定部と、回帰制御部とを備えたコンピュータを用いるものである。そして、座標入力部にて、上記座標値ｘ，ｙの入力を受け付ける。演算部に、直角座標の座標値ｘを与えて、地球楕円体の赤道から当該座標値ｘ，ｙが属する座標系原点までの子午線弧長に対し、当該座標値ｘを仮の縮尺係数ｍt で除した値を加算して、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとし、演算部にて、地球楕円体の長半径と第１離心率と緯度とから定まる子午線曲率半径について、地球楕円体の長半径ａと第１離心率ｅと共に任意の緯度φj を初期値として与え、別途に子午線曲率半径を以って赤道から上記緯度φj まで積分することにより、赤道からの基準子午線上の子午線弧長Ｓ（φj ）を算出し、演算部に、基準子午線上の緯度φj を与えることによって、式１を計算して基準子午線上の漸近緯度φj+1 漸近値として算出させる。
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）…式１
判定部にて、上記式１にて算出した基準子午線上の漸近緯度φj+1 と基準子午線上の緯度φj との差の絶対値と、基準値との比較を行わせるものである。そして、判定部にて、上記基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値を超えると判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線に下ろした垂線の緯度とするには十分収束していないと判定した場合、回帰制御部により、漸近緯度φj+1 を上記基準子午線上の緯度φj とし、演算部に基準子午線上の赤道からの子午線弧長Ｓ（φj ）を算出させると共に当該子午線弧長Ｓ（φj ）から上記の式１の計算にて再度漸近緯度φj+1 を算出させ、判定部に再度上記比較を行わせる。判定部にて、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線に下ろした垂線の緯度とするのに十分収束したと判定した場合、回帰制御部は、当該漸近緯度φj+1 を、入力座標から基準子午線上に下ろした垂線の緯度φk として、仮の縮尺係数ｍt と、上記垂線緯度φk と、当該垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式３を計算して緯度φk+1 を演算部に算出させる。
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …式３
そして、回帰制御部により、上記緯度φk+1 を入力座標に対応する緯度φn として、演算部に、縮尺係数ｍt+1 を算出させ、前記仮の縮尺係数ｍt に代え当該縮尺係数ｍt+1 を用いて求めた、基準子午線上の子午線弧長Ｗにて、上記式１の計算を行わせるものである。
【発明の効果】
【００２１】
上記本願の各発明は、直角座標上の座標点の入力により、自動的に、少なくとも当該緯度を算出するものであって、当該緯度計算中、必要とする垂線緯度の反復計算（収束）において、緯度に対応して縮尺係数の収束も反映させることが可能な装置・方法・プログラムやその記録媒体を提供し得た。これによって、極めて精度の高い、緯度計算を行うことができる。
即ち、この発明では、入力座標から基準子午線上に下ろした垂線の緯度（垂線緯度）を算出するに当たり、緯度の収束計算に適した式１をイテレーション式として用いて緯度の計算（反復計算）を行い、得た垂線緯度から式３を用いて当該点（入力座標値に対応する）緯度を算出し直すことにより、より精度の高い緯度を得ることができる。更に、このようにして求めた緯度から縮尺係数を算出して当該縮尺係数の収束を加味した後処理を行うことにより、極めて精度の高い当該点の緯度の算出を可能とした。
また、特に、本願第４の発明では、上記の計算を可能としつつも、基準子午線に近い場合など十分に誤差の小さい場合にまで、後処理として、式３を経て得た上記の緯度で算出した縮尺係数による、式１の計算を強要するものではなく、縮尺係数の収束状況で、後処理を行うか否かを決定し、無駄な反復計算を抑制し、迅速に計算を完了する手段を備えた装置を提供し得た。
更に本願第５の発明では、上記の当該点の緯度と共に、経度λを、また本願第６の発明では、子午線収差を、夫々算出し得る位置計算装置を提供した。
【発明を実施するための最良の形態】
【００２２】
以下、本願発明の好ましい実施の形態について、説明する。
【００２３】
１．本願発明に係る計算方法の概要
本願発明は、平面直角座標上の座標ｘ，ｙから、対応する地球楕円体上の緯度、更には、経度、子午線収差を精度よく算出する位置計算方法を提供するものである。即ち、ガウス・クリューゲルの等角投影法により地球楕円体を投影する平面直角座標の座標点の座標値から、当該地球楕円上の対応する点の緯度、経度及び子午線収差を、高精度に算出する方法を提供するものである。
この方法は、垂線緯度の算出のための収束計算において、縮尺係数の収束を加味して、計算結果の精度を向上させたものである。
【００２４】
即ち、この方法は、主工程において、下記の式１にて垂線緯度の反復計算を行い、反復計算後の値から下記の式３を計算して緯度を算出する。そして、式３を経て得た値を用いて下記の式２にて縮尺係数を算出し、算出した縮尺係数から導いた基準子午線上の子午線弧長を用いて再度式１にて垂線緯度の反復計算を、更に、再度の式１による反復計算後の値にて式３の計算を、後処理工程として行うものである。
また、上記後処理工程に先立ち、算出した縮尺係数が収束したかを調べ、その結果当該縮尺係数が収束していた場合、後処理工程を行わずに、計算を終了することもできる。
従って、主工程の処理に続いて後処理工程を行い上記にて算出した縮尺係数を加味した結果を得るものとしてもよく、また、主工程の結果に対して、算出した縮尺係数の収束の判定を行い、当該縮尺係数が収束していない場合にのみ、上記の後処理工程を行うものとしてもよい。
このように、主工程後必ず後処理工程を行う場合も、主工程後の縮尺係数の収束の判定により後処理工程を行う場合も、何れも、垂線緯度の収束計算において、縮尺係数の収束を加味したものとなり、緯度の計算結果について必要な精度を確保することができる。尚、縮尺係数の収束を調べて後処理工程を行うか否かの判断を行う場合、無駄な反復計算（後処理工程）を行う必要がない。例えば、基準子午線付近では、主工程のみにて緯度及び経度の十分な精度を確保している場合がある。このような場合にまで後処理工程を行う必要がないからである。
但し、近年のコンピュータの性能の向上により、主工程後、必ず後処理工程を行うものとしても、実用上、計算に要する時間の増加や、ハードウエアの負荷による問題はない。
また、後処理工程は、１回に限らない。更に、上記と同様、後処理工程を行う毎に縮尺係数の収束の判定を行い、その結果により更に後処理工程を行うか否か判断するものとしても実施可能である。
【００２５】
以下、主工程後の縮尺係数の収束の判定により後処理工程を行うものについてその概要を説明する。
この位置計算方法は、上記座標値ｘ，ｙの入力を受け付ける。そして当該座標値ｘが与えられることにより、地球楕円体の赤道から当該座標値ｘ，ｙが属する座標系原点までの子午線弧長に対し、当該座標値ｘを０以外の仮の縮尺係数ｍt で除した値を、加算して、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとする。また、地球楕円体の長半径と第１離心率と緯度とから定まる子午線曲率半径について、地球楕円体の長半径ａと第１離心率ｅと共に任意の数値（緯度φj ）を初期値として与え、別途に子午線曲率半径を以って赤道から上記緯度φj まで積分（近似計算）することにより、赤道からの基準子午線上の子午線弧長Ｓ（φj ）を算出する。そして、基準子午線上の上記漸近緯度φj から、式１を計算して基準子午線上の緯度φj+1 を漸近値（漸近緯度φj+1 ）として算出する。
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）…式１
上記式１にて算出した基準子午線上の緯度φj+1 と基準子午線上の緯度φj との差の絶対値と、基準値との比較を行い、上記基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値を超えると判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするには十分収束していないと判定した場合、この漸近緯度φj+1 を上記基準子午線上の緯度φj とし、基準子午線上の赤道からの子午線弧長Ｓ（φj ）を算出すると共に当該子午線弧長Ｓ（φj ）から上記式１の計算にて再度漸近緯度φj+1 を算出し、再度上記比較を行う。その結果、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、即ち、この漸近緯度φj+1 （以下必要に応じて最近似垂線緯度と呼ぶ。）を、上記入力座標から基準子午線上に下ろした垂線とするのに十分収束したと判定した場合、当該漸近緯度φj+1 を、入力座標から基準子午線上に下ろした垂線の緯度φk として、仮の縮尺係数ｍt と、当該垂線緯度φk と、当該垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式３を計算して緯度φk+1 を算出する。
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …式３
上記にて求めた緯度φk+1 を、当該（入力座標値に対応する）緯度φn として縮尺係数ｍt+1 を算出し、縮尺係数ｍt+1 と上記仮の縮尺係数ｍt との差の絶対値と、所定値との比較を行い（以上主工程）、当該縮尺係数ｍt+1 と上記仮の縮尺係数ｍt との差の絶対値が所定値を超えると判定した場合に、当該縮尺係数ｍt+1 を用いて求めた、基準子午線上の子午線弧長Ｗにて、上記式１の計算をやり直す（後処理工程）。
その結果得られた緯度を計算結果として出力する。また、緯度の計算の過程にて、経度及び子午線収差を計算する。
上記の式１をイテレーション式（反復式）として行う反復計算は、図３に示す平面座標において、基準子午線（Ｘ軸）上の漸近緯度φj+1 を、入力座標点（ｘ，ｙ）から基準子午線（Ｘ軸）に下ろした垂線ｈと、基準子午線（Ｘ軸）との交点Ｇに近似させる（近づける）計算である。基準子午線上において縮尺係数は0.9999である。式１による反復計算により、漸近緯度φj+1 を上記交点Ｇと取り扱うに十分近接した場合、当該漸近緯度φj+1 を最近似垂線緯度とし、式３の計算により、座標値（ｘ，ｙ）に対応する緯度φk+1 を算出するのである。
また、上記の仮の縮尺係数とは、上記直角座標の座標値ｙと、原点における縮尺係数ｍ0 又はその代替値と、当該座標が属する座標系の原点における平均曲率半径Ｒ0 とから、式２の計算により、算出した縮尺係数ｍt である。
ｍt ＝｛１＋（ｙ2 ／２ｍ02Ｒ02）＋（ｙ4 ／24ｍ04Ｒ04）｝×0.9999…式２
上記の原点における平均曲率半径Ｒ0 は、当該座標が属する座標系の原点の緯度φ0 又はその代替値と共に、地球楕円体の長半径ａ及び短半径ｂ、或いは当該半径ａ，ｂから導かれる第１離心率ｅ、第２離心率ｅ’を用いて定めたものである。
上記の仮の縮尺係数には、式２にて算出するものに限定するものではなく、０以外の任意の数値を採用することが可能である。但し、無駄な後処理工程を減ずる観点から、上記式２を用いて算出したものを採用するのが好ましい。
また上記において、式２を用いる場合も、原点における縮尺係数ｍ0 の代替値として０以外の任意の数値を用いることが可能であるが、上記と同様の観点から、このような代替値ではなく、当該原点における縮尺係数ｍ0 を用いて計算を行うのが好ましい。
また、上記の上記縮尺係数ｍt+1 は、原点縮尺係数ｍ0 に代え上記縮尺係数ｍt を用い、且つ、原点における平均曲率半径Ｒ0 に代え式３にて算出した緯度φn における平均曲率半径Ｒφn （後述のＲk2）を用いて、上記式２により算出するものである（言い換えると、式２において、ｍt をｍ0 に代入し、Ｒφn をＲ0 に代入して計算する）。
【００２６】
上記の通り、主工程は、式１にて垂線緯度を反復収束計算させ、その後式３、式４（式１９）、式５、及び式２を経ることにより、緯度、経度、子午線収差角及び縮尺係数を算出する工程である。垂線緯度の反復収束計算は、上記の式１にて、他の値は変更せずに、漸近緯度φj+1 を緯度φj とし直して新たに漸近緯度φj+1 を求め、当該計算を｜φj+1 −φj ｜が収束するまで行う。）。そして、主工程で｜φj+1 −φj ｜が収束した際、即ち当該漸近緯度φj+1 が最近似垂線緯度となった際、漸近緯度φj+1 を垂線緯度φk として用い、式３を計算する。式３の計算にて得られた値を緯度φk +1をφn として、上記の通り当該緯度φn における平均曲率半径から縮尺係数ｍt+1 を求める。
上記にて算出した縮尺係数ｍt+1 と、それまで各値を算出するのに用いた縮尺係数ｍt とについて、｜ｍt+1 −ｍt ｜が収束しているか判定する。その結果、収束していれば、計算を終了し、式３、式４（式１９）、式５及び式２を経て得た上記の緯度φn 、経度λ、子午線収差角γ、及び縮尺係数ｍt+1 を、目的の数値として出力する。
上記において｜ｍt+1 −ｍt ｜が収束していないと判定した場合、後処理工程を行う。後処理工程において、上記主工程で得た縮尺係数ｍt+1 を用い、他の値は、（主工程で用いたものと）変更することなく、主工程と同様の処理を行う。当該後処理工程において｜ｍt+1 −ｍt ｜が収束していなかった場合、当該後処理工程で得た縮尺係数ｍt+1 を用いて、更なる後処理工程を行う。
要するに、式２において、上記縮尺係数ｍt を原点縮尺係数ｍ0 の代わりに用いて、式３にて計算した緯度φn における平均曲率半径Ｒφn （後述のＲk2）を原点における平均曲率半径Ｒ0 の代わりに用いることにより、当該点の縮尺係数を算出する。
これを反復させると、当該点の縮尺係数は収束する性質を持ち、完全に収束した場合の縮尺係数は、当該点（入力座標）における縮尺係数である。
【００２７】
２．本願発明に係る装置の概要とその構成
図１は、この装置の説明図である。図２は、その処理のフローを示す説明図である。
２−１）位置計算装置の概要
この装置は、オペレータから平面直角座標上の座標値ｘ，ｙの入力を受付け、自動的に（出力まで人手を要することなく専ら装置のみにて）計算を行い、対応する地球楕円体（準拠楕円体を指す。ここではＷＧＳ８４楕円体を採用する。但し、ベッセル、ＧＲＳ８０等の他の楕円体を採用しても実施可能である。）上の緯度及び経度を出力するものである。
【００２８】
２−２）ハードウエア構成
この位置計算装置は、位置計算プログラムをコンピュータへ導入することにより、構築される。このコンピュータは、入力装置と、演算装置と、制御装置と、記憶装置（内部記憶装置及び外部記憶装置）とを備えた一般的なＰＣ（パーソナルコンピュータ）やＷＳ（ワークステーション）を採用して実施することができる。即ち、キーボード、マウスに代表される入力装置、ディスプレイやプリンタに代表される出力装置、ＣＰＵ（中央上方処理装置）に代表されるされる演算装置や制御装置（演算制御装置）、ＲＡＭ（ランダムアクセスメモリ）に代表される内部記憶装置、ハードディスクやフレキシブルディスクその他の磁気ディスク或いは光学ディスクに代表される外部記憶装置を備えた一般的なコンピュータを採用して実施することができる。
【００２９】
２−３）ソフトウエア
上記の位置計算プログラムは、コンピュータの資源を全て占有する専用のソフトウエアであっても実施可能であるが、ＭＳ−ＤＯＳ（登録商標）、Ｗｉｎｄｏｗｓ（登録商標）やＵＮＩＸ（登録商標）に代表される周知のＯＳ（オペレーティングシステム）上で動作するものとしても実施可能である。
また、位置計算プログラムは、インタプリンタ型のものであっても、アセンブラやコンパイラを使用して生成されたものであっても何れでもよい。
この実施の形態において、位置計算プログラムは、ＯＳが導入されたコンピュータ上に、導入されて位置計算装置を実現する。以下、ＯＳと、位置計算プログラムとを位置計算ソフトウエアと呼ぶ（必要に応じて単にソフトと呼ぶ）。
【００３０】
２−４）計算装置の構成
この実施の形態に示す位置計算装置は、等角投影法により地球楕円体（準拠楕円体）を投影する平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙを入力することにより、少なくとも当該座標値に対応する地球楕円体上の緯度φn 、経度λ、子午線収差γを出力するものである。
図１に示す通り、この装置は、データ置部１０を有する座標入力部１と、出力部２と、定数設定部３１と、初期値設定部３２と、条件設定部３３と、第１仮緯度収容部４１と、第２仮緯度収容部４２と、第１仮縮尺係数収容部５１と、第２仮縮尺係数収容部５２と、第１〜第７の７つの演算部６１，６２，６３，６４，６５，６６，６７を有する演算部６と、演算制御部７と、主判定部８１及び副判定部８２を有する判定部８と、回帰制御部９とを備える。
この実施の形態において、座標入力部１は、入力されたデータを一時的に保持するデータ置部１０を備える。また、定数設定部３１は、算出値置部３０を備える。
座標入力部１は、上記ソフトと、コンピュータが有する、上記の入力装置と、演算制御装置とにて構築される。出力部２は、上記ソフトと、コンピュータが有する上記出力装置と、演算制御装置とにて構築される。
定数設定部３１、初期値設定部３２及び条件設定部３３の夫々は、上記ソフトと、コンピュータが有する、上記の演算制御装置、内部記憶装置及び外部記憶装置とにて構築される。
第１仮緯度収容部４１、第２仮緯度収容部４２、第１仮縮尺係数収容部５１及び第２仮縮尺係数収容部５２は、夫々、上記ソフトと、コンピュータが有する、上記の演算制御装置と、内部記憶装置とにて構築される。
第１〜第７の７つの演算部６１，６２，６３，６４，６５，６６，６７と、演算制御部７と、主判定部８１と、副判定部８２と、回帰制御部９とは、夫々、上記ソフトと、コンピュータが有する、上記の演算制御装置と、内部記憶装置とにて構築される。
【００３１】
２−５）計算装置の構成の詳細
上記の座標入力部１は、オペレータからの上記座標値ｘ，ｙの入力を受け付け。データ置部１０に収容する。
定数設定部３１は、上記地球楕円体の、長半径ａ、短半径ｂ、第１離心率ｅ、第２離心率ｅ’の夫々、又はこれらに準じた値と、上記の座標ｘ，ｙが属する座標系の原点の、緯度φ0 又はその代替値、経度λ0 の夫々と、当該原点における縮尺係数ｍ0 又はその代替値と、後述する定数Ａ〜Ｆとが、既知の数値として設定されている。
ここで設定とは、この装置を用いた位置計算の度にオペレータがこれらの値（定数）を上記入力装置を通じて入力するもの、オペレータが予め外部記憶装置にこれらの値（定数）を収容しておくもの、位置計算プログラムの一部としてプログラムが備えプログラムの導入によって自動的に外部記憶装置に収容されるものの何れも含む。
ここでは、必要な定数は、プログラムに含まれるものとして外部記憶装置に導入され、導入後オペレータによってコンピュータが備える上記入力装置を通じて、書き換え可能とする（書き換えられるのは、主として、上記の、原点の緯度φ0 又はその代替値と、経度λ0 と、原点の縮尺係数ｍ0 又はその代替値である）。
定数設定部３１は、算出値置部３０を備える。算出値置部３０は、演算部６が算出した値を収容する。
【００３２】
演算部６について説明する。尚、ここでは（この２−５）の説明では）、演算部６を構成する第１〜第７の演算部６１〜６７夫々について、他の演算部との関係を無視し（演算部６１〜６７間のデータの流れについては無視し）、個々の演算部の内容を中心に説明する。従って、この２−５）において、演算部夫々の説明に用いる記号は、特に断りがない限り、他の演算部で用いるものと独立したものとする（他の演算部が、同じ記号で表された定数やパラメータを使うものであっても、特に断りがない限り、当該記号が示す数値は別のものである）。
夫々の演算にて、処理される加減乗除及び、正弦、余弦、正接は、周知の方法によって、コンピュータを用いたソフトウエアによる情報処理（数値計算）にて行われる（近似計算も含む）。
【００３３】
第１演算部６１は、少なくとも次の４つの演算１）〜４）を行うことが可能である。
演算１）
ｃ＝ａ2 ／ｂ＝ａ／（１−ｅ2 ）1/2 ＝ａ（１＋ｅ'2）1/2 …式６
を、地球楕円体の長半径ａ、地球楕円体の短半径ｂ、地球楕円体の第１離心率ｅ、地球楕円体の第２離心率ｅ’の各値から計算し、地球楕円体の極の曲率半径ｃを求めること。尚、上記第１及び第２の離心率は、地球楕円体の長半径ａと短半径ｂとから算定することができる。また、この地球楕円体の極の曲率半径ｃは、定数として算出によらずに、定数設定部３１が保持するものとしてもよい。
演算２）
Ｖ＝（１＋ｅ'2cos2φ）1/2 …式７
上記式７を、緯度φと第２離心率ｅ’の各値から計算し、Ｖを求めること。
演算３）
Ｒ＝ｃ／Ｖ2 …式８
上記式８を、地球楕円体の極の曲率半径ｃと、上記Ｖの値とから計算し、緯度φにおける平均曲率半径を求めること。
演算４）
ｍn+1 ＝（１＋ｙ2 ／２ｍn 2 Ｒ2 ＋ｙ4 ／２４ｍn 4 Ｒ4 ）×0.9999
…式９
上記式９を、平面直角座標上の座標値ｘ，ｙのｙの値と、仮の縮尺係数ｍn とから計算し、縮尺係数ｍn+1 を求めること。
この演算４）は、特許請求の範囲に示す式２を処理する機能である（式９の縮尺係数ｍn を原点縮尺係数ｍ0 とし、平均曲率半径Ｒを原点における平均曲率半径Ｒ0 としたものが、請求項２に示す式２である）。
【００３４】
第２演算部６２は、次の演算５）を行うことが可能である。
演算５）
Ｗ＝ａ（１−ｅ2 ）・｛Ａ・φ−（Ｂ／２）・（sin ２φ）
＋（Ｃ／４）・（sin ４φ）−（Ｄ／６）・（sin ６φ）
＋（Ｅ／８）・（sin ８φ）−（Ｆ／１０）・（sin 10φ）｝
＋ｘ／ｍ …式１０
上記式１０を、第２演算部６２は、地球楕円体の長半径ａの値と、定数Ａ〜Ｆの各値と、緯度φの値と、平面直角座標上の座標値ｘ，ｙのｘの値と、０以外の任意の実数ｍとから計算し、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗを求めること。尚、各定数Ａ〜Ｆは、第１離心率ｅによって定まるものである（後述のステップ１０７参照）。
【００３５】
第３演算部６３は、次の演算６）を行うことが可能である。
演算６）
Ｓ (φ )＝ａ（１−ｅ2 ）・｛Ａ・φ−（Ｂ／２）・（sin ２φ）
＋（Ｃ／４）・（sin ４φ）−（Ｄ／６）・（sin ６φ）
＋（Ｅ／８）・（sin ８φ）−（Ｆ／１０）・（sin10 φ）｝
…式１１
上記の式１１を、地球楕円体の長半径ａと、第１離心率ｅと、上記の各定数Ａ〜Ｆと、緯度φとから計算し、赤道から緯度φまでの子午線弧長Ｓ (φ )を求めること。
【００３６】
第４演算部６４は、次の演算７）を行うことが可能である。
演算７）
φn+1 ＝φn
＋（Ｗ−Ｓ (φn ) ）（１−ｅ2 sin2φn ）3/2 ／ａ（１−ｅ2 ）
…式１２
上記式１２を、仮の垂線緯度φn と、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗと、仮の垂線緯度φn における赤道からの子午線弧長Ｓ (φn ) と、地球楕円体の長半径ａと、第１離心率ｅの、各値から計算し垂線緯度φn+1 を求めること。
この式１２は、特許請求の範囲に示す式１と対応する。
【００３７】
第５演算部６５は、次の５つの演算８）〜１２）を行うことが可能である。
演算８）
ｔ＝tan φ …式１３
上記式１３を、与えられた緯度φ値から、計算し、ｔを求めること。
演算９）
η＝ｅ’cos φ …式１４
上記式１４を、地球楕円体の第２離心率ｅ’と、緯度φの、夫々の値から、計算し、ηを求めること。
演算１０）
Ｍ＝ｃ／Ｖ 3 …式１５
上記式１５を、地球楕円体の極の曲率半径ｃの値と、緯度φにて定まるＶの値とから計算し、当該緯度φについての子午線曲率半径Ｍを求めること。
演算１１）
Ｎ＝ｃ／Ｖ …式１６
上記式１６を、地球楕円体の極の曲率半径ｃの値と、緯度φにて定まるＶの値とから計算し、当該緯度φについての卯酉線曲率半径Ｎを求めること。
演算１２）
φn+1 ＝φn −ｔn （ｙ／ｍt ）2 ／２Ｍn Ｎn ＋ｔn （ｙ／ｍt ）4 （５＋３ｔn 2 ＋ηn 2 −９ｔn 2 ηn 2 −４ηn 4 ）／24Ｍn Ｎn 3 −ｔn （ｙ／ｍt ）6 （61＋90ｔn 2 ＋45ｔn 4 ）／720 Ｍn Ｎn 5 …式１７
上記式１７を、平面直角座標上の座標値ｘ，ｙのｙの値、最近似垂線緯度φn 、縮尺係数ｍt 、最近似垂線緯度φn についてのｔn 、最近似垂線緯度φn についてのηn 、最近似垂線緯度φn についての子午線曲率半径Ｍn 、及び、最近似垂線緯度φn についての卯酉線曲率半径Ｎn の、夫々の値から計算し、当該緯度φn+1 を求めること。
この式１７は、特許請求の範囲に示す式３に対応する。
【００３８】
第６演算部６６は、次の２つの演算１３）及び１４）を行うことが可能である。
演算１３）
Δλ＝（ｙ／ｍt ）／Ｎn cos φn
−（ｙ／ｍt ）3 （１＋２ｔn 2 ＋ηn 2 ）／６Ｎn 3cosφn
＋（ｙ／ｍt ）5 （５＋２８ｔn 2 ＋２４ｔn 4 ）／１２０Ｎn 5cosφn
…式１８
上記式１８を、平面直角座標上の座標値ｘ，ｙのｙの値、最近似垂線緯度φn 、縮尺係数ｍt 、最近似垂線緯度φn についてのｔn 、最近似垂線緯度φn についてのηn 、及び、最近似垂線緯度φn についての卯酉線曲率半径Ｎn の、夫々の値から計算し、経度差分Δλを求めること。
この式１８は、特許請求の範囲（請求項４）の式４に対応する。
演算１４）
λ＝λ0 ＋Δλ …式１９
上記式１９を、平面直角座標上の座標値ｘ，ｙの属する原点の経度λ0 と、経度差分の夫々の値とから、計算し、経度λを求めること。
【００３９】
第７演算部６７は、次の演算１５）を行うことが可能である。
演算１５）
γ＝（ｙ／ｍt ）／Ｎn
−（ｙ／ｍt ）3 ｔn （１＋ｔn 2 −ηn 2 ）／３Ｎn 3
＋（ｙ／ｍt ）5 ｔn （２＋５ｔn 2 ＋３ｔn 4 ）／１５Ｎn 5
…式２０
上記式２０を、平面直角座標上の座標値ｘ，ｙのｙの値、縮尺係数ｍt 、最近似垂線緯度φn についてのｔn 、最近似垂線緯度φn についてのηn 、及び、最近似垂線緯度φn についての卯酉線曲率半径Ｎn の、夫々の値から計算し、子午線収差γを求めること。
この式２０は、特許請求の範囲（請求項５）の式５に対応する。
【００４０】
演算制御部７は、演算部６の各演算部６１〜６７へ演算の指令を行う。各演算部６１〜６７から、演算及び演算に必要な数値の所在（置き場）の指令を受けて、演算制御部７の指令の順に各演算部６１〜６７は、指令された演算を行う。
また、演算制御部７は、演算部６１〜６７の演算の結果得た緯度及び縮尺係数の収束状況について、判定部８に判定を行わせる。
判定部８は、垂線緯度の収束状態の判定を行う主判定部８１と、縮尺係数の判定状態の判定を行う副判定部８２とを備える。
回帰制御部９は、判定部８の判定を受けて、演算制御部７に、判定に対応した指示を行う。当該指示に従って、演算制御部７は、各演算部６１〜６７に指令を出す。即ち、回帰演算部９は、判定部８による判定の結果に従って、演算制御部７を通じて、各演算部６１〜６７を制御する。
【００４１】
初期値設定部３２は、（演算部６が行う）垂線緯度計算の初期値を決定するための数値（以下、必要に応じて初期値決定用助数と呼ぶ。）が設定され、当該数値を保持する。
条件設定部３３は、判定部８の判定（比較）に必要な基準値が設定され、設定された基準値を保持する。
【００４２】
第１仮緯度収容部４１は、演算部６が算出した緯度の値を収容する。第１仮緯度収容部４１は、新たに演算部６が算出した緯度の値を受け取ることにより、それまで収容していた緯度の値を、第２仮緯度収容部４２へ渡し、それまで収容した緯度の値を新たに受け取った緯度の値に書き換えて、新たな緯度の値を収容する。
第２仮緯度収容部４２は、第１仮緯度収容部４１から、第１仮緯度収容部４１が収容していた緯度の値を受け取り、当該緯度の値を収容する。第２仮緯度収容部４２は、更に第１仮緯度収容部４１から緯度の値を受け取ることにより、それまで収容していた緯度の値を受け取った後の緯度の値に書き換えて、当該後の緯度の値を収容する。
【００４３】
第１仮縮尺係数収容部５１は、演算部６が算出した縮尺係数の値を収容する。第１仮縮尺係数収容部５１は、新たに演算部６が算出した縮尺係数の値を受け取ることにより、それまで収容していた縮尺係数の値を、第２仮縮尺係数収容部５２へ渡し、それまで収容した縮尺係数の値を新たに受け取った縮尺係数の値に書き換えて、新たな縮尺係数の値を収容する。
第２仮縮尺係数収容部５２は、第１仮縮尺係数収容部５１から、第１仮縮尺係数収容部５１が収容していた縮尺係数の値を受け取り、当該縮尺係数の値を収容する。第２仮緯度収容部４２は、更に第１仮縮尺係数収容部５１から縮尺係数の値を受け取ることにより、それまで収容していた縮尺係数の値を受け取った後の縮尺係数の値に書き換えて、当該後の縮尺係数の値を収容する。
出力部２は、演算制御部７と回帰演算部９の指令を受けて、演算部６の演算結果を出力する。
【００４４】
３．計算の手順
以下、この装置を利用した緯度・経度・子午線収差の算出の流れについて、図２を用いて、説明する。
この計算方法は、主工程１００を遂行し、更に縮尺係数の収束状況の判定結果により後処理工程を遂行するものである。
また、上記の主工程１００の遂行に先立ち、事前に、前処理工程を行っておく。
【００４５】
３−０）前処理工程
前処理工程は、オペレータにより、位置計算前に行われる種々の設定工程である。
詳しくは、前処理工程において、オペレータにより、初期値設定部３２にて緯度計算の初期値を決定するための数値（初期値決定用助数）が設定される。このような助数として、１”（１秒）をrad （ラジアン）に換算した値である 1／206265が、条件設定部３３において設定される。但しこのような数値は変更可能である。
又、前処理工程において、オペレータにより、条件設定部３３に、垂線緯度の収束条件ＩＰＳ（しきい値）が設定される。
この収束条件ＩＰＳは、計算の精度を確保する上で重要である。即ち、演算部における各式にて的確な計算がなされても、収束条件ＩＰＳが広いと、十分な精度を確保するのが困難になるからである。
ここでは、収束条件ＩＰＳは、２”× 1０-5／206265rad と設定する。
更に、前処理工程において、オペレータにより、初期値設定部３２にて、縮尺係数の収束条件ＩＰＳ２が設定される。ここでは、ＩＰＳ２は、座標値の有効数字（最大９桁）を考慮して、４× 1０-12 と設定してある。
但し、各収束条件ＩＰＳ、ＩＰＳ２の各数値は、上記の値に限定するものではなく、変更可能である。
主工程１００以降の工程は、オペレータから座標値ｘ，ｙを受け付けたこの位置計算装置が、位置計算を実際に行う工程である。
主工程は、次のステップ１０１〜１２５を順次遂行する。
【００４６】
３−１）ステップ１０１（座標値受付工程）
オペレータは、座標入力部１により、直角座標の座標値ｘ，ｙ、及び当該座標値が属する座標系Ｚを入力する。座標系Ｚの入力については、日本に設定されている１９の座標系のうち何れに属するかを特定するものである。
上記の入力作業は、オペレータが行う（これ以降の他の作業は、位置計算装置が行う。尚、この実施の形態において、後述する初期値設定部３２の設定、及び定数設定部３１内の各代替値の設定は、位置計算前に、オペレータにより行われる）。
座標入力部１は、入力された座標値ｘ，ｙを、データ置部１０に収容し、演算制御部７に座標値の入力があったことを通知する。
演算制御部７は、上記通知を受けて、第１演算部６１に、定数設定部３１を参照させる。即ち、定数設定部３１は、各座標系の原点の緯度及び経度の数値を保持しており、第１演算部６１は、入力された（座標値ｘ，ｙが属する）座標系Ｚの、原点の緯度φ0 及び経度λ0 を特定する。これにて、以下、各演算部が定数設定部３１にて参照する原点の緯度φ0 及び経度λ0 は、確定する。
【００４７】
３−２）ステップ１０２（極の曲率半径計算工程）
このステップ１０２において、演算制御部７は、第１演算部６１に、定数設定部３１を参照させて、定数設定部３１が保持している、地球楕円体の長半径ａと短半径ｂから、或いは、地球楕円体の長半径ａと第１離心率ｅから、或いは、地球楕円体の長半径ａと第２離心率ｅ’から、下記の計算１を行わせ、地球楕円体の極の曲率半径ｃを算出させる。第１演算部６１は、算出した極の曲率半径ｃを、算出値置部３０に収容する。
ｃ＝ａ2 ／ｂ＝ａ／（１−ｅ2 ）1/2 ＝ａ（１＋ｅ'2）1/2 …計算１
計算１は、第１演算部６１の演算１）の式６の計算機能にて実行される。
【００４８】
３−３）ステップ１０３（Ｖ0 計算工程）
ステップ１０３において、演算制御部７は、第１演算部６１に、定数設定部３１を参照させて、第２離心率ｅ’と当該座標が属する原点の緯度φ0 （以下単に原点緯度φ0 と呼ぶ。）とから、計算２を行わせ、Ｖ0 を算出させる。第１演算部６１は、算出したＶφ0 を、算出値置部３０に収容する。
Ｖ0 ＝（１＋ｅ'2cos2φ0 ）1/2 …計算２
計算２は、第１演算部６１の演算２）の式７の計算機能にて実行される。
【００４９】
３−４）ステップ１０４（原点の平均曲率半径計算工程）
ステップ１０４において、演算制御部７は、第１演算部６１に、定数設定部３１（の算出値置部３０）を参照させて、極の曲率半径ｃとＶ0 の値とから、計算３を行わせ、当該座標系原点の平均曲率半径Ｒ0 を算出させる。第１演算部６１は、算出した平均曲率半径Ｒ0 を算出値置部３０に収容する。
Ｒ0 ＝ｃ／Ｖ02…計算３
計算３は、第１演算部６１の演算３）の式８の計算機能にて実行される。
【００５０】
３−５）ステップ１０５（縮尺係数計算工程）
ステップ１０５において、演算制御部７は、定数設定部３１の原点における縮尺係数ｍ0 （以下原点縮尺係数ｍ0 と呼ぶ。）を、第１仮縮尺係数収容部５１に収容させる（定数設定部３１の上記原点縮尺係数ｍ0 を第１仮縮尺係数収容部５１へコピーする）。この原点縮尺係数ｍ0 は、0.9999という数値を採る。
そして、演算制御部７は、第１演算部６１に、座標入力部１のデータ置部１０及び定数設定部３１（又は第１仮縮尺係数収容部５１）を参照させて、入力された座標値ｙの値と、原点縮尺係数ｍ0 と、上記平均曲率半径Ｒ0 とから、計算４を行わせ、縮尺係数ｍt を算出させる。
ｍt ＝（１＋ｙ2 ／２ｍ02Ｒ02＋ｙ4 ／２４ｍ04Ｒ04）×0.9999…計算４
（右辺の0.9999は、数値として原点縮尺係数ｍ0 と同じであるが、これは基準子午線上の縮尺係数を0.9999とするための一定の縮尺係数である。）
計算４は、第１演算部６１の演算４）の式９の計算機能にて実行される。
第１演算部６１は、算出した縮尺係数ｍt を第１仮縮尺係数収容部５１に収容する。第１仮縮尺係数収容部５１は、収容していた原点縮尺係数ｍ0 を第２仮縮尺係数収容部５２へ渡し、原点縮尺係数ｍ0 に代え上記の通り上記縮尺係数ｍt を収容する。
【００５１】
３−６）ステップ１０６（垂線緯度の初期値設定工程）
ステップ１０６において、演算制御部７は、第１演算部６１に、定数設定部３１と、初期値設定部３２とを参照させて、原点緯度φ0 から助数を減じた値、即ち、原点緯度φ0 − 1／206265を計算させる。第１演算部６１は、算出した値を緯度φj として、第１仮緯度収容部４１に収容する。
【００５２】
３−７）ステップ１０７（赤道から当該座標値ｘまでの基準子午線上の子午線弧長計算工程）
このステップ１０７において、演算制御部７は、第２演算部６２に、定数設定部３１と、座標入力部１のデータ置部１０と、第１仮縮尺係数収容部５１とを参照させて、長半径ａと、第１離心率ｅと、原点緯度φ0 と、座標値ｘと、縮尺係数ｍt とから、計算５をさせ、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとして算出させる。
Ｗ＝ａ（１−ｅ2 ）・｛Ａ・φ0 −（Ｂ／２）・（sin ２φ0 ）
＋（Ｃ／４）・（sin ４φ0 ）−（Ｄ／６）・（sin ６φ0 ）
＋（Ｅ／８）・（sin ８φ0 ）−（Ｆ／10）・（sin10 φ0 ）｝
＋ｘ／ｍt …計算５
計算５の、Ａ〜Ｆは、次の通りである。
Ａ＝１＋（３／４）ｅ2 ＋（45／64）ｅ4 ＋（175 ／256 ）ｅ6
＋（11025 ／16384 ）ｅ8 ＋（43659 ／65536 ）ｅ10
Ｂ＝ （３／４）ｅ2 ＋（15／16）ｅ4 ＋（525 ／512 ）ｅ6
＋（2205／2048）ｅ8 ＋（72765 ／65536 ）ｅ10
Ｃ＝ （15／64）ｅ4 ＋（105 ／256 ）ｅ6
＋（2205／4096）ｅ8 ＋（10395 ／16384 ）ｅ10
Ｄ＝ （35／512 ）ｅ6
＋（315 ／2048）ｅ8 ＋（31185 ／131072）ｅ10
Ｅ＝ （315 ／16384 ）ｅ8 ＋（3465／65536 ）ｅ10
Ｆ＝ （693 ／131072）ｅ10
計算５は、第２演算部６２の演算５）の式１０の計算機能にて実行される。
前述の通り、この定数Ａ〜Ｆは、定数設定部３１に予め収容されており、第２演算部６２や他の演算部が参照して利用する。
ここで、計算５の右辺のｘ／ｍt 以外の部分が、地球楕円体の赤道から座標値ｘ，ｙが属する座標系原点までの子午線弧長を示している。また、この縮尺係数ｍt は、特許請求の範囲でいう、仮の縮尺係数である。
第２演算部６２は、算出した、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗを定数設定部３１の算出値置部３０に収容する。
【００５３】
３−８）ステップ１０８（緯度に対する赤道からの子午線弧長算出工程）
このステップ１０８において、演算制御部７は、第３演算部６３に、定数設定部３１と、第１仮緯度収容部４１とを参照させ、長半径ａ、第１離心率ｅ、定数Ａ〜Ｆと、緯度φj とから、計算６を行わせ、緯度φj に対する赤道からの子午線弧長Ｓ (φj ) を算出させる。
Ｓ (φj ) ＝ａ（１−ｅ2 ）・｛Ａ・φj −（Ｂ／２）・（sin ２φj ）
＋（Ｃ／４）・（sin ４φj ）−（Ｄ／６）・（sin ６φj ）
＋（Ｅ／８）・（sin ８φj ）−（Ｆ／１０）・（sin10 φj ）｝
…計算６
即ち、このステップ１０８は、子午線曲率半径Ｍを以って赤道から垂線緯度φj まで積分するステップであり、上記の通り計算６は、第３演算部６３にて当該積分の近似計算である（ここでいう積分とは、当該近似計算を指している）。
計算６は、第３演算部６３の演算６）の式１１の計算機能にて実行される。
第３演算部６３は、上記にて算出した子午線弧長Ｓ (φj ) を、上記算出値置部３０に収容する。
【００５４】
３−９）ステップ１０９（初期の垂線緯度の収束計算工程）
ステップ１０９において、演算制御部７は、第４演算部６４に、定数設定部３１と、第１仮緯度収容部４１とを参照させ、緯度φj 、長半径ａ、第１離心率ｅ、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗ、子午線弧長Ｓ (φj ) とから、計算７を行わせ、垂線緯度φj+1 を算出させる。
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ） …計算７
計算７は、第４演算部６４の演算７）の式１２の計算機能にて実行される。
第４演算部６４は、算出した緯度φj+1 を、第１仮緯度収容部４１へ収容する。第４演算部６４から上記の緯度φj+1 を受けた第１仮緯度収容部４１は、収容していた緯度φj を第２仮緯度収容部４２へ渡し、緯度φj に代え上記緯度φj+1 の値を収容する。
【００５５】
３−１０）ステップ１１０（垂線緯度の収束判定工程）
ステップ１１０において、判定部８の主判定部８１は、第１仮緯度収容部４１と第２仮緯度収容部４２と条件設定部３３とを参照し、緯度φj+1 と緯度φj の差の絶対値即ち、｜φj+1 −φj ｜と、緯度の収束条件ＩＰＳ（しきい値）とを比較する。そして、｜φj+1 −φj ｜＞ＩＰＳと判定した場合、判定部８（主判定部８１）が、その結果を回帰制御部９に通知する。
この通知を受けた回帰制御部９は、第１仮緯度収容部４１の緯度φj+1 の値を緯度φj の値として（第１仮緯度収容部４１の緯度φj+1 を緯度φj に書き換えて）、装置の処理を再びステップ１０８に移す。そして、回帰制御部９は、ステップ１０８からステップ１１０を、判定部８（主判定部８１）が｜φj+1 −φj ｜≦ＩＰＳと判定するまで繰り返させる（垂線緯度の反復計算をさせる）。この反復計算において、垂線緯度φj+1 を垂線緯度φj に置き換える（書き換える）のみであり、演算部６（の各演算部）が参照する他の数値を変更しない（当該他の数値については、当初の値を用いて最初から計算される）。
以上により、ステップ１１０において、判定部８（の主判定部８１）は、｜φj+1 −φj ｜≦ＩＰＳと判定した場合、その結果を回帰制御部９に通知する。
回帰制御部９は、当該収束の通知を受けて、次のステップ１１１へ処理を移行させる。即ち、回帰制御部９は、演算制御部７に次のステップ１１１の処理を行わせる。
【００５６】
３−１１）ステップ１１１（Ｖk 計算工程）
ステップ１１１において、演算制御部７は、第１演算部６１に、定数設定部３１と、第１仮緯度収容部４１とを参照させて、第２離心率ｅ’と垂線緯度φj+1 とから、計算８を行わせ、Ｖk を算出させる。尚、説明の便宜上、この垂線緯度φj+1 を緯度φk として説明する。
Ｖk ＝（１＋ｅ'2cos2φk ）1/2 …計算８
計算８は、第１演算部６１の演算２）の式７の計算機能にて実行される。
第１演算部６１は、算出したＶk を、算出値置部３０に収容する。
【００５７】
３−１２）ステップ１１２（第２の値計算工程）
ステップ１１２において、演算制御部７は、第５演算部６５に、第１仮緯度収容部４１を参照させ、緯度φk の値を用いて（即ち、緯度φj+1 の値を緯度φk の値として用いて）、計算９を行わせ、ｔk （説明の便宜上第２の値ｔk と呼ぶ。）を算出させる。
ｔk ＝tan φk …計算９
計算９は、第５演算部６５の演算８）の式１３の計算機能にて実行される。
第５演算部６５は、算出した第２の値ｔk を、算出値置部３０へ収容する。
【００５８】
３−１３）ステップ１１３（第１の値計算工程）
ステップ１１３において、演算制御部７は、第５演算部６５に、定数設定部３１と第１仮緯度収容部４１を参照させ、第２離心率ｅ’と共に、緯度φk の値を用いて、計算１０を行わせ、ηk （説明の便宜上第１の値ηk と呼ぶ。）を算出させる。
ηk ＝ｅ’cos φk …計算１０
計算１０は第５演算部６５の演算９）の式１４の計算機能にて実行される。
第５演算部６５は、算出した第１の値ηk を、算出値置部３０へ収容する。
【００５９】
３−１４）ステップ１１４（子午線曲率半径の計算工程）
ステップ１１４において、演算制御部７は、第５演算部６５に、定数設定部３１を参照させ、極の曲率半径ｃとＶk とから、計算１１にて、緯度φk における子午線曲率半径Ｍk を算出させる。
Ｍk ＝ｃ／Ｖk 3 …計算１１
計算１１は、第５演算部６５の演算１０）の式１５の計算機能にて実行される。
第５演算部６５は、算出した子午線曲率半径Ｍk の値を、算出値置部３０へ収容する。
【００６０】
３−１５）ステップ１１５（卯酉線曲率半径の計算工程）
ステップ１１５において、演算制御部７は、第５演算部６５に、定数設定部３１を参照させ、極の曲率半径ｃとＶk とから、計算１２にて、緯度φk における卯酉線曲率半径Ｎk を算出させる。
Ｎk ＝ｃ／Ｖk …計算１２
計算１２は、第５演算部６５の演算１１）の式１６の計算機能にて実行される。
第５演算部６５は、算出した卯酉線曲率半径Ｎk の値を、算出値置部３０へ収容する。
【００６１】
３−１６）ステップ１１６（入力座標に対応する緯度の計算工程）
ステップ１１６において、回帰制御部９は、演算部制御部７を通じて、第５演算部６５に、緯度φk を与え、計算１３を行わせる。具体的には、演算制御部７は、回帰制御部９からの指令を受けて、第５演算部６５に、座標入力部１のデータ置部１０、第１仮緯度収容部４１、第１仮縮尺収容部５１及び定数設定部３１を参照させて、上記にて、第４演算部６４に算出させた最近似垂線緯度φk と、当該点の縮尺係数ｍt と、最近似垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、計算１３の計算にて、緯度φk+1 を計算させる。
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …計算１３
計算１３は、第５演算部６５の演算１２）の式１７の計算機能にて実行される。
第５演算部６５は、算出した緯度φk+1 を、当該点（入力座標点）の緯度φn として、第１仮緯度収容部４１へ、それまで第１仮緯度収容部４１に収容されていた値と置き換えて収容する。このとき、第１仮緯度収容部４１は、それまで保持していた緯度φk の値を第２仮緯度収容部４２へ渡す。
【００６２】
３−１７）ステップ１１７（経度差分計算工程）
このステップ１１７は、上記のステップ１１６と平行として行われる。即ち、ステップ１１７は、ステップ１１６で得られた結果（緯度φk+1 ）を用いるのではなく、ステップ１１６の処理に用いたデータを用いて、実行される。
ステップ１１７において、演算制御部７は、第６演算部６６に、座標入力部１のデータ置部１０、第１仮縮尺係数収容部５１及び算出値置部３０を参照させて、最近似垂線緯度φk （ステップ１１１に至るまでの一連の処理で得られた緯度、即ちステップ１１１〜１１６の処理のために用いた緯度φk ）と、当該点の縮尺係数ｍt と、最近似垂線緯度φk における子午線曲率半径Ｍk （ステップ１１４の結果）及び卯酉線曲率半径Ｎk （ステップ１１５の結果）と、ｅ’cos φk にて定まる第１の値ηk （ステップ１１３の結果）と、tan φk にて定まる第２の値ｔk （ステップ１１２の結果）と、上記座標値ｙとから、計算１４を行わせ、経度差分Δλを算出させる。
Δλ＝（ｙ／ｍt ）／Ｎk cos φk
−（ｙ／ｍt ）3 （１＋２ｔk 2 ＋ηk 2 ）／６Ｎk 3cosφk
＋（ｙ／ｍt ）5 （５＋２８ｔk 2 ＋２４ｔk 4 ）／１２０Ｎk 5cosφk
…計算１４
第６演算部６６は、算出した経度差分Δλを、定数設定部３１の算出値置部３０に、収容する。
【００６３】
３−１８）ステップ１１８（経度計算工程）
このステップ１１８と当該ステップ１１８に続くステップ１１９は、ステップ１１７と同様、上記のステップ１１６と平行として行われる。
ステップ１１８において、演算制御部７は、第６演算部６６に、定数設定部３１を参照させて、原点経度λ0 と経度差分Δλの和の計算を即ち、計算１５をさせて、経度λを算出させる。
第６演算部６６は、算出した経度λの値を算出値置部３０に収容する。
λ＝λ0 ＋Δλ …計算１５
計算１５は、第６演算部６６の演算１４）の式１９の計算機能にて実行される。
【００６４】
３−１９）ステップ１１９（子午線収差の計算工程）
このステップ１１９は、上記の通り、ステップ１１８に続く工程であり、ステップ１１６と平行して行われる。即ち、ステップ１１６で得られた結果（緯度φk+1 ）を用いるのではなく、ステップ１１６の処理に用いたデータを用いて、実行される。
ステップ１１９において、演算制御部７は、第７演算部６７に、座標入力部１のデータ置部１０、第１仮縮尺係数収容部５１及び算出値置部３０を参照させて、最近似垂線緯度φk （ステップ１１１に至るまでの一連の処理で得られた緯度、即ちステップ１１１〜１１６の処理のために用いた緯度φk ）と、当該点の縮尺係数ｍt と、最近似垂線緯度φk における卯酉線曲率半径Ｎk （ステップ１１５の結果）と、ｅ’cos φk にて定まる第１の値ηk （ステップ１１３の結果）と、tan φk にて定まる第２の値ｔk （ステップ１１２の結果）と、上記座標値ｙとから、計算１６を行わせ、子午線収差γを算出させる。
γ＝（ｙ／ｍt ）／Ｎk
−（ｙ／ｍt ）3 ｔk （１＋ｔk 2 −ηk 2 ）／３Ｎk 3
＋（ｙ／ｍt ）5 ｔk （２＋５ｔk 2 ＋３ｔk 4 ）／１５Ｎk 5
…計算１６
計算１６は、第７演算部６７の演算１５）の式２０の計算機能にて実行される。
第７演算部６７は、算出した子午線収差γを算出値置部３０に収容する。
【００６５】
３−２０）ステップ１２０（Ｖk2計算工程）
ステップ１２０は、ステップ１１６に続く工程である。
ステップ１２０において、演算制御部７は、第１演算部６１に、定数設定部３１と、第１仮緯度収容部４１とを参照させて、第２離心率ｅ’とステップ１１６で得た当該点の緯度φk+1 とから、計算１７を行わせ、Ｖk2を算出させる。尚、説明の便宜上、この緯度φk+1 を緯度φk2として説明する。
Ｖk2＝（１＋ｅ'2cos2φk2）1/2 …計算１７
計算１７は、第１演算部６１の演算２）の式７の計算機能にて実行される。
第１演算部６１は、算出したＶk2を、算出値置部３０に収容する。
【００６６】
３−２１）ステップ１２１（子午線曲率半径の計算工程）
ステップ１２１において、演算制御部７は、第５演算部６５に、定数設定部３１を参照させ、極の曲率半径ｃとＶk2とから、計算１８にて、緯度φk2における子午線曲率半径Ｍk2を算出させる。
Ｍk2＝ｃ／Ｖk23 …計算１８
計算１８は、第５演算部６５の演算１０）の式１５の計算機能にて実行される。
第５演算部６５は、算出した子午線曲率半径Ｍk2の値を、算出値置部３０へ収容する。
【００６７】
３−２２）ステップ１２２（卯酉線曲率半径の計算工程）
ステップ１２２において、演算制御部７は、第５演算部６５に、定数設定部３１を参照させ、極の曲率半径ｃとＶk2とから、計算１９にて、緯度φk2における卯酉線曲率半径Ｎk2を算出させる。
Ｎk2＝ｃ／Ｖk2 …計算１９
計算１９は、第５演算部６５の演算１１）の式１６の計算機能にて実行される。
第５演算部６５は、算出した卯酉線曲率半径Ｎk2の値を、算出値置部３０へ収容する。
【００６８】
３−２３）ステップ１２３（平均曲率半径の計算工程）
ステップ１２３において、演算制御部７は、第１演算部６１に、定数設定部３１を参照させて、極の曲率半径ｃとＶk2の値とから、計算２０を行わせ、当該座標（入力座標）の平均曲率半径Ｒk2を算出させる。第１演算部６１は、算出した平均曲率半径Ｒk2を算出値置部３０に収容する。
Ｒk2＝ｃ／Ｖk22 …計算２０
計算２０は、第１演算部６１の演算３）の式８の計算機能にて実行される。
【００６９】
３−２４）ステップ１２４（縮尺係数計算工程）
ステップ１２４において、演算制御部７は、第１演算部６１に、座標入力部１のデータ置部１０及び第１仮縮尺係数収容部５１を参照させて、入力された座標値ｙの値と、縮尺係数ｍt1（第１仮縮尺係数収容部５１に収容されている前記縮尺係数ｍt を縮尺係数ｍt1として説明する。以下同じ。）と、上記平均曲率半径Ｒk2とから、計算２１を行わせ、縮尺係数ｍt2を算出させる。
ｍt2＝（１＋ｙ2 ／２ｍt12 Ｒk22 ＋ｙ4 ／２４ｍt14 Ｒk24 ）×0.9999
…計算２１
計算２１は、第１演算部６１の演算４）の式９の計算機能にて実行される。
第１演算部６１は、算出した縮尺係数ｍt2を第１仮縮尺係数収容部５１に収容する。第１仮縮尺係数収容部５１は、収容していた原点縮尺係数ｍt を第２仮縮尺係数収容部５２へ渡し、上記の通り縮尺係数ｍt に代え上記縮尺係数ｍt2を収容する。
【００７０】
３−２５）ステップ１２５（縮尺係数の収束判定工程）
ステップ１２５において、判定部８の副判定部８２は、第１仮縮尺係数収容部５１と、第２仮縮尺係数収容部５２と、条件設定部３３とを参照し、縮尺係数ｍt1と縮尺係数ｍt2の差の絶対値即ち｜ｍt2−ｍt1｜と、縮尺係数の収束条件ＩＰＳ２とを比較する。そして、縮尺係数ｍt2（と縮尺係数ｍt1の差の絶対値）が収束したと判定された場合、即ち、副判定部８２が｜ｍt2−ｍt1｜≦ＩＰＳ２と判定した場合、判定部８は、その結果を回帰制御部９に通知する。
【００７１】
３−２６）ステップ１２６（出力工程）
上記ステップ１２５にて、縮尺係数の収束の通知を受けた回帰制御部９は、出力部２により、少なくとも、第１仮緯度収容部４１が収容している当該点の緯度φk+1 の値（φn ) と、算出値置部３０が収容している経度λの値及び子午線収差γの値とを、出力させる。このとき、必要に応じて、座標入力部１のデータ置部１０の入力座標値ｘ，ｙや、定数設定部３１が保持しその算出値置部３０が収容している、上記入力座標が属する原点の緯度や経度の数値、地球楕円体の長半径や短半径などの定数、或いは計算工程中に算出した中間の値を、併せて、出力部２に出力させるものとしても、上記の結果がどのような座標値に対応するか、また、どのような経緯を経て算出されたか、オペレータに確認させることができ、便利である。
【００７２】
一方、ステップ１２５において、判定部８（副判定部８２）が｜ｍt2−ｍt1｜＞ＩＰＳ２と判定した場合、ステップ１２６に移行せずに、前述の後処理工程を遂行する。
即ち、上記にて、｜ｍt2−ｍt1｜＞ＩＰＳ２と判定された場合、判定部８は、その結果を回帰制御部９に通知する。この通知を受けた回帰制御部９は、演算制御部７を通じて、上記主工程１００の各ステップ（ステップ１０２〜１２５）を繰り返させる。この後処理工程において、主工程１００と異なるのは、ステップ１０５の処理において、前述の原点縮尺係数ｍ0 に代え、第１仮縮尺係数収容部５１に収容されている縮尺係数ｍt2を第１演算部６１に参照させて、計算４の計算を行わせるということである（原点縮尺係数ｍ0 として、上記主工程１００を経て得た縮尺係数ｍt2を用いる。即ち、後処理工程で原点縮尺係数ｍ0 として用いる値は、上記主工程１００を経て得た縮尺係数ｍt2である）。他の値については、上記主工程１００と同様のものを用いる。即ち、後処理工程は、縮尺係数ｍt2を用いて、主工程１００をやり直す工程である。従って、後処理工程において、主工程１００の結果として加味されるのは、主工程１００で算出された縮尺係数のみであり、主工程１００で算出した緯度等の値は用いず、再度主工程１００で用いたのと同様の値を用いて計算する（縮尺係数を使わない計算や、縮尺係数を用いる計算であっても、縮尺係数以外の数値については、主工程１００で用いたものと同様の数値を用いて計算を行う）。
この後処理工程でのステップ１２５において、｜ｍt2−ｍt1｜＞ＩＰＳ２と判定された場合、判定部８は、その結果を回帰制御部９に通知する。そして、当該後処理工程を経て得た縮尺係数ｍt2を原点縮尺係数ｍ0 として、再度後処理工程（主工程１００のステップ１０２からのやり直し）を行う。このように、１回の後処理工程を経ても、縮尺係数ｍt2（と縮尺係数ｍt1の差の絶対値）が収束したと判定されない場合、収束するまで、即ち、副判定部８２が｜ｍt2−ｍt1｜≦ＩＰＳ２となるまで、上記の後処理工程を繰り返す。
後処理工程におけるステップ１２５にて、縮尺係数ｍt2が収束した即ち｜ｍt2−ｍt1｜≦ＩＰＳ２と副判定部８２に判定された場合、上記の通り、回帰制御部９は、演算制御部７に、ステップ１２６を遂行させる。
【００７３】
上記の通り、後処理工程（主工程１００のやり直し工程）が繰り返されると、ステップ１０７における計算５のｘ／ｍt が変化し、座標値ｘまでの基準子午線上の子午線弧長Ｗの値が変化し、その結果、計算７の計算結果がより精度の高いものとなるのである。
上記の実施の形態において、後処理工程は、判定部８（副判定部８２）による縮尺係数の収束不十分、即ち、｜ｍt2−ｍt1｜＞ＩＰＳ２の判定により、行うものとした。
この他、前述の通り、このような縮尺係数の収束の判定を行わずに、設定した回数必ず後処理工程を実行するものとしても実施可能である。例えば、後処理工程を必ず１回（或いは２回以上）行うものとして実施することが可能である。
また、上記の実施の形態において、ステップ１０７で用いる仮の数値として、ステップ１０５で算出した縮尺係数ｍt を用いたが、初回のステップ１０７（後処理工程でのステップ１０７は除く。）については、０以外の任意の実数を当該縮尺係数ｍt の代替値として用いて計算を行うものとしても実施可能である。
但し、ステップ１０７で用いる仮の縮尺係数として、ステップ１０５で算出した縮尺係数ｍt を用いることにより、収束が効率良く行われるので、例えば、ステップ１１０における判定の結果ステップ１０８への反復（小ループ）を３回程度に抑えることができ、また、ステップ１２５における後処理工程の反復（大ループ）を２回以内に抑えることができる。従って、この場合、後処理工程をステップ１２５の判定無しに２回に固定しても、事実上良好な結果が得られる。
【００７４】
ステップ１２５における縮尺係数の判定を行わずに後処理工程を必ず行うものとした場合、通常、主工程１００のステップ１０５において、仮の数値として、原点縮尺係数ｍ0 を用いて計算４をした縮尺係数ｍt を採用する場合、入力した座標値ｘ，ｙがその座標系のどの位置にあっても、主工程１００において、ステップ１１０の判定による反復計算（ステップ１０８〜１１０）を繰り返し、後処理工程を１回（当該後処理工程中においてもテップ１１０の判定による反復計算（ステップ１０８〜１１０）を繰り返す。）行えば、十分な精度の結果が得られる。
また、ステップ１０６において、原点緯度φ0 に対して、１秒（−1/206265ラジアン）ずらすものとした。このずらす値については、原点緯度φ0 に対して±５秒を超えない範囲で行うのが好ましい。この範囲であれば、上記の通り、主工程１００のステップ１１０の判定による反復計算（ステップ１０８〜１１０）を最大３回とし後処理工程を１回として十分な精度が得られることが発明者によって確認されている。従って、このステップ１０８〜１１０の反復計算を３回を限度として実施するのが好ましい。但し、反復計算は、収束条件ＩＰＳに依存しているため、このような回数に限定するものではない。
【００７５】
平面上の点を楕円体上の点とを対応させるのには、近似値を求める計算（反復計算即ち収束計算）が必要であり、精度良く対応させるには本願発明のような工夫が必要となるが、逆に、楕円体上の点を平面上の点に正確に対応させるのは、対応する点を簡単に算出することができる（反復計算は不要である）。即ち、周知の計算方法や装置で、緯度・経度から平面直角座標上の座標値を逆算して正確な値を算出することができる。
従って、このような逆算によって、本願発明に係る方法で得られた計算結果の精度を検証することができる。
【００７６】
検証法の例について簡単に説明する。
本願発明に係る上記の方法にて得られた結果（出力値）を緯度φ（緯度φn ）、経度λとして、下記の計算を行う。
計算２２にて、経度差分Δλを逆算する。
Δλ＝λ−λ0 …計算２２
上記結果の緯度φから、
計算２３にて、第１の値ηを求める。
η＝ｅ’cos φ…計算２３
また、上記の結果の緯度φから、計算２４にて、第２の値ｔを求める。
ｔ＝tan φ…計算２４
更に、上記の結果の緯度φから、計算２５にて、Ｖの値を求める。
Ｖ＝（１＋ｅ'2cos2φ）1/2 …計算２５
求めた当該Ｖの値から、計算２６にて子午線曲率半径Ｍを、計算２７にて、卯酉線曲率半径Ｎを、更に計算２８から平均曲率半径Ｒを、夫々求める。
Ｍ＝ｃ／Ｖ3 …計算２６
Ｎ＝ｃ／Ｖ…計算２７
Ｒ＝ｃ／Ｖ2 …計算２８
上記結果の緯度φ及び経度λ、更に、上記にて求めた第１の値η、第２の値ｔとから、計算２９にて、当該緯度φ及び経度λにおける縮尺係数ｍを求める。
ｍ＝｛（１＋Δλ2cos2 φ（１＋η2 ）／２＋Δλ4cos4 φ（５−４ｔ2 ）／24）｝×0.9999…計算２９
そして、上記の結果の緯度φと、当該φについての上記子午線曲率半径Ｍ、卯酉線曲率半径Ｎ、平均曲率半径Ｒ、上記の経度差分Δλ、第１の値η、第２の値ｔ、縮尺係数ｍとから、計算３０にて、当該緯度φ及び経度λに対応する平面直角座標上の座標値ｙを求めることができる（計算２９の右辺の0.9999は、数値として原点縮尺係数ｍ0 と同じであるが、前述の通り、これは基準子午線上の縮尺係数を0.9999とするための一定の縮尺係数である）。
ｙ＝｛ΔλＮcos φ＋Δλ3 Ｎcos 3 φ（１−ｔ2 ＋η2 ）／６
＋Δλ5 Ｎcos 5 φ（５−１８ｔ2 ＋ｔ4 ＋１４η2 −５８η2 ｔ2 ）／120 ｝×ｍ…計算３０
【００７７】
そして、地球楕円体の長半径ａと、第１離心率ｅと、上記結果の緯度φと、座標値ｘと、前述の定数Ａ〜Ｆとから、計算３１にて、上記結果の緯度φに対する赤道からの子午線弧長Ｓ (φ )を求める。
Ｓ (φ )＝ａ（１−ｅ2 ）・｛Ａ・φ−（Ｂ／２）・（sin ２φ）
＋（Ｃ／４）・（sin ４φ）−（Ｄ／６）・（sin ６φ）
＋（Ｅ／８）・（sin ８φ）−（Ｆ／１０）・（sin10 φ）｝
…計算３１
上記で求めた緯度φに対する赤道からの子午線弧長Ｓ (φ )と、原点緯度φ0 に対する赤道からの子午線弧長Ｓ (φ0 ) とにより、計算３２にて、その差Ｂを求める。
Ｂ＝Ｓ (φ )−Ｓ (φ0 ) …計算３２
更に、この差Ｂと、上記の結果の緯度φと、当該φについての上記子午線曲率半径Ｍ、卯酉線曲率半径Ｎ、平均曲率半径Ｒ、上記の経度差分Δλ、第１の値η、第２の値ｔ、縮尺係数ｍとから、計算３３にて、にて、当該緯度φ及び経度λに対応する平面直角座標上の座標値ｘを求めることができる。
ｘ＝｛Ｂ＋Δλ2sinφcos φ／２＋Δλ4 Ｎsin φcos3φ（５−ｔ2 ＋９η2 ＋４η4 ）／24＋Δλ6 Ｎsin φcos 5 φ（61−58ｔ2 ＋ｔ4 ＋ 270η2 −330 η2 ｔ2 ）／720 ｝×ｍ…計算３３
また、上記の結果の緯度φと、上記の経度差分Δλ、第１の値η、第２の値ｔとから、計算３４にて、当該緯度φ及び経度λにおける子午線収差γを求めることができる。
γ＝Δλsin φ＋Δλ3 sin φcos2φ（１＋３η2 ＋２η4 ）／３
＋Δλ5 sin φcos4φ（２−ｔ2 ）／15…計算３４
【００７８】
本明細書の背景技術の欄で、国土地理院のプログラムによる計算に使用した値について、本願発明に係る計算装置に入力して計算し、その結果の値を上記の検算法で検算したところ、国土地理院のプログラムが設定してある有効数字の範囲内（算出値の桁数内）において、本願発明に係る計算装置に入力した値（ｘ，ｙ）と、検算で得た値（ｘ，ｙ）とは、完全に一致することが分った。
更に、具体的な実施例について、以下に掲げる。
【実施例】
【００７９】
本願発明に係る装置を用いて、座標系３に属する５箇所の座標値ｘ，ｙ（入力）から、夫々の緯度φ・経度λを算出（出力）した（結果データ１〜５）。ここでは、座標値ｘを固定し、結果データ１から結果データ５へかけて、座標値ｙが原点から遠ざかる点を調べた。結果データ５は、縮尺係数が、1.0001となる点である。
結果データ１〜５で得た緯度・経度から、上記の検算法（検算プログラム）にて、座標値についての検証データ１〜５を得た（検証データの番号は、検証した結果データの番号と対応している）。
各データの計算において、ＷＧＳ８４地球楕円体の長半径ａ、地球楕円体の極の曲率半径ｃ（＝ａ2 ／ｂ）、離心率（第１離心率）ｅの２乗、第２離心率ｅ’の２乗、座標系Ｚ、当該座標系Ｚの原点緯度φ0 、当該座標系Ｚの原点経度λ0 、原点縮尺係数ｍ0 について、次の数値を定数として用いた。
ａ＝6378137.000
ｃ＝6399593.62576
ｅ2 ＝ 0.006694379990141
ｅ’2 ＝0.006739496742276
Ｚ＝3
φ0 ＝36- 0- 0.000
λ0 ＝ 132-10- 0.000
ｍ0 ＝0.9999
また、以下において、
ZAHYOCHI X( 1)は、入力座標値ｘを示し、単位はメートルである。
ZAHYOCHI Y( 1)は、入力座標値ｙを示し、単位はメートルである。
PHAIは、最近似垂線緯度φj+1 を示す。
LATITUDEは、当該点の緯度φn （北緯）を示す。例えば 34- 7-45.93167 は、 34 度 7分45.93167秒（34°7'45".93167 ）を示している。
LONGITUDE は、当該点の経度λ（東経）を示す。
GAMMA ANGLEは、当該点の子午線収差（真北方向角）γを示す。
MERIDIAN RADIUSは、当該点の子午線曲率半径Ｍを示す。
ORTHOG. MERIDIAN RADIUS は、当該点の卯酉線曲率半径Ｎを示す。
MEAN RADIUS は、当該点の平均曲率半径Ｒを示す。
SCALE FACTOR は、当該点の縮尺係数を示す。
【００８０】
（結果データ１）
１ ZAHYOCHI X( 1) = -207504.143
ZAHYOCHI Y( 1) = 0.000
PHAI = 34- 7-45.93167
LATITUDE = 34- 7-45.93167
LONGITUDE =132-10- 0.00000
GAMMA ANGLE = 0- 0- 0.00000
MERIDIAN RADIUS = 6355518.63673
ORTHOG. MERIDIAN RADIUS =6384868.10519
MEAN RADIUS = 6370176.46816
SCALE FACTOR = 0.9999000000
【００８１】
（結果データ２）
2 ZAHYOCHI X( 2) = -207504.143
ZAHYOCHI Y( 2) = -60000.000
PHAI = 34- 7-46.23045
LATITUDE = 34- 7-40.02854
LONGITUDE = 131-30-58.32910
GAMMA ANGLE = - 0-21-53.81142
MERIDIAN RADIUS = 6355516.93662
ORTHOG. MERIDIAN RADIUS =6384867.53587
MEAN RADIUS = 6370175.33215
SCALE FACTOR = 0.9999443586
【００８２】
（結果データ３）
3 ZAHYOCHI X( 3) = -207504.143
ZAHYOCHI Y( 3) = -90091.612
PHAI = 34- 7-46.60518
LATITUDE = 34-7-32.62478
LONGITUDE = 131-11-24.23440
GAMMA ANGLE = -0-32-52.51513
MERIDIAN RADIUS = 6355514.80440
ORTHODG. MERIDIAN RADIUS = 6384866.82185
MEAN RADIUS = 6370173.90739
SCALE FACTOR = 1.0000000000
【００８３】
（結果データ４）
4 ZAHYOCHI X( 4) = -207504.143
ZAHYOCHI Y( 4) = -119361.021
PHAI = 34-7-47.111365
LATITUDE = 34-7-22.57909
LONGITUDE = 130-52-22.58341
GAMMA ANGLE = -0-43-32.98171
MERIDIAN RADIUS = 6355511.91142
ORTHOG. MERIDIAN RADIUS = 6384865.85306
MEAN RADIUS = 6370171.97429
SCALE FACTOR = 1.0000755078
【００８４】
（結果データ５）
5 ZAHYOCHI X( 5) = -207504.143
ZAHYOCHI Y( 5) = -127420.415
PHAI = 34-7-47.27856
LATITUDE = 34-7-19.32099
LONGITUDE = 130-47-8.30603
GAMMA ANGLE = -0-46-29.28380
MERIDIAN RADIUS = 6355510.97317
ORTHOG. MERIDIAN RADIUS = 6384865.53887
MEAN RADIUS = 6370171.34734
SCALE FACTOR = 1.0001000000
【００８５】
（検証データ１）
１ ZAHYOCHI X( 1) = -207504.143
ZAHYOCHI Y( 1) = 0.000
LATITUDE = 34-7-45.93167
LONGITUDE = 132-10-0.00000
GAMMA ANGLE = 0-0-0.00000
MERIDIAN RADIUS = 6355518.63672
ORTHODG. MERIDIAN RADIUS =6384868.10519
MEAN RADIUS = 6370176.46816
SCALE FACTOR = 0.9999000000
【００８６】
（検証データ２）
2 ZAHYOCHI X( 2) = -207504.143
ZAHYOCHI Y( 2) = -60000.000
LATITUDE = 34-7-40.02854
LONGITUDE = 131-30-58.32910
GAMMA ANGLE = -0-21-53.81142
MERIDIAN RADIUS = 6355516.93662
ORTHODG. MERIDIAN RADIUS =6384867.53587
MEAN RADIUS = 6370175.33215
SCALE FACTOR = 0.9999443586
【００８７】
（検証データ３）
3 ZAHYOCHI X( 3) = -207504.143
ZAHYOCHI Y( 3) = -90091.612
LATITUDE = 34-7-32.62478
LONGITUDE = 131-11-24.23440
GAMMA ANGLE = -0-32-52.51514
MERIDIAN RADIUS = 6355514.80440
ORTHODG. MERIDIAN RADIUS = 6384866.82185
MEAN RADIUS = 6370173.90739
SCALE FACTOR = 0.9999999999
【００８８】
（検証データ４）
4 ZAHYOCHI X( 4) = -207504.143
ZAHYOCHI Y( 4) = -119361.021
LATITUDE = 34-7-22.57909
LONGITUDE = 130-52-22.58341
GAMMA ANGLE = -0-43-32.98172
MERIDIAN RADIUS = 6355511.91142
ORTHOG. MERIDIAN RADIUS = 6384865.85306
MEAN RADIUS = 6370171.97429
SCALE FACTOR = 1.0000755076
【００８９】
（検証データ５）
5 ZAHYOCHI X( 5) = -207504.143
ZAHYOCHI Y( 5) = -127420.415
LATITUDE = 34-7-19.32099
LONGITUDE = 130-47-8.30603
GAMMA ANGLE = -0-46-29.28381
MERIDIAN RADIUS = 6355510.97317
ORTHOG. MERIDIAN RADIUS= 6384865.53887
MEAN RADIUS = 6370171.34734
SCALE FACTOR = 1.0000999998
【００９０】
上記の通り、何れも、結果データの入力値と、検証データ出力値は、完全に一致している。ここで求めた緯度・経度は、小数点以下第５位であり、この範囲で誤差が生じている国土地理院のプログラムとでは、その精度の差は歴然としている。
また、結果データ５を見れば、縮尺係数を示す SCALE FACTORも、1.0001に、精度よく収束しているのが分かる。
【図面の簡単な説明】
【００９１】
【図１】本願発明に係る位置計算装置の一実施の形態の説明図である。
【図２】本願発明に係る位置計算方法のフローを示す説明図である。
【図３】本願発明に係る位置計算方法の、垂線緯度の計算についての説明図である。
【符号の説明】
【００９２】
１ 座標入力部
２ 主力部
６ 演算部
７ 演算制御部
８ 判定部
９ 回帰制御部
【特許請求の範囲】
【請求項１】
等角投影法により地球楕円体を投影する平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙを入力することにより、少なくとも当該座標値に対応する地球楕円体上の緯度φn を出力することが可能な位置計算装置において、
座標入力部と、出力部と、演算部と、判定部と、回帰制御部とを備えるものであり、
座標入力部は、上記座標値ｘ，ｙの入力を受け付け、
演算部は、直角座標の座標値ｘが与えられることにより、地球楕円体の赤道から当該座標値ｘ，ｙが属する座標系原点までの子午線弧長に対し、当該座標値ｘを仮の縮尺係数ｍt で除した値を、加算して、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとし、
この子午線弧長Ｗに対応する基準子午線上の緯度、即ち、上記入力座標から基準子午線上に下ろした垂線の緯度を算出するため、
演算部は、地球楕円体の長半径と第１離心率と緯度とから定まる子午線曲率半径について、地球楕円体の長半径ａと第１離心率ｅと共に基準子午線上の任意の緯度φj を初期値として与え、別途に子午線曲率半径を以って赤道から上記緯度φj まで積分することにより、赤道からの基準子午線上の子午線弧長Ｓ（φj ）を算出し、
演算部は、基準子午線上の緯度φj を与えることによって、式１を計算して基準子午線上の緯度φj+1 を漸近値として算出し、
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）…式１
判定部は、上記式１にて算出した基準子午線上の上記漸近緯度φj+1 と基準子午線上の緯度φj との差の絶対値と、基準値との比較を行うものであり、
判定部にて、基準子午線上の、上記漸近緯度φj+1 と緯度φj との差の絶対値が基準値を超えると判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするには十分収束していないと判定した場合、回帰制御部は、漸近緯度φj+1 を上記基準子午線上の緯度φj とし、演算部に基準子午線上の赤道からの子午線弧長Ｓ（φj ）を算出させると共に当該子午線弧長Ｓ（φj ）から上記式１の計算にて再度漸近緯度φj+1 を算出させ、判定部に再度上記比較を行わせるものであり、
判定部が、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値を基準値以内と判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするのに十分収束したと判定した場合、回帰制御部は、当該漸近緯度φj+1 を、入力座標から基準子午線上に下ろした垂線の緯度φk として、仮の縮尺係数ｍt と、当該垂線緯度φk と、当該垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式３を計算して緯度φk+1 を演算部に算出させ、
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …式３
回帰制御部は、上記緯度φk+1 を入力座標に対応する緯度φn として、演算部に、縮尺係数ｍt+1 を算出させ、前記仮の縮尺係数ｍt に代え当該縮尺係数ｍt+1 を用いて求めた、基準子午線上の子午線弧長Ｗにて、上記式１の計算を行わせるものであることを特徴とする位置計算装置。
【請求項２】
上記の仮の縮尺係数とは、演算部に、直角座標の座標値ｙを与えることにより、原点における縮尺係数ｍ0 又はその代替値と、当該座標が属する座標系の原点における平均曲率半径Ｒ0 とから、式２の計算により、演算部が算出した縮尺係数ｍt であり、
ｍt ＝｛１＋（ｙ2 ／２ｍ02Ｒ02）＋（ｙ4 ／24ｍ04Ｒ04）｝×0.9999…式２
上記の原点における平均曲率半径Ｒ0 は、当該座標が属する座標系の原点の緯度φ0 又はその代替値と共に、地球楕円体の長半径ａ及び短半径ｂ、或いは当該半径ａ，ｂから導かれる第１離心率ｅ、第２離心率ｅ’を用いて定めたものであることを特徴とする請求項１記載の位置計算装置。
【請求項３】
上記縮尺係数ｍt+1 は、原点縮尺係数ｍ0 に代え上記縮尺係数ｍt を用い、且つ、原点における平均曲率半径Ｒ0 に代え演算部に式３にて算出させた緯度φn における平均曲率半径Ｒφn を用いて、演算部が、上記式２により算出するものであることを特徴とする請求項２記載の位置計算装置。
【請求項４】
上記にあって、判定部が、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合に、回帰制御部は、上記緯度φn に対応する縮尺係数ｍt+1 と上記仮の縮尺係数ｍt との差の絶対値と、所定値との比較を、更に判定部にさせ、判定部が当該縮尺係数ｍt+1 と上記仮の縮尺数値ｍt との差の絶対値が所定値を超えると判定した場合にのみ、上記の通り、前記の仮の縮尺係数ｍt に代え、当該縮尺係数ｍt+1 を用いて求めた基準子午線上の子午線弧長Ｗにて、上記式１の計算を演算部に行わせるものであることを特徴とする請求項１乃至３の何れかに記載の位置計算装置。
【請求項５】
上記にあって、判定部が、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、回帰制御部は、垂線緯度φk として上記漸近緯度φj+1 を演算部に与えるものであり、演算部は、上記の仮の縮尺係数ｍt と、上記垂線緯度φk における卯酉線曲率半径Ｎk と、当該垂線緯度φk と第２離心率とによりｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式４にて、経度差分Δλを算出することが可能であり、
Δλ＝（ｙ／ｍt ）／Ｎk cos φk
−（ｙ／ｍt ）3 （１＋２ｔk 2 ＋ηk 2 ）／６Ｎk 3cosφk
＋（ｙ／ｍt ）5 （５＋２８ｔk 2 ＋２４ｔk 4 ）／１２０Ｎk 5cosφk …式４
演算部は、当該座標系が属する原点における経度λ0 に対し当該経度差分Δλを加算することにより、上記座標値ｘ，ｙにおける経度λを算出するものであることを特徴とする請求項１乃至４の何れかに記載の位置計算装置。
【請求項６】
上記にあって、判定部が、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値を基準値以内と判定した場合、回帰制御部は、垂線緯度φk としてこの漸近緯度φj+1 を演算部に与えるものであり、演算部は、上記の縮尺係数ｍt と、当該垂線緯度φk における、子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、当該垂線緯度φk と第２離心率とによりｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、
γ＝（ｙ／ｍt ）ｔk ／Ｎk
−（ｙ／ｍt ）3 ｔk （１＋ｔk 2 −ηk 2 ）／３Ｎk 3
＋（ｙ／ｍt ）5 ｔk （２＋５ｔk 2 ＋３ｔk 4 ）／１５Ｎk 5 …式５
式５にて、座標ｘ，ｙにおける子午線収差γを算出することが可能なことを特徴とする請求項１乃至５の何れかに記載の位置計算装置。
【請求項７】
等角投影法により地球楕円体を投影する平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙを入力することにより、少なくとも当該座標値に対応する地球楕円体上の緯度φn を出力することが可能な位置計算プログラムにおいて、
座標入力部と、出力部と、演算部と、判定部と、回帰制御部とを備えたコンピュータに導入されるものであり、
座標入力部にて、上記座標値ｘ，ｙの入力を受け付け、
演算部に、直角座標の座標値ｘを与えて、地球楕円体の赤道から当該座標値ｘ，ｙが属する座標系原点までの子午線弧長に対し、当該座標値ｘを仮の縮尺係数ｍt で除した値を加算して、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとし、
この子午線弧長Ｗに対応する基準子午線上の緯度、即ち、上記入力座標から基準子午線上に下ろした垂線の緯度を算出するため、
演算部にて、地球楕円体の長半径と第１離心率と緯度とから定まる子午線曲率半径について、地球楕円体の長半径ａと第１離心率ｅと共に任意の緯度φj を初期値として与え、別途に子午線曲率半径を以って赤道から上記緯度φj まで積分することにより、赤道からの基準子午線上の子午線弧長Ｓ（φj ）を算出し、
演算部に、基準子午線上の緯度φj を与えることによって、式１を計算して基準子午線上の緯度φj+1 を漸近値として算出させ、
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）…式１
判定部に、上記式１にて算出した基準子午線上の緯度φj+1 と基準子午線上の緯度φj との差の絶対値と、基準値との比較を行わせ、
判定部にて、基準子午線上の、上記漸近緯度φj+1 と緯度φj との差の絶対値が基準値を超えると判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするには十分収束していないと判定した場合、回帰制御部により、漸近緯度φj+1 を上記基準子午線上の緯度φj とし、演算部に基準子午線上の赤道からの子午線弧長Ｓ（φj ）を算出させると共に当該子午線弧長Ｓ（φj ）から上記の式１の計算にて再度漸近緯度φj+1 を算出させ、判定部に再度上記比較を行わせるものであり、
判定部にて、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするのに十分収束したと判定した場合、
当該漸近緯度φj+1 を、入力座標から基準子午線上に下ろした垂線の緯度φk として、回帰制御部により、仮の縮尺係数ｍt と、当該垂線緯度φk と、当該垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式３を計算して緯度φk+1 を演算部に算出させ、
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …式３
回帰制御部により、上記緯度φk+1 を入力座標に対応する緯度φn として、演算部に、縮尺係数ｍt+1 を算出させ、前記仮の縮尺係数ｍt に代え当該縮尺係数ｍt+1 を用いて求めた、基準子午線上の子午線弧長Ｗにて、上記式１の計算を行わせるものであることを特徴とする位置計算プログラム。
【請求項８】
請求項７記載の位置計算プログラムの記録媒体。
【請求項９】
等角投影法により地球楕円体を投影する平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙを入力することにより、少なくとも当該座標値に対応する地球楕円体上の緯度φn を出力することが可能な位置計算方法において、
座標入力部と、出力部と、演算部と、判定部と、回帰制御部とを備えたコンピュータを用いるものであり、
座標入力部にて、上記座標値ｘ，ｙの入力を受け付け、
演算部に、直角座標の座標値ｘを与えて、地球楕円体の赤道から当該座標値ｘ，ｙが属する座標系原点までの子午線弧長に対し、当該座標値ｘを仮の縮尺係数ｍt で除した値を加算して、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとし、
演算部にて、地球楕円体の長半径と第１離心率と緯度とから定まる子午線曲率半径について、地球楕円体の長半径ａと第１離心率ｅと共に任意の緯度φj を初期値として与え、別途に子午線曲率半径を以って赤道から上記緯度φj まで積分することにより、赤道からの基準子午線上の子午線弧長Ｓ（φj ）を算出し、
演算部に、基準子午線上の緯度φj を与えることによって、式１を計算して基準子午線上の緯度φj+1 を漸近値として算出させ、
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）…式１
判定部にて、上記の式１にて算出した基準子午線上の上記漸近緯度φj+1 と基準子午線上の緯度φj との差の絶対値と、基準値との比較を行わせるものであり、
判定部にて、上記基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値を超えると判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするには十分収束していないと判定した場合、回帰制御部により、漸近緯度φj+1 を上記基準子午線上の緯度φj とし、演算部に基準子午線上の赤道からの子午線弧長Ｓ（φj ）を算出させると共に当該子午線弧長Ｓ（φj ）から上記の式１の計算にて再度漸近緯度φj+1 を算出させ、判定部に再度上記比較を行わせるものであり、
判定部にて、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするのに十分収束したと判定した場合、回帰制御部は、当該漸近緯度φj+1 を、入力座標から基準子午線上に下ろした垂線の緯度φk として、仮の縮尺係数ｍt と、当該垂線緯度φk と、当該垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式３を計算して緯度φk+1 を演算部に算出させ、
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …式３
回帰制御部により、上記緯度φk+1 を入力座標に対応する緯度φn として、演算部に、縮尺係数ｍt+1 を算出させ、前記仮の縮尺係数ｍt に代え当該縮尺係数ｍt+1 を用いて求めた、基準子午線上の子午線弧長Ｗにて、上記式１の計算を行わせるものであることを特徴とする位置計算方法。
【請求項１】
等角投影法により地球楕円体を投影する平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙを入力することにより、少なくとも当該座標値に対応する地球楕円体上の緯度φn を出力することが可能な位置計算装置において、
座標入力部と、出力部と、演算部と、判定部と、回帰制御部とを備えるものであり、
座標入力部は、上記座標値ｘ，ｙの入力を受け付け、
演算部は、直角座標の座標値ｘが与えられることにより、地球楕円体の赤道から当該座標値ｘ，ｙが属する座標系原点までの子午線弧長に対し、当該座標値ｘを仮の縮尺係数ｍt で除した値を、加算して、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとし、
この子午線弧長Ｗに対応する基準子午線上の緯度、即ち、上記入力座標から基準子午線上に下ろした垂線の緯度を算出するため、
演算部は、地球楕円体の長半径と第１離心率と緯度とから定まる子午線曲率半径について、地球楕円体の長半径ａと第１離心率ｅと共に基準子午線上の任意の緯度φj を初期値として与え、別途に子午線曲率半径を以って赤道から上記緯度φj まで積分することにより、赤道からの基準子午線上の子午線弧長Ｓ（φj ）を算出し、
演算部は、基準子午線上の緯度φj を与えることによって、式１を計算して基準子午線上の緯度φj+1 を漸近値として算出し、
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）…式１
判定部は、上記式１にて算出した基準子午線上の上記漸近緯度φj+1 と基準子午線上の緯度φj との差の絶対値と、基準値との比較を行うものであり、
判定部にて、基準子午線上の、上記漸近緯度φj+1 と緯度φj との差の絶対値が基準値を超えると判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするには十分収束していないと判定した場合、回帰制御部は、漸近緯度φj+1 を上記基準子午線上の緯度φj とし、演算部に基準子午線上の赤道からの子午線弧長Ｓ（φj ）を算出させると共に当該子午線弧長Ｓ（φj ）から上記式１の計算にて再度漸近緯度φj+1 を算出させ、判定部に再度上記比較を行わせるものであり、
判定部が、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値を基準値以内と判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするのに十分収束したと判定した場合、回帰制御部は、当該漸近緯度φj+1 を、入力座標から基準子午線上に下ろした垂線の緯度φk として、仮の縮尺係数ｍt と、当該垂線緯度φk と、当該垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式３を計算して緯度φk+1 を演算部に算出させ、
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …式３
回帰制御部は、上記緯度φk+1 を入力座標に対応する緯度φn として、演算部に、縮尺係数ｍt+1 を算出させ、前記仮の縮尺係数ｍt に代え当該縮尺係数ｍt+1 を用いて求めた、基準子午線上の子午線弧長Ｗにて、上記式１の計算を行わせるものであることを特徴とする位置計算装置。
【請求項２】
上記の仮の縮尺係数とは、演算部に、直角座標の座標値ｙを与えることにより、原点における縮尺係数ｍ0 又はその代替値と、当該座標が属する座標系の原点における平均曲率半径Ｒ0 とから、式２の計算により、演算部が算出した縮尺係数ｍt であり、
ｍt ＝｛１＋（ｙ2 ／２ｍ02Ｒ02）＋（ｙ4 ／24ｍ04Ｒ04）｝×0.9999…式２
上記の原点における平均曲率半径Ｒ0 は、当該座標が属する座標系の原点の緯度φ0 又はその代替値と共に、地球楕円体の長半径ａ及び短半径ｂ、或いは当該半径ａ，ｂから導かれる第１離心率ｅ、第２離心率ｅ’を用いて定めたものであることを特徴とする請求項１記載の位置計算装置。
【請求項３】
上記縮尺係数ｍt+1 は、原点縮尺係数ｍ0 に代え上記縮尺係数ｍt を用い、且つ、原点における平均曲率半径Ｒ0 に代え演算部に式３にて算出させた緯度φn における平均曲率半径Ｒφn を用いて、演算部が、上記式２により算出するものであることを特徴とする請求項２記載の位置計算装置。
【請求項４】
上記にあって、判定部が、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合に、回帰制御部は、上記緯度φn に対応する縮尺係数ｍt+1 と上記仮の縮尺係数ｍt との差の絶対値と、所定値との比較を、更に判定部にさせ、判定部が当該縮尺係数ｍt+1 と上記仮の縮尺数値ｍt との差の絶対値が所定値を超えると判定した場合にのみ、上記の通り、前記の仮の縮尺係数ｍt に代え、当該縮尺係数ｍt+1 を用いて求めた基準子午線上の子午線弧長Ｗにて、上記式１の計算を演算部に行わせるものであることを特徴とする請求項１乃至３の何れかに記載の位置計算装置。
【請求項５】
上記にあって、判定部が、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、回帰制御部は、垂線緯度φk として上記漸近緯度φj+1 を演算部に与えるものであり、演算部は、上記の仮の縮尺係数ｍt と、上記垂線緯度φk における卯酉線曲率半径Ｎk と、当該垂線緯度φk と第２離心率とによりｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式４にて、経度差分Δλを算出することが可能であり、
Δλ＝（ｙ／ｍt ）／Ｎk cos φk
−（ｙ／ｍt ）3 （１＋２ｔk 2 ＋ηk 2 ）／６Ｎk 3cosφk
＋（ｙ／ｍt ）5 （５＋２８ｔk 2 ＋２４ｔk 4 ）／１２０Ｎk 5cosφk …式４
演算部は、当該座標系が属する原点における経度λ0 に対し当該経度差分Δλを加算することにより、上記座標値ｘ，ｙにおける経度λを算出するものであることを特徴とする請求項１乃至４の何れかに記載の位置計算装置。
【請求項６】
上記にあって、判定部が、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値を基準値以内と判定した場合、回帰制御部は、垂線緯度φk としてこの漸近緯度φj+1 を演算部に与えるものであり、演算部は、上記の縮尺係数ｍt と、当該垂線緯度φk における、子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、当該垂線緯度φk と第２離心率とによりｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、
γ＝（ｙ／ｍt ）ｔk ／Ｎk
−（ｙ／ｍt ）3 ｔk （１＋ｔk 2 −ηk 2 ）／３Ｎk 3
＋（ｙ／ｍt ）5 ｔk （２＋５ｔk 2 ＋３ｔk 4 ）／１５Ｎk 5 …式５
式５にて、座標ｘ，ｙにおける子午線収差γを算出することが可能なことを特徴とする請求項１乃至５の何れかに記載の位置計算装置。
【請求項７】
等角投影法により地球楕円体を投影する平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙを入力することにより、少なくとも当該座標値に対応する地球楕円体上の緯度φn を出力することが可能な位置計算プログラムにおいて、
座標入力部と、出力部と、演算部と、判定部と、回帰制御部とを備えたコンピュータに導入されるものであり、
座標入力部にて、上記座標値ｘ，ｙの入力を受け付け、
演算部に、直角座標の座標値ｘを与えて、地球楕円体の赤道から当該座標値ｘ，ｙが属する座標系原点までの子午線弧長に対し、当該座標値ｘを仮の縮尺係数ｍt で除した値を加算して、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとし、
この子午線弧長Ｗに対応する基準子午線上の緯度、即ち、上記入力座標から基準子午線上に下ろした垂線の緯度を算出するため、
演算部にて、地球楕円体の長半径と第１離心率と緯度とから定まる子午線曲率半径について、地球楕円体の長半径ａと第１離心率ｅと共に任意の緯度φj を初期値として与え、別途に子午線曲率半径を以って赤道から上記緯度φj まで積分することにより、赤道からの基準子午線上の子午線弧長Ｓ（φj ）を算出し、
演算部に、基準子午線上の緯度φj を与えることによって、式１を計算して基準子午線上の緯度φj+1 を漸近値として算出させ、
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）…式１
判定部に、上記式１にて算出した基準子午線上の緯度φj+1 と基準子午線上の緯度φj との差の絶対値と、基準値との比較を行わせ、
判定部にて、基準子午線上の、上記漸近緯度φj+1 と緯度φj との差の絶対値が基準値を超えると判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするには十分収束していないと判定した場合、回帰制御部により、漸近緯度φj+1 を上記基準子午線上の緯度φj とし、演算部に基準子午線上の赤道からの子午線弧長Ｓ（φj ）を算出させると共に当該子午線弧長Ｓ（φj ）から上記の式１の計算にて再度漸近緯度φj+1 を算出させ、判定部に再度上記比較を行わせるものであり、
判定部にて、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするのに十分収束したと判定した場合、
当該漸近緯度φj+1 を、入力座標から基準子午線上に下ろした垂線の緯度φk として、回帰制御部により、仮の縮尺係数ｍt と、当該垂線緯度φk と、当該垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式３を計算して緯度φk+1 を演算部に算出させ、
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …式３
回帰制御部により、上記緯度φk+1 を入力座標に対応する緯度φn として、演算部に、縮尺係数ｍt+1 を算出させ、前記仮の縮尺係数ｍt に代え当該縮尺係数ｍt+1 を用いて求めた、基準子午線上の子午線弧長Ｗにて、上記式１の計算を行わせるものであることを特徴とする位置計算プログラム。
【請求項８】
請求項７記載の位置計算プログラムの記録媒体。
【請求項９】
等角投影法により地球楕円体を投影する平面直角座標Ｘ，Ｙ上の座標値ｘ，ｙを入力することにより、少なくとも当該座標値に対応する地球楕円体上の緯度φn を出力することが可能な位置計算方法において、
座標入力部と、出力部と、演算部と、判定部と、回帰制御部とを備えたコンピュータを用いるものであり、
座標入力部にて、上記座標値ｘ，ｙの入力を受け付け、
演算部に、直角座標の座標値ｘを与えて、地球楕円体の赤道から当該座標値ｘ，ｙが属する座標系原点までの子午線弧長に対し、当該座標値ｘを仮の縮尺係数ｍt で除した値を加算して、赤道から当該座標値ｘまでの基準子午線上の子午線弧長Ｗとし、
演算部にて、地球楕円体の長半径と第１離心率と緯度とから定まる子午線曲率半径について、地球楕円体の長半径ａと第１離心率ｅと共に任意の緯度φj を初期値として与え、別途に子午線曲率半径を以って赤道から上記緯度φj まで積分することにより、赤道からの基準子午線上の子午線弧長Ｓ（φj ）を算出し、
演算部に、基準子午線上の緯度φj を与えることによって、式１を計算して基準子午線上の緯度φj+1 を漸近値として算出させ、
φj+1 ＝φj ＋（Ｗ−Ｓ (φj ) ）（１−ｅ2 sin2φj ）3/2 ／ａ（１−ｅ2 ）…式１
判定部にて、上記の式１にて算出した基準子午線上の上記漸近緯度φj+1 と基準子午線上の緯度φj との差の絶対値と、基準値との比較を行わせるものであり、
判定部にて、上記基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値を超えると判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするには十分収束していないと判定した場合、回帰制御部により、漸近緯度φj+1 を上記基準子午線上の緯度φj とし、演算部に基準子午線上の赤道からの子午線弧長Ｓ（φj ）を算出させると共に当該子午線弧長Ｓ（φj ）から上記の式１の計算にて再度漸近緯度φj+1 を算出させ、判定部に再度上記比較を行わせるものであり、
判定部にて、基準子午線上の、漸近緯度φj+1 と緯度φj との差の絶対値が基準値以内と判定した場合、即ち、この漸近緯度φj+1 を、上記入力座標から基準子午線上に下ろした垂線の緯度とするのに十分収束したと判定した場合、回帰制御部は、当該漸近緯度φj+1 を、入力座標から基準子午線上に下ろした垂線の緯度φk として、仮の縮尺係数ｍt と、当該垂線緯度φk と、当該垂線緯度φk における子午線曲率半径Ｍk 及び卯酉線曲率半径Ｎk と、ｅ’cos φk にて定まる第１の値ηk と、tan φk にて定まる第２の値ｔk と、上記座標値ｙとから、式３を計算して緯度φk+1 を演算部に算出させ、
φk+1 ＝φk −ｔk （ｙ／ｍt ）2 ／２Ｍk Ｎk ＋ｔk （ｙ／ｍt ）4 （５＋３ｔk 2 ＋ηk 2 −９ｔk 2 ηk 2 −４ηk 4 ）／24Ｍk Ｎk 3 −ｔk （ｙ／ｍt ）6 （61＋90ｔk 2 ＋45ｔk 4 ）／720 Ｍk Ｎk 5 …式３
回帰制御部により、上記緯度φk+1 を入力座標に対応する緯度φn として、演算部に、縮尺係数ｍt+1 を算出させ、前記仮の縮尺係数ｍt に代え当該縮尺係数ｍt+1 を用いて求めた、基準子午線上の子午線弧長Ｗにて、上記式１の計算を行わせるものであることを特徴とする位置計算方法。
【図１】
【図２】
【図３】
【図２】
【図３】
【公開番号】特開２００６−６４８９２（Ｐ２００６−６４８９２Ａ）
【公開日】平成１８年３月９日（２００６．３．９）
【国際特許分類】
【出願番号】特願２００４−２４６０６４（Ｐ２００４−２４６０６４）
【出願日】平成１６年８月２６日（２００４．８．２６）
【出願人】（５０４３２５４５０）有限会社長谷川商事 (3)
【Ｆターム（参考）】
【公開日】平成１８年３月９日（２００６．３．９）
【国際特許分類】
【出願日】平成１６年８月２６日（２００４．８．２６）
【出願人】（５０４３２５４５０）有限会社長谷川商事 (3)
【Ｆターム（参考）】
[ Back to top ]