diff --git a/docs/source/images/thumbnail.png b/docs/source/images/thumbnail.png index 7d8b362..9417775 100644 Binary files a/docs/source/images/thumbnail.png and b/docs/source/images/thumbnail.png differ diff --git a/docs/source/index.rst b/docs/source/index.rst index f4e2f8f..2ae2fe4 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -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 .. _Installation: diff --git a/docs/source/notebooks/Graphics.ipynb b/docs/source/notebooks/Graphics.ipynb index af44108..6f16460 100644 --- a/docs/source/notebooks/Graphics.ipynb +++ b/docs/source/notebooks/Graphics.ipynb @@ -26,7 +26,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 16, "id": "7ae99533", "metadata": {}, "outputs": [ @@ -59,7 +59,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 2, "id": "f3351c95", "metadata": {}, "outputs": [ @@ -88,7 +88,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "id": "aa8f032c", "metadata": {}, "outputs": [ @@ -118,7 +118,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 4, "id": "82cd9e82", "metadata": {}, "outputs": [ @@ -171,7 +171,7 @@ "metadata": {}, "source": [ "\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." ] }, { @@ -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)" ] }, { @@ -231,7 +234,7 @@ "metadata": {}, "source": [ "\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": [ + "" ] } diff --git a/docs/source/notebooks/animation1.mp4 b/docs/source/notebooks/animation1.mp4 index 86ed0bd..912bf5c 100644 Binary files a/docs/source/notebooks/animation1.mp4 and b/docs/source/notebooks/animation1.mp4 differ diff --git a/docs/source/notebooks/animation2.mp4 b/docs/source/notebooks/animation2.mp4 index 7338a1e..e42cb6f 100644 Binary files a/docs/source/notebooks/animation2.mp4 and b/docs/source/notebooks/animation2.mp4 differ diff --git a/docs/source/notebooks/animation3.mp4 b/docs/source/notebooks/animation3.mp4 index fa431e9..e9e9584 100644 Binary files a/docs/source/notebooks/animation3.mp4 and b/docs/source/notebooks/animation3.mp4 differ diff --git a/docs/source/notebooks/animation4.mp4 b/docs/source/notebooks/animation4.mp4 index 6ffb965..64d89a3 100644 Binary files a/docs/source/notebooks/animation4.mp4 and b/docs/source/notebooks/animation4.mp4 differ diff --git a/docs/source/notebooks/animation5.mp4 b/docs/source/notebooks/animation5.mp4 new file mode 100644 index 0000000..3f2a9d9 Binary files /dev/null and b/docs/source/notebooks/animation5.mp4 differ diff --git a/kerrgeopy/constants.py b/kerrgeopy/constants.py index b9ad26b..abc4d19 100644 --- a/kerrgeopy/constants.py +++ b/kerrgeopy/constants.py @@ -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) @@ -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) @@ -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 diff --git a/kerrgeopy/orbit.py b/kerrgeopy/orbit.py index 33789d3..9f70497 100644 --- a/kerrgeopy/orbit.py +++ b/kerrgeopy/orbit.py @@ -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 @@ -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. @@ -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 @@ -367,10 +367,14 @@ 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 @@ -378,14 +382,38 @@ def animate(self,filename,lambda0=0, lambda1=10, elevation=30 ,azimuth=-60, init 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)) @@ -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) @@ -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: @@ -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() @@ -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) \ No newline at end of file diff --git a/kerrgeopy/plunge.py b/kerrgeopy/plunge.py index ce9e12d..9b9787f 100644 --- a/kerrgeopy/plunge.py +++ b/kerrgeopy/plunge.py @@ -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 diff --git a/kerrgeopy/stable.py b/kerrgeopy/stable.py index 7b813cc..48f13db 100644 --- a/kerrgeopy/stable.py +++ b/kerrgeopy/stable.py @@ -37,6 +37,7 @@ class StableOrbit(Orbit): :ivar E: dimensionless energy :ivar L: dimensionless angular momentum :ivar Q: dimensionless carter constant + :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 upsilon_r: dimensionless radial orbital frequency in Mino time @@ -92,7 +93,7 @@ def mino_frequencies(self, units="natural"): :param units: units to return the frequencies in (options are "natural", "mks" and "cgs"), defaults to "natural" :type units: str, optional - :return: tuple of orbital frequencies in the form :math:`(\Upsilon_r, \Upsilon_\theta, \Upsilon_\phi, \Gamma)` + :return: tuple of orbital frequencies :math:`(\Upsilon_r, \Upsilon_\theta, \Upsilon_\phi, \Gamma)` :rtype: tuple(double, double, double, double) """ upsilon_r, upsilon_theta, upsilon_phi, gamma = self.upsilon_r, self.upsilon_theta, self.upsilon_phi, self.gamma @@ -113,7 +114,7 @@ def fundamental_frequencies(self, units="natural"): :param units: units to return the frequencies in (options are "natural", "mks", "cgs" and "mHz"), defaults to "natural" :type units: str, optional - :return: tuple of orbital frequencies in the form :math:`(\Omega_r, \Omega_\theta, \Omega_\phi)` + :return: tuple of orbital frequencies :math:`(\Omega_r, \Omega_\theta, \Omega_\phi)` :rtype: tuple(double, double, double) """ upsilon_r, upsilon_theta, upsilon_phi, gamma = self.upsilon_r, self.upsilon_theta, self.upsilon_phi, self.gamma @@ -140,7 +141,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 @@ -161,7 +162,7 @@ def radial_solutions(a,constants,radial_roots): :param radial_roots: tuple of roots :math:`(r_1,r_2,r_3,r_4)`. Assumes that motion is between :math:`r_1` and :math:`r_2` and that roots are otherwise in decreasing order. :type radial_roots: tuple(double,double,double,double) - :return: tuple of functions in the form :math:`(r, t^{(r)}, \phi^{(r)})` + :return: tuple of functions :math:`(r, t^{(r)}, \phi^{(r)})` :rtype: tuple(function, function, function) """ E, L, Q = constants @@ -229,7 +230,7 @@ def polar_solutions(a,constants,polar_roots): :param polar_roots: tuple of roots :math:`(z_-,z_+)` :type polar_roots: tuple(double,double) - :return: tuple of functions in the form :math:`(\theta, t^{(\theta)}, \phi^{(\theta)})` + :return: tuple of functions :math:`(\theta, t^{(\theta)}, \phi^{(\theta)})` :rtype: tuple(function, function, function) """ E, L, Q = constants @@ -283,7 +284,7 @@ def stable_trajectory(a,p,e,x,initial_phases=(0,0,0,0),M=None,distance_units="na :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 ((distance_units != "natural") or (time_units != "natural")) and M is None: raise ValueError("M must be specified to convert to physical units") diff --git a/setup.py b/setup.py index d1eb288..c7c0ae4 100644 --- a/setup.py +++ b/setup.py @@ -13,5 +13,5 @@ "License :: OSI Approved :: GNU General Public License (GPL)", "Natural Language :: English", ), - install_requires=["scipy>=1.8","numpy","matplotlib>=3.3","tqdm"] + install_requires=["scipy>=1.8","numpy","matplotlib>=3.7","tqdm"] ) \ No newline at end of file diff --git a/tests/README.txt b/tests/README.txt index e38e925..fc11fc6 100644 --- a/tests/README.txt +++ b/tests/README.txt @@ -5,21 +5,21 @@ Below is a list of data files used by each test suite: test_constants.py: const_parameters.txt - orbital parameters (a,p,e,x) -mathematica_const_output.txt - (E,L,Q) computed using KerrGeoConstantsOfMotion[a,p,e,x] for each orbit in const_values.txt +mathematica_const_output.txt - (E,L,Q) computed using KerrGeoConstantsOfMotion[a,p,e,x] for each orbit in const_parameters.txt separatrix_parameters.txt - orbital parameters (a,p,e,x) -mathematica_separatrix_output.txt - output of KerrGeoSeparatrix[a,p,e,x] for each orbit in separatrix_values.txt +mathematica_separatrix_output.txt - output of KerrGeoSeparatrix[a,p,e,x] for each orbit in separatrix_parameters.txt test_frequencies.py: freq_parameters.txt - orbital parameters (a,p,e,x) -mathematica_freq_output.txt - (upsilon_r, upsilon_theta, upsilon_phi, gamma) computed using KerrGeoFrequencies[a,p,e,x] for each orbit in freq_values.txt +mathematica_freq_output.txt - (upsilon_r, upsilon_theta, upsilon_phi, gamma) computed using KerrGeoFrequencies[a,p,e,x] for each orbit in freq_parameters.txt test_stable_solutions.py: stable_orbit_parameters.txt - orbital parameters (a,p,e,x) stable_orbit_times.txt - list of mino time values to test -stable_orbits/trajectory{i}.txt - (t, r, theta, phi) evaluated at each time from stable_orbit_times.txt for the i-th orbit defined in stable_orbit_values.txt -stable_solutions/trajectory{i}.txt - (t_r, t_theta, phi_r, phi_theta) evaluated at each time from stable_orbit_times.txt for the i-th orbit defined in stable_orbit_values.txt +stable_orbits/trajectory{i}.txt - (t, r, theta, phi) evaluated at each time from stable_orbit_times.txt for the i-th orbit defined in stable_orbit_parameters.txt +stable_solutions/trajectory{i}.txt - (t_r, t_theta, phi_r, phi_theta) evaluated at each time from stable_orbit_times.txt for the i-th orbit defined in stable_orbit_parameters.txt test_plunging_solutions.py: @@ -37,4 +37,4 @@ stable_orbit_parameters.txt - orbital parameters (a,p,e,x) stable_orbit_times.txt - list of mino time values to test plunging_orbit_parameters_real.txt - list of orbital parameters (a,E,L,Q) for which the radial polynomial has all real roots plunging_orbit_parameters_complex.txt - list of orbital parameters (a,E,L,Q) for which the radial polynomial has complex roots -four_velocity/trajectory{i}.txt - (u_t, u_r, u_theta, u_phi) evaluated at each time from stable_orbit_times.txt for the i-th orbit in stable_orbit_values.txt \ No newline at end of file +four_velocity/trajectory{i}.txt - (u_t, u_r, u_theta, u_phi) evaluated at each time from stable_orbit_times.txt for the i-th orbit in stable_orbit_parameters.txt \ No newline at end of file