天体の位置計算 – 恒星の高度・方位角計算
「天球座標-赤道座標から地平座標への変換」28-30頁
恒星の位置（赤経・赤緯）から高度・方位角を求めます。いわゆる赤道座標から地平座標への変換になります。この変換はよく使われるもので、赤経・赤緯をいわれてもどこ？となりがちですが、高度はいくつ、方位角は北から測って何度というようにすれば大体の位置がわかります。恒星を例に取り上げていますが、赤経および赤緯がわかれば太陽、月、惑星や彗星、小惑星などの天体の高度・方位角を求めることができます。
赤道座標から地平座標の計算式
\[ \begin{array}{lr} \sin h &=& \sin\phi\sin\delta + \cos\phi\cos\delta\cos(\theta\,-\,\alpha)\\ \cos h\cos A &=& -\cos\phi\sin\delta + \sin\phi\cos\delta\cos(\theta\,-\,\alpha)\\ \cos h\sin A &=& \cos\delta\sin(\theta\,-\,\alpha) \end{array} \] \[ 方位角(\tan^{-1} A) = \frac{\cos\delta\sin(\theta\,-\,\alpha)}{-\cos\phi\sin\delta + \sin\phi\cos\delta\cos(\theta\,-\,\alpha)} \]象限について
赤経や方位角、黄経は0°から360°をとることになります。ですがtan-1は−90°から＋90°しかとることができないために象限の判断をするときに分母や分子の正負を判断して決める必要があります。分母が0より小さいときは＋180°、分子が0より小さいときは＋360°を求めた方位角の値に加算するという手順で求めます。分母および分子とも0より小さいときは分母で判断します。計算例では分母が0より小さいために＋180°を加算しています。
象限の判断を行うサブルーチンは、共通サブルーチンの中に含まれていますので、あえて掲載していません。
Fortranによるコード
▶ サブルーチン AltitudeAzimuth
function AltitudeAzimuth(lsdt, la, ra, dc) result(hc)
! 赤道座標から地平座標への変換
! 「天体の位置計算 増補版」28-31頁
double precision, intent(in) :: lsdt, la, ra, dc
double precision :: hg, al, hi, ss, cc, hc(2)
! 時角の計算
hg = lsdt - ra
! 方位角の計算
ss = cos(dc) * sin(hg)
cc = -cos(la) * sin(dc) + sin(la) * cos(dc) * cos(hg)
al = Quadrant(ss, cc)
! 高度の計算
hi = asin(sin(la) * sin(dc) + cos(la) * cos(dc) * cos(hg))
hc(1) = al
hc(2) = hi
end function AltitudeAzimuth▶ メインルーチン pcc_ec2hc.f90
program pcc_ec2hc
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 ! 年周視差
! 地点データ
character(21), parameter :: OBS_NAME = 'Tokyo Observatory PZT' ! 地点名
double precision, parameter :: OBS_LG(1:3) = (/9d0, 18d0, 9.936d0/) ! 経度
double precision, parameter :: OBS_LA(1:3) = (/35d0, 40d0, 20.707d0/) ! 緯度
double precision, parameter :: GAST(1:3) = (/17d0, 51d0, 24.267d0/)
! 検証用地点データ
! character(21), parameter :: OBS_NAME = 'Ngasawa home ' ! 地点名
! double precision, parameter :: OBS_LG(1:3) = (/9d0, 18d0, 7.573d0/) ! 経度
! double precision, parameter :: OBS_LA(1:3) = (/35d0, 47d0, 20.0d0/) ! 緯度
! double precision, parameter :: GAST(1:3) = (/17d0, 11d0, 58.714d0/)
double precision :: dy, dt, julian
double precision :: t
double precision :: year, month, day, hour, minute, second
double precision :: ra, dc, la, al, hi
double precision :: rh, rm, rs, dh, dm, ds
double precision :: gasdt, lsdt, long, cv
double precision :: sdt_hour, sdt_minute, sdt_second
double precision :: hc(2), sdt(3)
character(255) :: sfmt, dfmt
character(1) :: sg
! 日付と時刻の入力
write(*,*)
write(*, fmt='(a)', advance='no') 'DATE AND TIME(JST) ? '
read(*,*) dy, dt
! 日付と時刻を抽出
year = int(dy / 10000.0d0)
month = mod(int(dy / 100.0d0), 100)
day = mod(dy, 100.0d0)
hour = int(dt / 10000.0d0)
minute = mod(int(dt / 100.0d0), 100)
second = mod(dt, 100.0d0)
! ユリウス日の計算
julian = JulianDay(year, month, day)
julian = julian + (hour / 24d0 + minute / 1440d0 + second / 86400d0)
julian = julian - 0.375
! グリニッジ視恒星時を時の小数点に換算
gasdt = dms2deg(GAST(1), GAST(2), GAST(3))
! 経度を時の小数点に換算
long = dms2deg(OBS_LG(1), OBS_LG(2), OBS_LG(3))
! 世界0時からの経過時間を時の小数点に換算
hour = hour - 9d0;
t = dms2deg(hour, minute, second)
! 補正値の計算
cv = (hour * 3600d0 + minute * 60d0 + second) * 0.00273791d0
cv = cv / 3600d0
! 地方恒星時
lsdt = gasdt + long + t + cv
sdt = deg2dms(lsdt)
sdt_hour = sdt(1)
sdt_minute = sdt(2)
sdt_second = sdt(3)
! 恒星時が範囲を超えたときの処理
if (sdt_hour > 24d0) then
sdt_hour = sdt_hour - 24d0
end if
if (sdt_second >= 60d0) then
sdt_second = sdt_second - 60d0
sdt_minute = sdt_minute + 1d0
end if
if (sdt_minute == 60d0) then
sdt_minute = 0d0
sdt_hour = sdt_hour + 1d0
end if
if (sdt_hour == 24d0) then
sdt_hour = 0d0
end if
! 高度・方位角
ra = deg2rad(hms2deg(SRA(1), SRA(2), SRA(3)))
dc = deg2rad(dms2deg(SDC(1), SDC(2), SDC(3)))
la = deg2rad(dms2deg(OBS_LA(1), OBS_LA(2), OBS_LA(3)))
lsdt = deg2rad(hms2deg(sdt_hour, sdt_minute, sdt_second))
hc = AltitudeAzimuth(lsdt, la, ra, dc)
al = rad2deg(hc(1))
hi = rad2deg(hc(2))
! 赤経・赤緯をそれぞれの変数に代入
rh = SRA(1)
rm = SRA(2)
rs = SRA(3)
dh = SDC(1)
dm = SDC(2)
ds = SDC(3)
sg = '+'
if (dh < 0) then
dh = abs(dh)
sg = '-'
end if
! 出力
! 角度の譚雨を表す°記号は出力に合わせて調整すること
sfmt = '(i4, "年 ", i2, "月 ", i2, "日 ", i2, "時 ", i2, "分 ", i2, "秒 (JST)")'
dfmt = '(a16, " ", i2, " ", i2, " ", f6.3, " ", a1, i2, " ", i2, " ", f5.2, " ", f10.6, " ", f10.6)'
write(*, *)
write(*, sfmt) int(year), int(month), int(day), int(hour + 9), int(minute), int(second)
write(*, '("Julian Date = ", f14.5, " UTC")') julian
write(*, *)
write(*, *) '星名 赤経α (1950) 赤緯δ 方位角 高度'
write(*, *) '------------------------------------------------------------------------------'
write(*, *) ' h m s ° ′ ″ ° °'
write(*, dfmt) SNM, int(rh), int(rm), rs, sg, int(dh), int(dm), ds, al, hi
write(*, *)
stop
end program pcc_ec2hcコンパイル
gfortran pcc_lib.f90 pcc_ec2hc.f90実行
./a.out実行結果
DATE AND TIME(JST) ? 19780620,223217
1978年 6月 20日 22時 32分 17秒 (JST)
Julian Date = 2443680.06409 UTC
星名 赤経α (1950) 赤緯δ 方位角 高度
------------------------------------------------------------------------------
h m s ° ′ ″ ° °
α CMa Sirius 6 42 56.714 -16 38 46.36 117.999123 -57.459200求められた高度と方位角の値は、「天体の位置計算 増補版」245頁の結果と比較すると、微妙に違ってきています。これは演算精度によるものです。同書では、演算精度は基本的に10桁です。Fortranによる演算は17桁となります。たとえば、30頁の緯度の場合、35.788889として計算していますが、Fortranでは、35.788888888888884で計算します。小数点以下15桁目はおかしくなっていますが、とりあえず無視しておきます。この2つの数値を使って計算した場合、答えは違ってきます。たとえば、3600で割り算したとします。
前者は0.009941358055556、後者は0.009941358024691ですので、その後の計算によっては、おのずと答えも微妙に違ってきます。途中結果も同書と同じにしたい場合は、次に述べる方法で、すべての計算を倍精度ではなく、必要桁数にして計算することです。
もし、値を4.641にしたい場合は、次のような処理をサブルーチンAltitudeAzimuthを呼び出す前に入れておきます。こうすることで、秒の値が4.6410000000000000になり、実行すると245頁の値と同じになります。
write(sfmt, '(f6.3)') sdt_second ! 文字列に変換
read(sfmt, *) sdt_second ! 文字列から数値に変換実行結果は以下のとおりになります。
DATE AND TIME(JST) ? 19780620,223217
1978年 6月 20日 22時 32分 17秒 (JST)
Julian Date = 2443680.06409 UTC
星名 赤経α (1950) 赤緯δ 方位角 高度
------------------------------------------------------------------------------
h m s ° ′ ″ ° °
α CMa Sirius 6 42 56.714 -16 38 46.36 117.999125 -57.459201
