!
!satellite_orbital.f90 
!

!system of unit
!length 1AU (1AU = 1.495978707*10^{11} m)
!mass Msun = 1.0
!time 1yr

!Caltesian Coordinate
!The origin is sun.
!
!
!
!

Module constants_set
  implicit none

  real(kind=8),parameter :: pi=acos(-1.0d0)
  real(kind=8),parameter :: GMsun = 4.0d0*pi*pi


  !Mass & Orbital Elements
  !Earth
  real(kind=8), parameter:: M_E = 3.0D0*(10.0D0**(-6.0D0))  !Mass
  

  !Mars
  real(kind=8), parameter:: M_Ma = 3.0D0*(10.0D0**(-7.0D0))  !Mass
  real(kind=8), parameter:: A_Ma_ini = 1.52368D0  !semi-major axis
  real(kind=8), parameter:: e_Ma_ini = 0.09341233D0  !eccentricity
  real(kind=8), parameter:: i_Ma_ini = 1.8497D0*pi/180.0D0  !inclination
  real(kind=8), parameter:: omega_p_Ma_ini = 336.2075D0*pi/180.0D0  !pericentre
  real(kind=8), parameter:: Omega_a_Ma_ini = 49.6198D0*pi/180.0D0  !asceding node


end Module constants_set






program main
  use constants_set
  implicit none

  integer(kind=4),parameter:: neq = 18  !number of equations, (x, y, z, vx, vy, vz) * 3
  integer(kind=4),parameter:: nmax = 10000  !number of step

  real(kind=8),parameter:: eps = 1.0d-10

  external ex2
  integer(kind=4) ierr,i,n,iseed
  real(kind=8),dimension(neq,12):: work
  real(kind=8),dimension(neq):: y0,yn !position & velocity
  real(kind=8) x0,xn,dx !time

  integer(4) counter, n_acc

  real(8) a_stllt, e_stllt, cosf_stllt

  counter = 0
  n_acc = 0
  
  n = 0
  x0 = 0.0d0				!x = time
  dx = 0.001d0				!timestep
  !Timesteps need enough short to describe satellite. 
  
  
  !initial_condition
  !Caltesian coordinate
  !Earth (simplified)
  !position
  y0(4) = 1.0D0
  y0(5) = 0
  y0(6) = 0
  !velocity
  y0(1) = 0
  y0(2) = sqrt(GMsun)
  y0(3) = 0

  !2nd body
  !x
  y0(10) = A_Ma_ini* (1.0D0 - e_Ma_ini)* ( cos(omega_p_Ma_ini)* cos(Omega_a_Ma_ini) &
       &- sin(omega_p_Ma_ini)* sin(Omega_a_Ma_ini)* cos(i_Ma_ini)  )
  y0(11) = A_Ma_ini* (1.0D0 - e_Ma_ini)* ( cos(omega_p_Ma_ini)* sin(Omega_a_Ma_ini) &
       &+ sin(omega_p_Ma_ini)* cos(Omega_a_Ma_ini)* cos(i_Ma_ini)  )
  y0(12) = A_Ma_ini* (1.0D0 - e_Ma_ini)* ( sin(omega_p_Ma_ini)*  sin(i_Ma_ini)  )
  !v
  y0(7) = sqrt(GMsun)* (1.0D0 + e_Ma_ini)* (-cos(Omega_a_Ma_ini)*sin(omega_p_Ma_ini) &
       &- sin(Omega_a_Ma_ini)* cos(omega_p_Ma_ini)* cos(i_Ma_ini)  ) /sqrt(A_Ma_ini* (1.0D0 - e_Ma_ini* e_Ma_ini) )
  y0(8) = sqrt(GMsun)* (1.0D0 + e_Ma_ini)* (-sin(Omega_a_Ma_ini)*sin(omega_p_Ma_ini) &
       &+ cos(Omega_a_Ma_ini)* cos(omega_p_Ma_ini)* cos(i_Ma_ini)  ) /sqrt(A_Ma_ini* (1.0D0 - e_Ma_ini* e_Ma_ini) )
  y0(9) = sqrt(GMsun)* (1.0D0 + e_Ma_ini)* ( cos(omega_p_Ma_ini)* sin(i_Ma_ini)  )/sqrt(A_Ma_ini* (1.0D0 - e_Ma_ini* e_Ma_ini) )

  !satellite
  !Initiallly on 
  !x
  y0(16) = y0(4) + 0.0025D0
  y0(17) = y0(5)
  y0(18) = y0(6)
  !v
  y0(13) = 0
  y0(14) = sqrt(GMsun/1.0D0) + sqrt((GMsun*M_E)/0.0025D0)
  y0(15) = 0
  
  
  a_stllt = 1.0D0/( (2.0D0/ sqrt( y0(16)*y0(16) + y0(17)*y0(17) &
       &+ y0(18)*y0(18)  ) )  - (( y0(13)*y0(13) + y0(14)*y0(14) + y0(15)*y0(15) )/GMsun ) )  
  e_stllt = sqrt(1.0D0 - ( y0(16)*y0(14)-y0(17)*y0(13) )**2.0D0/(GMsun*a_stllt) )
  
  open(10,file='orbit.dat')
  write(10,'(10f20.10)') x0,  y0(4), y0(5), y0(6), y0(10), y0(11), y0(12)&
       &, y0(16), y0(17), y0(18)

  
  
  !Calculation
  do n = 1,nmax
     
     
     if(e_stllt==0) then
        cosf_stllt = 0
     else        
        cosf_stllt = ( a_stllt* (1.0D0 - e_stllt* e_stllt)/sqrt(y0(16)*y0(16) + y0(17)*y0(17) &
             &+ y0(18)*y0(18) ) - 1.0D0)/e_stllt
     end if

     !Kick
     if(counter < n_acc ) then
        if(cosf_stllt<1.000001D0 .and. cosf_stllt>0.999999D0) then
           y0(14) = y0(14) + 0.010D0*sqrt(GMsun/a_stllt)  !Kick
           counter = counter+1
        end if
     end if
     
        xn = x0 + dx
     
     call runge(neq,ex2,x0,xn,y0,yn,eps,ierr,work)	!runge.f90
     x0 = xn  !time evolution
     
     do i = 1,neq
        y0(i) = yn(i)  !update y0
     end do
     
   a_stllt = 1.0D0/( (2.0D0/ sqrt( y0(16)*y0(16) + y0(17)*y0(17) &
       &+ y0(18)*y0(18)  ) )  - (( y0(13)*y0(13) + y0(14)*y0(14) + y0(15)*y0(15) )/GMsun ) )  
  e_stllt = sqrt(1.0D0 - ( y0(16)*y0(14)-y0(17)*y0(13) )**2.0D0/(GMsun*a_stllt) )
    
     if( mod(n, 10)==0 ) then
        
        write(10,'(10f20.10)') x0,  y0(4), y0(5), y0(6), y0(10), y0(11), y0(12)&
             &, y0(16), y0(17), y0(18)
     end if
  end do
  
  close(10)  
end program main



real(kind=8) function ex2(i,x,y)
  use constants_set
  implicit none
  
  integer(kind=4),parameter:: neq=18
  
  integer(kind=4) i
  real(kind=8) x
  real(kind=8),dimension(neq):: y
  real(kind=8) indirect_x, indirect_y, indirect_z
  real(kind=8) r_SE, r_S2, r_Ss
  real(kind=8) r_E2, r_2s, r_Es 
  real(kind=8) M_2

  M_2 = M_Ma  !

  
  r_SE = sqrt( y(4)*y(4) +  y(5)*y(5) +  y(6)*y(6) )
  r_S2 = sqrt( y(10)*y(10) +  y(11)*y(11) +  y(12)*y(12) )
  r_Ss = sqrt( y(16)*y(16) +  y(17)*y(17) +  y(18)*y(18) )
  
  r_E2 = sqrt( (y(4)-y(10))*(y(4)-y(10)) + (y(5)-y(11))*(y(5)-y(11)) + (y(6)-y(12))*(y(6)-y(12)) )
  r_2s = sqrt( (y(10)-y(16))*(y(10)-y(16)) + (y(11)-y(17))*(y(11)-y(17)) + (y(12)-y(18))*(y(12)-y(18)) )
  r_Es = sqrt( (y(4)-y(16))*(y(4)-y(16)) + (y(5)-y(17))*(y(5)-y(17)) + (y(6)-y(18))*(y(6)-y(18)) )
  
  
  
  indirect_x = M_E*y(4)/(r_SE**(3.0D0)) + M_2*y(10)/(r_S2**(3.0D0))  
  indirect_y = M_E*y(5)/(r_SE**(3.0D0)) + M_2*y(11)/(r_S2**(3.0D0))  
  indirect_z = M_E*y(6)/(r_SE**(3.0D0)) + M_2*y(12)/(r_S2**(3.0D0))  
  
  
  
  if (i==1) then
     !Earth
     !solar gravity + gravity of 2nd body + indirect tetrm
     ex2= -GMsun* ( y(4)/(r_SE**(3.0D0)) + M_2*( y(4)-y(10) )/(r_E2**(3.0D0)) + indirect_x )  !dvx/dt = 
  else if (i==2) then
     ex2= -GMsun* ( y(5)/(r_SE**(3.0D0)) + M_2*( y(5)-y(11) )/(r_E2**(3.0D0)) + indirect_y )		!dvy/dt = 
  else if (i==3) then
     ex2= -GMsun* ( y(6)/(r_SE**(3.0D0)) + M_2*( y(6)-y(12) )/(r_E2**(3.0D0)) + indirect_z )	!dvz/dt =
     
  else if (i==4) then
     ex2= y(1)  !dx/dt
  else if (i==5) then
     ex2= y(2)  !dy/dt =
  else if (i==6) then
     ex2= y(3)  !dz/dt = 
     
     !2
  else if (i==7) then
     ex2= -GMsun* ( y(10)/(r_S2**(3.0D0)) + M_E*( y(10)-y(4) )/(r_E2**(3.0D0)) + indirect_x )
  else if (i==8) then
     ex2= -GMsun* ( y(11)/(r_S2**(3.0D0)) + M_E*( y(11)-y(5) )/(r_E2**(3.0D0)) + indirect_y )
  else if (i==9) then
     ex2= -GMsun* ( y(12)/(r_S2**(3.0D0)) + M_E*( y(12)-y(6) )/(r_E2**(3.0D0)) + indirect_z )
     
  else if (i==10) then
     ex2= y(7)
  else if (i==11) then
     ex2= y(8)
  else if (i==12) then
     ex2= y(9)
     
     !satellite
  else if (i==13) then
     ex2= -GMsun* ( y(16)/(r_Ss**(3.0D0)) + M_E*( y(16)-y(4) )/(r_Es**3.0D0) + M_2*( y(16)-y(10) )/(r_2s**(3.0D0)) + indirect_x )
  else if (i==14) then
     ex2= -GMsun* ( y(17)/(r_Ss**(3.0D0)) + M_E*( y(17)-y(5) )/(r_Es**3.0D0) + M_2*( y(17)-y(11) )/(r_2s**(3.0D0)) + indirect_y )
  else if (i==15) then
     ex2= -GMsun* ( y(18)/(r_Ss**(3.0D0)) + M_E*( y(18)-y(6) )/(r_Es**3.0D0) + M_2*( y(18)-y(12) )/(r_2s**(3.0D0)) + indirect_z )
     
  else if (i==16) then
     ex2=y(13)
  else if (i==17) then
     ex2=y(14)
  else
     ex2=y(15)
  end if
  
  return
end function ex2

