Skip to content

Commit

Permalink
more improvements
Browse files Browse the repository at this point in the history
  • Loading branch information
@syp2001
syp2001 committed Nov 5, 2023
1 parent 654a4c0 commit 9ebc8dc
Show file tree
Hide file tree
Showing 14 changed files with 155 additions and 55 deletions.
Binary file modified docs/source/images/thumbnail.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
9 changes: 3 additions & 6 deletions docs/source/index.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,15 +14,12 @@ The library also provides a set of methods for computing constants of motion and
:align: left
:width: 45%

.. https://github.com/sphinx-doc/sphinx/issues/823
.. raw:: html

<video width="45%" autoplay loop>
<source src="https://kerrgeopy.readthedocs.io/en/latest/_downloads/e0f5e701241a3dc4ff27d031b7edad24/animation4.mp4" type="video/mp4">
<p>

[](animation4.mp4)

</p>
<source src="https://raw.githubusercontent.com/BlackHolePerturbationToolkit/KerrGeoPy/main/docs/source/notebooks/animation4.mp4" type="video/mp4">
</video>

.. _Installation:
Expand Down
64 changes: 53 additions & 11 deletions docs/source/notebooks/Graphics.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
},
{
"cell_type": "code",
"execution_count": 1,
"execution_count": 16,
"id": "7ae99533",
"metadata": {},
"outputs": [
Expand Down Expand Up @@ -59,7 +59,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 2,
"id": "f3351c95",
"metadata": {},
"outputs": [
Expand Down Expand Up @@ -88,7 +88,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 3,
"id": "aa8f032c",
"metadata": {},
"outputs": [
Expand Down Expand Up @@ -118,7 +118,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 4,
"id": "82cd9e82",
"metadata": {},
"outputs": [
Expand Down Expand Up @@ -171,7 +171,7 @@
"metadata": {},
"source": [
"<video width=\"75%\" controls autoplay loop>\n",
" <source src=\"https://kerrgeopy.readthedocs.io/en/latest/_downloads/f073aec72d9326b380458508c6f93966/animation1.mp4\" type=\"video/mp4\">\n",
" <source src=\"https://raw.githubusercontent.com/BlackHolePerturbationToolkit/KerrGeoPy/main/docs/source/notebooks/animation1.mp4\" type=\"video/mp4\">\n",
" <p>\n",
"\n",
" [](animation1.mp4)\n",
Expand Down Expand Up @@ -203,15 +203,15 @@
"metadata": {},
"source": [
"<video width=\"75%\" controls autoplay loop>\n",
" <source src=\"https://kerrgeopy.readthedocs.io/en/latest/_downloads/97d0b4212f1e160a4e6b68b2313aeada/animation2.mp4\" type=\"video/mp4\">\n",
" <source src=\"https://raw.githubusercontent.com/BlackHolePerturbationToolkit/KerrGeoPy/main/docs/source/notebooks/animation2.mp4\" type=\"video/mp4\">\n",
" <p>\n",
"\n",
" [](animation2.mp4)\n",
" \n",
" </p>\n",
"</video>\n",
"\n",
"Use the `background_color` option to change the background color of the animation."
"Use the `background_color` option to change the background color of the animation. If necesssary, use `plt.style.use('dark_background')` or `with plt.style.context('dark_background'):` to change the color of the grid and axes to white."
]
},
{
Expand All @@ -221,8 +221,11 @@
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"orbit = kg.StableOrbit(0.999,4,0.1,cos(pi/20))\n",
"orbit.animate(\"animation3.mp4\", 0, 10, background_color=\"black\", elevation=0)"
"with plt.style.context('dark_background'):\n",
" orbit.animate(\"animation3.mp4\",0,10, background_color=\"black\", elevation=10)"
]
},
{
Expand All @@ -231,7 +234,7 @@
"metadata": {},
"source": [
"<video width=\"75%\" controls autoplay loop>\n",
" <source src=\"https://kerrgeopy.readthedocs.io/en/latest/_downloads/3ee6be26a849b6eb69a4edcef89e4207/animation3.mp4\" type=\"video/mp4\">\n",
" <source src=\"https://raw.githubusercontent.com/BlackHolePerturbationToolkit/KerrGeoPy/main/docs/source/notebooks/animation3.mp4\" type=\"video/mp4\">\n",
" <p>\n",
"\n",
" [](animation3.mp4)\n",
Expand All @@ -254,7 +257,8 @@
"orbit = kg.StableOrbit(0.999,3,0.1,cos(pi/20))\n",
"orbit.animate(\"animation4.mp4\", \n",
" elevation_pan = lambda x: 180*x/10,\n",
" roll = lambda x: 180/(1+exp(10-2*x)))"
" roll = lambda x: 180/(1+exp(10-2*x)),\n",
" )"
]
},
{
Expand All @@ -263,12 +267,50 @@
"metadata": {},
"source": [
"<video width=\"75%\" controls autoplay loop>\n",
" <source src=\"https://kerrgeopy.readthedocs.io/en/latest/_downloads/e0f5e701241a3dc4ff27d031b7edad24/animation4.mp4\" type=\"video/mp4\">\n",
" <source src=\"https://raw.githubusercontent.com/BlackHolePerturbationToolkit/KerrGeoPy/main/docs/source/notebooks/animation4.mp4\" type=\"video/mp4\">\n",
" <p>\n",
"\n",
" [](animation4.mp4)\n",
" \n",
" </p>\n",
"</video>\n",
"\n",
"The primary black hole is always placed in the center of the plot, so the axis limits are symmetric about the origin and are the same along all three axes. To zoom in and out of the orbit, use the `axis_limit` option to set the axis limit as a function of Mino time. Set `plot_components` to `True` in order to plot each of the individual components of the trajectory alongside the orbit itself."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f5f218ee",
"metadata": {},
"outputs": [],
"source": [
"orbit = kg.StableOrbit(0.998,2.1,0.8,cos(pi/20))\n",
"t, r, theta, phi = orbit.trajectory()\n",
"\n",
"with plt.style.context('dark_background'):\n",
" orbit.animate(\"animation5.mp4\",0,30,\n",
" elevation = 10,\n",
" azimuthal_pan = lambda x: 360*x/20,\n",
" axis_limit = lambda x: r(x)+1,\n",
" plot_components=True,\n",
" grid=False,\n",
" background_color=\"black\"\n",
" )"
]
},
{
"cell_type": "markdown",
"id": "ab760a07",
"metadata": {},
"source": [
"<video width=\"100%\" controls autoplay loop>\n",
" <source src=\"https://raw.githubusercontent.com/BlackHolePerturbationToolkit/KerrGeoPy/main/docs/source/notebooks/animation5.mp4\" type=\"video/mp4\">\n",
" <p>\n",
"\n",
" [](animation5.mp4)\n",
" \n",
" </p>\n",
"</video>"
]
}
Expand Down
Binary file modified docs/source/notebooks/animation1.mp4
View file
Binary file not shown.
Binary file modified docs/source/notebooks/animation2.mp4
View file
Binary file not shown.
Binary file modified docs/source/notebooks/animation3.mp4
View file
Binary file not shown.
Binary file modified docs/source/notebooks/animation4.mp4
View file
Binary file not shown.
Binary file added docs/source/notebooks/animation5.mp4
View file
Binary file not shown.
6 changes: 3 additions & 3 deletions kerrgeopy/constants.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,7 @@ def stable_polar_roots(a,p,e,x,constants=None):
:param constants: dimensionless constants of motion for the orbit
:type constants: tuple(double, double, double)
:return: tuple of roots in the form :math:`(z_-, z_+)`
:return: tuple of roots :math:`(z_-, z_+)`
:rtype: tuple(double, double, double, double)
"""
if constants is None: constants = constants_of_motion(a,p,e,x)
Expand Down Expand Up @@ -326,7 +326,7 @@ def constants_of_motion(a,p,e,x):
:param x: cosine of the orbital inclination (must satisfy 0 <= x^2 <= 1)
:type x: double
:return: tuple of constants in the form :math:`(E, L, Q)`
:return: tuple of constants of motion :math:`(E, L, Q)`
:rtype: tuple(double, double, double)
"""
E = energy(a,p,e,x)
Expand All @@ -350,7 +350,7 @@ def apex_from_constants(a,E,L,Q):
:return: tuple of orbital parameters :math:`(a,p,e,x)`
:rtype: tuple(double,double,double,double)
"""
# radial polynomial written in terms of z = cos^2(theta)
# Radial polynomial
R = Polynomial([-a**2*Q, 2*L**2+2*Q+2*a**2*E**2-4*a*E*L, a**2*E**2-L**2-Q-a**2, 2, E**2-1])
radial_roots = R.roots()
# numpy returns roots in increasing order
Expand Down
102 changes: 80 additions & 22 deletions kerrgeopy/orbit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@ def trajectory(self,initial_phases=None,distance_units="natural",time_units="nat
:param time_units: units to compute the time component of the trajectory in (options are "natural", "mks", "cgs", and "days"), defaults to "natural"
:type time_units: str, optional
:return: tuple of functions in the form :math:`(t(\lambda), r(\lambda), \theta(\lambda), \phi(\lambda))`
:return: tuple of functions :math:`(t(\lambda), r(\lambda), \theta(\lambda), \phi(\lambda))`
:rtype: tuple(function, function, function, function)
"""
if initial_phases is None: initial_phases = self.initial_phases
Expand Down Expand Up @@ -332,7 +332,7 @@ def is_visible(self,points,elevation,azimuth):
return ((normal_component >= 0) | (np.linalg.norm(projection,axis=1) > event_horizon)) & (np.linalg.norm(points) > event_horizon)

def animate(self,filename,lambda0=0, lambda1=10, elevation=30 ,azimuth=-60, initial_phases=None, grid=True, axes=True, color="red", tau=2, tail_length=5,
lw=2, azimuthal_pan=None, elevation_pan=None, roll=None, speed=1, background_color=None):
lw=2, azimuthal_pan=None, elevation_pan=None, roll=None, speed=1, background_color=None, axis_limit=None, plot_components=False):
r"""
Saves an animation of the orbit as an mp4 file.
Note that this function requires ffmpeg to be installed and may take several minutes to run depending on the length of the animation.
Expand All @@ -355,7 +355,7 @@ def animate(self,filename,lambda0=0, lambda1=10, elevation=30 ,azimuth=-60, init
:type axes: bool, optional
:param color: color of the orbital tail, defaults to "red"
:type color: str, optional
:param tau: time constant for the exponential decay in the opacity of the tail, defaults to infinity
:param tau: time constant for the exponential decay in the opacity of the tail, defaults to 2
:type tau: double, optional
:param tail_length: length of the tail in units of mino time, defaults to 5
:type tail_length: double, optional
Expand All @@ -367,25 +367,53 @@ def animate(self,filename,lambda0=0, lambda1=10, elevation=30 ,azimuth=-60, init
:type elevation_pan: function, optional
:param roll: function defining the roll angle of the camera in degrees as a function of mino time, defaults to None
:type roll: function, optional
:param axis_limit: sets the axis limit as a function of mino time, defaults to None
:type axis_limit: function, optional
:param speed: playback speed of the animation in units of mino time per second (must be a multiple of 1/8), defaults to 1
:type speed: double, optional
:param background_color: color of the background, defaults to None
:type background_color: str, optional
:param plot_components: if true, plots the components of the trajectory in addition to the trajectory itself, defaults to False
:type plot_components: bool, optional
"""
lambda_range = lambda1 - lambda0
point_density = 240 # number of points per unit of mino time
num_pts = int(lambda_range*point_density) # total number of points
time = np.linspace(lambda0,lambda1,num_pts)
speed_multiplier = int(speed*8)
num_frames = int(num_pts/speed_multiplier)
# compute trajectory
t, r, theta, phi = self.trajectory(initial_phases)

fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(projection='3d')
eh = 1+sqrt(1-self.a**2)
fig = plt.figure(figsize=((18,12) if plot_components else (12,12)))
if plot_components:
ax_dict = fig.subplot_mosaic(
"""
OOOOTT
OOOORR
OOOOΘΘ
OOOOΦΦ
""",
per_subplot_kw={"O":{"projection":"3d"},"T":{"facecolor":"none"},"R":{"facecolor":"none"},"Θ":{"facecolor":"none"},"Φ":{"facecolor":"none"}}
)
ax = ax_dict["O"]
ax_dict["T"].set_ylabel("$t$")
ax_dict["R"].set_ylabel("$r$")
ax_dict["Θ"].set_ylabel(r"$\theta$")
ax_dict["Φ"].set_ylabel(r"$\phi$")
t_plot, = ax_dict["T"].plot(time,t(time))
r_plot, = ax_dict["R"].plot(time,r(time))
theta_plot, = ax_dict["Θ"].plot(time,theta(time))
phi_plot, = ax_dict["Φ"].plot(time,phi(time))
# add text with parameters and time
text = ax.text2D(0.05, 0.95, "", transform=ax.transAxes, fontsize=20, bbox=dict(facecolor='none', pad=10.0))

else:
ax = fig.add_subplot(projection='3d')

t, r, theta, phi = self.trajectory(initial_phases)
eh = 1+sqrt(1-self.a**2) # event horizon radius

ax.view_init(elevation,azimuth)
# compute trajectory in cartesian coordinates
trajectory_x = r(time)*sin(theta(time))*cos(phi(time))
trajectory_y = r(time)*sin(theta(time))*sin(phi(time))
trajectory_z = r(time)*cos(theta(time))
Expand All @@ -401,7 +429,7 @@ def animate(self,filename,lambda0=0, lambda1=10, elevation=30 ,azimuth=-60, init
# plot black hole
black_hole_color = "#333" if background_color == "black" else "black"
ax.plot_surface(x_sphere, y_sphere, z_sphere, color=black_hole_color,shade=(background_color == "black"), zorder=0)
# create tail
# create orbital tail
decay = np.flip(0.1+0.9*np.exp(-(time-time[0])/tau)) # exponential decay
tail = Line3DCollection([], color=color, linewidths=lw, zorder=1)
ax.add_collection(tail)
Expand All @@ -420,13 +448,22 @@ def animate(self,filename,lambda0=0, lambda1=10, elevation=30 ,azimuth=-60, init
ax.set_aspect("equal")
# https://matplotlib.org/stable/gallery/mplot3d/projections.html
ax.set_proj_type('ortho')
# set viewing angle
ax.view_init(elevation,azimuth)
# turn off grid and axes if specified
if not grid or background_color is not None: ax.grid(False)
if not axes or background_color is not None: ax.axis("off")

# turn off grid and axes if specified
if not grid: ax.grid(False)
if not axes: ax.axis("off")

# remove margins
fig.tight_layout()

# set background color if specified
if background_color is not None: ax.set_facecolor(background_color)
if background_color is not None:
fig.set_facecolor(background_color)
ax.set_facecolor(background_color)
# make the panes transparent so that the background color shows through the grid
ax.xaxis.pane.fill = False
ax.yaxis.pane.fill = False
ax.zaxis.pane.fill = False

# start progress bar
with tqdm(total=num_frames,ncols=80) as pbar:
Expand All @@ -436,15 +473,21 @@ def draw_frame(i,body,tail):

j = speed_multiplier*i
j0 = max(0,j-tail_length*point_density)
mino_time = time[j]
current_time = time[j]

# update camera angles
updated_azimuth = azimuthal_pan(mino_time) if azimuthal_pan is not None else azimuth
updated_elevation = elevation_pan(mino_time) if elevation_pan is not None else elevation
updated_roll = roll(mino_time) if roll is not None else 0

updated_azimuth = azimuthal_pan(current_time) if azimuthal_pan is not None else azimuth
updated_elevation = elevation_pan(current_time) if elevation_pan is not None else elevation
updated_roll = roll(current_time) if roll is not None else 0
ax.view_init(updated_elevation,updated_azimuth,updated_roll)

# update axis limits
if axis_limit is not None:
updated_limit = axis_limit(current_time)
ax.set_xlim(-updated_limit,updated_limit)
ax.set_ylim(-updated_limit,updated_limit)
ax.set_zlim(-updated_limit,updated_limit)

# filter out points behind the black hole
visible = self.is_visible(trajectory[j0:j],updated_elevation,updated_azimuth)
trajectory_z_visible = trajectory_z[j0:j].copy()
Expand All @@ -455,13 +498,28 @@ def draw_frame(i,body,tail):
# update tail
tail.set_segments(segments)
tail.set_alpha(decay[-(j-j0):])
#tail.set_array(decay[-(j-j0):])
# update body
body._offsets3d = ([trajectory_x[j]],[trajectory_y[j]],[trajectory_z[j]])

# update plots
if plot_components:
t_plot.set_data(time[:j],t(time[:j]))
r_plot.set_data(time[:j],r(time[:j]))
theta_plot.set_data(time[:j],theta(time[:j]))
phi_plot.set_data(time[:j],phi(time[:j]))
# set text
if self.stable:
text.set_text(f"$a = {self.a}\quad p = {self.p}\quad e = {self.e}\quad x = {self.x:.3f}\quad \lambda = {current_time:.2f}$")
else:
text.set_text(f"$a = {self.a}\quad E = {self.E:.3f}\quad L = {self.L:.3f}\quad Q = {self.Q:.3f}\quad \lambda = {current_time:.2f}$")

# save to file
ani = FuncAnimation(fig,draw_frame,num_frames,fargs=(body,tail))
FFwriter = FFMpegWriter(fps=30)
ani.save(filename, writer = FFwriter)
# savefig overrides the facecolor so we need to set it again
if background_color is not None:
ani.save(filename,savefig_kwargs={"facecolor":background_color}, writer = FFwriter)
else:
ani.save(filename, writer = FFwriter)
# close figure so it doesn't show up in notebook
plt.close(fig)
2 changes: 2 additions & 0 deletions kerrgeopy/plunge.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,8 @@ class PlungingOrbit(Orbit):
:ivar E: dimensionless energy
:ivar L: dimensionless angular momentum
:ivar Q: dimensionless Carter constant
:ivar initial_phases: tuple of initial phases :math:`(q^t_0, q^r_0, q^\theta_0, q^\phi_0)`
:ivar stable: boolean indicating whether the orbit is stable
:ivar initial_position: tuple of initial position coordinates :math:`(t_0, r_0, \theta_0, \phi_0)`
:ivar initial_velocity: tuple of initial four-velocity components :math:`(u^t_0, u^r_0, u^\theta_0, u^\phi_0)`
:ivar M: mass of the primary in solar masses
Expand Down
Loading

0 comments on commit 9ebc8dc

Please sign in to comment.