天体の位置計算 - 赤道座標の方向余弦の表示

天体の位置計算 – 赤道座標の方向余弦表示

「天球座標-赤道座標の方向余弦表示」34頁

「天体の位置計算増補版」では、ここから行列計算をおこないます。そこで、今回からFortranの行列を利用して恒星の位置計算を行っていきます。行列を使わない式もでてきますが、とくに断りのない限り行列を使用することにします。Fortranの行列計算で積だけはmatmulという関数を使いますが、和や差については、配列同士の計算で通常の＋やーで計算できます。また、Fortranの行列は配列として計算しますが、配列の中に直接数値を記述しなくても、変数や数式が使えます。

赤道座標系から地平座標の方向余弦変換
天球上の位置を方向余弦で表すとき。赤道座標に対して、X軸は春分手方向、Y軸は赤経6h、赤緯０°、Z軸は天の北極方向へとるとすると、赤経αと赤緯δの方向余弦（L, M,N）は次の関係式で表すことができます。

\[ \begin{align} L &= \cos a = \cos\delta\cos\alpha\\ M &= \cos b = \cos\delta\sin\alpha\\ N &= \cos c = \sin\delta\\ \end{align} \] \[ 引用元: 「天体の位置計算増補版」34頁式2-18 \]

Fortranによるコード
次回の行列の計算の前に、赤道座標の方向余弦を求めてみます。サブルーチンはありません。

▶ メインルーチン pcc_matrix.f90

program matrix
    use pcclib

    implicit none
    
    ! 恒星データ
    ! character(16), parameter :: SNM = 'α CMa  Sirius  '              ! 恒星名
    ! double precision, parameter :: SRA(1:3) = (/6d0, 42d0, 56.714d0/)   ! 赤経
    ! double precision, parameter :: SDC(1:3) = (/-16d0, 38d0, 46.36d0/)  ! 赤緯
    ! double precision, parameter :: SRA_PROPER_MOTION = -0.03791d0     ! 赤経方向の固有運動量
    ! double precision, parameter :: SDC_PROPER_MOTION = -1.2114d0      ! 赤緯方向の固有運動量
    ! double precision, parameter :: SRV = -7.6d0                       ! 視線速度
    ! double precision, parameter :: SAP = 0.377d0                      ! 年周視差

    !検証用恒星データ
    character(16), parameter :: SNM = '61 Cyg          ';           ! 恒星名
    double precision, parameter :: SRA(1:3) = (/21d0, 4d0, 39.935d0/) ! 赤経
    double precision, parameter :: SDC(1:3) = (/38d0, 29d0, 59.10d0/) ! 赤緯
    double precision, parameter :: SRA_PROPER_MOTION = 0.35227d0    ! 赤経方向の固有運動量
    double precision, parameter :: SDC_PROPER_MOTION = 3.1847d0     ! 赤緯方向の固有運動量
    double precision, parameter :: SRV = -64.3d0                    ! 視線速度
    double precision, parameter :: SAP = 0.296d0                    ! 年周視差

    double precision :: ra_rad, dc_rad
    double precision :: l, m, n
    character(255) :: dfmt

    ra_rad = deg2rad(hms2deg(SRA(1), SRA(2), SRA(3)))
    dc_rad = deg2rad(dms2deg(SDC(1), SDC(2), SDC(3)))

    ! 方向余弦
    l = cos(dc_rad) * cos(ra_rad)
    m = cos(dc_rad) * sin(ra_rad)
    n = sin(dc_rad)

    ! 出力
    dfmt = '(a16, "      ", f11.8, "     ",  f11.8, "     ", f11.8, "    ", f11.8)'
    write(*, *)
    write(*, *) '星名                      L               M               N            二乗和'
    write(*, *) '-------------------------------------------------------------------------------'
    write(*, dfmt) SNM, l, m, n, l**2 + m**2 + n**2
    write(*, *)

    stop
end program matrix

コンパイル

gfortran pcc_lib.f90 pcc_matrix.f90

実行

./a.out

実行結果

 星名                      L               M               N            二乗和
 -------------------------------------------------------------------------------
α CMa  Sirius        -0.17848222      0.94132039     -0.28646120     1.00000000

老猫こてつ