どせいたんさき。

ナスダヨー

銀河座標から赤道座標に変換するプログラムを書いてみた

精度はいらないけれど銀河座標 (l,b) から赤道座標 (RA.,Dec.) の値を計算する必要性が出てきたのでちょいちょいっと作ってみた.必要性がありそうだから変換ツールとかどこかに落ちていると思ったのだけどなぁ...全然見つからなかった.仕方ないので自作.

使用した言語は octave なので, matlab 使っている人でも利用できるはず. octave ユーザーの人であれば,以下にあるソースをコピペして gc2ec.m という名前でパスの通ったところに保存すれば使える.


やっていることは至極簡単.

  • 銀河座標 (l,b) から極座標 (x,y,z) に変換
  • 赤道座標系とマッチするように変換行列をかけて座標を回転させる
  • 変換後の極座標 (x',y',z') から赤道座標 (RA.,Dec.) に変換
  • フォーマットに合わせて出力

エラー処理とかはまったくしていないので変な値を代入するととんでもないことになるかも.

つかいかた
[RA Dec] = gc2ec (l,b);

RA,Dec は出力結果を代入する変数. l,b は変換したい銀河座標.銀河座標は degree で入力すること. radian じゃないよ. RA は hh:mm:ss.sss で, Dec は dd:mm:ss.sss のフォーマットで出力される.どちらも文字列として出力される.数値としてさらに計算に使いたい場合はプログラム内の deg2hms とか deg2dms の処理をはさまなければよい.

そーす

出てくる関数とそれぞれの役割は次のとおり.ただし xyz2lb とかとてもひどい関数(レポートで提出したら突き返されるかも)なのであまり参考にはしない方がいいかも.

  • gc2ec: メインプログラム
  • deg2hms: 角度を hh:mm:ss.sss のフォーマットで表示
  • deg2dms: 角度を dd:mm:ss.sss のフォーマットで表示
  • lb2xyz: 角度 l,b を x,y,z の座標で表示
  • xyz2lb: 座標 x,y,z を l,b の角度にて表示
  • vecprod: 外積を計算.ただし絶対値が 1 のベクトルを返す
## Copyright (C) 1995, 1996, 1997, 1999, 2000, 2002, 2005, 2006, 2007
##               John W. Eaton
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## Author: Ryou


function [RA DEC] = gc2ec(l,b)
  Gzero.l = 266.4051000; Gzero.b = -28.9167889;
  Gpole.l = 192.8595083; Gpole.b =  27.1283361;

  G.zero = lb2xyz(Gzero.l,Gzero.b);
  G.pole = lb2xyz(Gpole.l,Gpole.b);
  G.perp = vecprod(G.pole,G.zero);

  Cmat = [G.zero G.perp G.pole];

  GCxyz = lb2xyz(l,b);
  ECxyz = Cmat*GCxyz;

  EClb = xyz2lb(ECxyz);

  RA = deg2hms(EClb(1));
  DEC = deg2dms(EClb(2));
endfunction

function dms = deg2dms(deg)
  sgn = sign(deg); deg = abs(deg);
  d = floor(deg);
  m = floor((deg-d)*60);
  s = (((deg-d)*60-m)*60);
  dms = [sprintf("%02d:",sgn*d)  \
         sprintf("%02d:",m)      \
         sprintf("%06.3f",s) ];
endfunction

function hms = deg2hms(deg)
  h = floor(deg/15);
  m = floor((deg/15-h)*60);
  s = (((deg/15-h)*60-m)*60);
  hms = [sprintf("%02d:",h)    \
         sprintf("%02d:",m)    \
         sprintf("%06.3f",s) ];
endfunction

function X = lb2xyz(l,b)
  theta = pi/180*(90-b);
  phi = pi/180*l;
  X = [ sin(theta)*cos(phi); \
        sin(theta)*sin(phi); \
        cos(theta) ];
endfunction

function X = xyz2lb(v)
  phi = atan(v(2)/v(1));
  theta = atan(norm(v(3)/[v(1);v(2)]));
  phi = phi/pi*180 + \
        45*(1-sign(v(1)))*(1+sign(v(2))) + \
        45*(1-sign(v(1)))*(1-sign(v(2))) + \
        90*(1+sign(v(1)))*(1-sign(v(2)));
  theta = theta/pi*180;
  X = [ phi theta ];
endfunction

function c = vecprod(a,b)
  x = [a(2)*b(3) - a(3)*b(2);  \
       a(3)*b(1) - a(1)*b(3);  \
       a(1)*b(2) - a(2)*b(1) ];
  c = x./norm(x);
endfunction