!======================================================================
!＜orbit.f90＞
!・ルンゲクッタステップは、サブルーチン"runge"としてrunge.f90に記述されている。
!・コンパイルするときは、runge.f90も一緒にすることを忘れずに。
!======================================================================
program main
  implicit none

  !---- 変数の宣言 ----
  integer(kind=4),parameter:: neq=4 !方程式の数
  integer(kind=4),parameter:: nmax=10000 !ステップ数
  real(kind=8),parameter:: eps=1.0d-10

  external ex2 !外部関数の利用を宣言

  real(kind=8),dimension(neq):: y0,yn,y_r
  real(kind=8) x0,xn,dx
  !y0,ynは慣性系での座標、y_rは回転系での座標
  !y0(1)=x座標,y0(2)=y座標,y0(3)=速度のx成分,y0(4)=速度のy成分 
  !x0=角速度*時間、角速度=1
  !x0=2*piで惑星の１公転周期
  !mu:太陽と惑星の質量比
  real(kind=8),parameter:: mu=0.001d0 !(A)tadpole
  !real(kind=8),parameter:: mu=0.000953875d0 !(B)horseshoe

  real(kind=8),dimension(neq,12):: work
  integer(kind=4) ierr,i,n,iseed
  !---- ここまで宣言 ----
  
  n=0
  x0=0.0d0 !ステップの初期値
  dx=1.0d-2 !ステップの刻み幅

  !初期値：(A)(B)のうち一方
  !(A)tadpole_orbit
  y_r(1)=0.5d0-mu+0.0065d0
  y_r(2)=0.5d0*sqrt(3.0d0)+0.0065d0
  y_r(3)=0.0d0
  y_r(4)=0.0d0

  !(B)horseshoe_orbit
  !y_r(1)=-0.97668d0
  !y_r(2)=0.0d0
  !y_r(3)=0.0d0
  !y_r(4)=-0.06118d0

  !慣性系へ変換
  y0(1) = cos(x0)*y_r(1) - sin(x0)*y_r(2)
  y0(2) = sin(x0)*y_r(1) + cos(x0)*y_r(2)
  y0(3) = cos(x0)*(y_r(3)-y_r(2)) - sin(x0)*(y_r(4)+y_r(1))
  y0(4) = sin(x0)*(y_r(3)-y_r(2)) + cos(x0)*(y_r(4)+y_r(1))

  open(10,file='tadpole.dat')
  !データの書き出し
  write(10,'(i5,5f20.10)') n,x0,y0(1),y0(2),y_r(1),y_r(2)
       
  do n=1,nmax
     !dx後の値を計算：慣性系
     !rungeを呼び出すことで計算された値がynに入る
     xn=x0+dx
     call runge(neq,ex2,x0,xn,y0,yn,eps,ierr,work)

     !x0,y0を更新
     x0=xn
     do i=1,neq
        y0(i)=yn(i)
     end do

     !回転系へ変換
     y_r(1) = cos(x0)*y0(1) + sin(x0)*y0(2)
     y_r(2) = -sin(x0)*y0(1) + cos(x0)*y0(2)

     !データの書き出し
     write(10,'(i5,5f20.10)') n,x0,y0(1),y0(2),y_r(1),y_r(2)
  end do
  
  close(10)
end program main


real(kind=8) function ex2(i,x,y)
  implicit none
  
  !変数の宣言
  integer(kind=4),parameter:: neq=4  
  integer(kind=4) i
  real(kind=8) x
  real(kind=8),dimension(neq):: y, ys, yp !y:粒子の座標、ys:太陽の座標、yp:惑星の座標
  real(kind=8) rs, rp !rs:太陽と粒子の距離、rp:惑星と粒子の距離
  !mu:太陽と惑星の質量比
  real(kind=8),parameter:: mu=0.001d0 !(A)tadpole
  !real(kind=8),parameter:: mu=0.000953875d0 !(B)horseshoe

  !準備
  ys(1)=-mu*cos(x) !太陽のx座標
  ys(2)=-mu*sin(x) !太陽のy座標
  yp(1)=(1.0d0-mu)*cos(x) !惑星のx座標
  yp(2)=(1.0d0-mu)*sin(x) !惑星のy座標
  rs=sqrt((y(1)-ys(1))**2+(y(2)-ys(2))**2) !太陽との距離
  rp=sqrt((y(1)-yp(1))**2+(y(2)-yp(2))**2) !惑星との距離
  
  !dy/dx=の式の右辺を書く
  if (i==1) then
     ex2=y(3) !dy(1)/dx（＝速度のx成分）
  else if (i==2) then
     ex2=y(4) !dy(2)/dx（＝速度のy成分）     
  else if (i==3) then
     ex2=(1.0d0-mu)*(ys(1)-y(1))/(rs**3) + mu*(yp(1)-y(1))/(rp**3)
     ! dy(3)/dx（＝加速度のx成分：太陽からの重力と惑星からの重力のx成分）
  else 
     ex2=(1.0d0-mu)*(ys(2)-y(2))/(rs**3) + mu*(yp(2)-y(2))/(rp**3)
     ! dy(4)/dx（＝加速度のy成分：太陽からの重力と惑星からの重力のy成分）
  end if
  
  return
end function ex2









