-
Notifications
You must be signed in to change notification settings - Fork 0
/
frames.Rmd
212 lines (158 loc) · 9.11 KB
/
frames.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
---
title: "Reference frames"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
if(!require("r2d3")) { install.packages('r2d3', repos = "http://cran.us.r-project.org"); library("r2d3") }
if(!require("ggplot2")) { install.packages('ggplot2', repos = "http://cran.us.r-project.org"); library("ggplot2") }
if(!require("gifski")) { install.packages('gifski', repos = "http://cran.us.r-project.org"); library("gifski") }
if(!require("gganimate")) { install.packages(c('stringi', 'gganimate'), repos = "http://cran.us.r-project.org"); library("gganimate") }
```
***
A very good place to start analyzing any motion is choosing a frame of reference and setting up coordinates. Orbital mechanics systems are often considered with respect to the *fixed stars*.
We usually try to find a reference frame that is not accelerating. This is hard! Everything turns out to be wagging a little, even what we can't see or feel. But supposing we do find an *inertial frame*, we would see that Newton's [laws of motion](https://en.wikipedia.org/wiki/Newton's_laws_of_motion) apply:
1. A body remains at rest, or in motion at a constant speed in a straight line, unless acted upon by a force.
2. When a body is acted upon by a force, the time rate of change of its momentum equals the force.
3. If two bodies exert forces on each other, these forces have the same magnitude but opposite directions.
In this way we avoid having to think through [fictitious forces](https://en.wikipedia.org/wiki/Fictitious_force) such as the centrifugal or Coriolis.
## A system of two bodies
We first look at the mutual gravitation of two masses, $m_0$ and $m_1$ with positions ${\bf r}_0$ and ${\bf r}_1$ described with their [Cartesian coordinates](https://en.wikipedia.org/wiki/Cartesian_coordinate_system):
$$
{\bf r}_0 = x_0 {\bf \hat{x}} + y_0 {\bf \hat{y}} + z_0 {\bf \hat{z}} \\
{\bf r}_1 = x_1 {\bf \hat{x}} + y_1 {\bf \hat{y}} + z_1 {\bf \hat{z}}
$$
The relative position, ${\bf r}$ and its unit vector, $\hat{\bf u}_\mathrm{r}$ can be calculated.
$$
{\bf r} = {\bf r}_1 - {\bf r}_0 = (x_1 - x_0) {\bf \hat{i}} + (y_1 - y_0) {\bf \hat{j}} + (z_1 - z_0) {\bf \hat{k}} \\
\hat{\bf u}_\mathrm{r} = \frac{\bf r}{r}
$$
***
**Example**
```{r}
df <- data.frame(m = 1.0e26, x = 100, y = 120)
df[2, ] <- data.frame(m = 5.0e25, x = 150, y = 50)
```
```{r echo=FALSE, fig.height=2, out.width="50%"}
r2d3(script="../scripts/rcog.js", data=df, options=list(relative = "true", lines = "true", text = "true", masses = "true"))
```
***
Here we are imagining the only two bodies in the universe! (Or that everyone else is really, really, far away.) So there are no other forces acting and the force on $m_0$ is opposite that on $m_1$ by Newtons' third law, i.e. ${\bf F}_0 = -{\bf F}_1$. Using the [law of gravitation](https://en.wikipedia.org/wiki/Newton%27s_law_of_universal_gravitation) ($F = G \frac{m_0 m_1}{r^2}$) we can find the force acting on the line between the masses.
$$
\begin{align}
{\bf F}_0 &= \frac{G m_0 m_1}{r^2} \hat{\bf u}_\mathrm{r} \\
{\bf F}_1 &= - \frac{G m_0 m_1}{r^2} \hat{\bf u}_\mathrm{r}
\end{align}
$$
Now applying Newton's second law (${\bf F}_i = m {\bf \ddot{r}}_i$) and the definition of $\hat{\bf u}_\mathrm{r}$ above we get the accelerations of the two masses, which depend only on the relative position.
$$
\begin{align}
\ddot{\bf r}_0 &= G m_1 \frac{\bf r}{r^3} \\
\ddot{\bf r}_1 &= - G m_0 \frac{\bf r}{r^3}
\end{align}
$$
We can integrate the accelerations to get equations for the velocities and positions of each mass $m_\mathrm{i}$. Keeping in mind that in general the positions, ${\bf r}_\mathrm{i}$, velocities, $\dot{\bf r}_\mathrm{i}$, and accelerations, $\ddot{\bf r}_\mathrm{i}$, depend on time.
$$
\dot{\bf r}(t)_\mathrm{i} = \dot{\bf r}(0)_\mathrm{i} + \int_{0}^{t} \ddot{\bf r}(t)_\mathrm{i} \mathrm{d}t \\
{\bf r}(t)_\mathrm{i} = {\bf r}(0)_\mathrm{i} + \int_{0}^{t} \dot{\bf r}(t)_\mathrm{i} \mathrm{d}t
$$
We see that we will need the initial positions, ${\bf r}(0)_\mathrm{i}$, and velocities, $\dot{\bf r}(0)_\mathrm{i}$, to solve the motion. Together we call these 12 quantities the *state vector* of the system: $S = [{\bf r}_0, {\bf r}_1, {\bf \dot{r}}_0, {\bf \dot{r}}_1]$.
### Numerically in R
#### Initial conditions
We input the [gravitational constant](https://en.wikipedia.org/wiki/Gravitational_constant), $G = 6.674 \times 10^{-2} \text{km}^3 \text{kg}^{-1} \text{s}^{-2}$, and an equal mass for the particles, $m_0 = m_1 = 10^{26} \text{kg}$. For reference Earth's mass is about $6 \times 10^{24} \text{kg}$, so these are somewhat larger.
```{r}
G <- 6.674e-2
m_0 <- 1.0e26
m_1 <- m_0
```
We represent the state vector as a `data.frame` where each record contains the state vector of a single mass and a timestamp.
```{r}
s_vector <- data.frame(n_mass=0, ts=0, x=0, y=0, z=0, dot.x=0, dot.y=0, dot.z=0)
s_vector[2, ] <- data.frame(n_mass=1, ts=0, x=2000, y=0, z=0, dot.x=-2000, dot.y=500, dot.z=0)
```
We have $m_0$ starting out stationary at $(0, 0, 0)$ and $m_1$ starting out at $(2000, 0, 0)$ with some initially counterclockwise velocity.
#### Stepwise solutions
We can advance the state vector by an arbitrary time step $\Delta t$ if we assume the velocities and accelerations are constant over that time slice. This can be a good approximation if the time slice is small compared to the characteristic time scale of the motion.
$$
\dot{\bf r}(t_0 + \Delta t)_\mathrm{i} = \dot{\bf r}(t_0)_\mathrm{i} + \ddot{\bf r}(t_0)_\mathrm{i} \Delta t \\
{\bf r}(t_0 + \Delta t)_\mathrm{i} = {\bf r}(t_0)_\mathrm{i} + \dot{\bf r}(t_0)_\mathrm{i} \Delta t.
$$
We develop a function that takes this state vector and returns a new state vector including the current relative distance and accelerations as well as the next positions and velocities.
<details>
<summary>Show code</summary>
```{r}
#' Get relative distance and accelerations from initial positions
#' @return A data frame containing the relative distance and accelerations
get_derived_initial_values <- function(pos_0, pos_1) {
# initial distance
r0 <- sqrt(sum((pos_1 - pos_0)^2))
# initial accelerations
ddot_r0_0 <- m_1 * G * (pos_1 - pos_0) / (r0^3)
ddot_r0_1 <- m_0 * G * (pos_0 - pos_1) / (r0^3)
init_vals <- data.frame(r=r0, ddot_r0=ddot_r0_0, ddot_r1=ddot_r0_1)
init_vals
}
#' Step forward linearly
step <- function(initial, dot_initial, delta) {
final <- initial + dot_initial * delta
final
}
#' Return new state vector after advancing by stepsize
#'
#' @return A data frame containing the original state vector with new values
#' appended after stepsize as well as the derived initial quantities
get_next_s_vector <- function(s_vector, stepsize=1) {
if (nrow(s_vector) < 2) {
return(666)
}
# extract current positions (last rows)
pos_0 <- s_vector[nrow(s_vector) - 1, 3:5]
pos_1 <- s_vector[nrow(s_vector), 3:5]
# extract current velocities (last rows)
vel_0 <- s_vector[nrow(s_vector) - 1, 6:8]
vel_1 <- s_vector[nrow(s_vector), 6:8]
# calculate current distance and accelerations
init_vals = get_derived_initial_values(pos_0, pos_1)
# annotate current state vector with derived initial values
s_vector[nrow(s_vector) - 1, 9:12] <- data.frame(init_vals[1], init_vals[2:4])
s_vector[nrow(s_vector), 9:12] <- data.frame(init_vals[1], init_vals[5:7])
# step forward
new_pos_0 <- step(pos_0, vel_0, stepsize)
new_pos_1 <- step(pos_1, vel_1, stepsize)
new_vel_0 <- step(vel_0, init_vals[2:4], stepsize)
new_vel_1 <- step(vel_1, init_vals[5:7], stepsize)
s_vector[nrow(s_vector) + 1, 1:8] <- data.frame(n_mass=0, ts=s_vector[nrow(s_vector) - 1, 2] + stepsize,
new_pos_0, new_vel_0)
s_vector[nrow(s_vector) + 1, 1:8] <- data.frame(n_mass=1, ts=s_vector[nrow(s_vector) - 1, 2] + stepsize,
new_pos_1, new_vel_1)
s_vector
}
```
</details>
And now we can see how the motion of the two masses evolves in time.
```{r}
stepsize <- 5e-10
new_s_vector <- s_vector
for(timestep in 1:58) { new_s_vector <- get_next_s_vector(new_s_vector, stepsize) }
```
Here are the first few time slices for $m_0$:
```{r echo = FALSE}
knitr::kable(head(new_s_vector[new_s_vector$n_mass == 0, ], 4), row.names=FALSE, digits=9)
```
and for $m_1$:
```{r echo = FALSE}
knitr::kable(head(new_s_vector[new_s_vector$n_mass == 1, ], 4), row.names=FALSE, digits=9)
```
The derived quantities $r$, and ${\bf \ddot{r}}$ are not populated for the last record as they are generated during the next time step.
All the motion has been set in the $x$-$y$ plane through the initial conditions.
```{r echo=FALSE, message=FALSE, warning=FALSE, out.width="50%"}
g = ggplot(new_s_vector, aes(y = y, x = x))
g = g + geom_point(data = new_s_vector, size=4)
g = g + transition_time(ts) + ease_aes('linear')
g
```
We can see the mutual attraction pull $m_0$ off of the origin as the two bodies approach each other, eventually getting within $80 \text{km}$ where the maximum gravitational acceleration slingshots them back away from each other at high speed.
***
[Center of Gravity ->](./center.html)
***
[Orbits!](../index.html) | [code](https://github.com/je-miralles/orbits/blob/main/rmd_sources/physics/frames.Rmd)