<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Yarkovsky Effect
This example shows how to add the Yarkovsky effect in a Rebound simulation. There are two versions, which we call the 'Full Version' and the 'Simple Version.' A special parameter called 'ye_flag' is used to switch between the two. The difference between the versions and what situations they're better suited for is discussed in more detail below. 
For more information on this effect, please visit
Step1:  As with all REBOUNDx effects, the parameters must be inputed with the same units as the simulation (in this case it's AU/Msun/yr). We'll use the astropy units module to help avoid errors
Step2:  We then add the Yarkovsky effect and the required parameters for this version. Importantly, we must set 'ye_flag' to 0 to get the Full Version. Physical constants and the stellar luminosity get added to the effect yark
Step3:  Other parameters need to be added to each particle feeling the Yarkovsky effect
Step4:  We integrate this system for 100,000 years and print out the difference between the particle's semi-major axis before and after the simulation.
Step5:  Simple Version
This version of the effect is based off of equations from Veras et al. (2019). Once again, a link to this paper is provided below. This version simplifies the equations by placing constant values in a rotation matrix that in general is time-dependent. It requires fewer parameters than the full version and takes less computational time. However, it is mostly useful only to get a general idea on how much the effect can push bodies inwards or outwards. This version of the effect is better for simulating large groups of asteroids or trying to see general trends in the behavior of a body. 
Link to paper
Step6:  We then add the Yarkovsky effect from Reboundx and the necesary parameters for this version. This time, we must make sure that 'ye_flag' is set to 1 or -1 to get the Simple Version of the effect. Setting it to 1 will push the asteroid outwards, while setting it to -1 will push it inwards. We'll push out our original asteroid and push in our new one. We use the same physical properties as in the example above
Step7:  Now we run the sim for 100,000 years and print out the results for both asteroids. Note the difference in simulation times between the versions. Even with an extra particle, the simple version was faster than the full version.
<ASSISTANT_TASK:>
Python Code:
import rebound
import reboundx
import numpy as np
import astropy.units as u
import astropy.constants as constants
import matplotlib.pyplot as plt
%matplotlib inline
#Simulation begins here
sim = rebound.Simulation()
sim.units = ('yr', 'AU', 'Msun') #changes simulation and G to units of solar masses, years, and AU  
sim.integrator = "whfast" #integrator for sim
sim.dt = .05 #timestep for sim
sim.add(m=1) #Adds Sun 
sim.add(a=.5, f=0, Omega=0, omega=0, e=0, inc=0, m=0) #adds test particle 
#Moves all particles to center of momentum frame
sim.move_to_com()
#Gives orbital information before the simulation begins
print("\n***INITIAL ORBITS:***")
for orbit in sim.calculate_orbits():
    print(orbit)
Explanation: Yarkovsky Effect
This example shows how to add the Yarkovsky effect in a Rebound simulation. There are two versions, which we call the 'Full Version' and the 'Simple Version.' A special parameter called 'ye_flag' is used to switch between the two. The difference between the versions and what situations they're better suited for is discussed in more detail below. 
For more information on this effect, please visit: (implementation paper in progress) 
We'll start with the Full Version.
Full Version
This version of the effect is based off of the equations found in Veras et al. (2015). A link to the paper is provided below. The Full Version can be used to get detailed calculations of the Yarkovsky effect on a particular body. However, it requires a large amount of parameters that may be difficult to find for a particular object. It also takes more computational time due to the large amount of equations that must be calaculated between each time step of the simulation. This version of the effect can be used to get accurate calculations on how a body is perturbed by the Yarkovsky effect.
Link to paper: https://ui.adsabs.harvard.edu/abs/2015MNRAS.451.2814V/abstract
Below is a simple example to show how to add the effect to a simulation. First, we create a Rebound simulation with the Sun and a test particle (which will be considered an asteroid) at .5 AU.
End of explanation
density = (3000.0*u.kg/u.m**3).to(u.Msun/u.AU**3)
c = (constants.c).to(u.AU/u.yr) #speed of light
lstar = (3.828e26*u.kg*u.m**2/u.s**3).to(u.Msun*u.AU**2/u.yr**3) #luminosity of star
radius = (1000*u.m).to(u.AU) #radius of object
albedo = .017 #albedo of object
stef_boltz = constants.sigma_sb.to(u.Msun/u.yr**3/u.K**4) #Stefan-Boltzmann constant
emissivity = .9 #emissivity of object
k = .25 #constant between
Gamma = (310*u.kg/u.s**(5/2)).to(u.Msun/u.yr**(5/2)) #thermal inertia of object
rotation_period = (15470.9*u.s).to(u.yr) #rotation period of object
Explanation: As with all REBOUNDx effects, the parameters must be inputed with the same units as the simulation (in this case it's AU/Msun/yr). We'll use the astropy units module to help avoid errors
End of explanation
#Loads the effect into Rebound
rebx = reboundx.Extras(sim)
yark = rebx.load_force("yarkovsky_effect")
#Sets the parameters for the effect
yark.params["ye_c"] = c.value #set on the sim and not a particular particle
yark.params["ye_lstar"] = lstar.value #set on the sim and not a particular particle
yark.params["ye_stef_boltz"] = stef_boltz.value #set on the sim and not a particular particle
Explanation: We then add the Yarkovsky effect and the required parameters for this version. Importantly, we must set 'ye_flag' to 0 to get the Full Version. Physical constants and the stellar luminosity get added to the effect yark
End of explanation
# Sets parameters for the particle
ps = sim.particles
ps[1].r = radius.value #remember radius is not inputed as a Rebx parameter - it's inputed on the particle in the Rebound sim
ps[1].params["ye_flag"] = 0 #setting this flag to 0 will give us the full version of the effect
ps[1].params["ye_body_density"] = density.value
ps[1].params["ye_albedo"] = albedo
ps[1].params["ye_emissivity"] = emissivity
ps[1].params["ye_k"] = k
ps[1].params["ye_thermal_inertia"] = Gamma.value
ps[1].params["ye_rotation_period"] = rotation_period.value
# For this example we assume the object has a spin axis perpendicular to the orbital plane: unit vector = (0,0,1)
ps[1].params["ye_spin_axis_x"] = 0
ps[1].params["ye_spin_axis_y"] = 0
ps[1].params["ye_spin_axis_z"] = 1
rebx.add_force(yark) #adds the force to the simulation
Explanation: Other parameters need to be added to each particle feeling the Yarkovsky effect
End of explanation
%%time
tmax=100000 # in yrs
Nout = 1000
times = np.linspace(0, tmax, Nout)
a_start = .5 #starting semi-major axis for the asteroid
a = np.zeros(Nout)
for i, time in enumerate(times):
    a[i] = ps[1].a
    sim.integrate(time)
a_final = ps[1].a #semi-major axis of asteroid after the sim    
                      
print("CHANGE IN SEMI-MAJOR AXIS:", a_final-a_start, "AU\n") #prints difference between the initial and final semi-major axes of asteroid
fig, ax = plt.subplots()
ax.plot(times, a-a_start, '.')
ax.set_xlabel('Time (yrs)')
ax.set_ylabel('Change in semimajor axis (AU)')
Explanation: We integrate this system for 100,000 years and print out the difference between the particle's semi-major axis before and after the simulation.
End of explanation
sim = rebound.Simulation()
sim.units = ('yr', 'AU', 'Msun') #changes simulation and G to units of solar masses, years, and AU  
sim.integrator = "whfast" #integrator for sim
sim.dt = .05 #timestep for sim
sim.add(m=1) #Adds Sun 
sim.add(a=.5, f=0, Omega=0, omega=0, e=0, inc=0, m=0) #adds test particle 
sim.add(a=.75, f=0, Omega=0, omega=0, e=0, inc=0, m=0) #adds a second test particle
#Moves all particles to center of momentum frame
sim.move_to_com()
#Gives orbital information before the simulation begins
print("\n***INITIAL ORBITS:***")
for orbit in sim.calculate_orbits():
    print(orbit)
Explanation: Simple Version
This version of the effect is based off of equations from Veras et al. (2019). Once again, a link to this paper is provided below. This version simplifies the equations by placing constant values in a rotation matrix that in general is time-dependent. It requires fewer parameters than the full version and takes less computational time. However, it is mostly useful only to get a general idea on how much the effect can push bodies inwards or outwards. This version of the effect is better for simulating large groups of asteroids or trying to see general trends in the behavior of a body. 
Link to paper: https://ui.adsabs.harvard.edu/abs/2019MNRAS.485..708V/abstract
We'll use the same setup as before, but we'll also add another asteroid at .75 AU with identical physical properties. Let's start by creating a Rebound simulation again.
End of explanation
#Loads the effect into Rebound
rebx = reboundx.Extras(sim)
yark = rebx.load_force("yarkovsky_effect")
#Sets the parameters for the effect
yark.params["ye_c"] = c.value
yark.params["ye_lstar"] = lstar.value
ps = sim.particles #simplifies way to access particles parameters 
ps[1].params["ye_flag"] = 1 #setting this flag to 1 will give us the outward version of the effect 
ps[1].params["ye_body_density"] = density.value
ps[1].params["ye_albedo"] = albedo
ps[1].r = radius.value #remember radius is not inputed as a Rebx parameter - it's inputed on the particle in the Rebound sim
ps[2].params["ye_flag"] = -1 #setting this flag to -1 will give us the inward version of the effect 
ps[2].params["ye_body_density"] = density.value
ps[2].params["ye_albedo"] = albedo
ps[2].r = radius.value 
rebx.add_force(yark) #adds the force to the simulation
Explanation: We then add the Yarkovsky effect from Reboundx and the necesary parameters for this version. This time, we must make sure that 'ye_flag' is set to 1 or -1 to get the Simple Version of the effect. Setting it to 1 will push the asteroid outwards, while setting it to -1 will push it inwards. We'll push out our original asteroid and push in our new one. We use the same physical properties as in the example above:
End of explanation
%%time
tmax=100000 # in yrs
a_start_1 = .5 #starting semi-major axis for the 1st asteroid
a_start_2 = .75 #starting semi-major axis for the 2nd asteroid
a1, a2 = np.zeros(Nout), np.zeros(Nout)
for i, time in enumerate(times):
    a1[i] = ps[1].a
    a2[i] = ps[2].a
    sim.integrate(time)
a_final_1 = ps[1].a #semi-major axis of 1st asteroid after the sim
a_final_2 = ps[2].a #semi-major axis of 2nd asteroid after the sim
                      
print("CHANGE IN SEMI-MAJOR AXIS(Asteroid 1):", a_final_1-a_start_1, "AU\n")
print("CHANGE IN SEMI-MAJOR AXIS(Asteroid 2):", a_final_2-a_start_2, "AU\n")
fig, ax = plt.subplots()
ax.plot(times, a1-a_start_1, '.', label='Asteroid 1')
ax.plot(times, a2-a_start_2, '.', label='Asteroid 2')
ax.set_xlabel('Time (yrs)')
ax.set_ylabel('Change in semimajor axis (AU)')
ax.legend()
Explanation: Now we run the sim for 100,000 years and print out the results for both asteroids. Note the difference in simulation times between the versions. Even with an extra particle, the simple version was faster than the full version.
End of explanation
<END_TASK>