program kepler
  implicit none
  integer,parameter :: kosuu=256
  real,dimension(0:kosuu) :: t
  real,dimension(2,0:kosuu) :: pos0,pos1,vel1,pos2,vel2,pos4,vel4
  real,dimension(2) :: k1p,k1v,k2p,k2v,k3p,k3v,k4p,k4v
  real :: h,epot,ekin,etot,l,e,a,b,tperiod,alpha
  real :: r,dt,theta,theta_old
  integer :: i
  real :: x0,y0,vx0,vy0
  real,parameter :: t_kaishi=0.0,t_shuryo=40.0
  real,parameter :: pi=4.0*atan(1.0)
  !
  dt = (t_shuryo-t_kaishi)/real(kosuu)
  do i = 0,kosuu
     t(i) = t_kaishi+dt*real(i)
  end do
  ! 初期条件 ! ここを変更すると、異なった軌道が得られる
  x0 = 1.0 ! 位置 (x座標)
  y0 = 0.0 !      (y座標)
  vx0 = 0.2/sqrt(sqrt(x0**2+y0**2)) ! 速度 (x方向成分)
  vy0 = 1.2/sqrt(sqrt(x0**2+y0**2)) !      (y方向成分)
  ! 全エネルギーの計算
  epot = -1.0/sqrt(x0**2+y0**2) ! 重力ポテンシャル (太陽は原点にあるとしている)
  ekin = 0.5*(vx0**2+vy0**2)    ! 運動エネルギー
  etot = epot+ekin              ! 全力学的エネルギー
  if (etot >= 0.0) then ! 全エネルギーが0以上だと楕円軌道にならない
     write(*,*)"楕円軌道になりません!! 全エネルギーを小さくしてください!!"
     stop
  endif
  ! 与えた初期条件のもとに得られる楕円軌道を理論的に求める
  h = abs(x0*vy0-y0*vx0)   ! 2×面積速度 (太陽は原点にあるとしている)
  l = h**2                 ! 半直弦
  e = sqrt(1.0+2.0*l*etot) ! 離心率
  a = l/(1.0-e**2)         ! 長半径
  b = sqrt(a*l)            ! 短半径
  tperiod = 2.0*a*b*pi/h   ! 公転周期
  alpha = atan2(y0,x0) &                  ! 楕円軌道の長軸と
       -acos((l/sqrt(x0**2+y0**2)-1.0)/e) ! x軸のなす角度
  ! 計算結果を出力するファイル名を指定する
  open(10,file="kepler.txt")
  write(10,'("# 離心率 = ",f10.8)')e
  write(10,'("# 長半径 = ",f11.8)')a
  write(10,'("# 短半径 = ",f11.8)')b
  write(10,'("# 公転周期 = ",f12.8)')tperiod
  write(10,'("# 長軸とx軸の角度 = ",f12.8)')alpha
  ! 理論的な計算
  pos0(1,0) = x0
  pos0(2,0) = y0
  theta_old = atan2(y0,x0)
  do i = 1,kosuu
     theta = theoretically_updated_theta(theta_old,tperiod,dt,e,alpha)
     r = l/(1.0+e*cos(theta-alpha))
     pos0(1,i) = r*cos(theta)
     pos0(2,i) = r*sin(theta)
     theta_old = theta
  end do
  ! オイラー法による計算
  pos1(1,0) = x0
  pos1(2,0) = y0
  vel1(1,0) = vx0
  vel1(2,0) = vy0
  do i = 0,kosuu-1
     vel1(:,i+1) = vel1(:,i)+dt*dvdt(pos1(1,i),pos1(2,i))
     pos1(:,i+1) = pos1(:,i)+dt*dxdt(vel1(1,i),vel1(2,i))
  end do
  ! ホイン法による計算
  pos2(1,0) = x0
  pos2(2,0) = y0
  vel2(1,0) = vx0
  vel2(2,0) = vy0
  do i = 0,kosuu-1
     k1v(:) = dt*dvdt(pos2(1,i),pos2(2,i))
     k1p(:) = dt*dxdt(vel2(1,i),vel2(2,i))
     !
     k2v(:) = dt*dvdt(pos2(1,i)+k1p(1),pos2(2,i)+k1p(2))
     k2p(:) = dt*dxdt(vel2(1,i)+k1v(1),vel2(2,i)+k1v(2))
     !
     vel2(:,i+1) = vel2(:,i)+0.5*(k1v(:)+k2v(:))
     pos2(:,i+1) = pos2(:,i)+0.5*(k1p(:)+k2p(:))
  end do
  ! (4次の) ルンゲ-クッタ法による計算
  pos4(1,0) = x0
  pos4(2,0) = y0
  vel4(1,0) = vx0
  vel4(2,0) = vy0
  do i = 0,kosuu-1
     k1v(:) = dt*dvdt(pos4(1,i),pos4(2,i))
     k1p(:) = dt*dxdt(vel4(1,i),vel4(2,i))
     !
     k2v(:) = dt*dvdt(pos4(1,i)+0.5*k1p(1),pos4(2,i)+0.5*k1p(2))
     k2p(:) = dt*dxdt(vel4(1,i)+0.5*k1v(1),vel4(2,i)+0.5*k1v(2))
     !
     k3v(:) = dt*dvdt(pos4(1,i)+0.5*k2p(1),pos4(2,i)+0.5*k2p(2))
     k3p(:) = dt*dxdt(vel4(1,i)+0.5*k2v(1),vel4(2,i)+0.5*k2v(2))
     !
     k4v(:) = dt*dvdt(pos4(1,i)+k3p(1),pos4(2,i)+k3p(2))
     k4p(:) = dt*dxdt(vel4(1,i)+k3v(1),vel4(2,i)+k3v(2))
     !
     vel4(:,i+1) = vel4(:,i)+(k1v(:)+2.0*k2v(:)+2.0*k3v(:)+k4v(:))/6.0
     pos4(:,i+1) = pos4(:,i)+(k1p(:)+2.0*k2p(:)+2.0*k3p(:)+k4p(:))/6.0
  end do
  ! 結果の表示
  do i = 0,kosuu
     write(10,'(9f16.8)')t(i),pos0(1:2,i),pos1(1:2,i),pos2(1:2,i),pos4(1:2,i)
  end do
contains
  function dvdt(x,y)
    real,intent(in) :: x,y
    real,dimension(2) :: dvdt
    real :: r3i
    r3i = 1.0/sqrt(x**2+y**2)**3
    dvdt(1) = -x*r3i
    dvdt(2) = -y*r3i
  end function dvdt
  function dxdt(vx,vy)
    real,intent(in) :: vx,vy
    real,dimension(2) :: dxdt
    dxdt(1) = vx
    dxdt(2) = vy
  end function dxdt
  function theoretically_updated_theta(theta_old,tperiod,dt,e,alpha)
    real,intent(in) :: theta_old,tperiod,dt,e,alpha
    real :: theoretically_updated_theta
    real :: beta,tan_work
    real :: dfdb,f,uhen,deltabeta
    real,parameter :: very_small=1.e-6
    integer :: n
    integer,parameter :: nmax=999 ! ニュートン法の最大反復回数
    tan_work = sqrt((1.0-e)/(1.0+e))*tan(0.5*(theta_old-alpha))
    beta = atan(tan_work)
    uhen = beta-0.5*e*sin(2.0*beta)+pi*dt/tperiod
    ! ここからニュートン法
    newton: do n = 1,nmax
       f = beta-0.5*e*sin(2.0*beta)-uhen
       dfdb = 1.0-e*cos(2.0*beta)
       deltabeta = f/dfdb
       if (abs(deltabeta) <= very_small) exit newton
       beta = beta-deltabeta
    end do newton
    tan_work = sqrt((1.0+e)/(1.0-e))*tan(beta)
    theoretically_updated_theta = 2.0*atan(tan_work)+alpha
  end function theoretically_updated_theta
end program kepler