Copyright 2016, Pedro Belin Castellucci,
This work is licensed under a Creative Commons Attribution 4.0 International License.</i>
Gravity is a familiar concept. Every object near the Earth is subject to Earth gravitational force, which attract the object to its center. This gravitational force is associated with gravitational acceleration $g$. Near the surface of the Earth $g \approx 9.8\ m/s²$. However, assuming the acceleration is constant may be unrealistic for some applications. In fact, different parameters affect $g$ such as latitude and altitude. In this Notebook, we will play a little bit with the dependence of the gravitational acceleration and altitude. Our goal is not to do rigorous Math or Physics, but to use the context to learn some Julia programming basics (for some Math and Physics, the reader is refered to this link).
Instead of assuming a constant acceleration, a more realistic model for the gravity acceleration $g_e$ is:
$$g_e = g_0 \Big(\frac{r_e}{r_e + h}\Big)^2,$$in which $g_0$ is called the standard acceleration of gravity and is defined as $g_0 = 9.80665\ m/s^2$. The actual acceleration is $g_e$ and depends on $r_e$ and $h$, the average radius of Earth and the height (above sea level) of the point from which we are computing $g_e$, respectively.
Let us define Julia variables to keep $g_0 = 9.80665\ m/s^2$ and the average radius of earth $r_e = 6 378 000\ m$.
In [1]:
g0 = 9.80665 # g (m/s²)
re = 6378000 # Average radius of Earth in meters (m).
Out[1]:
Note that we do not need to declare the type of our variables. However they do have a type, which we can check with typeof.
In [2]:
typeof(g0), typeof(re), typeof("Earth")
Out[2]:
Julia supports Unicode variable names. This allows the definition of variables with more mathy names, e.g. with subscripts or with greek symbols. To declare a variable $\alpha$, for example, we can type \alpha{TAB}. Thus, let us redefine our variables with the g\_0{TAB} and r\_e{TAB}.
In [3]:
g₀ = g0
rₑ = re
Out[3]:
In [4]:
function localGravity(g₀, rₑ, h)
g₀*(rₑ/(rₑ + h))^2 # This is the return of the function
end
Out[4]:
Note that the return of a function is the statement in its last line. Another way of defining the same function is using a one line definition.
In [5]:
localGravity(g₀, rₑ, h) = g₀*(rₑ/(rₑ + h))^2
Out[5]:
In [6]:
Δt = 0.5 # Time step.
t = 0 # Initial time.
T = 4 # Final time.
while t <= T
println(t, "s")
t += Δt
end
Using a while-loop, we can simulate a vertical launch of a particle. We will assume that the variation of position and velocity of our particle is given by:
$$h_k = h_{k-1} + v_{k-1} \Delta t,$$$$v_k = v_{k-1} - g_e^{k-1} \Delta t,$$in which $h_k$ and $v_k$ are the height and velocity at time $k$, respectively. $g_e^{k-1}$ is given by the localGravity function and $\Delta t$ is the discrete time step we choose for our simulation. As a first example, we will consider $(h_0, v_0) = (0, 10)$ and simulate 2.2 seconds with $\Delta t = 0.1s.$
In [7]:
# Initial conditions:
v0 = 10.0 # m/s
h0 = 0.0 # m
time = 0.0 # s
# Simulation time:
endOfSimulation = 2.2
# Time step:
Δt = 0.1 # s
while time <= endOfSimulation
# Printing a formatted output:
toPrint = @sprintf("t = %.2fs \t h = %.2fm", time, h0)
println(toPrint)
# Updating height and velocity:
h0, v0 = h0 + v0*Δt, v0 - localGravity(g₀, rₑ, h0)*Δt
# Incrementing time:
time += Δt
end
We printed out the height of the particle as time evolves. To make it a little prettier, we used the macro @sprintf. We are not going to discuss macros here, but superficially they are like functions.
However, we can do a better job with the output. Julia has different plotting libraries, we are going to use Plots to plot the trajectory of the particle over time. To install it you can use the following:
In [8]:
Pkg.add("Plots");
In [9]:
myRange = 1:10
Out[9]:
We can check the type of the object.
In [10]:
typeof(myRange)
Out[10]:
To check the elements inside the range object, we can use collect.
In [11]:
collect(myRange)
Out[11]:
Note that the syntax a:b creates a collection of the integer number in the interval $[a, b]$. By default, each element is separated by a unit but we can also specify the step of the range to be created.
In [12]:
collect(1:0.1:2)
Out[12]:
Exercise: Use ranges to create the following arrays:
In [ ]:
For-loops are suitable to iterate over range objects.
In [14]:
for i in myRange
print("$i ")
end
In Julia, true and false are Bool types.
In [15]:
5 == 2, 5 >= 2, 1 < 2 && 2 < 3, 1 < 2 < 3, typeof(5 == 2)
Out[15]:
The Bool types can be used in conditional statements like the one that follows.
In [16]:
a, b = 5, 2
if a == b
println("The numbers are equal.")
else
println("The numbers are different.")
end
Exercise: Write an algorithm to tell which integer numbers between 1 and 10 are even and odd.
In [ ]:
In [17]:
a = [1, 3, 4, 2]
Out[17]:
We have an array with four Int64 elements. Arrays are indexed starting at one, knowing that, try to predict the outcome of each following statement.
In [18]:
println(a[1])
println(a[2])
println(a[end])
println(a[1:3])
println(a[1:end])
Exercise: try to create an array with values of different types.
In [19]:
b = [1, 2.0, 3, "String"]
Out[19]:
We can add elements to an array using the push! method.
In [20]:
push!(a, 10)
Out[20]:
As mentioned, arrays can also be used in a mathematical context.
In [21]:
println(2a)
println(2 + a)
println([1, 2, 3] + [1, 2, 3])
Julia also provides different mathematical operations. Like the dot product.
In [22]:
dot([1, 2, 3], [1, 2, 3])
Out[22]:
And the cross product.
In [23]:
cross([1, 2, 3], [1, 1, 1])
Out[23]:
We can check all the operations available issuing the following command.
In [24]:
methodswith(Array)
Out[24]:
Getting back to our example, we are ready to plot the trajectory of the particle. We will use an array to keep track of the positions of the particle over time.
In [25]:
using Plots # Importing the library we will use.
h0 = 0.0 # Initial height (m).
v0 = 10.0 # Initial velocity (m/s).
endOfSimulation = 2.2 # We will stop the simulation at 2.2s.
posList = [] # The array to store the positions.
Δt = 0.1 # Time step
for i in 0:Δt:endOfSimulation
push!(posList, h0) # Storing the position in posList array.
# Updating position and velocity:
h0, v0 = h0 + v0*Δt, v0 - localGravity(g₀, rₑ, h0)*Δt
if h0 < 0
# If the particle returns to the ground, we stop the simulation:
endOfSimulation = i
break
end
end
# To plot the results we use the plot function from Plots:
plot(0:Δt:endOfSimulation, posList, xlabel="Time", ylabel="Height", w=3, label="")
Out[25]:
The objects reaches a height of nearly $5.5\ m$ when launched with a velocity of $10\ m/s$.
There is the concept of Escape Velocity, which is the velocity something needs to be launched from a planet's surface to escape its gravitational influence (and getting lost in space forever). So, what is the escape velocity of our model? We can try to estimate it, empirically. First, let us convert our simulation code into a function.
In [26]:
function simul(h0, v0, endOfSimulation)
posList = [] # The array to store the positions.
Δt = 0.1 # Time step
for i = 0:Δt:endOfSimulation
push!(posList, h0) # Storing the position in posList array.
# Updating position and velocity:
h0, v0 = h0 + v0*Δt, v0 - localGravity(g₀, rₑ, h0)*Δt
if h0 < 0
# If the particle returns to the ground we stop the simulation
endOfSimulation = i
break
end
end
0:Δt:endOfSimulation, posList # This is the return of our function
end
Out[26]:
Now, we plot the trajectory of the particle for different initial velocities.
In [67]:
fig = plot()
for v0 in 11500:500:11500
time, toPlot = simul(0, v0, 400000)
plot!(time, toPlot, label=string(v0) * " m/s", w=3)
end
fig
Out[67]:
We can see from the plot that the escape velocity is bigger than 10500 $m/s$. The actual Earth escape velocity is estimated to be 11200 $m/s$. You can play with the values to see if you can get close to the actual value from our model.
In this Notebook, we explored some basic features of Julia, namely, variables and function declarations, while and for-loops and used the Plots library to visualize some output.