Saved in:
| Main Authors: | , |
|---|---|
| Format: | Recurso digital |
| Language: | |
| Published: |
Zenodo
2026
|
| Online Access: | https://doi.org/10.5281/zenodo.19369572 |
| Tags: |
Add Tag
No Tags, Be the first to tag this record!
|
Table of Contents:
- <p>This dataset contains angular momentum deposition and torque excitation profiles from a planet embedded in an astrophysical disc as a function of planetary mass, cooling timescale and viscosity.</p> <p>Any use of this data set should cite both the zenodo record and the associated paper Cordwell & Rafikov (2026).</p> <p>Questions or issues using this dataset should be sent to ajc356@cam.ac.uk</p> <h3>Data Format</h3> <p>Physical descriptions of each term and calculation method are included section 2.1 of the associated paper. </p> <p>Each .zip file contains data from our full set of simulations using either the Bessel-type potential or the second order (plummer) potential. Inside the .zip file there are 375 python pickle (.p) files named for the simulation parameters used in their calculation. Each pickle file contains a single dictionary which records both simulation setup details (keys: 'time', 'h_p', 'm_p', 'sig_slope', 't_slope', 'alpha', 'beta', 'final_time') as floats, or simulation radial profile outputs as numpy arrays (keys: 'F_dep', 'F_wave', 'G_background', 'G_total', 'G_wave', 'R', 'T', 'dF_dep', 'dF_wave', 'dG', 'dTdR', 'dsigmadt', 'nu', 'sigma', 'u_phi', 'u_r'). All radial profile outputs have the same shape as 'R'.</p> <p>See the example code below for further usage detail. </p> <p> </p> <h2>Example Code</h2> <p>The below example code will reprodice a combination of Figures 1 and 2 from the associated paper. It shows both how to get the filename for selected physical parameters and how to extract angular momentum deposition and torque excitation profiles.</p> <p> </p> <pre><code>import numpy as np import matplotlib.pyplot as plt import scipy.integrate import pickle import string from matplotlib import cm from scipy.signal import savgol_filter INNER_DAMP_R = 0.28 OUTER_DAMP_R = 3.6 ## define the grid MP_GRID = [2.5, 1, 2.5e-1, 0.1, 2.5e-2] # BETA_GRID = [np.inf, 100, 10, 3, 1, 0.3, 0.1, 0.01, 0] BETA_DISPLAY_GRID = [np.inf, 10, 1, 0.1, 0] # 0 = strict iso, inf = adiabatic ALPHA_GRID = [0, 1e-4, 1e-3,1e-2, 1e-1] HP_GRID = [0.05] GAMMA = 1.4 FILE_FORMAT = 'output_second/abrun_mp_{mp_mth:.2e}_hp_{hp:.2f}_alpha_{alpha:.1e}_beta_{beta:.1e}_orbit_11.p' column_width = 4 def get_color_alpha(alpha, max_value=-1, min_value=-5): if alpha == 0: return 'black' if alpha == np.inf: alpha = 1e4 alpha = np.log10(alpha) return cm.viridis((alpha - min_value)/(max_value - min_value)) def get_color(beta, max_value=2.1, min_value=-3.5): if beta == 0: return 'black' if beta == np.inf: beta = 1e8 beta = np.log10(beta) return cm.plasma((beta - min_value)/(max_value - min_value)) filter_window = 12 # Half a scale height filter_number = 3 def moving_average(a): return savgol_filter(a, filter_window, filter_number) for mp in [2.5, 1, 2.5e-2]: alpha = 0 hp = 0.05 h_p = 0.05 dt_scaling = mp**-2 * h_p**-3 fig, axs = plt.subplots( 2, 2, sharex=True, sharey='row', figsize=(column_width, column_width)) h_p = 0.05 hp = 0.05 orbit_number = 11 for beta in BETA_DISPLAY_GRID: d = FILE_FORMAT.format( mp_mth=mp, hp = hp, alpha=alpha, beta=beta) try: file = open(d, 'rb') details = pickle.load(file) file.close() axs[0, 1].plot(details['R'][300:-100], details['dTdR'][300:-100] * dt_scaling, color=get_color(beta), label="$\\beta = {}$".format(beta)) axs[1, 1].plot(details['R'][300:-100], moving_average(np.gradient(details['F_dep'], details['R'])[300:-100]) * dt_scaling, color=get_color(beta), label="$\\beta = {}$".format(beta)) except FileNotFoundError: print('File not found,', d) continue beta = 0 for alpha in ALPHA_GRID: d = FILE_FORMAT.format( mp_mth=mp, hp = hp, alpha=alpha, beta=beta) try: file = open(d, 'rb') details = pickle.load(file) file.close() if alpha != 0: axs[0, 0].plot(details['R'][300:-100], details['dTdR'][300:-100] * dt_scaling, color=get_color_alpha(alpha), label="$\\alpha = 10^{" + "{}".format(np.log10(alpha)) + "}$") axs[1, 0].plot(details['R'][300:-100], moving_average(np.gradient(details['F_dep'], details['R'])[300:-100]) * dt_scaling, color=get_color_alpha(alpha), label="$\\alpha = 10^{" + "{}".format(np.log10(alpha)) + "}$") else: axs[0, 0].plot(details['R'][300:-100], details['dTdR'][300:-100] * dt_scaling, color=get_color_alpha(alpha), label="$\\alpha = 0$") axs[1, 0].plot(details['R'][300:-100], moving_average(np.gradient(details['F_dep'], details['R'])[300:-100]) * dt_scaling, color=get_color_alpha(alpha), label="$\\alpha = 0$") except FileNotFoundError: print('File not found,', d) continue axs[0, 1].text(0.02, 0.92, "(b) $\\alpha = 0$", transform=axs[0, 1].transAxes) axs[0, 0].text(0.02, 0.92, "(a) Isothermal", transform=axs[0, 0].transAxes) axs[1, 1].text(0.02, 0.92, "(d) $\\alpha = 0$", transform=axs[1, 1].transAxes) axs[1, 0].text(0.02, 0.92, "(c) Isothermal", transform=axs[1, 0].transAxes) axs[0, 0].set_ylabel('$dT/dR / F_{j, 0}$') axs[1, 0].set_ylabel('$dF_{\mathrm{dep}}/dR / F_{j, 0}$') axs[-1, 0].set_xlabel('$R/R_{\mathrm{p}}$') axs[-1, 1].set_xlabel('$R/R_{\mathrm{p}}$') plt.subplots_adjust(hspace=0, wspace=0) plt.xlim([0.7, 1.3]) axs[0, 0].set_ylim([-5, 6]) axs[1, 0].set_ylim([-3, 3]) axs[1, 0].legend(loc='lower right', ncol=2, bbox_to_anchor=(1, -0.6)) axs[1, 1].legend(loc='lower right', ncol=2, bbox_to_anchor=(1, -0.6)) plt.show()<br><br></code></pre>