天体の位置計算 – 共通サブルーチン
「天体の位置計算 増補版」で使用する共通サブルーチンを掲載します。グリニッジ平均恒星時や固有運動といった部分については、その都度紹介していきますが、時分秒から角度、象限、ベッセル年初からの経過年数などは、よく使われるので共通サブルーチンの中に組み入れました。一部のサブルーチンを「マイコン宇宙講座」からFortranに書き換えたものを使用しています。
サブルーチンは、関数として宣言し直しましたので、コンパイル時には、注意が必要です。また、今後は共通サブルーチンに追加していく形になります。
コンパイル
コンパイルすると、pcc_lib.oとpcclib.modファイルが作成されます。コンパイルは新しく追加するたびに行ってください。エラーなどでコンパイル失敗すると、すでにpcclib.modがあると削除されますので、注意してください。削除された場合は、再度コンパイルして作成します。
gfortran -c pcc_lib.f90▶ 共通サブルーチン pcc_lib.f90
module pcclib
implicit none
contains
function Frc(x) result(val)
! 小数部を返す
double precision, intent(in) :: x
double precision :: val
val = x - int(x);
end Function Frc
function Mla(x) result(val)
! 360°の整数倍を返す
double precision, parameter :: PI = 3.141592653589793238462643d0
double precision, parameter :: PI2 = 6.283185307179586476925287d0
double precision, intent(in) :: x
double precision :: val
val = PI2 * (x / 360d0 - int(x / 360d0));
end function Mla
function Mla2(x) result(val)
! 360°の整数倍を返す
double precision, intent(in) :: x
double precision :: val
val = 360d0 * frc(x / 360d0)
if (val < 0) then
val = val + 360d0;
end if
end function Mla2
function Sgn(x) result(flg)
! 正負の符号を返す
double precision, intent(in) :: x
integer :: flg
flg = 1
if (x < 0) then
flg = -1
end if
end function Sgn
function hms2deg(rh, rm, rs) result(deg)
! 時分秒を角度に換算して返す
double precision, intent(in) :: rh, rm, rs
double precision :: deg
deg = 15.0d0 * (abs(rh) + rm / 60d0 +rs / 3600d0)
end function hms2deg
function dms2deg(dh, dm, ds) result(deg)
! 度分秒を角度に換算して返す
double precision, intent(in) :: dh, dm, ds
double precision :: deg
deg = Sgn(dh) * (abs(dh) + dm / 60d0 + ds / 3600d0)
end function dms2deg
function deg2hms(deg) result(ra)
! 角度から時分秒を返す
double precision, intent(in) :: deg
double precision :: degree
double precision :: rh, rm, rs, ra(3)
degree = abs(deg) / 15.0d0
rh = int(degree)
rs = 60.0d0 * (degree - int(degree))
rm = int(rs)
rs = 60.0d0 * (rs-rm)
ra(1) = rh
ra(2) = rm
ra(3) = rs
end function deg2hms
function deg2dms(degree) result(dc)
! 角度から度分秒を返す
double precision, intent(in) :: degree
double precision :: dh, dm, ds, dc(3)
integer :: sg
sg = Sgn(degree)
dh = int(abs(degree))
ds = 60.0d0 * (abs(degree) - int(abs(degree)))
dm = int(ds)
ds = 60.0d0 * (ds-dm)
dc(1) = sg * dh
dc(2) = dm
dc(3) = ds
end function deg2dms
function deg2rad(x) result(radian)
! 角度をラジアンに換算する
double precision, parameter :: PI = 3.141592653589793238462643d0
double precision, parameter :: RAD = 180d0 / PI
double precision, intent(in) :: x
double precision :: radian
radian = x / RAD
end function deg2rad
function rad2deg(x) result(degree)
! ラジアンを角度に換算する
double precision, parameter :: PI = 3.141592653589793238462643d0
double precision, parameter :: RAD = 180d0 / PI
double precision, intent(in) :: x
double precision :: degree
degree = x * RAD
end function rad2deg
function Quadrant(ss, cc) result(tt)
! 象限の判断
!「マイコン宇宙講座」79頁 リスト3-13
double precision, parameter :: PI = 3.141592653589793238462643d0
double precision, parameter :: PI2 = 6.283185307179586476925287d0
double precision, intent(in) :: ss, cc
double precision :: tt
tt = atan(ss / cc)
if (cc < 0) then
tt = tt + PI
else if (ss < 0) then
tt = tt + PI2
end if
end function Quadrant
function ModifiedJulianDay(year, month, day, hour, minute, second) result(mjd)
! 修正ユリウス日
double precision, intent(in) :: year, month, day, hour, minute, second
double precision :: jd, mjd
jd = JulianDay(year, month, day)
jd = jd + (hour / 24d0 + minute / 1440d0 + second / 86400d0)
jd = jd - 0.375d0
mjd = jd - 2400000.5d0
end function ModifiedJulianDay
function JulianDay(year, month, day) result(julian)
! ユリウス日(0 UTC)
!「マイコン天文学Ⅰ」20-22頁
double precision, intent(in) :: year, month, day
double precision :: julian
double precision :: yy, mm, dd, branch
yy = year
mm = month
dd = day
branch = yy + (mm - 1d0) / 12d0 + dd / 365.25d0
if (int(mm) < 3d0) then
mm = mm + 12d0
yy = yy - 1d0
end if
if (branch >= 1582.78d0) then
julian = int(yy * 365.25d0) + int(yy / 400d0) - int(yy / 100d0) + int(30.59d0 * (mm - 2d0)) + dd + 1721088.5d0
else if (branch >= 0d0) then
julian = int(yy * 365.25d0) + int(30.59d0 * (mm - 2d0)) + dd + 1721086.5d0
else if (year < 0d0) then
julian = Sgn(year) * int(abs(yy) * 365.25d0) + int(30.59d0 * (mm - 2d0)) + dd + 1721085.5d0
end if
end function JulianDay
function Julian2date(julian) result(dt)
! ユリウス日から年月日
! 「マイコン宇宙講座」サブルーチン JDATE リスト2-2
integer, dimension(1:12) :: T = (/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
double precision, intent(in) :: julian
double precision :: dt(3), year, month, day
double precision :: jj, jd, r1, r2
integer :: m
jj = julian - 2400000.5d0
year = int(2.7379093d-03 * (julian - 2400000.5d0) + 1858.877d0)
month = 1.0d0
day = 0.0d0
jd = JulianDay(year, month, day)
r2 = jj - (jd - 2400000.5d0)
if (mod(year, 4d0) == 0d0 .and. mod(year, 100d0) /= 0d0 .or. mod(year, 400d0) == 0d0) then
T(2) = 29
end if
r1 = 0d0
m = 1
do while (m < 13)
if (int(r2) - r1 - T(m) <= 0) then
exit
end if
r1 = r1 + T(m)
m = m + 1
end do
month = m
day = r2 - r1
T(2) = 28
if (month == 13d0) then
year = year + 1d0
month = month - 12d0
end if
dt(1) = year
dt(2) = month
dt(3) = day
end function Julian2date
function BesselianYear(julian) result(t)
! 1950のベッセル年初からの経過日数
! 「マイコン宇宙講座」31頁 式（2.3.1）
double precision, parameter :: JDBEPOCH = 2433282.42346d0
double precision, intent(in) :: julian
double precision :: t
t = julian - JDBEPOCH
t = t * (0.00002737909288d0 + 1.25013d-17 * t)
end function BesselianYear
! ここに新しい関数（サブルーチン）を記述する
! end module pcclibより後ろに記述しないこと！
end module pcclib
