! -*- 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
