1
$\begingroup$

If Earth and Sun were in a isolated system, will the Earth's motion around Sun will be similar? What will be Sun's and Earth's velocity when Earth is at its aphelion?

Please note that it's not a homework question. I am creating a simulation of two-body system and I just wanted to apply it to an idealised Sun-Earth system. I am having trouble choosing the initial conditions. The simulation runs but I don't get desired closed orbit.

module var
  implicit none
  real, parameter :: mScale = 1.988544e30   ! Mass of Sun
  real, parameter :: lScale = 1.49597870700e11   ! 1 AU
  real, parameter :: tScale = 86400    ! Mean Solar day
  real, parameter :: G=0.0002959122083
  real, parameter :: mEarth = 3.0024584e-6
  real, parameter :: mSun = 1.0
  integer, parameter :: n=2
  real, dimension(3,n) :: xyz,vel,acc
  real, parameter, dimension(n) :: m=[mSun, mEarth]
  real, parameter :: tot_t=100.0
  real, parameter :: dt=0.001
  real, parameter, dimension(n) :: rad=[0.1, 0.1]
  real :: ke,pe
end module var
module initial
  use var, only: xyz,vel,acc,m,orig,angmom,G,tScale,lScale,mEarth,mSun
  implicit none
  private
  public :: init
contains
  subroutine init()
    real :: temp
    xyz(:,1)=[(0.025e11)/lScale,0.0,0.0]
    xyz(:,2)=[(-1.471e11)/lScale,0.0,0.0]
    vel(:,1)=[0.0,(mEarth/mSun)*30300*tScale/lScale,0.0]
    vel(:,2)=[0.0,-30300*tScale/lScale,0.0]
    acc(:,1)=[-G*m(2)/(1.496e11/lScale)**2,0.0,0.0]
    acc(:,2)=[G*m(1)/(1.496e11/lScale)**2,0.0,0.0]
  end subroutine init
end module initial
module update
  use var, only: xyz,vel,acc,m,orig,angmom,G,n,dt
  implicit none
  private
  integer :: i,j
  real, dimension(3,n) :: tempacc
  real, dimension(3) :: dist
  public :: posUpd,velUpd,accUpd
contains
  subroutine posUpd()
    implicit none
    integer :: i,k
    do i=1,n
      do k=1,3
        xyz(k,i)=xyz(k,i) + vel(k,i)*dt + (acc(k,i)*(dt**2))/2
      enddo
    enddo
  end subroutine posUpd
  subroutine accUpd()
    implicit none
    integer :: i,j,k
    real :: r,temp
    do i=1,n
      tempacc(:,i)=acc(:,i)
      acc(:,i)=0.0
    enddo
    do i=1,n
      do j=i+1,n
        dist=xyz(:,i)-xyz(:,j)
        r=sqrt(sum(dist**2))
        temp=(G*m(i)*m(j))/(r**3)
        do k=1,3
          acc(k,i)=acc(k,i) - temp*dist(k)
          acc(k,j)=acc(k,j) + temp*dist(k)
        enddo
      enddo
    enddo
  end subroutine accUpd
  subroutine velUpd()
    implicit none
    integer :: i,k
    do i=1,n
      do k=1,3
        vel(k,i)=vel(k,i) + 0.5*(acc(k,i)+tempacc(k,i))*dt
      enddo
    enddo
  end subroutine velUpd
end module update
module energy
  use var, only: xyz,vel,m,ke,pe,n,G
  implicit none
  private
  real, dimension(3) ::dist
  real :: r
  integer :: i,j
  public :: kinetic,potential
contains
  subroutine kinetic()
    implicit none
    do i=1,n
      ke=0.5*m(i)*sum(vel(:,i)*vel(:,i))
    enddo
  end subroutine kinetic
  subroutine potential()
    real :: temp
    pe=0.0
    do i=1,n-1
      do j=i+1,n
        dist=xyz(:,i)-xyz(:,j)
        r=sqrt(sum(dist**2))
        pe= pe + (G*m(i)*m(j))/r       !r is relative distance.
      enddo
    enddo
  end subroutine potential
end module energy
program planet
  use var, only : xyz,vel,m,n,rad,ke,pe,tot_t,dt
  use energy, only : kinetic,potential
  use update, only : posUpd,velUpd,accUpd
  use initial, only: init
  implicit none
  integer :: i,t,iter
  character(len=30) :: fmt
  open(100,file="xyz.dat",status="replace")
  open(200,file="vel.dat",status="replace")
  open(300,file="acc.dat",status="replace")
  open(400,file="energy.dat",status="replace")
  open(500,file="params.dat",status="replace")
  call init()
  print*, "The Simulation is running."
  iter=int(tot_t/dt)
  do t=1,iter
    if(t==1) then
      call kinetic()
      call potential()
      write(400,*) t,ke,pe,ke+pe
    endif
    if(mod(t,100)==0) then
      call kinetic()
      call potential()
      write(400,*) t,ke,pe,ke+pe
    endif
    do i=1,n
      write(100,*) xyz(:,i)
    enddo
    call posUpd()
    call accUpd()
    call velUpd()
  enddo
  write(500,*) n
  write(500,*) rad
  call kinetic()
  print*, "The Kinetic Energy is ", ke
  call potential()
  print*, "The Potential Energy is ", pe
  call execute_command_line("start /b python show.py")
end program planet
$\endgroup$
  • $\begingroup$ The average speed of the earth is easy to figure out from knowing what the length of a year is. If you're using realistic numbers and units, the motion of the sun will be very slow, and possibly chasing numerical error depending on how careful you're being with your doubles and floats. $\endgroup$ – Jerry Schirmer Jun 30 '15 at 16:12
  • $\begingroup$ en.wikipedia.org/wiki/Orbital_speed#Precise_orbital_speed $\endgroup$ – lemon Jun 30 '15 at 16:42
  • $\begingroup$ Note in particular that in the real world the Sun responds much more strongly to the influence of Jupiter. $\endgroup$ – Emilio Pisanty Jun 30 '15 at 17:32
  • 1
    $\begingroup$ It is far more likely that your simulation is running into problems ("I don't get the desired closed orbit") because of your method of integration, than because you are not accounting for the finite velocity of the Sun. I have certainly run into that problem myself. You need some higher order method (e.g. leapfrog or Runge-Kutta type integration) - a simple Newton integration will not lead to a closed / stable orbit regardless of time step. $\endgroup$ – Floris Jun 30 '15 at 19:40
  • 1
    $\begingroup$ @Floris I found the mistake. While updating the acceleration, I was not dividing by the object's mass and thus using the force in place of acceleration. It works now. Thanks anyways. $\endgroup$ – Yogesh Yadav Jul 1 '15 at 5:43
1
$\begingroup$

In a reference frame where the center of mass of the Earth-Sun system is at rest, we will have $$ m_\text{Sun} \vec{v}_\text{Sun} + m_\text{Earth} \vec{v}_\text{Earth} = 0 \quad \Rightarrow \quad \vec{v}_\text{Sun} = - \frac{m_\text{Earth}}{m_\text{Sun}} \vec{v}_\text{Earth}. $$ In particular, if this is true at some initial time, then it's true at all later times as well. So you would probably want to set up your initial conditions so that this was true at $t = 0$, and then integrate from there. Note that (as pointed out by Jerry Schirmer in the comments) this velocity will be incredibly slow compared to the velocity of the Earth, and you may run into numerical errors.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.