lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
8aa83adbdcb1544bf9916c0cc0d4de0049dd4668
FYS-STK4155/doc/src/LectureNotes/_build/html/testbook/chapter6.html
T
mhjensen 8aa83adbdc update on book
2020-12-19 22:47:17 +01:00

1979 lines
156 KiB
HTML
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>Two-body Problems &#8212; Applied Machine Learning and Data Analysis</title>
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.11.2/css/all.min.css" integrity="sha384-KA6wR/X5RY4zFAHpv/CnoG2UW1uogYfdnP67Uv7eULvTveboZJg0qUpmJZb5VqzN" crossorigin="anonymous">
<link href="../_static/css/index.css" rel="stylesheet">
<link rel="stylesheet" href="../_static/sphinx-book-theme.css" type="text/css" />
<link rel="stylesheet" href="../_static/pygments.css" type="text/css" />
<link rel="stylesheet" type="text/css" href="../_static/togglebutton.css" />
<link rel="stylesheet" type="text/css" href="../_static/copybutton.css" />
<link rel="stylesheet" type="text/css" href="../_static/mystnb.css" />
<link rel="stylesheet" type="text/css" href="../_static/sphinx-thebe.css" />
<link rel="stylesheet" type="text/css" href="../_static/panels-main.c949a650a448cc0ae9fd3441c0e17fb0.css" />
<link rel="stylesheet" type="text/css" href="../_static/panels-variables.06eb56fa6e07937060861dad626602ad.css" />
<script id="documentation_options" data-url_root="../" src="../_static/documentation_options.js"></script>
<script src="../_static/jquery.js"></script>
<script src="../_static/underscore.js"></script>
<script src="../_static/doctools.js"></script>
<script src="../_static/language_data.js"></script>
<script src="../_static/togglebutton.js"></script>
<script src="../_static/clipboard.min.js"></script>
<script src="../_static/copybutton.js"></script>
<script src="../_static/sphinx-book-theme.js"></script>
<script >var togglebuttonSelector = '.toggle, .admonition.dropdown, .tag_hide_input div.cell_input, .tag_hide-input div.cell_input, .tag_hide_output div.cell_output, .tag_hide-output div.cell_output, .tag_hide_cell.cell, .tag_hide-cell.cell';</script>
<script async="async" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.7/latest.js?config=TeX-AMS-MML_HTMLorMML"></script>
<script type="text/x-mathjax-config">MathJax.Hub.Config({"tex2jax": {"inlineMath": [["\\(", "\\)"]], "displayMath": [["\\[", "\\]"]], "processRefs": false, "processEnvironments": false}})</script>
<script async="async" src="https://unpkg.com/thebelab@latest/lib/index.js"></script>
<script >
const thebe_selector = ".thebe"
const thebe_selector_input = "pre"
const thebe_selector_output = ".output"
</script>
<script async="async" src="../_static/sphinx-thebe.js"></script>
<link rel="index" title="Index" href="../genindex.html" />
<link rel="search" title="Search" href="../search.html" />
<meta name="viewport" content="width=device-width, initial-scale=1">
<meta name="docsearch:language" content="en">
</head>
<body data-spy="scroll" data-target="#bd-toc-nav" data-offset="80">
<div class="container-xl">
<div class="row">
<div class="col-12 col-md-3 bd-sidebar site-navigation show" id="site-navigation">
<div class="navbar-brand-box">
<a class="navbar-brand text-wrap" href="../index.html">
<img src="../_static/Picture.png" class="logo" alt="logo">
<h1 class="site-logo" id="site-title">Applied Machine Learning and Data Analysis</h1>
</a>
</div>
<form class="bd-search d-flex align-items-center" action="../search.html" method="get">
<i class="icon fas fa-search"></i>
<input type="search" class="form-control" name="q" id="search-input" placeholder="Search this book..." aria-label="Search this book..." autocomplete="off" >
</form>
<nav class="bd-links" id="bd-docs-nav" aria-label="Main navigation">
<ul class="nav sidenav_l1">
<li class="toctree-l1">
<a class="reference internal" href="../chapter1.html">
Introduction to Applied Data Analysis and Machine Learning
</a>
</li>
</ul>
<p class="caption">
<span class="caption-text">
Supervised Learning
</span>
</p>
<ul class="nav sidenav_l1">
<li class="toctree-l1">
<a class="reference internal" href="../chapter2.html">
1. Elements of Probability Theory and Statistical Data Analysis
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter2.html#random-numbers">
2. Random Numbers
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter2.html#random-numbers-better-name-pseudo-random-numbers">
3. Random Numbers, better name: pseudo random numbers
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter2.html#random-number-generator-rng">
4. Random number generator RNG
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter2.html#random-number-generator-rng-and-periodic-outputs">
5. Random number generator RNG and periodic outputs
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter2.html#random-number-generator-rng-and-its-period">
6. Random number generator RNG and its period
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter2.html#random-number-generator-rng-other-examples">
7. Random number generator RNG, other examples
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter2.html#id9">
8. Random number generator RNG, other examples
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter2.html#random-number-generator-rng-ran0">
9. Random number generator RNG, RAN0
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter2.html#id10">
10. Random number generator RNG, RAN0
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter2.html#id11">
11. Random number generator RNG, RAN0
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter2.html#id12">
12. Random number generator RNG, RAN0
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter2.html#random-number-generator-rng-ran0-code">
13. Random number generator RNG, RAN0 code
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter2.html#which-rng-should-i-use">
14. Which RNG should I use?
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter3.html">
15. Getting started, our first data and Machine Learning encounters
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter4.html">
16. Linear Regression and more Advanced Regression Analysis
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter5.html">
17. Logistic Regression
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter6.html">
18. Neural networks, from the simple perceptron to deep learning
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter7.html">
19. Support Vector Machines, overarching aims
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="../chapter8.html">
20. Dimensionality Reduction
</a>
</li>
</ul>
</nav>
<!-- To handle the deprecated key -->
<div class="navbar_extra_footer">
Powered by <a href="https://jupyterbook.org">Jupyter Book</a>
</div>
</div>
<main class="col py-md-3 pl-md-4 bd-content overflow-auto" role="main">
<div class="row topbar fixed-top container-xl">
<div class="col-12 col-md-3 bd-topbar-whitespace site-navigation show">
</div>
<div class="col pl-2 topbar-main">
<button id="navbar-toggler" class="navbar-toggler ml-0" type="button" data-toggle="collapse"
data-toggle="tooltip" data-placement="bottom" data-target=".site-navigation" aria-controls="navbar-menu"
aria-expanded="true" aria-label="Toggle navigation" aria-controls="site-navigation"
title="Toggle navigation" data-toggle="tooltip" data-placement="left">
<i class="fas fa-bars"></i>
<i class="fas fa-arrow-left"></i>
<i class="fas fa-arrow-up"></i>
</button>
<div class="dropdown-buttons-trigger">
<button id="dropdown-buttons-trigger" class="btn btn-secondary topbarbtn" aria-label="Download this page"><i
class="fas fa-download"></i></button>
<div class="dropdown-buttons">
<!-- ipynb file if we had a myst markdown file -->
<!-- Download raw file -->
<a class="dropdown-buttons" href="../_sources/testbook/chapter6.ipynb"><button type="button"
class="btn btn-secondary topbarbtn" title="Download source file" data-toggle="tooltip"
data-placement="left">.ipynb</button></a>
<!-- Download PDF via print -->
<button type="button" id="download-print" class="btn btn-secondary topbarbtn" title="Print to PDF"
onClick="window.print()" data-toggle="tooltip" data-placement="left">.pdf</button>
</div>
</div>
<!-- Source interaction buttons -->
<!-- Full screen (wrap in <a> to have style consistency -->
<a class="full-screen-button"><button type="button" class="btn btn-secondary topbarbtn" data-toggle="tooltip"
data-placement="bottom" onclick="toggleFullScreen()" title="Fullscreen mode"><i
class="fas fa-expand"></i></button></a>
<!-- Launch buttons -->
</div>
<!-- Table of contents -->
<div class="d-none d-md-block col-md-2 bd-toc show">
<div class="tocsection onthispage pt-5 pb-3">
<i class="fas fa-list"></i> Contents
</div>
<nav id="bd-toc-nav">
<ul class="nav section-nav flex-column">
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#relative-and-center-of-mass-motion">
Relative and Center of Mass Motion
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#deriving-elliptical-orbits">
Deriving Elliptical Orbits
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#effective-or-centrifugal-potential">
Effective or Centrifugal Potential
</a>
<ul class="nav section-nav flex-column">
<li class="toc-h3 nav-item toc-entry">
<a class="reference internal nav-link" href="#gravitational-force-example">
Gravitational force example
</a>
</li>
<li class="toc-h3 nav-item toc-entry">
<a class="reference internal nav-link" href="#harmonic-oscillator-in-two-dimensions">
Harmonic Oscillator in two dimensions
</a>
</li>
</ul>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#stability-of-orbits">
Stability of Orbits
</a>
<ul class="nav section-nav flex-column">
<li class="toc-h3 nav-item toc-entry">
<a class="reference internal nav-link" href="#code-example-with-gravitional-force">
Code example with gravitional force
</a>
</li>
</ul>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#scattering-and-cross-sections">
Scattering and Cross Sections
</a>
<ul class="nav section-nav flex-column">
<li class="toc-h3 nav-item toc-entry">
<a class="reference internal nav-link" href="#solution">
Solution
</a>
</li>
</ul>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#rutherford-scattering">
Rutherford Scattering
</a>
<ul class="nav section-nav flex-column">
<li class="toc-h3 nav-item toc-entry">
<a class="reference internal nav-link" href="#id1">
Solution
</a>
</li>
</ul>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#tidal-forces">
Tidal Forces
</a>
</li>
</ul>
</nav>
</div>
</div>
<div id="main-content" class="row">
<div class="col-12 col-md-9 pl-md-3 pr-md-0">
<div>
<div class="section" id="two-body-problems">
<h1>Two-body Problems<a class="headerlink" href="#two-body-problems" title="Permalink to this headline"></a></h1>
<p>The gravitational potential energy and forces involving two masses <span class="math notranslate nohighlight">\(a\)</span> and <span class="math notranslate nohighlight">\(b\)</span> are</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
U_{ab}&amp;=&amp;-\frac{Gm_am_b}{|\boldsymbol{r}_a-\boldsymbol{r}_b|},\\
\nonumber
F_{ba}&amp;=&amp;-\frac{Gm_am_b}{|\boldsymbol{r}_a-\boldsymbol{r}_b|^2}\hat{r}_{ab},\\
\nonumber
\hat{r}_{ab}&amp;=&amp;\frac{\boldsymbol{r}_b-\boldsymbol{r}_a}{|\boldsymbol{r}_a-\boldsymbol{r}_b|}.
\end{eqnarray}
\end{split}\]</div>
<p>Here <span class="math notranslate nohighlight">\(G=6.67\times 10^{-11}\)</span> Nm<span class="math notranslate nohighlight">\(^2\)</span>/kg<span class="math notranslate nohighlight">\(^2\)</span>, and <span class="math notranslate nohighlight">\(F_{ba}\)</span> is the force
on <span class="math notranslate nohighlight">\(b\)</span> due to <span class="math notranslate nohighlight">\(a\)</span>. By inspection, one can see that the force on <span class="math notranslate nohighlight">\(b\)</span>
due to <span class="math notranslate nohighlight">\(a\)</span> and the force on <span class="math notranslate nohighlight">\(a\)</span> due to <span class="math notranslate nohighlight">\(b\)</span> are equal and opposite. The
net potential energy for a large number of masses would be</p>
<!-- Equation labels as ordinary links -->
<div id="_auto1"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
U=\sum_{a&lt;b}U_{ab}=\frac{1}{2}\sum_{a\ne b}U_{ab}.
\label{_auto1} \tag{1}
\end{equation}
\]</div>
<div class="section" id="relative-and-center-of-mass-motion">
<h2>Relative and Center of Mass Motion<a class="headerlink" href="#relative-and-center-of-mass-motion" title="Permalink to this headline"></a></h2>
<p>Thus far, we have considered the trajectory as if the force is
centered around a fixed point. For two bodies interacting only with
one another, both masses circulate around the center of mass. One
might think that solutions would become more complex when both
particles move, but we will see here that the problem can be reduced
to one with a single body moving according to a fixed force by
expressing the trajectories for <span class="math notranslate nohighlight">\(\boldsymbol{r}_1\)</span> and <span class="math notranslate nohighlight">\(\boldsymbol{r}_2\)</span> into the
center-of-mass coordinate <span class="math notranslate nohighlight">\(\boldsymbol{R}_{\rm cm}\)</span> and the relative
coordinate <span class="math notranslate nohighlight">\(\boldsymbol{r}\)</span>,</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
\boldsymbol{R}_{\rm cm}&amp;\equiv&amp;\frac{m_1\boldsymbol{r}_1+m_2\boldsymbol{r}_2}{m_1+m_2},\\
\nonumber
\boldsymbol{r}&amp;\equiv&amp;\boldsymbol{r}_1-\boldsymbol{r_2}.
\end{eqnarray}
\end{split}\]</div>
<p>Here, we assume the two particles interact only with one another, so
<span class="math notranslate nohighlight">\(\boldsymbol{F}_{12}=-\boldsymbol{F}_{21}\)</span> (where <span class="math notranslate nohighlight">\(\boldsymbol{F}_{ij}\)</span> is the force on <span class="math notranslate nohighlight">\(i\)</span>
due to <span class="math notranslate nohighlight">\(j\)</span>. The equations of motion then become</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
\ddot{\boldsymbol{R}}_{\rm cm}&amp;=&amp;\frac{1}{m_1+m_2}\left\{m_1\ddot{\boldsymbol{r}}_1+m_2\ddot{\boldsymbol{r}}_2\right\}\\
\nonumber
&amp;=&amp;\frac{1}{m_1+m_2}\left\{\boldsymbol{F}_{12}+\boldsymbol{F}_{21}\right\}=0.\\
\ddot{\boldsymbol{r}}&amp;=&amp;\ddot{\boldsymbol{r}}_1-\ddot{\boldsymbol{r}}_2=\left(\frac{\boldsymbol{F}_{12}}{m_1}-\frac{\boldsymbol{F}_{21}}{m_2}\right)\\
\nonumber
&amp;=&amp;\left(\frac{1}{m_1}+\frac{1}{m_2}\right)\boldsymbol{F}_{12}.
\end{eqnarray}
\end{split}\]</div>
<p>The first expression simply states that the center of mass coordinate
<span class="math notranslate nohighlight">\(\boldsymbol{R}_{\rm cm}\)</span> moves at a fixed velocity. The second expression
can be rewritten in terms of the reduced mass <span class="math notranslate nohighlight">\(\mu\)</span>.</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
\mu \ddot{\boldsymbol{r}}&amp;=&amp;\boldsymbol{F}_{12},\\
\frac{1}{\mu}&amp;=&amp;\frac{1}{m_1}+\frac{1}{m_2},~~~~\mu=\frac{m_1m_2}{m_1+m_2}.
\end{eqnarray}
\end{split}\]</div>
<p>Thus, one can treat the trajectory as a one-body problem where the
reduced mass is <span class="math notranslate nohighlight">\(\mu\)</span>, and a second trivial problem for the center of
mass. The reduced mass is especially convenient when one is
considering gravitational problems because then</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
\mu \ddot{r}&amp;=&amp;-\frac{Gm_1m_2}{r^2}\hat{r}\\
\nonumber
&amp;=&amp;-\frac{GM\mu}{r^2}\hat{r},~~~M\equiv m_1+m_2.
\end{eqnarray}
\end{split}\]</div>
<p>For the gravitational problem, the reduced mass then falls out and the
trajectory depends only on the total mass <span class="math notranslate nohighlight">\(M\)</span>.</p>
<p>The kinetic energy and momenta also have analogues in center-of-mass
coordinates. The total and relative momenta are</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
\boldsymbol{P}&amp;\equiv&amp;\boldsymbol{p}_1+\boldsymbol{p}_2=M\dot{\boldsymbol{R}}_{\rm cm},\\
\nonumber
\boldsymbol{q}&amp;\equiv&amp;\mu\dot{\boldsymbol{r}}.
\end{eqnarray}
\end{split}\]</div>
<p>With these definitions, a little algebra shows that the kinetic energy becomes</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
T&amp;=&amp;\frac{1}{2}m_1|\boldsymbol{v}_1|^2+\frac{1}{2}m_2|\boldsymbol{v}_2|^2\\
\nonumber
&amp;=&amp;\frac{1}{2}M|\dot{\boldsymbol{R}}_{\rm cm}|^2
+\frac{1}{2}\mu|\dot{\boldsymbol{r}}|^2\\
\nonumber
&amp;=&amp;\frac{P^2}{2M}+\frac{q^2}{2\mu}.
\end{eqnarray}
\end{split}\]</div>
<p>The standard strategy is to transform into the center of mass frame,
then treat the problem as one of a single particle of mass <span class="math notranslate nohighlight">\(\mu\)</span>
undergoing a force <span class="math notranslate nohighlight">\(\boldsymbol{F}_{12}\)</span>. Scattering angles can also be
expressed in this frame, then transformed into the lab frame. In
practice, one sees examples in the literature where <span class="math notranslate nohighlight">\(d\sigma/d\Omega\)</span>
expressed in both the “center-of-mass” and in the “laboratory”
frame.</p>
</div>
<div class="section" id="deriving-elliptical-orbits">
<h2>Deriving Elliptical Orbits<a class="headerlink" href="#deriving-elliptical-orbits" title="Permalink to this headline"></a></h2>
<p>Kepler’s laws state that a gravitational orbit should be an ellipse
with the source of the gravitational field at one focus. Deriving this
is surprisingly messy. To do this, we first use angular momentum
conservation to transform the equations of motion so that it is in
terms of <span class="math notranslate nohighlight">\(r\)</span> and <span class="math notranslate nohighlight">\(\theta\)</span> instead of <span class="math notranslate nohighlight">\(r\)</span> and <span class="math notranslate nohighlight">\(t\)</span>. The overall strategy
is to</p>
<ol class="simple">
<li><p>Find equations of motion for <span class="math notranslate nohighlight">\(r\)</span> and <span class="math notranslate nohighlight">\(t\)</span> with no angle (<span class="math notranslate nohighlight">\(\theta\)</span>) mentioned, i.e. <span class="math notranslate nohighlight">\(d^2r/dt^2=\cdots\)</span>. Angular momentum conservation will be used, and the equation will involve the angular momentum <span class="math notranslate nohighlight">\(L\)</span>.</p></li>
<li><p>Use angular momentum conservation to find an expression for <span class="math notranslate nohighlight">\(\dot{\theta}\)</span> in terms of <span class="math notranslate nohighlight">\(r\)</span>.</p></li>
<li><p>Use the chain rule to convert the equations of motions for <span class="math notranslate nohighlight">\(r\)</span>, an expression involving <span class="math notranslate nohighlight">\(r,\dot{r}\)</span> and <span class="math notranslate nohighlight">\(\ddot{r}\)</span>, to one involving <span class="math notranslate nohighlight">\(r,dr/d\theta\)</span> and <span class="math notranslate nohighlight">\(d^2r/d\theta^2\)</span>. This is quitecomplicated because the expressions will also involve a substitution <span class="math notranslate nohighlight">\(u=1/r\)</span> so that one finds an expression in terms of <span class="math notranslate nohighlight">\(u\)</span> and <span class="math notranslate nohighlight">\(\theta\)</span>.</p></li>
<li><p>Once <span class="math notranslate nohighlight">\(u(\theta)\)</span> is found, you need to show that this can be converted to the familiar form for an ellipse.</p></li>
</ol>
<p>The equations of motion give</p>
<!-- Equation labels as ordinary links -->
<div id="eq:radialeqofmotion"></div>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
\label{eq:radialeqofmotion} \tag{2}
\frac{d}{dt}r^2&amp;=&amp;\frac{d}{dt}(x^2+y^2)=2x\dot{x}+2y\dot{y}=2r\dot{r},\\
\nonumber
\dot{r}&amp;=&amp;\frac{x}{r}\dot{x}+\frac{y}{r}\dot{y},\\
\nonumber
\ddot{r}&amp;=&amp;\frac{x}{r}\ddot{x}+\frac{y}{r}\ddot{y}
+\frac{\dot{x}^2+\dot{y}^2}{r}
-\frac{\dot{r}^2}{r}.
\end{eqnarray}
\end{split}\]</div>
<p>Recognizing that the numerator of the third term is the velocity squared, and that it can be written in polar coordinates,</p>
<!-- Equation labels as ordinary links -->
<div id="_auto2"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
v^2=\dot{x}^2+\dot{y}^2=\dot{r}^2+r^2\dot{\theta}^2,
\label{_auto2} \tag{3}
\end{equation}
\]</div>
<p>one can write <span class="math notranslate nohighlight">\(\ddot{r}\)</span> as</p>
<!-- Equation labels as ordinary links -->
<div id="eq:radialeqofmotion2"></div>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
\label{eq:radialeqofmotion2} \tag{4}
\ddot{r}&amp;=&amp;\frac{F_x\cos\theta+F_y\sin\theta}{m}+\frac{\dot{r}^2+r^2\dot{\theta}^2}{r}-\frac{\dot{r}^2}{r}\\
\nonumber
&amp;=&amp;\frac{F}{m}+\frac{r^2\dot{\theta}^2}{r}\\
\nonumber
m\ddot{r}&amp;=&amp;F+\frac{L^2}{mr^3}.
\end{eqnarray}
\end{split}\]</div>
<p>This derivation used the fact that the force was radial,
<span class="math notranslate nohighlight">\(F=F_r=F_x\cos\theta+F_y\sin\theta\)</span>, and that angular momentum is
<span class="math notranslate nohighlight">\(L=mrv_{\theta}=mr^2\dot{\theta}\)</span>. The term <span class="math notranslate nohighlight">\(L^2/mr^3=mv^2/r\)</span> behaves
like an additional force. Sometimes this is referred to as a
centrifugal force, but it is not a force. Instead, it is the
consequence of considering the motion in a rotating (and therefore
accelerating) frame.</p>
<p>Now, we switch to the particular case of an attractive inverse square
force, <span class="math notranslate nohighlight">\(F=-\alpha/r^2\)</span>, and show that the trajectory, <span class="math notranslate nohighlight">\(r(\theta)\)</span>, is
an ellipse. To do this we transform derivatives w.r.t. time to
derivatives w.r.t. <span class="math notranslate nohighlight">\(\theta\)</span> using the chain rule combined with angular
momentum conservation, <span class="math notranslate nohighlight">\(\dot{\theta}=L/mr^2\)</span>.</p>
<!-- Equation labels as ordinary links -->
<div id="eq:rtotheta"></div>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
\label{eq:rtotheta} \tag{5}
\dot{r}&amp;=&amp;\frac{dr}{d\theta}\dot{\theta}=\frac{dr}{d\theta}\frac{L}{mr^2},\\
\nonumber
\ddot{r}&amp;=&amp;\frac{d^2r}{d\theta^2}\dot{\theta}^2
+\frac{dr}{d\theta}\left(\frac{d}{dr}\frac{L}{mr^2}\right)\dot{r}\\
\nonumber
&amp;=&amp;\frac{d^2r}{d\theta^2}\left(\frac{L}{mr^2}\right)^2
-2\frac{dr}{d\theta}\frac{L}{mr^3}\dot{r}\\
\nonumber
&amp;=&amp;\frac{d^2r}{d\theta^2}\left(\frac{L}{mr^2}\right)^2
-\frac{2}{r}\left(\frac{dr}{d\theta}\right)^2\left(\frac{L}{mr^2}\right)^2
\end{eqnarray}
\end{split}\]</div>
<p>Equating the two expressions for <span class="math notranslate nohighlight">\(\ddot{r}\)</span> in Eq.s (<a class="reference external" href="#eq:radialeqofmotion2">4</a>) and (<a class="reference external" href="#eq:rtotheta">5</a>) eliminates all the derivatives w.r.t. time, and provides a differential equation with only derivatives w.r.t. <span class="math notranslate nohighlight">\(\theta\)</span>,</p>
<!-- Equation labels as ordinary links -->
<div id="eq:rdotdot"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\label{eq:rdotdot} \tag{6}
\frac{d^2r}{d\theta^2}\left(\frac{L}{mr^2}\right)^2
-\frac{2}{r}\left(\frac{dr}{d\theta}\right)^2\left(\frac{L}{mr^2}\right)^2
=\frac{F}{m}+\frac{L^2}{m^2r^3},
\end{equation}
\]</div>
<p>that when solved yields the trajectory, i.e. <span class="math notranslate nohighlight">\(r(\theta)\)</span>. Up to this
point the expressions work for any radial force, not just forces that
fall as <span class="math notranslate nohighlight">\(1/r^2\)</span>.</p>
<p>The trick to simplifying this differential equation for the inverse
square problems is to make a substitution, <span class="math notranslate nohighlight">\(u\equiv 1/r\)</span>, and rewrite
the differential equation for <span class="math notranslate nohighlight">\(u(\theta)\)</span>.</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
r&amp;=&amp;1/u,\\
\nonumber
\frac{dr}{d\theta}&amp;=&amp;-\frac{1}{u^2}\frac{du}{d\theta},\\
\nonumber
\frac{d^2r}{d\theta^2}&amp;=&amp;\frac{2}{u^3}\left(\frac{du}{d\theta}\right)^2-\frac{1}{u^2}\frac{d^2u}{d\theta^2}.
\end{eqnarray}
\end{split}\]</div>
<p>Plugging these expressions into Eq. (<a class="reference external" href="#eq:rdotdot">6</a>) gives an
expression in terms of <span class="math notranslate nohighlight">\(u\)</span>, <span class="math notranslate nohighlight">\(du/d\theta\)</span>, and <span class="math notranslate nohighlight">\(d^2u/d\theta^2\)</span>. After
some tedious algebra,</p>
<!-- Equation labels as ordinary links -->
<div id="_auto3"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\frac{d^2u}{d\theta^2}=-u-\frac{F m}{L^2u^2}.
\label{_auto3} \tag{7}
\end{equation}
\]</div>
<p>For the attractive inverse square law force, <span class="math notranslate nohighlight">\(F=-\alpha u^2\)</span>,</p>
<!-- Equation labels as ordinary links -->
<div id="_auto4"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\frac{d^2u}{d\theta^2}=-u+\frac{m\alpha}{L^2}.
\label{_auto4} \tag{8}
\end{equation}
\]</div>
<p>The solution has two arbitrary constants, <span class="math notranslate nohighlight">\(A\)</span> and <span class="math notranslate nohighlight">\(\theta_0\)</span>,</p>
<!-- Equation labels as ordinary links -->
<div id="eq:Ctrajectory"></div>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
\label{eq:Ctrajectory} \tag{9}
u&amp;=&amp;\frac{m\alpha}{L^2}+A\cos(\theta-\theta_0),\\
\nonumber
r&amp;=&amp;\frac{1}{(m\alpha/L^2)+A\cos(\theta-\theta_0)}.
\end{eqnarray}
\end{split}\]</div>
<p>The radius will be at a minimum when <span class="math notranslate nohighlight">\(\theta=\theta_0\)</span> and at a
maximum when <span class="math notranslate nohighlight">\(\theta=\theta_0+\pi\)</span>. The constant <span class="math notranslate nohighlight">\(A\)</span> is related to the
eccentricity of the orbit. When <span class="math notranslate nohighlight">\(A=0\)</span> the radius is a constant
<span class="math notranslate nohighlight">\(r=L^2/(m\alpha)\)</span>, and the motion is circular. If one solved the
expression <span class="math notranslate nohighlight">\(mv^2/r=-\alpha/r^2\)</span> for a circular orbit, using the
substitution <span class="math notranslate nohighlight">\(v=L/(mr)\)</span>, one would reproduce the expression
<span class="math notranslate nohighlight">\(r=L^2/(m\alpha)\)</span>.</p>
<p>The form describing the elliptical trajectory in
Eq. (<a class="reference external" href="#eq:Ctrajectory">9</a>) can be identified as an ellipse with one
focus being the center of the ellipse by considering the definition of
an ellipse as being the points such that the sum of the two distances
between the two foci are a constant. Making that distance <span class="math notranslate nohighlight">\(2D\)</span>, the
distance between the two foci as <span class="math notranslate nohighlight">\(2a\)</span>, and putting one focus at the
origin,</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
2D&amp;=&amp;r+\sqrt{(r\cos\theta-2a)^2+r^2\sin^2\theta},\\
\nonumber
4D^2+r^2-4Dr&amp;=&amp;r^2+4a^2-4ar\cos\theta,\\
\nonumber
r&amp;=&amp;\frac{D^2-a^2}{D+a\cos\theta}=\frac{1}{D/(D^2-a^2)-a\cos\theta/(D^2-a^2)}.
\end{eqnarray}
\end{split}\]</div>
<p>By inspection, this is the same form as Eq. (<a class="reference external" href="#eq:Ctrajectory">9</a>) with <span class="math notranslate nohighlight">\(D/(D^2-a^2)=m\alpha/L^2\)</span> and <span class="math notranslate nohighlight">\(a/(D^2-a^2)=A\)</span>.</p>
<p>Let us remind ourselves about what an ellipse is before we proceed.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="o">%</span><span class="k">matplotlib</span> inline
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">from</span> <span class="nn">matplotlib</span> <span class="kn">import</span> <span class="n">pyplot</span> <span class="k">as</span> <span class="n">plt</span>
<span class="kn">from</span> <span class="nn">math</span> <span class="kn">import</span> <span class="n">pi</span>
<span class="n">u</span><span class="o">=</span><span class="mf">1.</span> <span class="c1">#x-position of the center</span>
<span class="n">v</span><span class="o">=</span><span class="mf">0.5</span> <span class="c1">#y-position of the center</span>
<span class="n">a</span><span class="o">=</span><span class="mf">2.</span> <span class="c1">#radius on the x-axis</span>
<span class="n">b</span><span class="o">=</span><span class="mf">1.5</span> <span class="c1">#radius on the y-axis</span>
<span class="n">t</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="p">,</span> <span class="mi">100</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span> <span class="n">u</span><span class="o">+</span><span class="n">a</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">t</span><span class="p">)</span> <span class="p">,</span> <span class="n">v</span><span class="o">+</span><span class="n">b</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">t</span><span class="p">)</span> <span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">grid</span><span class="p">(</span><span class="n">color</span><span class="o">=</span><span class="s1">&#39;lightgray&#39;</span><span class="p">,</span><span class="n">linestyle</span><span class="o">=</span><span class="s1">&#39;--&#39;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<img alt="../_images/chapter6_37_01.png" src="../_images/chapter6_37_01.png" />
</div>
</div>
</div>
<div class="section" id="effective-or-centrifugal-potential">
<h2>Effective or Centrifugal Potential<a class="headerlink" href="#effective-or-centrifugal-potential" title="Permalink to this headline"></a></h2>
<p>The total energy of a particle is</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
E&amp;=&amp;U(r)+\frac{1}{2}mv_\theta^2+\frac{1}{2}m\dot{r}^2\\
\nonumber
&amp;=&amp;U(r)+\frac{1}{2}mr^2\dot{\theta}^2+\frac{1}{2}m\dot{r}^2\\
\nonumber
&amp;=&amp;U(r)+\frac{L^2}{2mr^2}+\frac{1}{2}m\dot{r}^2.
\end{eqnarray}
\end{split}\]</div>
<p>The second term then contributes to the energy like an additional
repulsive potential. The term is sometimes referred to as the
“centrifugal” potential, even though it is actually the kinetic energy
of the angular motion. Combined with <span class="math notranslate nohighlight">\(U(r)\)</span>, it is sometimes referred
to as the “effective” potential,</p>
<div class="math notranslate nohighlight">
\[
\begin{eqnarray}
U_{\rm eff}(r)&amp;=&amp;U(r)+\frac{L^2}{2mr^2}.
\end{eqnarray}
\]</div>
<p>Note that if one treats the effective potential like a real potential, one would expect to be able to generate an effective force,</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
F_{\rm eff}&amp;=&amp;-\frac{d}{dr}U(r) -\frac{d}{dr}\frac{L^2}{2mr^2}\\
\nonumber
&amp;=&amp;F(r)+\frac{L^2}{mr^3}=F(r)+m\frac{v_\perp^2}{r},
\end{eqnarray}
\end{split}\]</div>
<p>which is indeed matches the form for <span class="math notranslate nohighlight">\(m\ddot{r}\)</span> in Eq. (<a class="reference external" href="#eq:radialeqofmotion2">4</a>), which included the <strong>centrifugal</strong> force.</p>
<p>The following code plots this effective potential for a simple choice of parameters, with a standard gravitational potential <span class="math notranslate nohighlight">\(-\alpha/r\)</span>. Here we have chosen <span class="math notranslate nohighlight">\(L=m=\alpha=1\)</span>.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1"># Common imports</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">from</span> <span class="nn">math</span> <span class="kn">import</span> <span class="o">*</span>
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="n">Deltax</span> <span class="o">=</span> <span class="mf">0.01</span>
<span class="c1">#set up arrays</span>
<span class="n">xinitial</span> <span class="o">=</span> <span class="mf">0.3</span>
<span class="n">xfinal</span> <span class="o">=</span> <span class="mf">5.0</span>
<span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.0</span> <span class="c1"># spring constant</span>
<span class="n">m</span> <span class="o">=</span> <span class="mf">1.0</span> <span class="c1"># mass, you can change these</span>
<span class="n">AngMom</span> <span class="o">=</span> <span class="mf">1.0</span> <span class="c1"># The angular momentum</span>
<span class="n">n</span> <span class="o">=</span> <span class="n">ceil</span><span class="p">((</span><span class="n">xfinal</span><span class="o">-</span><span class="n">xinitial</span><span class="p">)</span><span class="o">/</span><span class="n">Deltax</span><span class="p">)</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
<span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">xinitial</span><span class="o">+</span><span class="n">i</span><span class="o">*</span><span class="n">Deltax</span>
<span class="n">V</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">V</span> <span class="o">=</span> <span class="o">-</span><span class="n">alpha</span><span class="o">/</span><span class="n">x</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="n">AngMom</span><span class="o">*</span><span class="n">AngMom</span><span class="o">/</span><span class="p">(</span><span class="n">m</span><span class="o">*</span><span class="n">x</span><span class="o">*</span><span class="n">x</span><span class="p">)</span>
<span class="c1"># Plot potential</span>
<span class="n">fig</span><span class="p">,</span> <span class="n">ax</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">subplots</span><span class="p">()</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;r[m]&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;V[J]&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">V</span><span class="p">)</span>
<span class="n">fig</span><span class="o">.</span><span class="n">tight_layout</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<img alt="../_images/chapter6_45_01.png" src="../_images/chapter6_45_01.png" />
</div>
</div>
<div class="section" id="gravitational-force-example">
<h3>Gravitational force example<a class="headerlink" href="#gravitational-force-example" title="Permalink to this headline"></a></h3>
<p>Using the above parameters, we can now study the evolution of the system using for example the velocity Verlet method.
This is done in the code here for an initial radius equal to the minimum of the potential well. We seen then that the radius is always the same and corresponds to a circle (the radius is always constant).</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1"># Common imports</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">import</span> <span class="nn">pandas</span> <span class="k">as</span> <span class="nn">pd</span>
<span class="kn">from</span> <span class="nn">math</span> <span class="kn">import</span> <span class="o">*</span>
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="kn">import</span> <span class="nn">os</span>
<span class="c1"># Where to save the figures and data files</span>
<span class="n">PROJECT_ROOT_DIR</span> <span class="o">=</span> <span class="s2">&quot;Results&quot;</span>
<span class="n">FIGURE_ID</span> <span class="o">=</span> <span class="s2">&quot;Results/FigureFiles&quot;</span>
<span class="n">DATA_ID</span> <span class="o">=</span> <span class="s2">&quot;DataFiles/&quot;</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">os</span><span class="o">.</span><span class="n">path</span><span class="o">.</span><span class="n">exists</span><span class="p">(</span><span class="n">PROJECT_ROOT_DIR</span><span class="p">):</span>
<span class="n">os</span><span class="o">.</span><span class="n">mkdir</span><span class="p">(</span><span class="n">PROJECT_ROOT_DIR</span><span class="p">)</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">os</span><span class="o">.</span><span class="n">path</span><span class="o">.</span><span class="n">exists</span><span class="p">(</span><span class="n">FIGURE_ID</span><span class="p">):</span>
<span class="n">os</span><span class="o">.</span><span class="n">makedirs</span><span class="p">(</span><span class="n">FIGURE_ID</span><span class="p">)</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">os</span><span class="o">.</span><span class="n">path</span><span class="o">.</span><span class="n">exists</span><span class="p">(</span><span class="n">DATA_ID</span><span class="p">):</span>
<span class="n">os</span><span class="o">.</span><span class="n">makedirs</span><span class="p">(</span><span class="n">DATA_ID</span><span class="p">)</span>
<span class="k">def</span> <span class="nf">image_path</span><span class="p">(</span><span class="n">fig_id</span><span class="p">):</span>
<span class="k">return</span> <span class="n">os</span><span class="o">.</span><span class="n">path</span><span class="o">.</span><span class="n">join</span><span class="p">(</span><span class="n">FIGURE_ID</span><span class="p">,</span> <span class="n">fig_id</span><span class="p">)</span>
<span class="k">def</span> <span class="nf">data_path</span><span class="p">(</span><span class="n">dat_id</span><span class="p">):</span>
<span class="k">return</span> <span class="n">os</span><span class="o">.</span><span class="n">path</span><span class="o">.</span><span class="n">join</span><span class="p">(</span><span class="n">DATA_ID</span><span class="p">,</span> <span class="n">dat_id</span><span class="p">)</span>
<span class="k">def</span> <span class="nf">save_fig</span><span class="p">(</span><span class="n">fig_id</span><span class="p">):</span>
<span class="n">plt</span><span class="o">.</span><span class="n">savefig</span><span class="p">(</span><span class="n">image_path</span><span class="p">(</span><span class="n">fig_id</span><span class="p">)</span> <span class="o">+</span> <span class="s2">&quot;.png&quot;</span><span class="p">,</span> <span class="nb">format</span><span class="o">=</span><span class="s1">&#39;png&#39;</span><span class="p">)</span>
<span class="c1"># Simple Gravitational Force -alpha/r</span>
<span class="n">DeltaT</span> <span class="o">=</span> <span class="mf">0.01</span>
<span class="c1">#set up arrays </span>
<span class="n">tfinal</span> <span class="o">=</span> <span class="mf">100.0</span>
<span class="n">n</span> <span class="o">=</span> <span class="n">ceil</span><span class="p">(</span><span class="n">tfinal</span><span class="o">/</span><span class="n">DeltaT</span><span class="p">)</span>
<span class="c1"># set up arrays for t, v and r</span>
<span class="n">t</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">v</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">r</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="c1"># Constants of the model, setting all variables to one for simplicity</span>
<span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="n">AngMom</span> <span class="o">=</span> <span class="mf">1.0</span> <span class="c1"># The angular momentum</span>
<span class="n">m</span> <span class="o">=</span> <span class="mf">1.0</span> <span class="c1"># scale mass to one</span>
<span class="n">c1</span> <span class="o">=</span> <span class="n">AngMom</span><span class="o">*</span><span class="n">AngMom</span><span class="o">/</span><span class="p">(</span><span class="n">m</span><span class="o">*</span><span class="n">m</span><span class="p">)</span>
<span class="n">c2</span> <span class="o">=</span> <span class="n">AngMom</span><span class="o">*</span><span class="n">AngMom</span><span class="o">/</span><span class="n">m</span>
<span class="n">rmin</span> <span class="o">=</span> <span class="p">(</span><span class="n">AngMom</span><span class="o">*</span><span class="n">AngMom</span><span class="o">/</span><span class="n">m</span><span class="o">/</span><span class="n">alpha</span><span class="p">)</span>
<span class="c1"># Initial conditions</span>
<span class="n">r0</span> <span class="o">=</span> <span class="n">rmin</span>
<span class="n">v0</span> <span class="o">=</span> <span class="mf">0.0</span>
<span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">r0</span>
<span class="n">v</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">v0</span>
<span class="c1"># Start integrating using the Velocity-Verlet method</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
<span class="c1"># Set up acceleration</span>
<span class="n">a</span> <span class="o">=</span> <span class="o">-</span><span class="n">alpha</span><span class="o">/</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">+</span><span class="n">c1</span><span class="o">/</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">**</span><span class="mi">3</span><span class="p">)</span>
<span class="c1"># update velocity, time and position using the Velocity-Verlet method</span>
<span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="n">DeltaT</span><span class="o">*</span><span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="p">(</span><span class="n">DeltaT</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">*</span><span class="n">a</span>
<span class="n">anew</span> <span class="o">=</span> <span class="o">-</span><span class="n">alpha</span><span class="o">/</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">+</span><span class="n">c1</span><span class="o">/</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">**</span><span class="mi">3</span><span class="p">)</span>
<span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">DeltaT</span><span class="o">*</span><span class="p">(</span><span class="n">a</span><span class="o">+</span><span class="n">anew</span><span class="p">)</span>
<span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="n">DeltaT</span>
<span class="c1"># Plot position as function of time</span>
<span class="n">fig</span><span class="p">,</span> <span class="n">ax</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">subplots</span><span class="p">(</span><span class="mi">2</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;time&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;radius&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span><span class="n">r</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;time&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;Velocity&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span><span class="n">v</span><span class="p">)</span>
<span class="n">save_fig</span><span class="p">(</span><span class="s2">&quot;RadialGVV&quot;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<img alt="../_images/chapter6_47_01.png" src="../_images/chapter6_47_01.png" />
</div>
</div>
<p>Changing the value of the initial position to a value where the energy is positive, leads to an increasing radius with time, a so-called unbound orbit. Choosing on the other hand an initial radius that corresponds to a negative energy and different from the minimum value leads to a radius that oscillates back and forth between two values.</p>
</div>
<div class="section" id="harmonic-oscillator-in-two-dimensions">
<h3>Harmonic Oscillator in two dimensions<a class="headerlink" href="#harmonic-oscillator-in-two-dimensions" title="Permalink to this headline"></a></h3>
<p>Consider a particle of mass <span class="math notranslate nohighlight">\(m\)</span> in a 2-dimensional harmonic oscillator with potential</p>
<div class="math notranslate nohighlight">
\[
U=\frac{1}{2}kr^2=\frac{1}{2}k(x^2+y^2).
\]</div>
<p>If the orbit has angular momentum <span class="math notranslate nohighlight">\(L\)</span>, we can find the radius and angular velocity of the circular orbit as well as the b) the angular frequency of small radial perturbations.</p>
<p>We consider the effective potential. The radius of a circular orbit is at the minimum of the potential (where the effective force is zero).
The potential is plotted here with the parameters <span class="math notranslate nohighlight">\(k=m=0.1\)</span> and <span class="math notranslate nohighlight">\(L=1.0\)</span>.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1"># Common imports</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">from</span> <span class="nn">math</span> <span class="kn">import</span> <span class="o">*</span>
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="n">Deltax</span> <span class="o">=</span> <span class="mf">0.01</span>
<span class="c1">#set up arrays</span>
<span class="n">xinitial</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="n">xfinal</span> <span class="o">=</span> <span class="mf">5.0</span>
<span class="n">k</span> <span class="o">=</span> <span class="mf">0.1</span> <span class="c1"># spring constant</span>
<span class="n">m</span> <span class="o">=</span> <span class="mf">0.1</span> <span class="c1"># mass, you can change these</span>
<span class="n">AngMom</span> <span class="o">=</span> <span class="mf">1.0</span> <span class="c1"># The angular momentum</span>
<span class="n">n</span> <span class="o">=</span> <span class="n">ceil</span><span class="p">((</span><span class="n">xfinal</span><span class="o">-</span><span class="n">xinitial</span><span class="p">)</span><span class="o">/</span><span class="n">Deltax</span><span class="p">)</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
<span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">xinitial</span><span class="o">+</span><span class="n">i</span><span class="o">*</span><span class="n">Deltax</span>
<span class="n">V</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">V</span> <span class="o">=</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">k</span><span class="o">*</span><span class="n">x</span><span class="o">*</span><span class="n">x</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="n">AngMom</span><span class="o">*</span><span class="n">AngMom</span><span class="o">/</span><span class="p">(</span><span class="n">m</span><span class="o">*</span><span class="n">x</span><span class="o">*</span><span class="n">x</span><span class="p">)</span>
<span class="c1"># Plot potential</span>
<span class="n">fig</span><span class="p">,</span> <span class="n">ax</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">subplots</span><span class="p">()</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;r[m]&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;V[J]&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">V</span><span class="p">)</span>
<span class="n">fig</span><span class="o">.</span><span class="n">tight_layout</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<img alt="../_images/chapter6_51_01.png" src="../_images/chapter6_51_01.png" />
</div>
</div>
<div class="math notranslate nohighlight">
\[
\begin{eqnarray*}
U_{\rm eff}&amp;=&amp;\frac{1}{2}kr^2+\frac{L^2}{2mr^2}
\end{eqnarray*}
\]</div>
<p>The effective potential looks like that of a harmonic oscillator for
large <span class="math notranslate nohighlight">\(r\)</span>, but for small <span class="math notranslate nohighlight">\(r\)</span>, the centrifugal potential repels the
particle from the origin. The combination of the two potentials has a
minimum for at some radius <span class="math notranslate nohighlight">\(r_{\rm min}\)</span>.</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray*}
0&amp;=&amp;kr_{\rm min}-\frac{L^2}{mr_{\rm min}^3},\\
r_{\rm min}&amp;=&amp;\left(\frac{L^2}{mk}\right)^{1/4},\\
\dot{\theta}&amp;=&amp;\frac{L}{mr_{\rm min}^2}=\sqrt{k/m}.
\end{eqnarray*}
\end{split}\]</div>
<p>For particles at <span class="math notranslate nohighlight">\(r_{\rm min}\)</span> with <span class="math notranslate nohighlight">\(\dot{r}=0\)</span>, the particle does not
accelerate and <span class="math notranslate nohighlight">\(r\)</span> stays constant, i.e. a circular orbit. The radius
of the circular orbit can be adjusted by changing the angular momentum
<span class="math notranslate nohighlight">\(L\)</span>.</p>
<p>For the above parameters this minimum is at <span class="math notranslate nohighlight">\(r_{\rm min}=1\)</span>.</p>
<p>Now consider small vibrations about <span class="math notranslate nohighlight">\(r_{\rm min}\)</span>. The effective spring constant is the curvature of the effective potential.</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray*}
k_{\rm eff}&amp;=&amp;\left.\frac{d^2}{dr^2}U_{\rm eff}(r)\right|_{r=r_{\rm min}}=k+\frac{3L^2}{mr_{\rm min}^4}\\
&amp;=&amp;4k,\\
\omega&amp;=&amp;\sqrt{k_{\rm eff}/m}=2\sqrt{k/m}=2\dot{\theta}.
\end{eqnarray*}
\end{split}\]</div>
<p>Here, the second step used the result of the last step from part
(a). Because the radius oscillates with twice the angular frequency,
the orbit has two places where <span class="math notranslate nohighlight">\(r\)</span> reaches a minimum in one
cycle. This differs from the inverse-square force where there is one
minimum in an orbit. One can show that the orbit for the harmonic
oscillator is also elliptical, but in this case the center of the
potential is at the center of the ellipse, not at one of the foci.</p>
<p>The solution is also simple to write down exactly in Cartesian coordinates. The <span class="math notranslate nohighlight">\(x\)</span> and <span class="math notranslate nohighlight">\(y\)</span> equations of motion separate,</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray*}
\ddot{x}&amp;=&amp;-kx,\\
\ddot{y}&amp;=&amp;-ky.
\end{eqnarray*}
\end{split}\]</div>
<p>So the general solution can be expressed as</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray*}
x&amp;=&amp;A\cos\omega_0 t+B\sin\omega_0 t,\\
y&amp;=&amp;C\cos\omega_0 t+D\sin\omega_0 t.
\end{eqnarray*}
\end{split}\]</div>
<p>The code here finds the solution for <span class="math notranslate nohighlight">\(x\)</span> and <span class="math notranslate nohighlight">\(y\)</span> using the code we developed in homework 4.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="n">DeltaT</span> <span class="o">=</span> <span class="mf">0.01</span>
<span class="c1">#set up arrays </span>
<span class="n">tfinal</span> <span class="o">=</span> <span class="mf">10.0</span>
<span class="n">n</span> <span class="o">=</span> <span class="n">ceil</span><span class="p">(</span><span class="n">tfinal</span><span class="o">/</span><span class="n">DeltaT</span><span class="p">)</span>
<span class="c1"># set up arrays</span>
<span class="n">t</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">v</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">n</span><span class="p">,</span><span class="mi">2</span><span class="p">))</span>
<span class="n">r</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">n</span><span class="p">,</span><span class="mi">2</span><span class="p">))</span>
<span class="n">radius</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="c1"># Constants of the model</span>
<span class="n">k</span> <span class="o">=</span> <span class="mf">0.1</span> <span class="c1"># spring constant</span>
<span class="n">m</span> <span class="o">=</span> <span class="mf">0.1</span> <span class="c1"># mass, you can change these</span>
<span class="n">omega02</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">k</span><span class="o">/</span><span class="n">m</span><span class="p">)</span> <span class="c1"># Frequency</span>
<span class="n">AngMom</span> <span class="o">=</span> <span class="mf">1.0</span> <span class="c1"># The angular momentum</span>
<span class="n">rmin</span> <span class="o">=</span> <span class="p">(</span><span class="n">AngMom</span><span class="o">*</span><span class="n">AngMom</span><span class="o">/</span><span class="n">k</span><span class="o">/</span><span class="n">m</span><span class="p">)</span><span class="o">**</span><span class="mf">0.25</span>
<span class="c1"># Initial conditions as compact 2-dimensional arrays</span>
<span class="c1">#x0 =rmin*0.5; y0 = sqrt(rmin*rmin-x0*x0)</span>
<span class="n">x0</span> <span class="o">=</span> <span class="mf">1.0</span><span class="p">;</span> <span class="n">y0</span><span class="o">=</span> <span class="mf">1.0</span>
<span class="n">r0</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="n">x0</span><span class="p">,</span><span class="n">y0</span><span class="p">])</span>
<span class="n">v0</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mf">0.0</span><span class="p">,</span><span class="mf">0.0</span><span class="p">])</span>
<span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">r0</span>
<span class="n">v</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">v0</span>
<span class="c1"># Start integrating using the Velocity-Verlet method</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
<span class="c1"># Set up the acceleration</span>
<span class="n">a</span> <span class="o">=</span> <span class="o">-</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">*</span><span class="n">omega02</span>
<span class="c1"># update velocity, time and position using the Velocity-Verlet method</span>
<span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="n">DeltaT</span><span class="o">*</span><span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="p">(</span><span class="n">DeltaT</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">*</span><span class="n">a</span>
<span class="n">anew</span> <span class="o">=</span> <span class="o">-</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">*</span><span class="n">omega02</span>
<span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">DeltaT</span><span class="o">*</span><span class="p">(</span><span class="n">a</span><span class="o">+</span><span class="n">anew</span><span class="p">)</span>
<span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="n">DeltaT</span>
<span class="c1"># Plot position as function of time</span>
<span class="n">radius</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">r</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="o">+</span><span class="n">r</span><span class="p">[:,</span><span class="mi">1</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
<span class="n">fig</span><span class="p">,</span> <span class="n">ax</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">subplots</span><span class="p">(</span><span class="mi">3</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;time&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;radius squared&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span><span class="n">r</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="o">+</span><span class="n">r</span><span class="p">[:,</span><span class="mi">1</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;time&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;x position&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span><span class="n">r</span><span class="p">[:,</span><span class="mi">0</span><span class="p">])</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;time&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;y position&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span><span class="n">r</span><span class="p">[:,</span><span class="mi">1</span><span class="p">])</span>
<span class="n">fig</span><span class="o">.</span><span class="n">tight_layout</span><span class="p">()</span>
<span class="n">save_fig</span><span class="p">(</span><span class="s2">&quot;2DimHOVV&quot;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<img alt="../_images/chapter6_62_01.png" src="../_images/chapter6_62_01.png" />
</div>
</div>
<p>With some work using double angle formulas, one can calculate</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray*}
r^2&amp;=&amp;x^2+y^2\\
\nonumber
&amp;=&amp;(A^2+C^2)\cos^2(\omega_0t)+(B^2+D^2)\sin^2\omega_0t+(AB+CD)\cos(\omega_0t)\sin(\omega_0t)\\
\nonumber
&amp;=&amp;\alpha+\beta\cos 2\omega_0 t+\gamma\sin 2\omega_0 t,\\
\alpha&amp;=&amp;\frac{A^2+B^2+C^2+D^2}{2},~~\beta=\frac{A^2-B^2+C^2-D^2}{2},~~\gamma=AB+CD,\\
r^2&amp;=&amp;\alpha+(\beta^2+\gamma^2)^{1/2}\cos(2\omega_0 t-\delta),~~~\delta=\arctan(\gamma/\beta),
\end{eqnarray*}
\end{split}\]</div>
<p>and see that radius oscillates with frequency <span class="math notranslate nohighlight">\(2\omega_0\)</span>. The
factor of two comes because the oscillation <span class="math notranslate nohighlight">\(x=A\cos\omega_0t\)</span> has two
maxima for <span class="math notranslate nohighlight">\(x^2\)</span>, one at <span class="math notranslate nohighlight">\(t=0\)</span> and one a half period later.</p>
<p>The following code shows first how we can solve this problem using the radial degrees of freedom only.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="n">DeltaT</span> <span class="o">=</span> <span class="mf">0.01</span>
<span class="c1">#set up arrays </span>
<span class="n">tfinal</span> <span class="o">=</span> <span class="mf">10.0</span>
<span class="n">n</span> <span class="o">=</span> <span class="n">ceil</span><span class="p">(</span><span class="n">tfinal</span><span class="o">/</span><span class="n">DeltaT</span><span class="p">)</span>
<span class="c1"># set up arrays for t, v and r</span>
<span class="n">t</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">v</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">r</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">E</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="c1"># Constants of the model</span>
<span class="n">AngMom</span> <span class="o">=</span> <span class="mf">1.0</span> <span class="c1"># The angular momentum</span>
<span class="n">m</span> <span class="o">=</span> <span class="mf">0.1</span>
<span class="n">k</span> <span class="o">=</span> <span class="mf">0.1</span>
<span class="n">omega02</span> <span class="o">=</span> <span class="n">k</span><span class="o">/</span><span class="n">m</span>
<span class="n">c1</span> <span class="o">=</span> <span class="n">AngMom</span><span class="o">*</span><span class="n">AngMom</span><span class="o">/</span><span class="p">(</span><span class="n">m</span><span class="o">*</span><span class="n">m</span><span class="p">)</span>
<span class="n">c2</span> <span class="o">=</span> <span class="n">AngMom</span><span class="o">*</span><span class="n">AngMom</span><span class="o">/</span><span class="n">m</span>
<span class="n">rmin</span> <span class="o">=</span> <span class="p">(</span><span class="n">AngMom</span><span class="o">*</span><span class="n">AngMom</span><span class="o">/</span><span class="n">k</span><span class="o">/</span><span class="n">m</span><span class="p">)</span><span class="o">**</span><span class="mf">0.25</span>
<span class="c1"># Initial conditions</span>
<span class="n">r0</span> <span class="o">=</span> <span class="n">rmin</span>
<span class="n">v0</span> <span class="o">=</span> <span class="mf">0.0</span>
<span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">r0</span>
<span class="n">v</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">v0</span>
<span class="n">E</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">m</span><span class="o">*</span><span class="n">v0</span><span class="o">*</span><span class="n">v0</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="n">k</span><span class="o">*</span><span class="n">r0</span><span class="o">*</span><span class="n">r0</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="n">c2</span><span class="o">/</span><span class="p">(</span><span class="n">r0</span><span class="o">*</span><span class="n">r0</span><span class="p">)</span>
<span class="c1"># Start integrating using the Velocity-Verlet method</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
<span class="c1"># Set up acceleration</span>
<span class="n">a</span> <span class="o">=</span> <span class="o">-</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">*</span><span class="n">omega02</span><span class="o">+</span><span class="n">c1</span><span class="o">/</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">**</span><span class="mi">3</span><span class="p">)</span>
<span class="c1"># update velocity, time and position using the Velocity-Verlet method</span>
<span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="n">DeltaT</span><span class="o">*</span><span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="p">(</span><span class="n">DeltaT</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">*</span><span class="n">a</span>
<span class="n">anew</span> <span class="o">=</span> <span class="o">-</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">*</span><span class="n">omega02</span><span class="o">+</span><span class="n">c1</span><span class="o">/</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">**</span><span class="mi">3</span><span class="p">)</span>
<span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">DeltaT</span><span class="o">*</span><span class="p">(</span><span class="n">a</span><span class="o">+</span><span class="n">anew</span><span class="p">)</span>
<span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="n">DeltaT</span>
<span class="n">E</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">m</span><span class="o">*</span><span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">*</span><span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="n">k</span><span class="o">*</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">*</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="n">c2</span><span class="o">/</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">*</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">])</span>
<span class="c1"># Plot position as function of time</span>
<span class="n">fig</span><span class="p">,</span> <span class="n">ax</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">subplots</span><span class="p">(</span><span class="mi">2</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;time&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;radius&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span><span class="n">r</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;time&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;Energy&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span><span class="n">E</span><span class="p">)</span>
<span class="n">save_fig</span><span class="p">(</span><span class="s2">&quot;RadialHOVV&quot;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<img alt="../_images/chapter6_66_01.png" src="../_images/chapter6_66_01.png" />
</div>
</div>
</div>
</div>
<div class="section" id="stability-of-orbits">
<h2>Stability of Orbits<a class="headerlink" href="#stability-of-orbits" title="Permalink to this headline"></a></h2>
<p>The effective force can be extracted from the effective potential, <span class="math notranslate nohighlight">\(U_{\rm eff}\)</span>. Beginning from the equations of motion, Eq. (<a class="reference external" href="#eq:radialeqofmotion">2</a>), for <span class="math notranslate nohighlight">\(r\)</span>,</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
m\ddot{r}&amp;=&amp;F+\frac{L^2}{mr^3}\\
\nonumber
&amp;=&amp;F_{\rm eff}\\
\nonumber
&amp;=&amp;-\partial_rU_{\rm eff},\\
\nonumber
F_{\rm eff}&amp;=&amp;-\partial_r\left[U(r)+(L^2/2mr^2)\right].
\end{eqnarray}
\end{split}\]</div>
<p>For a circular orbit, the radius must be fixed as a function of time,
so one must be at a maximum or a minimum of the effective
potential. However, if one is at a maximum of the effective potential
the radius will be unstable. For the attractive Coulomb force the
effective potential will be dominated by the <span class="math notranslate nohighlight">\(-\alpha/r\)</span> term for
large <span class="math notranslate nohighlight">\(r\)</span> because the centrifugal part falls off more quickly, <span class="math notranslate nohighlight">\(\sim
1/r^2\)</span>. At low <span class="math notranslate nohighlight">\(r\)</span> the centrifugal piece wins and the effective
potential is repulsive. Thus, the potential must have a minimum
somewhere with negative potential. The circular orbits are then stable
to perturbation.</p>
<p>The effective potential is sketched for two cases, a <span class="math notranslate nohighlight">\(1/r\)</span> attractive
potential and a <span class="math notranslate nohighlight">\(1/r^3\)</span> attractive potential. The <span class="math notranslate nohighlight">\(1/r\)</span> case has a
stable minimum, whereas the circular orbit in the <span class="math notranslate nohighlight">\(1/r^3\)</span> case is
unstable.</p>
<p>If one considers a potential that falls as <span class="math notranslate nohighlight">\(1/r^3\)</span>, the situation is
reversed and the point where <span class="math notranslate nohighlight">\(\partial_rU\)</span> disappears will be a local
maximum rather than a local minimum. <strong>Fig to come here with code</strong></p>
<p>The repulsive centrifugal piece dominates at large <span class="math notranslate nohighlight">\(r\)</span> and the attractive
Coulomb piece wins out at small <span class="math notranslate nohighlight">\(r\)</span>. The circular orbit is then at a
maximum of the effective potential and the orbits are unstable. It is
the clear that for potentials that fall as <span class="math notranslate nohighlight">\(r^n\)</span>, that one must have
<span class="math notranslate nohighlight">\(n&gt;-2\)</span> for the orbits to be stable.</p>
<p>Consider a potential <span class="math notranslate nohighlight">\(U(r)=\beta r\)</span>. For a particle of mass <span class="math notranslate nohighlight">\(m\)</span> with
angular momentum <span class="math notranslate nohighlight">\(L\)</span>, find the angular frequency of a circular
orbit. Then find the angular frequency for small radial perturbations.</p>
<p>For the circular orbit you search for the position <span class="math notranslate nohighlight">\(r_{\rm min}\)</span> where the effective potential is minimized,</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray*}
\partial_r\left\{\beta r+\frac{L^2}{2mr^2}\right\}&amp;=&amp;0,\\
\beta&amp;=&amp;\frac{L^2}{mr_{\rm min}^3},\\
r_{\rm min}&amp;=&amp;\left(\frac{L^2}{\beta m}\right)^{1/3},\\
\dot{\theta}&amp;=&amp;\frac{L}{mr_{\rm min}^2}=\frac{\beta^{2/3}}{(mL)^{1/3}}
\end{eqnarray*}
\end{split}\]</div>
<p>Now, we can find the angular frequency of small perturbations about the circular orbit. To do this we find the effective spring constant for the effective potential,</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray*}
k_{\rm eff}&amp;=&amp;\partial_r^2 \left.U_{\rm eff}\right|_{r_{\rm min}}\\
&amp;=&amp;\frac{3L^2}{mr_{\rm min}^4},\\
\omega&amp;=&amp;\sqrt{\frac{k_{\rm eff}}{m}}\\
&amp;=&amp;\frac{\beta^{2/3}}{(mL)^{1/3}}\sqrt{3}.
\end{eqnarray*}
\end{split}\]</div>
<p>If the two frequencies, <span class="math notranslate nohighlight">\(\dot{\theta}\)</span> and <span class="math notranslate nohighlight">\(\omega\)</span>, differ by an
integer factor, the orbit’s trajectory will repeat itself each time
around. This is the case for the inverse-square force,
<span class="math notranslate nohighlight">\(\omega=\dot{\theta}\)</span>, and for the harmonic oscillator,
<span class="math notranslate nohighlight">\(\omega=2\dot{\theta}\)</span>. In this case, <span class="math notranslate nohighlight">\(\omega=\sqrt{3}\dot{\theta}\)</span>,
and the angles at which the maxima and minima occur change with each
orbit.</p>
<div class="section" id="code-example-with-gravitional-force">
<h3>Code example with gravitional force<a class="headerlink" href="#code-example-with-gravitional-force" title="Permalink to this headline"></a></h3>
<p>The code example here is meant to illustrate how we can make a plot of the final orbit. We solve the equations in polar coordinates (the example here uses the minimum of the potential as initial value) and then we transform back to cartesian coordinates and plot <span class="math notranslate nohighlight">\(x\)</span> versus <span class="math notranslate nohighlight">\(y\)</span>. We see that we get a perfect circle when we place ourselves at the minimum of the potential energy, as expected.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1"># Simple Gravitational Force -alpha/r</span>
<span class="n">DeltaT</span> <span class="o">=</span> <span class="mf">0.01</span>
<span class="c1">#set up arrays </span>
<span class="n">tfinal</span> <span class="o">=</span> <span class="mf">8.0</span>
<span class="n">n</span> <span class="o">=</span> <span class="n">ceil</span><span class="p">(</span><span class="n">tfinal</span><span class="o">/</span><span class="n">DeltaT</span><span class="p">)</span>
<span class="c1"># set up arrays for t, v and r</span>
<span class="n">t</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">v</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">r</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">phi</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="c1"># Constants of the model, setting all variables to one for simplicity</span>
<span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="n">AngMom</span> <span class="o">=</span> <span class="mf">1.0</span> <span class="c1"># The angular momentum</span>
<span class="n">m</span> <span class="o">=</span> <span class="mf">1.0</span> <span class="c1"># scale mass to one</span>
<span class="n">c1</span> <span class="o">=</span> <span class="n">AngMom</span><span class="o">*</span><span class="n">AngMom</span><span class="o">/</span><span class="p">(</span><span class="n">m</span><span class="o">*</span><span class="n">m</span><span class="p">)</span>
<span class="n">c2</span> <span class="o">=</span> <span class="n">AngMom</span><span class="o">*</span><span class="n">AngMom</span><span class="o">/</span><span class="n">m</span>
<span class="n">rmin</span> <span class="o">=</span> <span class="p">(</span><span class="n">AngMom</span><span class="o">*</span><span class="n">AngMom</span><span class="o">/</span><span class="n">m</span><span class="o">/</span><span class="n">alpha</span><span class="p">)</span>
<span class="c1"># Initial conditions, place yourself at the potential min</span>
<span class="n">r0</span> <span class="o">=</span> <span class="n">rmin</span>
<span class="n">v0</span> <span class="o">=</span> <span class="mf">0.0</span> <span class="c1"># starts at rest</span>
<span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">r0</span>
<span class="n">v</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">v0</span>
<span class="n">phi</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.0</span>
<span class="c1"># Start integrating using the Velocity-Verlet method</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
<span class="c1"># Set up acceleration</span>
<span class="n">a</span> <span class="o">=</span> <span class="o">-</span><span class="n">alpha</span><span class="o">/</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">+</span><span class="n">c1</span><span class="o">/</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">**</span><span class="mi">3</span><span class="p">)</span>
<span class="c1"># update velocity, time and position using the Velocity-Verlet method</span>
<span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="n">DeltaT</span><span class="o">*</span><span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="p">(</span><span class="n">DeltaT</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">*</span><span class="n">a</span>
<span class="n">anew</span> <span class="o">=</span> <span class="o">-</span><span class="n">alpha</span><span class="o">/</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">+</span><span class="n">c1</span><span class="o">/</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">**</span><span class="mi">3</span><span class="p">)</span>
<span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">DeltaT</span><span class="o">*</span><span class="p">(</span><span class="n">a</span><span class="o">+</span><span class="n">anew</span><span class="p">)</span>
<span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="n">DeltaT</span>
<span class="n">phi</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">*</span><span class="n">c2</span><span class="o">/</span><span class="p">(</span><span class="n">r0</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
<span class="c1"># Find cartesian coordinates for easy plot </span>
<span class="n">x</span> <span class="o">=</span> <span class="n">r</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">phi</span><span class="p">)</span>
<span class="n">y</span> <span class="o">=</span> <span class="n">r</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">phi</span><span class="p">)</span>
<span class="n">fig</span><span class="p">,</span> <span class="n">ax</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">subplots</span><span class="p">(</span><span class="mi">3</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;time&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;radius&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span><span class="n">r</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;time&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;Angle $\cos{\phi}$&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span><span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">phi</span><span class="p">))</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;y&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;x&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span><span class="n">y</span><span class="p">)</span>
<span class="n">save_fig</span><span class="p">(</span><span class="s2">&quot;Phasespace&quot;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<img alt="../_images/chapter6_74_01.png" src="../_images/chapter6_74_01.png" />
</div>
</div>
<p>Try to change the initial value for <span class="math notranslate nohighlight">\(r\)</span> and see what kind of orbits you get.
In order to test different energies, it can be useful to look at the plot of the effective potential discussed above.</p>
<p>However, for orbits different from a circle the above code would need modifications in order to allow us to display say an ellipse. For the latter, it is much easier to run our code in cartesian coordinates, as done here. In this code we test also energy conservation and see that it is conserved to numerical precision. The code here is a simple extension of the code we developed for homework 4.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1"># Common imports</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">import</span> <span class="nn">pandas</span> <span class="k">as</span> <span class="nn">pd</span>
<span class="kn">from</span> <span class="nn">math</span> <span class="kn">import</span> <span class="o">*</span>
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="n">DeltaT</span> <span class="o">=</span> <span class="mf">0.01</span>
<span class="c1">#set up arrays </span>
<span class="n">tfinal</span> <span class="o">=</span> <span class="mf">10.0</span>
<span class="n">n</span> <span class="o">=</span> <span class="n">ceil</span><span class="p">(</span><span class="n">tfinal</span><span class="o">/</span><span class="n">DeltaT</span><span class="p">)</span>
<span class="c1"># set up arrays</span>
<span class="n">t</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">v</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">n</span><span class="p">,</span><span class="mi">2</span><span class="p">))</span>
<span class="n">r</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">n</span><span class="p">,</span><span class="mi">2</span><span class="p">))</span>
<span class="n">E</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="c1"># Constants of the model</span>
<span class="n">m</span> <span class="o">=</span> <span class="mf">1.0</span> <span class="c1"># mass, you can change these</span>
<span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="c1"># Initial conditions as compact 2-dimensional arrays</span>
<span class="n">x0</span> <span class="o">=</span> <span class="mf">0.5</span><span class="p">;</span> <span class="n">y0</span><span class="o">=</span> <span class="mf">0.</span>
<span class="n">r0</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="n">x0</span><span class="p">,</span><span class="n">y0</span><span class="p">])</span>
<span class="n">v0</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mf">0.0</span><span class="p">,</span><span class="mf">1.0</span><span class="p">])</span>
<span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">r0</span>
<span class="n">v</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">v0</span>
<span class="n">rabs</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="nb">sum</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">*</span><span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">]))</span>
<span class="n">E</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">m</span><span class="o">*</span><span class="p">(</span><span class="n">v</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="o">+</span><span class="n">v</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">-</span><span class="n">alpha</span><span class="o">/</span><span class="n">rabs</span>
<span class="c1"># Start integrating using the Velocity-Verlet method</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
<span class="c1"># Set up the acceleration</span>
<span class="n">rabs</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="nb">sum</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">*</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">]))</span>
<span class="n">a</span> <span class="o">=</span> <span class="o">-</span><span class="n">alpha</span><span class="o">*</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">/</span><span class="p">(</span><span class="n">rabs</span><span class="o">**</span><span class="mi">3</span><span class="p">)</span>
<span class="c1"># update velocity, time and position using the Velocity-Verlet method</span>
<span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="n">DeltaT</span><span class="o">*</span><span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="p">(</span><span class="n">DeltaT</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">*</span><span class="n">a</span>
<span class="n">rabs</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="nb">sum</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">*</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]))</span>
<span class="n">anew</span> <span class="o">=</span> <span class="o">-</span><span class="n">alpha</span><span class="o">*</span><span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">/</span><span class="p">(</span><span class="n">rabs</span><span class="o">**</span><span class="mi">3</span><span class="p">)</span>
<span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">DeltaT</span><span class="o">*</span><span class="p">(</span><span class="n">a</span><span class="o">+</span><span class="n">anew</span><span class="p">)</span>
<span class="n">E</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">m</span><span class="o">*</span><span class="p">(</span><span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="o">+</span><span class="n">v</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">-</span><span class="n">alpha</span><span class="o">/</span><span class="n">rabs</span>
<span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="n">DeltaT</span>
<span class="c1"># Plot position as function of time</span>
<span class="n">fig</span><span class="p">,</span> <span class="n">ax</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">subplots</span><span class="p">(</span><span class="mi">3</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;y&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;x&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">r</span><span class="p">[:,</span><span class="mi">0</span><span class="p">],</span><span class="n">r</span><span class="p">[:,</span><span class="mi">1</span><span class="p">])</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;time&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;y position&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span><span class="n">r</span><span class="p">[:,</span><span class="mi">0</span><span class="p">])</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;time&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;y position&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span><span class="n">r</span><span class="p">[:,</span><span class="mi">1</span><span class="p">])</span>
<span class="n">fig</span><span class="o">.</span><span class="n">tight_layout</span><span class="p">()</span>
<span class="n">save_fig</span><span class="p">(</span><span class="s2">&quot;2DimGravity&quot;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
<span class="nb">print</span><span class="p">(</span><span class="n">E</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<img alt="../_images/chapter6_76_01.png" src="../_images/chapter6_76_01.png" />
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[-1.5 -1.49999996 -1.49999984 -1.49999964 -1.49999935 -1.49999898
-1.49999853 -1.49999798 -1.49999734 -1.49999659 -1.49999574 -1.49999478
-1.49999369 -1.49999247 -1.49999112 -1.49998961 -1.49998794 -1.49998609
-1.49998404 -1.49998178 -1.49997929 -1.49997654 -1.49997351 -1.49997017
-1.49996649 -1.49996243 -1.49995796 -1.49995304 -1.49994762 -1.49994165
-1.49993508 -1.49992786 -1.49991994 -1.49991127 -1.4999018 -1.49989151
-1.49988041 -1.49986852 -1.49985595 -1.49984292 -1.49982978 -1.49981712
-1.49980588 -1.49979754 -1.49979431 -1.49979958 -1.4998184 -1.4998583
-1.49993028 -1.50005034 -1.50024115 -1.500534 -1.50097003 -1.50159939
-1.50247522 -1.50363822 -1.5050881 -1.50674417 -1.50841251 -1.50979454
-1.51056842 -1.5105255 -1.50967788 -1.5082518 -1.50657253 -1.50493047
-1.50350752 -1.50237442 -1.50152565 -1.50091821 -1.50049876 -1.50021789
-1.50003547 -1.49992116 -1.49985304 -1.4998157 -1.49979853 -1.49979431
-1.49979819 -1.49980692 -1.49981835 -1.49983109 -1.49984425 -1.49985725
-1.49986975 -1.49988156 -1.49989259 -1.49990279 -1.49991218 -1.49992077
-1.49992862 -1.49993577 -1.49994227 -1.49994819 -1.49995356 -1.49995843
-1.49996286 -1.49996688 -1.49997052 -1.49997383 -1.49997683 -1.49997955
-1.49998202 -1.49998426 -1.49998628 -1.49998812 -1.49998977 -1.49999126
-1.4999926 -1.49999381 -1.49999488 -1.49999583 -1.49999667 -1.49999741
-1.49999804 -1.49999858 -1.49999903 -1.49999939 -1.49999966 -1.49999986
-1.49999997 -1.5 -1.49999995 -1.49999982 -1.49999961 -1.49999932
-1.49999894 -1.49999848 -1.49999792 -1.49999727 -1.49999651 -1.49999565
-1.49999467 -1.49999357 -1.49999234 -1.49999097 -1.49998945 -1.49998776
-1.49998589 -1.49998382 -1.49998154 -1.49997902 -1.49997625 -1.49997319
-1.49996981 -1.4999661 -1.499962 -1.49995749 -1.49995252 -1.49994704
-1.49994101 -1.49993438 -1.49992709 -1.4999191 -1.49991035 -1.4999008
-1.49989043 -1.49987924 -1.49986728 -1.49985466 -1.49984159 -1.49982847
-1.4998159 -1.49980488 -1.49979693 -1.49979438 -1.49980076 -1.49982132
-1.49986388 -1.49993989 -1.50006592 -1.50026544 -1.50057068 -1.5010238
-1.50167564 -1.50257898 -1.50377191 -1.50524791 -1.50691602 -1.50857033
-1.50990493 -1.51060284 -1.51047421 -1.50955526 -1.50808857 -1.50640137
-1.50477517 -1.50337982 -1.50227653 -1.50145437 -1.5008683 -1.50046492
-1.50019563 -1.50002129 -1.49991251 -1.49984809 -1.4998132 -1.49979762
-1.4997944 -1.49979889 -1.49980797 -1.49981959 -1.49983241 -1.49984557
-1.49985853 -1.49987097 -1.49988271 -1.49989365 -1.49990377 -1.49991308
-1.49992159 -1.49992937 -1.49993645 -1.4999429 -1.49994875 -1.49995407
-1.4999589 -1.49996328 -1.49996726 -1.49997087 -1.49997414 -1.49997712
-1.49997981 -1.49998226 -1.49998447 -1.49998648 -1.49998829 -1.49998993
-1.4999914 -1.49999273 -1.49999392 -1.49999498 -1.49999592 -1.49999675
-1.49999747 -1.4999981 -1.49999863 -1.49999907 -1.49999942 -1.49999969
-1.49999987 -1.49999997 -1.5 -1.49999994 -1.49999981 -1.49999959
-1.49999929 -1.4999989 -1.49999842 -1.49999786 -1.49999719 -1.49999643
-1.49999556 -1.49999457 -1.49999346 -1.49999221 -1.49999083 -1.49998929
-1.49998758 -1.49998569 -1.4999836 -1.4999813 -1.49997876 -1.49997595
-1.49997286 -1.49996946 -1.4999657 -1.49996157 -1.49995701 -1.49995199
-1.49994646 -1.49994037 -1.49993368 -1.49992632 -1.49991825 -1.49990942
-1.49989979 -1.49988934 -1.49987807 -1.49986604 -1.49985336 -1.49984027
-1.49982716 -1.4998147 -1.49980391 -1.49979638 -1.49979455 -1.49980208
-1.49982445 -1.49986979 -1.49994999 -1.50008223 -1.50029079 -1.50060886
-1.50107959 -1.50175446 -1.50268571 -1.50390854 -1.50540978 -1.50708783
-1.50872493 -1.51000875 -1.51062866 -1.51041471 -1.509427 -1.50792313
-1.50623094 -1.50462231 -1.50325515 -1.50218152 -1.5013855 -1.50082025
-1.50043244 -1.50017432 -1.50000777 -1.4999043 -1.49984344 -1.49981089
-1.49979682 -1.49979457 -1.49979964 -1.49980906 -1.49982085 -1.49983373
-1.49984689 -1.49985981 -1.49987219 -1.49988385 -1.49989471 -1.49990475
-1.49991397 -1.49992241 -1.49993011 -1.49993713 -1.49994351 -1.49994931
-1.49995458 -1.49995936 -1.4999637 -1.49996764 -1.49997121 -1.49997446
-1.4999774 -1.49998007 -1.49998249 -1.49998468 -1.49998667 -1.49998846
-1.49999008 -1.49999154 -1.49999286 -1.49999403 -1.49999508 -1.49999601
-1.49999683 -1.49999754 -1.49999816 -1.49999868 -1.49999911 -1.49999945
-1.49999971 -1.49999988 -1.49999998 -1.5 -1.49999993 -1.49999979
-1.49999956 -1.49999925 -1.49999886 -1.49999837 -1.4999978 -1.49999712
-1.49999635 -1.49999546 -1.49999446 -1.49999334 -1.49999208 -1.49999068
-1.49998912 -1.4999874 -1.49998549 -1.49998338 -1.49998105 -1.49997849
-1.49997566 -1.49997253 -1.49996909 -1.4999653 -1.49996113 -1.49995652
-1.49995145 -1.49994587 -1.49993973 -1.49993297 -1.49992554 -1.4999174
-1.49990849 -1.49989878 -1.49988824 -1.49987689 -1.49986478 -1.49985205
-1.49983894 -1.49982587 -1.49981352 -1.49980297 -1.49979589 -1.49979481
-1.49980355 -1.49982782 -1.49987606 -1.49996062 -1.50009932 -1.50031724
-1.50064858 -1.50113745 -1.50183589 -1.50279545 -1.50404807 -1.50557355
-1.50725931 -1.50887594 -1.51010569 -1.51064581 -1.51034718 -1.50929344
-1.50775585 -1.50606146 -1.504472 -1.50313352 -1.50208936 -1.50131897
-1.50077398 -1.50040126 -1.50015393 -1.49999488 -1.49989651 -1.49983907
-1.49980877 -1.49979615 -1.49979481 -1.49980043 -1.49981016 -1.49982211
-1.49983506 -1.49984821 -1.49986109 -1.4998734 -1.49988498 -1.49989576
-1.49990571 -1.49991486 -1.49992322 -1.49993085 -1.4999378 -1.49994412
-1.49994987 -1.49995508 -1.49995982 -1.49996411 -1.49996802 -1.49997156
-1.49997477 -1.49997768 -1.49998032 -1.49998272 -1.49998489 -1.49998686
-1.49998863 -1.49999024 -1.49999168 -1.49999298 -1.49999414 -1.49999518
-1.4999961 -1.4999969 -1.49999761 -1.49999821 -1.49999872 -1.49999914
-1.49999948 -1.49999973 -1.4999999 -1.49999999 -1.49999999 -1.49999992
-1.49999977 -1.49999953 -1.49999922 -1.49999881 -1.49999832 -1.49999773
-1.49999705 -1.49999626 -1.49999537 -1.49999435 -1.49999322 -1.49999195
-1.49999053 -1.49998896 -1.49998722 -1.49998529 -1.49998316 -1.49998081
-1.49997821 -1.49997535 -1.4999722 -1.49996873 -1.4999649 -1.49996068
-1.49995603 -1.49995092 -1.49994528 -1.49993907 -1.49993225 -1.49992475
-1.49991654 -1.49990755 -1.49989775 -1.49988713 -1.4998757 -1.49986353
-1.49985074 -1.49983761 -1.49982457 -1.49981235 -1.49980207 -1.49979546
-1.49979517 -1.49980518 -1.49983142 -1.49988269 -1.49997179 -1.5001172
-1.50034484 -1.50068989 -1.50119744 -1.50191999 -1.5029082 -1.50419046
-1.50573905 -1.50743018 -1.50902301 -1.5101955 -1.51065423 -1.51027181
-1.50915493 -1.50758703 -1.50589314 -1.50432432 -1.50301493 -1.502
-1.50125473 -1.50072947 -1.50037135 -1.50013443 -1.4999826 -1.49988914
-1.49983497 -1.49980682 -1.49979559 -1.49979512 -1.49980126 -1.49981129
-1.49982339 -1.49983638 -1.49984953 -1.49986236 -1.4998746 -1.49988611
-1.4998968 -1.49990667 -1.49991573 -1.49992402 -1.49993158 -1.49993847
-1.49994473 -1.49995041 -1.49995558 -1.49996027 -1.49996453 -1.49996839
-1.49997189 -1.49997507 -1.49997796 -1.49998058 -1.49998295 -1.4999851
-1.49998704 -1.4999888 -1.49999039 -1.49999182 -1.4999931 -1.49999425
-1.49999528 -1.49999618 -1.49999698 -1.49999767 -1.49999827 -1.49999877
-1.49999918 -1.49999951 -1.49999975 -1.49999991 -1.49999999 -1.49999999
-1.49999991 -1.49999975 -1.49999951 -1.49999918 -1.49999877 -1.49999826
-1.49999767 -1.49999697 -1.49999618 -1.49999527 -1.49999424 -1.49999309
-1.49999181 -1.49999038 -1.49998879 -1.49998703 -1.49998508 -1.49998293
-1.49998056 -1.49997794 -1.49997505 -1.49997187 -1.49996836 -1.49996449
-1.49996023 -1.49995554 -1.49995037 -1.49994468 -1.49993841 -1.49993153
-1.49992396 -1.49991567 -1.4999066 -1.49989672 -1.49988602 -1.49987451
-1.49986226 -1.49984943 -1.49983628 -1.49982329 -1.4998112 -1.4998012
-1.49979509 -1.49979563 -1.49980697 -1.49983528 -1.49988969 -1.49998353
-1.50013591 -1.50037362 -1.50073284 -1.50125961 -1.5020068 -1.50302398
-1.50433563 -1.5059061 -1.50760012 -1.50916579 -1.5102779 -1.51065389
-1.51018882 -1.5090118 -1.50741701 -1.5057262 -1.50417936 -1.50289938
-1.5019134 -1.50119273 -1.50068664 -1.50034266 -1.50011579 -1.49997091
-1.49988216 -1.49983114 -1.49980505 -1.49979514 -1.49979549 -1.49980214
-1.49981244 -1.49982467 -1.49983771 -1.49985084 -1.49986362 -1.4998758
-1.49988722 -1.49989783 -1.49990762 -1.4999166 -1.49992482 -1.49993231
-1.49993912 -1.49994532 -1.49995096 -1.49995607 -1.49996072 -1.49996493
-1.49996876 -1.49997223 -1.49997538 -1.49997824 -1.49998083 -1.49998318
-1.4999853 -1.49998723 -1.49998897 -1.49999054 -1.49999196 -1.49999323
-1.49999436 -1.49999537 -1.49999627 -1.49999705 -1.49999774 -1.49999832
-1.49999881 -1.49999922 -1.49999954 -1.49999977 -1.49999992 -1.49999999
-1.49999999 -1.4999999 -1.49999973 -1.49999948 -1.49999914 -1.49999872
-1.49999821 -1.4999976 -1.4999969 -1.49999609 -1.49999517 -1.49999413
-1.49999297 -1.49999167 -1.49999023 -1.49998862 -1.49998684 -1.49998488
-1.4999827 -1.4999803 -1.49997766 -1.49997474 -1.49997153 -1.49996799
-1.49996408 -1.49995978 -1.49995504 -1.49994982 -1.49994408 -1.49993775
-1.49993079 -1.49992316 -1.49991479 -1.49990564 -1.49989568 -1.49988489
-1.49987331 -1.49986099 -1.49984811 -1.49983496 -1.49982202 -1.49981008
-1.49980037 -1.49979479 -1.4997962 -1.49980893 -1.4998394 -1.4998971
-1.49999586 -1.50015547 -1.50040363 -1.5007775 -1.50132402 -1.50209638
-1.50314281 -1.50448352 -1.50607451 -1.50776883 -1.50930394 -1.51035268
-1.5106448 -1.51009845 -1.50886441 -1.50724608 -1.50556083 -1.50403719
-1.50278686 -1.50182951 -1.50113291 -1.50064545 -1.50031516 -1.50009797
-1.49995978 -1.49987556 -1.49982755 -1.49980343 -1.49979479 -1.49979593
-1.49980304 -1.49981361 -1.49982597 -1.49983904 -1.49985215 -1.49986488
-1.49987698 -1.49988833 -1.49989886 -1.49990856 -1.49991746 -1.4999256
-1.49993302 -1.49993978 -1.49994592 -1.4999515 -1.49995656 -1.49996116
-1.49996533 -1.49996912 -1.49997256 -1.49997568 -1.49997851 -1.49998107
-1.4999834 -1.49998551 -1.49998741 -1.49998914 -1.49999069 -1.49999209
-1.49999335 -1.49999447 -1.49999547 -1.49999635 -1.49999713 -1.4999978
-1.49999838 -1.49999886 -1.49999925 -1.49999956 -1.49999979 -1.49999993
-1.5 -1.49999998 -1.49999988 -1.49999971 -1.49999945 -1.4999991
-1.49999867 -1.49999815 -1.49999754 -1.49999682 -1.499996 -1.49999507
-1.49999402 -1.49999285 -1.49999153 -1.49999007 -1.49998845 -1.49998665
-1.49998467 -1.49998247 -1.49998005 -1.49997738 -1.49997443 -1.49997119
-1.49996761 -1.49996367 -1.49995932 -1.49995454 -1.49994927 -1.49994346
-1.49993708 -1.49993006 -1.49992235 -1.4999139 -1.49990467 -1.49989463
-1.49988376 -1.4998721 -1.49985972 -1.49984679 -1.49983363 -1.49982075
-1.49980897 -1.49979958 -1.49979455 -1.49979688 -1.49981106 -1.49984379
-1.49990491 -1.50000879 -1.50017593 -1.5004349 -1.50082389 -1.50139073
-1.50218876 -1.50326468 -1.50463403 -1.50624407 -1.50793598 -1.5094371
-1.5104196 -1.51062698 -1.51000097 -1.50871311 -1.50707457 -1.50539721
-1.50389788 -1.50267736 -1.50174828 -1.50107521 -1.50060585 -1.50028879
-1.50008095 -1.49994919 -1.49986932 -1.4998242 -1.49980197 -1.49979453
-1.49979642 -1.49980398 -1.49981479 -1.49982726 -1.49984037 -1.49985346
-1.49986613 -1.49987816 -1.49988943 -1.49989987 -1.49990949 -1.49991832
-1.49992638 -1.49993373 -1.49994042 -1.4999465 -1.49995203 -1.49995705
-1.4999616 -1.49996573 -1.49996948 -1.49997289 -1.49997598 -1.49997878
-1.49998132 -1.49998362 -1.49998571 -1.4999876 -1.4999893 -1.49999084
-1.49999222 -1.49999347 -1.49999458 -1.49999556 -1.49999644 -1.4999972
-1.49999786 -1.49999843 -1.4999989 -1.49999929 -1.49999959 -1.49999981
-1.49999994 -1.5 -1.49999997 -1.49999987 -1.49999968 -1.49999942
-1.49999906 -1.49999862 -1.49999809 -1.49999747 -1.49999674 -1.49999591
-1.49999497 -1.49999391 -1.49999272 -1.49999139 -1.49998992 -1.49998828
-1.49998646 -1.49998446 -1.49998224 -1.49997979 -1.49997709 -1.49997412
-1.49997084 -1.49996723 -1.49996325 -1.49995886 -1.49995403 -1.49994871
-1.49994285 -1.4999364 -1.49992931 -1.49992153]
</pre></div>
</div>
</div>
</div>
</div>
</div>
<div class="section" id="scattering-and-cross-sections">
<h2>Scattering and Cross Sections<a class="headerlink" href="#scattering-and-cross-sections" title="Permalink to this headline"></a></h2>
<p>Scattering experiments don’t measure entire trajectories. For elastic
collisions, they measure the distribution of final scattering angles
at best. Most experiments use targets thin enough so that the number
of scatterings is typically zero or one. The cross section, <span class="math notranslate nohighlight">\(\sigma\)</span>,
describes the cross-sectional area for particles to scatter with an
individual target atom or nucleus. Cross section measurements form the
basis for MANY fields of physics. BThe cross section, and the
differential cross section, encapsulates everything measurable for a
collision where all that is measured is the final state, e.g. the
outgoing particle had momentum <span class="math notranslate nohighlight">\(\boldsymbol{p}_f\)</span>. y studying cross sections,
one can infer information about the potential interaction between the
two particles. Inferring, or constraining, the potential from the
cross section is a classic {\it inverse} problem. Collisions are
either elastic or inelastic. Elastic collisions are those for which
the two bodies are in the same internal state before and after the
collision. If the collision excites one of the participants into a
higher state, or transforms the particles into different species, or
creates additional particles, the collision is inelastic. Here, we
consider only elastic collisions.</p>
<p>For Coulomb forces, the cross section is infinite because the range of
the Coulomb force is infinite, but for interactions such as the strong
interaction in nuclear or particle physics, there is no long-range
force and cross-sections are finite. Even for Coulomb forces, the part
of the cross section that corresponds to a specific scattering angle,
<span class="math notranslate nohighlight">\(d\sigma/d\Omega\)</span>, which is a function of the scattering angle
<span class="math notranslate nohighlight">\(\theta_s\)</span> is still finite.</p>
<p>If a particle travels through a thin target, the chance the particle
scatters is <span class="math notranslate nohighlight">\(P_{\rm scatt}=\sigma dN/dA\)</span>, where <span class="math notranslate nohighlight">\(dN/dA\)</span> is the number
of scattering centers per area the particle encounters. If the density
of the target is <span class="math notranslate nohighlight">\(\rho\)</span> particles per volume, and if the thickness of
the target is <span class="math notranslate nohighlight">\(t\)</span>, the areal density (number of target scatterers per
area) is <span class="math notranslate nohighlight">\(dN/dA=\rho t\)</span>. Because one wishes to quantify the collisions
independently of the target, experimentalists measure scattering
probabilities, then divide by the areal density to obtain
cross-sections,</p>
<div class="math notranslate nohighlight">
\[
\begin{eqnarray}
\sigma=\frac{P_{\rm scatt}}{dN/dA}.
\end{eqnarray}
\]</div>
<p>Instead of merely stating that a particle collided, one can measure
the probability the particle scattered by a given angle. The
scattering angle <span class="math notranslate nohighlight">\(\theta_s\)</span> is defined so that at zero the particle is
unscattered and at <span class="math notranslate nohighlight">\(\theta_s=\pi\)</span> the particle is scattered directly
backward. Scattering angles are often described in the center-of-mass
frame, but that is a detail we will neglect for this first discussion,
where we will consider the scattering of particles moving classically
under the influence of fixed potentials <span class="math notranslate nohighlight">\(U(\boldsymbol{r})\)</span>. Because the
distribution of scattering angles can be measured, one expresses the
differential cross section,</p>
<!-- Equation labels as ordinary links -->
<div id="_auto5"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\frac{d^2\sigma}{d\cos\theta_s~d\phi}.
\label{_auto5} \tag{10}
\end{equation}
\]</div>
<p>Usually, the literature expresses differential cross sections as</p>
<!-- Equation labels as ordinary links -->
<div id="_auto6"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
d\sigma/d\Omega=\frac{d\sigma}{d\cos\theta d\phi}=\frac{1}{2\pi}\frac{d\sigma}{d\cos\theta},
\label{_auto6} \tag{11}
\end{equation}
\]</div>
<p>where the last equivalency is true when the scattering does not depend
on the azimuthal angle <span class="math notranslate nohighlight">\(\phi\)</span>, as is the case for spherically
symmetric potentials.</p>
<p>The differential solid angle <span class="math notranslate nohighlight">\(d\Omega\)</span> can be thought of as the area
subtended by a measurement, <span class="math notranslate nohighlight">\(dA_d\)</span>, divided by <span class="math notranslate nohighlight">\(r^2\)</span>, where <span class="math notranslate nohighlight">\(r\)</span> is the
distance to the detector,</p>
<div class="math notranslate nohighlight">
\[
\begin{eqnarray}
dA_d=r^2 d\Omega.
\end{eqnarray}
\]</div>
<p>With this definition <span class="math notranslate nohighlight">\(d\sigma/d\Omega\)</span> is independent of the distance
from which one places the detector, or the size of the detector (as
long as it is small).</p>
<p>Differential scattering cross sections are calculated by assuming a
random distribution of impact parameters <span class="math notranslate nohighlight">\(b\)</span>. These represent the
distance in the <span class="math notranslate nohighlight">\(xy\)</span> plane for particles moving in the <span class="math notranslate nohighlight">\(z\)</span> direction
relative to the scattering center. An impact parameter <span class="math notranslate nohighlight">\(b=0\)</span> refers to
being aimed directly at the target’s center. The impact parameter
describes the transverse distance from the <span class="math notranslate nohighlight">\(z=0\)</span> axis for the
trajectory when it is still far away from the scattering center and
has not yet passed it. The differential cross section can be expressed
in terms of the impact parameter,</p>
<!-- Equation labels as ordinary links -->
<div id="_auto7"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
d\sigma=2\pi bdb,
\label{_auto7} \tag{12}
\end{equation}
\]</div>
<p>which is the area of a thin ring of radius <span class="math notranslate nohighlight">\(b\)</span> and thickness <span class="math notranslate nohighlight">\(db\)</span>. In
classical physics, one can calculate the trajectory given the incoming
kinetic energy <span class="math notranslate nohighlight">\(E\)</span> and the impact parameter if one knows the mass and
potential. From the trajectory, one then finds the scattering angle
<span class="math notranslate nohighlight">\(\theta_s(b)\)</span>. The differential cross section is then</p>
<!-- Equation labels as ordinary links -->
<div id="_auto8"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\frac{d\sigma}{d\Omega}=\frac{1}{2\pi}\frac{d\sigma}{d\cos\theta_s}=b\frac{db}{d\cos\theta_s}=\frac{b}{(d/db)\cos\theta_s(b)}.
\label{_auto8} \tag{13}
\end{equation}
\]</div>
<p>Typically, one would calculate <span class="math notranslate nohighlight">\(\cos\theta_s\)</span> and <span class="math notranslate nohighlight">\((d/db)\cos\theta_s\)</span>
as functions of <span class="math notranslate nohighlight">\(b\)</span>. This is sufficient to plot the differential cross
section as a function of <span class="math notranslate nohighlight">\(\theta_s\)</span>.</p>
<p>The total cross section is</p>
<!-- Equation labels as ordinary links -->
<div id="_auto9"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\sigma_{\rm tot}=\int d\Omega\frac{d\sigma}{d\Omega}=2\pi\int d\cos\theta_s~\frac{d\sigma}{d\Omega}.
\label{_auto9} \tag{14}
\end{equation}
\]</div>
<p>Even if the total cross section is infinite, e.g. Coulomb forces, one
can still have a finite differential cross section as we will see
later on.</p>
<p>An asteroid of mass <span class="math notranslate nohighlight">\(m\)</span> and kinetic energy <span class="math notranslate nohighlight">\(E\)</span> approaches a planet of
radius <span class="math notranslate nohighlight">\(R\)</span> and mass <span class="math notranslate nohighlight">\(M\)</span>. What is the cross section for the asteroid to
impact the planet?</p>
<div class="section" id="solution">
<h3>Solution<a class="headerlink" href="#solution" title="Permalink to this headline"></a></h3>
<p>Calculate the maximum impact parameter, <span class="math notranslate nohighlight">\(b_{\rm max}\)</span>, for which the asteroid will hit the planet. The total cross section for impact is <span class="math notranslate nohighlight">\(\sigma_{\rm impact}=\pi b_{\rm max}^2\)</span>. The maximum cross-section can be found with the help of angular momentum conservation. The asteroid’s incoming momentum is <span class="math notranslate nohighlight">\(p_0=\sqrt{2mE}\)</span> and the angular momentum is <span class="math notranslate nohighlight">\(L=p_0b\)</span>. If the asteroid just grazes the planet, it is moving with zero radial kinetic energy at impact. Combining energy and angular momentum conservation and having <span class="math notranslate nohighlight">\(p_f\)</span> refer to the momentum of the asteroid at a distance <span class="math notranslate nohighlight">\(R\)</span>,</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray*}
\frac{p_f^2}{2m}-\frac{GMm}{R}&amp;=&amp;E,\\
p_fR&amp;=&amp;p_0b_{\rm max},
\end{eqnarray*}
\end{split}\]</div>
<p>allows one to solve for <span class="math notranslate nohighlight">\(b_{\rm max}\)</span>,</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray*}
b_{\rm max}&amp;=&amp;R\frac{p_f}{p_0}\\
&amp;=&amp;R\frac{\sqrt{2m(E+GMm/R)}}{\sqrt{2mE}}\\
\sigma_{\rm impact}&amp;=&amp;\pi R^2\frac{E+GMm/R}{E}.
\end{eqnarray*}
\end{split}\]</div>
</div>
</div>
<div class="section" id="rutherford-scattering">
<h2>Rutherford Scattering<a class="headerlink" href="#rutherford-scattering" title="Permalink to this headline"></a></h2>
<p>This refers to the calculation of <span class="math notranslate nohighlight">\(d\sigma/d\Omega\)</span> due to an inverse
square force, <span class="math notranslate nohighlight">\(F_{12}=\pm\alpha/r^2\)</span> for repulsive/attractive
interaction. Rutherford compared the scattering of <span class="math notranslate nohighlight">\(\alpha\)</span> particles
(<span class="math notranslate nohighlight">\(^4\)</span>He nuclei) off of a nucleus and found the scattering angle at
which the formula began to fail. This corresponded to the impact
parameter for which the trajectories would strike the nucleus. This
provided the first measure of the size of the atomic nucleus. At the
time, the distribution of the positive charge (the protons) was
considered to be just as spread out amongst the atomic volume as the
electrons. After Rutherford’s experiment, it was clear that the radius
of the nucleus tended to be roughly 4 orders of magnitude smaller than
that of the atom, which is less than the size of a football relative
to Spartan Stadium.</p>
<p>The incoming and outgoing angles of the trajectory are at
<span class="math notranslate nohighlight">\(\pm\theta'\)</span>. They are related to the scattering angle by
<span class="math notranslate nohighlight">\(2\theta'=\pi+\theta_s\)</span>.</p>
<p>In order to calculate differential cross section, we must find how the
impact parameter is related to the scattering angle. This requires
analysis of the trajectory. We consider our previous expression for
the trajectory where we derived the elliptic form for the trajectory,
Eq. (<a class="reference external" href="#eq:Ctrajectory">9</a>). For that case we considered an attractive
force with the particle’s energy being negative, i.e. it was
bound. However, the same form will work for positive energy, and
repulsive forces can be considered by simple flipping the sign of
<span class="math notranslate nohighlight">\(\alpha\)</span>. For positive energies, the trajectories will be hyperbolas,
rather than ellipses, with the asymptotes of the trajectories
representing the directions of the incoming and outgoing
tracks. Rewriting Eq. (<a class="reference external" href="#eq:Ctrajectory">9</a>),</p>
<!-- Equation labels as ordinary links -->
<div id="eq:ruthtraj"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}\label{eq:ruthtraj} \tag{15}
r=\frac{1}{\frac{m\alpha}{L^2}+A\cos\theta}.
\end{equation}
\]</div>
<p>Once <span class="math notranslate nohighlight">\(A\)</span> is large enough, which will happen when the energy is
positive, the denominator will become negative for a range of
<span class="math notranslate nohighlight">\(\theta\)</span>. This is because the scattered particle will never reach
certain angles. The asymptotic angles <span class="math notranslate nohighlight">\(\theta'\)</span> are those for which
the denominator goes to zero,</p>
<!-- Equation labels as ordinary links -->
<div id="_auto10"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\cos\theta'=-\frac{m\alpha}{AL^2}.
\label{_auto10} \tag{16}
\end{equation}
\]</div>
<p>The trajectory’s point of closest approach is at <span class="math notranslate nohighlight">\(\theta=0\)</span> and the
two angles <span class="math notranslate nohighlight">\(\theta'\)</span>, which have this value of <span class="math notranslate nohighlight">\(\cos\theta'\)</span>, are the
angles of the incoming and outgoing particles. From
Fig (<strong>to come</strong>), one can see that the scattering angle
<span class="math notranslate nohighlight">\(\theta_s\)</span> is given by,</p>
<!-- Equation labels as ordinary links -->
<div id="eq:sthetover2"></div>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
\label{eq:sthetover2} \tag{17}
2\theta'-\pi&amp;=&amp;\theta_s,~~~\theta'=\frac{\pi}{2}+\frac{\theta_s}{2},\\
\nonumber
\sin(\theta_s/2)&amp;=&amp;-\cos\theta'\\
\nonumber
&amp;=&amp;\frac{m\alpha}{AL^2}.
\end{eqnarray}
\end{split}\]</div>
<p>Now that we have <span class="math notranslate nohighlight">\(\theta_s\)</span> in terms of <span class="math notranslate nohighlight">\(m,\alpha,L\)</span> and <span class="math notranslate nohighlight">\(A\)</span>, we wish
to re-express <span class="math notranslate nohighlight">\(L\)</span> and <span class="math notranslate nohighlight">\(A\)</span> in terms of the impact parameter <span class="math notranslate nohighlight">\(b\)</span> and the
energy <span class="math notranslate nohighlight">\(E\)</span>. This will set us up to calculate the differential cross
section, which requires knowing <span class="math notranslate nohighlight">\(db/d\theta_s\)</span>. It is easy to write
the angular momentum as</p>
<!-- Equation labels as ordinary links -->
<div id="_auto11"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
L^2=p_0^2b^2=2mEb^2.
\label{_auto11} \tag{18}
\end{equation}
\]</div>
<p>Finding <span class="math notranslate nohighlight">\(A\)</span> is more complicated. To accomplish this we realize that
the point of closest approach occurs at <span class="math notranslate nohighlight">\(\theta=0\)</span>, so from
Eq. (<a class="reference external" href="#eq:ruthtraj">15</a>)</p>
<!-- Equation labels as ordinary links -->
<div id="eq:rminofA"></div>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
\label{eq:rminofA} \tag{19}
\frac{1}{r_{\rm min}}&amp;=&amp;\frac{m\alpha}{L^2}+A,\\
\nonumber
A&amp;=&amp;\frac{1}{r_{\rm min}}-\frac{m\alpha}{L^2}.
\end{eqnarray}
\end{split}\]</div>
<p>Next, <span class="math notranslate nohighlight">\(r_{\rm min}\)</span> can be found in terms of the energy because at the
point of closest approach the kinetic energy is due purely to the
motion perpendicular to <span class="math notranslate nohighlight">\(\hat{r}\)</span> and</p>
<!-- Equation labels as ordinary links -->
<div id="_auto12"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
E=-\frac{\alpha}{r_{\rm min}}+\frac{L^2}{2mr_{\rm min}^2}.
\label{_auto12} \tag{20}
\end{equation}
\]</div>
<p>One can solve the quadratic equation for <span class="math notranslate nohighlight">\(1/r_{\rm min}\)</span>,</p>
<!-- Equation labels as ordinary links -->
<div id="_auto13"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\frac{1}{r_{\rm min}}=\frac{m\alpha}{L^2}+\sqrt{(m\alpha/L^2)^2+2mE/L^2}.
\label{_auto13} \tag{21}
\end{equation}
\]</div>
<p>We can plug the expression for <span class="math notranslate nohighlight">\(r_{\rm min}\)</span> into the expression for <span class="math notranslate nohighlight">\(A\)</span>, Eq. (<a class="reference external" href="#eq:rminofA">19</a>),</p>
<!-- Equation labels as ordinary links -->
<div id="_auto14"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
A=\sqrt{(m\alpha/L^2)^2+2mE/L^2}=\sqrt{(\alpha^2/(4E^2b^4)+1/b^2}
\label{_auto14} \tag{22}
\end{equation}
\]</div>
<p>Finally, we insert the expression for <span class="math notranslate nohighlight">\(A\)</span> into that for the scattering angle, Eq. (<a class="reference external" href="#eq:sthetover2">17</a>),</p>
<!-- Equation labels as ordinary links -->
<div id="eq:scattangle"></div>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
\label{eq:scattangle} \tag{23}
\sin(\theta_s/2)&amp;=&amp;\frac{m\alpha}{AL^2}\\
\nonumber
&amp;=&amp;\frac{a}{\sqrt{a^2+b^2}}, ~~a\equiv \frac{\alpha}{2E}
\end{eqnarray}
\end{split}\]</div>
<p>The differential cross section can now be found by differentiating the
expression for <span class="math notranslate nohighlight">\(\theta_s\)</span> with <span class="math notranslate nohighlight">\(b\)</span>,</p>
<!-- Equation labels as ordinary links -->
<div id="eq:rutherford"></div>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray}
\label{eq:rutherford} \tag{24}
\frac{1}{2}\cos(\theta_s/2)d\theta_s&amp;=&amp;\frac{ab~db}{(a^2+b^2)^{3/2}}=\frac{bdb}{a^2}\sin^3(\theta_s/2),\\
\nonumber
d\sigma&amp;=&amp;2\pi bdb=\frac{\pi a^2}{\sin^3(\theta_s/2)}\cos(\theta_s/2)d\theta_s\\
\nonumber
&amp;=&amp;\frac{\pi a^2}{2\sin^4(\theta_s/2)}\sin\theta_s d\theta_s\\
\nonumber
\frac{d\sigma}{d\cos\theta_s}&amp;=&amp;\frac{\pi a^2}{2\sin^4(\theta_s/2)},\\
\nonumber
\frac{d\sigma}{d\Omega}&amp;=&amp;\frac{a^2}{4\sin^4(\theta_s/2)}.
\end{eqnarray}
\end{split}\]</div>
<p>where <span class="math notranslate nohighlight">\(a= \alpha/2E\)</span>. This the Rutherford formula for the differential
cross section. It diverges as <span class="math notranslate nohighlight">\(\theta_s\rightarrow 0\)</span> because
scatterings with arbitrarily large impact parameters still scatter to
arbitrarily small scattering angles. The expression for
<span class="math notranslate nohighlight">\(d\sigma/d\Omega\)</span> is the same whether the interaction is positive or
negative.</p>
<p>Consider a particle of mass <span class="math notranslate nohighlight">\(m\)</span> and charge <span class="math notranslate nohighlight">\(z\)</span> with kinetic energy <span class="math notranslate nohighlight">\(E\)</span>
(Let it be the center-of-mass energy) incident on a heavy nucleus of
mass <span class="math notranslate nohighlight">\(M\)</span> and charge <span class="math notranslate nohighlight">\(Z\)</span> and radius <span class="math notranslate nohighlight">\(R\)</span>. Find the angle at which the
Rutherford scattering formula breaks down.</p>
<div class="section" id="id1">
<h3>Solution<a class="headerlink" href="#id1" title="Permalink to this headline"></a></h3>
<p>Let <span class="math notranslate nohighlight">\(\alpha=Zze^2/(4\pi\epsilon_0)\)</span>. The scattering angle in Eq. (<a class="reference external" href="#eq:scattangle">23</a>) is</p>
<div class="math notranslate nohighlight">
\[
\sin(\theta_s/2)=\frac{a}{\sqrt{a^2+b^2}}, ~~a\equiv \frac{\alpha}{2E}.
\]</div>
<p>The impact parameter <span class="math notranslate nohighlight">\(b\)</span> for which the point of closest approach
equals <span class="math notranslate nohighlight">\(R\)</span> can be found by using angular momentum conservation,</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray*}
p_0b&amp;=&amp;b\sqrt{2mE}=Rp_f=R\sqrt{2m(E-\alpha/R)},\\
b&amp;=&amp;R\frac{\sqrt{2m(E-\alpha/R)}}{\sqrt{2mE}}\\
&amp;=&amp;R\sqrt{1-\frac{\alpha}{ER}}.
\end{eqnarray*}
\end{split}\]</div>
<p>Putting these together</p>
<div class="math notranslate nohighlight">
\[
\theta_s=2\sin^{-1}\left\{
\frac{a}{\sqrt{a^2+R^2(1-\alpha/(RE))}}
\right\},~~~a=\frac{\alpha}{2E}.
\]</div>
<p>It was from this departure of the experimentally measured
<span class="math notranslate nohighlight">\(d\sigma/d\Omega\)</span> from the Rutherford formula that allowed Rutherford
to infer the radius of the gold nucleus, <span class="math notranslate nohighlight">\(R\)</span>.</p>
<p>Just like electrodynamics, one can define “fields”, which for a small
additional mass <span class="math notranslate nohighlight">\(m\)</span> are the force per mass and the additional
potential energy per mass. The {\it gravitational field} related to
the force has dimensions of force per mass, or acceleration, and can
be labeled <span class="math notranslate nohighlight">\(\boldsymbol{g}(\boldsymbol{r})\)</span>. The potential energy per mass has
dimensions of energy per mass. This is analogous to the
electromagnetic potential, which is the potential energy per charge,
and the electric field which is the force per charge.</p>
<p>Because the field <span class="math notranslate nohighlight">\(\boldsymbol{g}\)</span> obeys the same inverse square law for a
point mass as the electric field does for a point charge, the
gravitational field also satisfies a version of Gauss’s law,</p>
<!-- Equation labels as ordinary links -->
<div id="eq:GravGauss"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\label{eq:GravGauss} \tag{25}
\oint d\boldsymbol{A}\cdot\boldsymbol{g}=-4\pi GM_{\rm inside}.
\end{equation}
\]</div>
<p>Here, <span class="math notranslate nohighlight">\(M_{\rm inside}\)</span> is the net mass inside a closed area.</p>
<p>Gauss’s law can be understood by considering a nozzle that sprays
paint in all directions uniformly from a point source. Let <span class="math notranslate nohighlight">\(B\)</span> be the
number of gallons per minute of paint leaving the nozzle. If the
nozzle is at the center of a sphere of radius <span class="math notranslate nohighlight">\(r\)</span>, the paint per
square meter per minute that is deposited on some part of the sphere
is</p>
<div class="math notranslate nohighlight">
\[
\begin{eqnarray}
F(r)&amp;=&amp;\frac{B}{4\pi r^2}.
\end{eqnarray}
\]</div>
<p>Now, let <span class="math notranslate nohighlight">\(F\)</span> also be assigned a direction, so that it becomes a vector
pointing along the direction of the flying paint. For any surface that
surrounds the nozzle, not necessarily a sphere, one can state that</p>
<!-- Equation labels as ordinary links -->
<div id="eq:paint"></div>
<div class="math notranslate nohighlight">
\[
\begin{eqnarray}
\label{eq:paint} \tag{26}
\oint \boldsymbol{dA}\cdot\boldsymbol{F}&amp;=&amp;B,
\end{eqnarray}
\]</div>
<p>regardless of the shape of the surface. This follows because the rate
at which paint is deposited on the surface should equal the rate at
which it leaves the nozzle. The dot product ensures that only the
component of <span class="math notranslate nohighlight">\(\boldsymbol{F}\)</span> into the surface contributes to the deposition
of paint. Similarly, if <span class="math notranslate nohighlight">\(\boldsymbol{F}\)</span> is any radial inverse-square forces,
that falls as <span class="math notranslate nohighlight">\(B/(4\pi r^2)\)</span>, then one can apply
Eq. (<a class="reference external" href="#eq:paint">26</a>). For gravitational fields, <span class="math notranslate nohighlight">\(B/(4\pi)\)</span> is replaced
by <span class="math notranslate nohighlight">\(GM\)</span>, and one quickly “derives” Gauss’s law for gravity,
Eq. (<a class="reference external" href="#eq:GravGauss">25</a>).</p>
<p>Consider Earth to have its mass <span class="math notranslate nohighlight">\(M\)</span> uniformly distributed in a sphere
of radius <span class="math notranslate nohighlight">\(R\)</span>. Find the magnitude of the gravitational acceleration as
a function of the radius <span class="math notranslate nohighlight">\(r\)</span> in terms of the acceleration of gravity
at the surface <span class="math notranslate nohighlight">\(g(R)\)</span>. Assume <span class="math notranslate nohighlight">\(r&lt;R\)</span>, i.e. you are inside the surface.</p>
<p>{\bf Solution}: Take the ratio of Eq. (<a class="reference external" href="#eq:GravGauss">25</a>) for two radii, <span class="math notranslate nohighlight">\(R\)</span> and <span class="math notranslate nohighlight">\(r&lt;R\)</span>,</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray*}
\frac{4\pi r^2 g(r)}{4\pi R^2 g(R)}&amp;=&amp;\frac{4\pi GM_{\rm inside~r}}{4\pi GM_{\rm inside~R}}\\
\nonumber
&amp;=&amp;\frac{r^3}{R^3}\\
\nonumber
g(r)&amp;=&amp;g(R)\frac{r}{R}~.
\end{eqnarray*}
\end{split}\]</div>
<p>The potential energy per mass is similar conceptually to the voltage, or electric potential energy per charge, that was studied in electromagnetism, if <span class="math notranslate nohighlight">\(V\equiv U/m\)</span>, <span class="math notranslate nohighlight">\(\boldsymbol{g}=-\nabla V\)</span>.</p>
</div>
</div>
<div class="section" id="tidal-forces">
<h2>Tidal Forces<a class="headerlink" href="#tidal-forces" title="Permalink to this headline"></a></h2>
<p>Consider a spherical planet of radius <span class="math notranslate nohighlight">\(r\)</span> a distance <span class="math notranslate nohighlight">\(D\)</span> from another
body of mass <span class="math notranslate nohighlight">\(M\)</span>. The magnitude of the force due to <span class="math notranslate nohighlight">\(M\)</span> on an small
object of mass <span class="math notranslate nohighlight">\(\delta m\)</span> on surface of the planet can be calculated
by performing a Taylor expansion about the center of the spherical
planet.</p>
<!-- Equation labels as ordinary links -->
<div id="_auto15"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
F=-\frac{GM\delta m}{D^2}+2\frac{GM\delta m}{D^3}\Delta D+\cdots
\label{_auto15} \tag{27}
\end{equation}
\]</div>
<p>If the <span class="math notranslate nohighlight">\(z\)</span> direction points toward the large object, <span class="math notranslate nohighlight">\(\Delta D\)</span> can be
referred to as <span class="math notranslate nohighlight">\(z\)</span>. In the accelerating frame of an observer at the
center of the planet,</p>
<!-- Equation labels as ordinary links -->
<div id="_auto16"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\delta m\frac{d^2 z}{dt^2}=F-\delta ma'+{\rm other~forces~acting~on~} \delta m,
\label{_auto16} \tag{28}
\end{equation}
\]</div>
<p>where <span class="math notranslate nohighlight">\(a'\)</span> is the acceleration of the observer. Because <span class="math notranslate nohighlight">\(\delta ma'\)</span>
equals the gravitational force on <span class="math notranslate nohighlight">\(\delta m\)</span> if it were located at the
planet’s center, one can write</p>
<!-- Equation labels as ordinary links -->
<div id="_auto17"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
m\frac{d^2z}{dt^2}=2\frac{GM\delta m}{D^3}z+{\rm other~forces~acting~on~}\delta m.
\label{_auto17} \tag{29}
\end{equation}
\]</div>
<p>Here the other forces could represent the forces acting on <span class="math notranslate nohighlight">\(\delta m\)</span>
from the spherical planet such as the gravitational force or the
contact force with the surface. If <span class="math notranslate nohighlight">\(\theta\)</span> is the angle w.r.t. the
<span class="math notranslate nohighlight">\(z\)</span> axis, the effective force acting on <span class="math notranslate nohighlight">\(\delta m\)</span> is</p>
<!-- Equation labels as ordinary links -->
<div id="_auto18"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
F_{\rm eff}\approx 2\frac{GM\delta m}{D^3}r\cos\theta\hat{z}+{\rm other~forces~acting~on~}\delta m.
\label{_auto18} \tag{30}
\end{equation}
\]</div>
<p>This first force is the “tidal” force. It pulls objects outward from the center of the object. If the object were covered with water, it would distort the objects shape so that the shape would be elliptical, stretched out along the axis pointing toward the large mass <span class="math notranslate nohighlight">\(M\)</span>. The force is always along (either parallel or antiparallel to) the <span class="math notranslate nohighlight">\(\hat{z}\)</span> direction.</p>
<p>Consider the Earth to be a sphere of radius <span class="math notranslate nohighlight">\(R\)</span> covered with water,
with the gravitational acceleration at the surface noted by <span class="math notranslate nohighlight">\(g\)</span>. Now
assume that a distant body provides an additional constant
gravitational acceleration <span class="math notranslate nohighlight">\(\boldsymbol{a}\)</span> pointed along the <span class="math notranslate nohighlight">\(z\)</span> axis. Find
the distortion of the radius as a function of <span class="math notranslate nohighlight">\(\theta\)</span>. Ignore
planetary rotation and assume <span class="math notranslate nohighlight">\(a&lt;&lt;g\)</span>.</p>
<p>{\bf Solution}: Because Earth would then accelerate with <span class="math notranslate nohighlight">\(a\)</span>, the
field <span class="math notranslate nohighlight">\(a\)</span> would seem invisible in the accelerating frame. A tidal
force would only appear if <span class="math notranslate nohighlight">\(a\)</span> depended on position, i.e. <span class="math notranslate nohighlight">\(\nabla
\boldsymbol{a}\ne 0\)</span>.</p>
<p>Now consider that the field is no longer constant, but that instead <span class="math notranslate nohighlight">\(a=-kz\)</span> with <span class="math notranslate nohighlight">\(|kR|&lt;&lt;g\)</span>.</p>
<p>{\bf Solution}: The surface of the planet needs to be at constant
potential (if the planet is not accelerating). The force per mass,
<span class="math notranslate nohighlight">\(-kz\)</span> is like a spring, and the potential per mass is
<span class="math notranslate nohighlight">\(kz^2/2\)</span>. Otherwise water would move to a point of lower
potential. Thus, the potential energy for a sample mass <span class="math notranslate nohighlight">\(\delta m\)</span> is</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray*}
V(R)+\delta m gh(\theta)-\frac{\delta m}{2}kr^2\cos^2\theta={\rm Constant}\\
V(R)+\delta mgh(\theta)-\frac{\delta m}{2}kR^2\cos^2\theta-\delta m kRh(\theta)\cos^2\theta-\frac{\delta m}{2}kh^2(\theta)\cos^2\theta={\rm Constant}.
\end{eqnarray*}
\end{split}\]</div>
<p>Here, the potential due to the external field is <span class="math notranslate nohighlight">\((1/2)kz^2\)</span> so that <span class="math notranslate nohighlight">\(-\nabla U=-kz\)</span>. One now needs to solve for <span class="math notranslate nohighlight">\(h(\theta)\)</span>. Absorbing all the constant terms from both sides of the equation into one constant <span class="math notranslate nohighlight">\(C\)</span>, and because both <span class="math notranslate nohighlight">\(h\)</span> and <span class="math notranslate nohighlight">\(kR\)</span> are small, we can through away terms of order <span class="math notranslate nohighlight">\(h^2\)</span> or <span class="math notranslate nohighlight">\(kRh\)</span>. This gives</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray*}
gh(\theta)-\frac{1}{2}kR^2\cos^2\theta&amp;=&amp;C,\\
h(\theta)&amp;=&amp;\frac{C}{g}+\frac{1}{2g}kR^2\cos^2\theta,\\
h(\theta)&amp;=&amp;\frac{1}{2g}kR^2(\cos^2\theta-1/3).
\end{eqnarray*}
\end{split}\]</div>
<p>The term with the factor of <span class="math notranslate nohighlight">\(1/3\)</span> replaced the constant and was chosen so that the average height of the water would be zero.</p>
<p>The Sun’s mass is <span class="math notranslate nohighlight">\(27\times 10^6\)</span> the Moon’s mass, but the Sun is 390 times further away from Earth as the Sun. What is ratio of the tidal force of the Sun to that of the Moon.</p>
<p>{\bf Solution}: The gravitational force due to an object <span class="math notranslate nohighlight">\(M\)</span> a distance <span class="math notranslate nohighlight">\(D\)</span> away goes as <span class="math notranslate nohighlight">\(M/D^2\)</span>, but the tidal force is only the difference of that force over a distance <span class="math notranslate nohighlight">\(R\)</span>,</p>
<div class="math notranslate nohighlight">
\[
F_{\rm tidal}\propto \frac{M}{D^3}R.
\]</div>
<p>Therefore the ratio of force is</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray*}
\frac{F_{\rm Sun's~tidal~force}}{F_{\rm Moon's~tidal~force}}
&amp;=&amp;\frac{M_{\rm sun}/D_{\rm sun}^3}{M_{\rm moon}/D_{\rm moon}^3}\\
&amp;=&amp;\frac{27\times 10^6}{390^3}=0.46.
\end{eqnarray*}
\end{split}\]</div>
<p>The Moon more strongly affects tides than the Sun.</p>
</div>
</div>
<script type="text/x-thebe-config">
{
requestKernel: true,
binderOptions: {
repo: "binder-examples/jupyter-stacks-datascience",
ref: "master",
},
codeMirrorConfig: {
theme: "abcdef",
mode: "python"
},
kernelOptions: {
kernelName: "python3",
path: "./testbook"
},
predefinedOutput: true
}
</script>
<script>kernelName = 'python3'</script>
</div>
</div>
</div>
<div class='prev-next-bottom'>
</div>
<footer class="footer mt-5 mt-md-0">
<div class="container">
<p>
By Morten Hjorth-Jensen<br/>
&copy; Copyright 2020.<br/>
</p>
</div>
</footer>
</main>
</div>
</div>
<script src="../_static/js/index.js"></script>
</body>
</html>
Reference in New Issue View Git Blame Copy Permalink