天体の位置計算 – 地方恒星時の計算
「天球座標-恒星時」25頁
恒星時には真春分点を基準にした視恒星時、平均春分点を基準にした平均恒星時、本初子午線を基準としたグリニッジ恒星時に観測地点の経度を加えた地方恒星時があります。それぞれに計算方法がありますが、ここでは、地方恒星時を求めてみることにします。
計算例では、必要な数値をすべて時分秒に換算してから計算しています。また、グリニッジ視恒星時\(\theta_{0}\)が必要なので天体位置表から計算日における世界時0時のグリニッジ視恒星時：17h11m58s.714を求めています。現在は天体位置表や天文観測年表は絶版になっていますので、国立天文台天文情報センター暦計算室にある暦象年表や理科年表や天文年鑑を利用するとよいでしょう。おすすめは国立天文台です。
地方恒星時は以下の式で求めることができます。
\[\theta = \theta_{0} + \lambda + t_{1} + 補正値 \] \[\theta: 地方恒星時 \theta_{0}: グリニッジ視恒星時 \lambda: 観測値の経度 t_{1}: 経過時間\]Fortranによるコード
本来ならグリニッジ視恒星時を計算する必要があります。今回は、計算例や問題に示されているグリニッジ視恒星時を使用して計算することにします。グリニッジ視恒星時については、章動のところで求め、後にメインルーチンに反映させることにします。
▶メインルーチン pcc_lsdt.f90
program pcc_lsdt
use pcclib
implicit none
! 地点
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 = "Tokyo Ngasawa" ! 地点名
! 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, t, julian
double precision :: year, month, day, hour, minute, second
double precision :: gsdt, lsdt, long, cv
double precision :: sdt_hour, sdt_minute, sdt_second
double precision :: sdt(3), dtime(3)
character(255) :: dfmt
! 日付の入力
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
! ユリウス日から年月日を求める
dtime = Julian2date(julian)
year = dtime(1)
month = dtime(2)
day = dtime(3)
! グリニッジ視恒星時を時の小数点に換算
gsdt = 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 = gsdt + 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
! 出力
dfmt = '(" ", i4, " ", i2, " ", f8.5, " ", f14.5, " ", i2, " ", i2, " ", f6.3)'
write(*, *)
write(*, *) 'Date/ Julian Date Local Sidereal Time'
write(*, *) '--------------------------------------------------------------'
write(*, *) ' h m s'
write(*, dfmt) int(year), int(month), day, julian, int(sdt_hour), int(sdt_minute), sdt_second
write(*, *)
stop
end program pcc_lsdtコンパイル
gfortran pcc_lib.f90 pcc_lsdt.f90実行
./a.out実行結果
DATE AND TIME(JST)? 19780620,223217
Date/ Julian Date Local Sidereal Time
--------------------------------------------------------------
h m s
1978 6 20.56409 2443680.06409 16 44 4.641グリニッジ地方恒星時から地方恒星時を求める
「マイコン天文学Ⅰ」49頁に掲載されている式2-4-18を使っても求めることができます。
▶ サブルーチン localSDT
function localSDT(julian) result(lsdt)
! グリニッジ地方恒星時
! 「マイコン天文学Ⅰ」49頁 式2-4-18
double precision, parameter :: PI = 3.141592653589793238462643d0
double precision, parameter :: PI2 = 6.283185307179586476925287d0
double precision, intent(in) :: julian
double precision :: lsdt
lsdt = PI2 * Frc(0.7769194d0 + 1.002737909265d0 * (julian - 2415020.0d0))
end function LocalSDT戻り値はラジアンですので、角度に換算し、経度を加算すれば、その地点の地方恒星時を求めることができます。こちらの方法では、グリニッジ視恒星時を使用しません。
▶ メインルーチン pcc_lsdt_2.f90
program pcc_lsdt_2
use pcclib
implicit none
! 地点
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 = "Tokyo Ngasawa" ! 地点名
! 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, t, julian
double precision :: year, month, day, hour, minute, second
double precision :: lsdt, lg
double precision :: sdt_hour, sdt_minute, sdt_second
double precision :: sdt(3), dtime(3)
character(255) :: dfmt
! 日付の入力
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
dtime = Julian2date(julian)
year = dtime(1)
month = dtime(2)
day = dtime(3)
! 地方恒星時の計算
lg = hms2deg(OBS_LG(1), OBS_LG(2), OBS_LG(3))
lsdt = rad2deg(localSDT(julian)) + lg
sdt = deg2hms(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
! 出力
dfmt = '(" ", i4, " ", i2, " ", f8.5, " ", f14.5, " ", i2, " ", i2, " ", f6.3)'
write(*, *)
write(*, *) 'Date/ Julian Date Local Sidereal Time'
write(*, *) '--------------------------------------------------------------'
write(*, *) ' h m s'
write(*, dfmt) int(year), int(month), day, julian, int(sdt_hour), int(sdt_minute), sdt_second
write(*, *)
stop
end program pcc_lsdt_2実行結果
DATE AND TIME(JST)? 19780620,223217
Date/ Julian Date Local Sidereal Time
--------------------------------------------------------------
h m s
1978 6 20.56409 2443680.06409 16 44 4.564