# What is the velocity of Sun due to Earth alone?

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

• 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. – Jerry Schirmer Jun 30 '15 at 16:12
• en.wikipedia.org/wiki/Orbital_speed#Precise_orbital_speed – lemon Jun 30 '15 at 16:42
• Note in particular that in the real world the Sun responds much more strongly to the influence of Jupiter. – Emilio Pisanty Jun 30 '15 at 17:32
• 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. – Floris Jun 30 '15 at 19:40
• @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. – Yogesh Yadav Jul 1 '15 at 5:43

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.