#======================================================================
#  R program to calculate and plot Earth’s orbit around the sun, using
#  Newtonian gravity.  Sun is at the origin, x and y coordinates are
#  measured in astronomical units.  Time is measured in sidereal years.
#======================================================================

#-------------------------------------------------------------
# 0. Set initial location, velocity and direction of the Earth here.
#-------------------------------------------------------------
x0=.983291     # Periapsis distance in astronomical units.
y0=0           # Start on x-axis, at periapsis.
phi=pi/2       # Initial direction angle of earth's orbit.
v0=30.3        # Earth's initial velocity at periapsis, km/s.

#-------------------------------------------------------------
# 1. Set number of time increments and end time for trajectory
#   calculation.
#-------------------------------------------------------------
n.int=10000    # Number of time increments per Earth year
               #   to be used. 
t.end=1        # Number of years duration of trajectory.

#-------------------------------------------------------------
# 2. Various constants gleaned from astronomy books or online.
#-------------------------------------------------------------
G=6.67428e-11        # Gravitational constant, m^3/(kg*s^2).
m=6.0e24             # Mass of the earth, kg.
M=m*3.33e5           # Mass of the sun.
km.per.AU=149.6e6            # One astronomical unit, km.
sec.per.year=31558149.540    # Seconds in a sidereal year.

#-------------------------------------------------------------
# 3. Do all calculations in units convenient for plotting.
#-------------------------------------------------------------
me=m/m      # Mass of earth expressed in earth units.
Me=M/m      # Sun's mass in earth units.
G=G*m*sec.per.year^2/((km.per.AU^3)*(1000^3)) 
                                 # New units of G are AU^2/(me*yr^2).
v0=v0*sec.per.year/km.per.AU     # New units of velocity are AU/yr.
 
#-------------------------------------------------------------
# 4. Initialize various quantities.
#-------------------------------------------------------------
v.x0=v0*cos(phi)      # Earth's velocity, x-component, at periapsis.
v.y0=v0*sin(phi)      # Earth's velocity, y-component, at periapsis.
dt=1/n.int            # Duration of one tiny time interval (delta t).
tot.int=round(t.end*n.int);  # Total number of tiny intervals
                             #   (must be an integer).

#--------------------------------------------------------------
# 5.  Pre-allocate vectors which will store all the values.
#   Note:  x is x-position, y is y-position, v.x is x-direction 
#   velocity, v.y is y-direction velocity, t is time.
#--------------------------------------------------------------
x=numeric(tot.int+1)      #  x starts as a vector of ninc+1 zeros.
                          #    The "plus one" is for the initial
                          #    condition.
y=x           # Allocate y the same as x.
v.x=x         # Allocate v.x the same as x.
v.y=x         # Allocate v.y the same as x.
t=x           # Allocate t the same as x;

#--------------------------------------------------------------
# 6.  Insert the initial conditions into the Vectors.
#--------------------------------------------------------------
x[1]=x0          # Initial value of x-position.
y[1]=y0          # Initial value of y-position.
v.x[1]=v.x0      # Initial value of x-direction velocity.
v.y[1]=v.y0      # Initial value of y-direction velocity.
t[1]=0           # Initial value of time.
c=G*Me           # Pre-calculate a constant that appears repeatedly
                 #   in the equations.

#---------------------------------------------------------------
# 7.  Loop to calculate trajectory.
#---------------------------------------------------------------
for (i in 1:tot.int) {  
  dx=v.x[i]*dt     # Change in x.
  dy=v.y[i]*dt     # Change in y.
  dv.x=-c*x[i]/(x[i]^2+y[i]^2)^(3/2)*dt    # Change in x-velocity.
  dv.y=-c*y[i]/(x[i]^2+y[i]^2)^(3/2)*dt    # Change in y-velocity.
  x[i+1]=x[i]+dx     # New value of x.
  y[i+1]=y[i]+dy     # New value of y.
  v.x[i+1]=v.x[i]+dv.x    # New value of x-velocity.
  v.y[i+1]=v.y[i]+dv.y    # New value of y-velocity.
  t[i+1]=t[i]+dt     # New value of time.
}

#---------------------------------------------------------------
# 8.  Plot the trajectory.
#---------------------------------------------------------------
par(pin=c(4,4))      # Equal screen size in x and y directions
plot(x,y,type="l",xlim=c(-1.1,1.1),ylim=c(-1.1,1.1)) # Line plot of y vs x.