! -*- coding: utf-8 -*- ! ! 軌道要素計算 ! Berger (1978) の方法と係数を用いる ! 係数は外部ファイルから読み込む ! ! 2013/09/20 Satoshi NODA 作成 ! program calc_orbit implicit none ! temporary integer, parameter:: DP = 8 integer, parameter:: TOKEN = 80 real(DP), parameter:: PI = 3.141592653589793 integer, parameter:: N_DEVICE = 1 ! device number for fileopen (5 is stdin) character(TOKEN), parameter:: DATAFILE = './CoefsForOrbElem_Earth_B1978.dat' ! filename (necessary if N_DEVICE is not 5) integer, save:: N_Obl ! Number of terms for obliquity integer, save:: N_EPi ! Number of terms for eccentricity * sin/cos (Periherion) integer, save:: N_GPre ! Number of terms for general precession real(DP), allocatable:: obl_ampl (:) ! Amplitude real(DP), allocatable:: obl_mrate (:) ! Mean rate real(DP), allocatable:: obl_phase (:) ! Phase real(DP), allocatable:: epi_ampl (:) ! Amplitude real(DP), allocatable:: epi_mrate (:) ! Mean rate real(DP), allocatable:: epi_phase (:) ! Phase real(DP), allocatable:: gpre_ampl (:) ! Amplitude real(DP), allocatable:: gpre_mrate (:) ! Mean rate real(DP), allocatable:: gpre_phase (:) ! Phase real(DP):: obl0 real(DP):: gpre_mrate0 real(DP):: gpre0 real(DP):: Obliquity real(DP):: Eccentricity real(DP):: LonOfPeri real(DP):: yearAD real(DP):: y0, y1, interval real(DP):: w integer:: iy call OrbitElemInit ! 計算する年の始点, 終点, 間隔. 西暦で指定. y0 = -1000000.0 y1 = 100000.0 interval = 100.0 do iy = 0, int((y1-y0)/interval) yearAD = y0 + iy*interval call OrbitElem( & & yearAD, & ! (in) year (A.D.) & Obliquity, & ! (out) & Eccentricity, & ! (out) & LonOfPeri & ! (out) & ) ! output ! radian は degree に変換して出力 w = 180.0 / PI print *, yearAD, w*Obliquity, Eccentricity, w*LonOfPeri end do call OrbitElemFinalize contains subroutine OrbitElem( & & yearAD, & ! (in) year (A.D.) & Obliquity, & ! (out) & Eccentricity, & ! (out) & LonOfPeri & ! (out) & ) implicit none real(DP), intent(in):: yearAD real(DP), intent(out):: Obliquity real(DP), intent(out):: Eccentricity real(DP), intent(out):: LonOfPeri ! Work variable ! real(DP):: y real(DP):: obl real(DP):: episin real(DP):: epicos real(DP):: gpre real(DP):: FixedEquinox integer:: i ! 計算に使用する年. 西暦 1950 年起点. y = yearAD - 1950.0 Obliquity = 0.0_DP episin = 0.0_DP epicos = 0.0_DP gpre = 0.0_DP ! 桁落ちを防ぐため, 振幅の小さい項から足す do i = N_Obl, 1, -1 Obliquity = Obliquity + obl_ampl(i) * cos(obl_mrate(i)*y + obl_phase(i)) end do Obliquity = Obliquity + obl0 do i = N_EPi, 1, -1 episin = episin + epi_ampl(i) * sin(epi_mrate(i)*y + epi_phase(i)) epicos = epicos + epi_ampl(i) * cos(epi_mrate(i)*y + epi_phase(i)) end do do i = N_GPre, 1, -1 gpre = gpre + gpre_ampl(i) * sin(gpre_mrate(i)*y + gpre_phase(i)) end do gpre = gpre + gpre0 + gpre_mrate0 * y Eccentricity = sqrt(episin**2 + epicos**2) FixedEquinox = acos(epicos / Eccentricity) if (asin(episin) < 0.0) then FixedEquinox = - FixedEquinox end if LonOfPeri = mod((fixedEquinox + gpre), 2.0*PI) if (LonOfPeri < 0.0) then LonOfPeri = 2.0*PI + lonOfPeri end if end subroutine OrbitElem subroutine OrbitElemInit implicit none integer:: i if (N_DEVICE .ne. 5) then open(N_DEVICE, file=trim(DATAFILE), status="old") end if read(N_DEVICE,*) N_Obl, N_EPi, N_GPre allocate( obl_ampl (1:N_Obl) ) allocate( obl_mrate (1:N_Obl) ) allocate( obl_phase (1:N_Obl) ) allocate( epi_ampl (1:N_EPi) ) allocate( epi_mrate (1:N_EPi) ) allocate( epi_phase (1:N_EPi) ) allocate( gpre_ampl (1:N_GPre) ) allocate( gpre_mrate (1:N_GPre) ) allocate( gpre_phase (1:N_GPre) ) do i=1,N_Obl read(N_DEVICE,*) obl_ampl(i), obl_mrate(i), obl_phase(i) end do do i=1,N_EPi read(N_DEVICE,*) epi_ampl(i), epi_mrate(i), epi_phase(i) end do do i=1,N_GPre read(N_DEVICE,*) gpre_ampl(i), gpre_mrate(i), gpre_phase(i) end do read(N_DEVICE,*) obl0, gpre_mrate0, gpre0 if (N_DEVICE .ne. 5) then close(N_DEVICE) end if ! 角度の単位を radian に直す obl_ampl = obl_ampl * PI / (180.0 * 3600.0) gpre_ampl = gpre_ampl * PI / (180.0 * 3600.0) obl_mrate = obl_mrate * PI / (180.0 * 3600.0) epi_mrate = epi_mrate * PI / (180.0 * 3600.0) gpre_mrate = gpre_mrate * PI / (180.0 * 3600.0) obl_phase = obl_phase * PI / 180.0 epi_phase = epi_phase * PI / 180.0 gpre_phase = gpre_phase * PI / 180.0 obl0 = obl0 * PI / 180.0 gpre_mrate0 = gpre_mrate0 * PI / (180.0 * 3600.0) gpre0 = gpre0 * PI / 180.0 end subroutine OrbitElemInit subroutine OrbitElemFinalize deallocate( obl_ampl ) deallocate( obl_mrate ) deallocate( obl_phase ) deallocate( epi_ampl ) deallocate( epi_mrate ) deallocate( epi_phase ) deallocate( gpre_ampl ) deallocate( gpre_mrate ) deallocate( gpre_phase ) end subroutine OrbitElemFinalize end program calc_orbit