diff --git a/docs/settings_examples/plot_LBV_flag.py b/docs/settings_examples/plot_LBV_flag.py index 8b37b8411..016c58f9d 100644 --- a/docs/settings_examples/plot_LBV_flag.py +++ b/docs/settings_examples/plot_LBV_flag.py @@ -4,7 +4,12 @@ This example shows the effect of the ``LBV_flag`` on the evolution of massive stars. -The ``LBV_flag`` controls which prescription to use for LBV-like mass loss for stars that are above the Humphreys-Davidson limit. The plots below show the evolution of a few single stars and highlight the LBV regime in which we'd expect to see differences with this flag. +The ``LBV_flag`` controls which prescription to use for LBV-like mass loss for +stars that are above the Humphreys-Davidson limit. The first figures show the +evolution of a few single stars on the HR diagram and highlight the LBV regime +where we'd expect to see differences with this flag. The final figure compares +the final primary mass over a wider initial-mass grid, showing how the LBV wind +prescription changes the amount of mass retained by massive single stars. """ #---------------------------------------------------------------------------------- @@ -22,6 +27,7 @@ #---------------------------------------------------------------------------------- import numpy as np +import pandas as pd import astropy.units as u import astropy.constants as consts @@ -35,25 +41,26 @@ def LBV_limit(T_eff): # create a small grid of single stars -masses = [30, 33, 36, 40, 45, 50] +hr_masses = [30, 33, 36, 40, 45, 50] grid = InitialBinaryTable.InitialBinaries( - m1=masses, - m2=np.zeros(len(masses)), - porb=np.ones(len(masses))*-1.0, - ecc=np.ones(len(masses))*0.0, - tphysf=np.ones(len(masses))*13700.0, - kstar1=np.ones(len(masses)), - kstar2=np.ones(len(masses)), - metallicity=np.ones(len(masses))*0.014*0.01 + m1=hr_masses, + m2=np.zeros(len(hr_masses)), + porb=np.ones(len(hr_masses))*-1.0, + ecc=np.ones(len(hr_masses))*0.0, + tphysf=np.ones(len(hr_masses))*13700.0, + kstar1=np.ones(len(hr_masses)), + kstar2=np.ones(len(hr_masses)), + metallicity=np.ones(len(hr_masses))*0.014*0.01 ) # loop over different flag choices for LBV_flag, label in zip([0, 1, 2], ["No LBV mass loss", "LBV mass loss following Hurley+2000", "LBV mass loss following Belczynski+2008"]): # evolve with updated BSEDict - BSEDict["LBV_flag"] = LBV_flag + settings = BSEDict.copy() + settings["LBV_flag"] = LBV_flag bpp, bcm, initC, kick_info = Evolve.evolve( - initialbinarytable=grid, BSEDict=BSEDict, + initialbinarytable=grid, BSEDict=settings, dtp=0 ) @@ -72,7 +79,7 @@ def LBV_limit(T_eff): ax.plot(log_teff[mask], log_lum[mask], color='grey', alpha=0.5, lw=0.5, zorder=-1) # annotate the first point with the initial mass - ax.annotate(f'{bcm[bcm["bin_num"] == bin_num]["mass_1"].values[0]:.0f} M$_\odot$', (log_teff[0], log_lum[0]), + ax.annotate(f'{bcm[bcm["bin_num"] == bin_num]["mass_1"].values[0]:.0f} M$_\\odot$', (log_teff[0], log_lum[0]), fontsize=14, ha="right", color='black', va='top') # annotate the various kstar types @@ -101,4 +108,82 @@ def LBV_limit(T_eff): ax.invert_xaxis() - plt.show() \ No newline at end of file + plt.show() + + +# compare the final masses over a wider initial-mass grid +LBV_flag_labels = { + 0: "0: no LBV winds", + 1: "1: Hurley+2000", + 2: "2: Belczynski+2008", +} +LBV_flag_values = list(LBV_flag_labels) +mass_grid = np.arange(25.0, 81.0, 1.0) + + +def make_final_mass_grid(): + """Build the massive single-star grid used for the final-mass comparison.""" + n_systems = len(mass_grid) + return InitialBinaryTable.InitialBinaries( + m1=mass_grid, + m2=np.zeros(n_systems), + porb=np.ones(n_systems) * -1.0, + ecc=np.zeros(n_systems), + tphysf=np.ones(n_systems) * 13700.0, + kstar1=np.ones(n_systems), + kstar2=np.zeros(n_systems), + metallicity=np.ones(n_systems) * 0.014, + ) + + +def summarize_final_masses(bcm, LBV_flag): + """Collect final primary masses for one LBV_flag value.""" + final_bcm = bcm.sort_values(["bin_num", "tphys"]).groupby("bin_num").last() + + summary = pd.DataFrame(index=final_bcm.index) + summary["LBV_flag"] = LBV_flag + summary["initial_mass_1"] = mass_grid[summary.index.astype(int)] + summary["final_mass_1"] = final_bcm["mass_1"] + return summary.reset_index(drop=True) + + +def evolve_final_mass_grid(final_mass_grid): + """Run the same single-star grid for each LBV_flag value.""" + summaries = [] + + for LBV_flag in LBV_flag_values: + settings = BSEDict.copy() + settings["LBV_flag"] = LBV_flag + settings["random_seed"] = 1 + + bpp, bcm, initC, kick_info = Evolve.evolve( + initialbinarytable=final_mass_grid, + BSEDict=settings, + dtp=0, + ) + summaries.append(summarize_final_masses(bcm, LBV_flag)) + + return pd.concat(summaries, ignore_index=True) + + +final_mass_grid = make_final_mass_grid() +final_mass_summaries = evolve_final_mass_grid(final_mass_grid) + +fig, ax = plt.subplots(figsize=(8.2, 6.0)) + +for LBV_flag in LBV_flag_values: + subset = final_mass_summaries[final_mass_summaries["LBV_flag"].eq(LBV_flag)] + ax.plot( + subset["initial_mass_1"], + subset["final_mass_1"], + label=LBV_flag_labels[LBV_flag], + ) + +ax.set_xlabel("Initial primary mass [$M_\\odot$]") +ax.set_ylabel("Final primary mass [$M_\\odot$]", fontsize=14) +ax.set_title("Effect of LBV_flag on final primary mass", fontsize=15) +ax.legend(fontsize=8) +ax.grid(True, alpha=0.3) + +plt.tight_layout() +plt.show() diff --git a/docs/settings_examples/plot_alpha1_lambdaf.py b/docs/settings_examples/plot_alpha1_lambdaf.py new file mode 100644 index 000000000..8a5e976ed --- /dev/null +++ b/docs/settings_examples/plot_alpha1_lambdaf.py @@ -0,0 +1,340 @@ +""" +alpha1 and lambdaf +================== + +This example compares two common-envelope parameters using real COSMIC +separation histories. ``alpha1`` changes the efficiency with which orbital +energy ejects the common envelope. ``lambdaf`` changes the envelope +binding-energy lambda prescription; the negative values used here select fixed +lambda values. + +Each panel holds the initial binary fixed and changes only the parameter shown +in that panel. Open circles mark surviving binaries, x markers mark mergers, +and triangles mark disruptions at the final real COSMIC output point shown. +Tracks that overlap exactly are given tiny horizontal plotting offsets for +visibility only; the COSMIC output times and separations are unchanged. +""" + +import copy +import sys +import warnings + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +from cosmic.evolve import Evolve +from cosmic.sample import InitialBinaryTable + +sys.path.append("..") +import generate_default_bsedict + +BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True) + + +EXAMPLES = [ + { + "title": "alpha1: CE efficiency", + "flag_key": "alpha1", + "values": [0.2, 1.0, 5.0], + "labels": { + 0.2: r"$\alpha$ = 0.2", + 1.0: r"$\alpha$ = 1.0", + 5.0: r"$\alpha$ = 5.0", + }, + "initial": { + "m1": 8.0, + "m2": 2.8, + "porb": 1584.893192461114, + "tphysf": 70.0, + }, + "plot_end_pad": 0.35, + }, + { + "title": "lambdaf: fixed envelope lambda", + "flag_key": "lambdaf", + "values": [-0.1, -0.5, -1.0], + "labels": { + -0.1: r"fixed $\lambda$ = 0.1", + -0.5: r"fixed $\lambda$ = 0.5", + -1.0: r"fixed $\lambda$ = 1.0", + }, + "initial": { + "m1": 13.0, + "m2": 3.25, + "porb": 2238.72113856834, + "tphysf": 50.0, + }, + "plot_end_pad": 0.5, + }, +] + +LINE_WIDTHS = { + "alpha1": { + 0.2: 3.0, + 1.0: 2.2, + 5.0: 1.4, + }, + "lambdaf": { + -0.1: 3.0, + -0.5: 2.2, + -1.0: 1.4, + }, +} + +LINE_ZORDERS = { + "alpha1": { + 0.2: 3, + 1.0: 4, + 5.0: 5, + }, + "lambdaf": { + -0.1: 3, + -0.5: 4, + -1.0: 5, + }, +} + +DISPLAY_X_OFFSETS = { + "alpha1": { + 0.2: 0.00, + 1.0: 0.15, + 5.0: 0.30, + }, + "lambdaf": { + -0.1: 0.00, + -0.5: 0.15, + -1.0: 0.30, + }, +} + + +def make_binary(initial): + """Build a one-row InitialBinaryTable.""" + return InitialBinaryTable.InitialBinaries( + m1=np.array([initial["m1"]]), + m2=np.array([initial["m2"]]), + porb=np.array([initial["porb"]]), + ecc=np.zeros(1), + tphysf=np.ones(1) * initial["tphysf"], + kstar1=np.ones(1), + kstar2=np.ones(1), + metallicity=np.ones(1) * 0.014, + ) + + +def set_parameter(settings, key, value): + """Set one CE parameter while preserving COSMIC's expected type.""" + if key == "alpha1": + settings[key] = [value, value] + else: + settings[key] = value + + +def evolve_one(example, value): + """Evolve one binary while changing only the selected CE parameter.""" + settings = copy.deepcopy(BSEDict) + set_parameter(settings, example["flag_key"], value) + + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + message="At least one of your initial binaries is starting in Roche Lobe Overflow:*", + category=UserWarning, + ) + bpp, bcm, initC, kick_info = Evolve.evolve( + initialbinarytable=make_binary(example["initial"]), + BSEDict=settings, + dtp=0.5, + randomseed=1, + ) + + return bpp, bcm + + +def first_ce_time(bpp): + """Return the first confirmed common-envelope event time.""" + ce_rows = bpp[bpp["evol_type"].eq(7)] + if ce_rows.empty: + raise RuntimeError("Selected binary does not enter common-envelope evolution.") + return ce_rows["tphys"].iloc[0] + + +def terminal_time(bpp): + """Return the first terminal-event time, or NaN if none is present.""" + terminal_rows = bpp[bpp["evol_type"].eq(6)] + if terminal_rows.empty: + return np.nan + return terminal_rows["tphys"].iloc[0] + + +def final_status(bpp, bcm): + """Read the final outcome from COSMIC outputs.""" + if np.isfinite(terminal_time(bpp)): + return "merges" + + final_row = bcm.sort_values("tphys").iloc[-1] + if final_row["bin_state"] == 0 and final_row["sep"] > 0: + return "survives" + if final_row["bin_state"] == 1: + return "merges" + return "disrupts" + + +def separation_track(bpp, bcm): + """Combine sampled bcm rows with bpp event rows into a real separation track.""" + sampled_rows = bcm[["tphys", "sep"]].copy() + sampled_rows["source_order"] = 0 + + event_rows = bpp[["tphys", "sep"]].copy() + event_rows["source_order"] = np.arange(len(event_rows)) + 1 + + track = pd.concat([sampled_rows, event_rows], ignore_index=True) + track = track[np.isfinite(track["tphys"]) & np.isfinite(track["sep"])] + track = track[track["sep"].gt(0)] + track = track.drop_duplicates(subset=["tphys", "sep"], keep="first") + track = track.sort_values(["tphys", "source_order"]).reset_index(drop=True) + + stop_time = terminal_time(bpp) + if np.isfinite(stop_time): + track = track[track["tphys"].le(stop_time)].copy() + + assert len(track) > 0 + assert np.all(np.isfinite(track["tphys"])) + assert np.all(np.isfinite(track["sep"])) + assert np.all(track["sep"] > 0) + assert np.all(np.diff(track["tphys"]) >= 0) + + return track + + +def mark_outcome(ax, track, status, color, x_offset=0.0, marker_scale=1.0, zorder=12): + """Place the endpoint marker at the final real plotted point.""" + final = track.iloc[-1] + marker_time = final["tphys"] + x_offset + if status == "survives": + ax.scatter( + marker_time, + final["sep"], + marker="o", + s=52 * marker_scale, + facecolor="white", + edgecolor=color, + linewidths=1.8, + zorder=zorder, + ) + elif status == "merges": + ax.scatter( + marker_time, + final["sep"], + marker="x", + s=62 * marker_scale, + color=color, + linewidths=2.1, + zorder=zorder, + ) + else: + ax.scatter( + marker_time, + final["sep"], + marker="^", + s=62 * marker_scale, + facecolor="white", + edgecolor=color, + linewidths=1.8, + zorder=zorder, + ) + + +fig, axes = plt.subplots( + 1, + 2, + figsize=(12, 5), + constrained_layout=True, +) + +for ax, example in zip(axes, EXAMPLES): + colors = plt.get_cmap("tab10").colors[: len(example["values"])] + records = [] + + for value, color in zip(example["values"], colors): + bpp, bcm = evolve_one(example, value) + ce_time = first_ce_time(bpp) + track = separation_track(bpp, bcm) + status = final_status(bpp, bcm) + records.append( + { + "value": value, + "label": example["labels"][value], + "color": color, + "track": track, + "status": status, + "ce_time": ce_time, + } + ) + + ce_times = np.array([record["ce_time"] for record in records]) + if not np.allclose(ce_times, ce_times[0]): + raise RuntimeError("Selected runs do not share the same CE onset time.") + ce_time = ce_times[0] + + for record in records: + linewidth = LINE_WIDTHS[example["flag_key"]][record["value"]] + zorder = LINE_ZORDERS[example["flag_key"]][record["value"]] + display_offset = DISPLAY_X_OFFSETS[example["flag_key"]][record["value"]] + marker_offset = display_offset + marker_scale = 1.0 + marker_zorder = 12 + zorder + if example["flag_key"] == "lambdaf" and record["status"] == "merges": + marker_offset += {-0.1: -0.06, -0.5: 0.06}.get(record["value"], 0.0) + marker_scale = {-0.1: 1.15, -0.5: 0.95}.get(record["value"], 1.0) + ax.plot( + record["track"]["tphys"] + display_offset, + record["track"]["sep"], + color=record["color"], + lw=linewidth, + linestyle="-", + marker=None, + label=record["label"], + zorder=zorder, + ) + mark_outcome( + ax, + record["track"], + record["status"], + record["color"], + x_offset=marker_offset, + marker_scale=marker_scale, + zorder=marker_zorder, + ) + + combined = pd.concat([record["track"] for record in records], ignore_index=True) + max_time = max(record["track"]["tphys"].max() for record in records) + ax.set_xlim(0, max_time + example["plot_end_pad"]) + ax.set_ylim(combined["sep"].min() * 0.85, combined["sep"].max() * 1.35) + ax.set_yscale("log") + ax.axvline(ce_time, color="0.45", lw=0.9, zorder=0) + x0, x1 = ax.get_xlim() + ax.text( + ce_time + 0.028 * (x1 - x0), + 0.93, + "CE onset", + transform=ax.get_xaxis_transform(), + ha="left", + va="top", + fontsize=7, + color="0.25", + ) + ax.set_title(example["title"], fontsize=11) + ax.set_xlabel("Time [Myr]") + ax.set_ylabel(r"Orbital separation [R$_\odot$]") + ax.grid(True, which="both", alpha=0.25) + ax.legend( + loc="best", + fontsize=8, + frameon=True, + title=None, + handlelength=1.8, + ) + +plt.show() diff --git a/docs/settings_examples/plot_ce_flags.py b/docs/settings_examples/plot_ce_flags.py new file mode 100644 index 000000000..d729d1852 --- /dev/null +++ b/docs/settings_examples/plot_ce_flags.py @@ -0,0 +1,348 @@ +""" +ceflag and cemergeflag +====================== + +This example compares two common-envelope-related COSMIC flags using one +representative binary per flag. Each panel evolves the same initial binary +multiple times while changing only the flag shown in that panel, then plots +orbital period as a function of time. + +Single-binary examples are useful for these flags because the important +behavior often happens at a boundary between common-envelope survival and +merger. + +Each panel varies one COSMIC flag while holding the initial binary fixed. +Tracks show orbital-period evolution; endpoint markers indicate whether the +binary ultimately survives or terminates through merger/disruption. The time +axis is limited to the interaction window where the flag differences are most +visible. +""" + +import copy +import sys +import warnings + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from matplotlib.lines import Line2D + +from cosmic.evolve import Evolve +from cosmic.sample import InitialBinaryTable + +sys.path.append("..") +import generate_default_bsedict + +BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True) + + +EXAMPLES = [ + { + "title": "ceflag: CE orbital energy", + "flag_key": "ceflag", + "flags": { + 0: "0: core", + 1: "1: total (default)", + }, + "initial": { + "m1": 32.0, + "m2": 8.0, + "porb": 200.0, + "kstar1": 1, + "kstar2": 1, + "tphysf": 300.0, + }, + "annotation": "same binary, different CE energy prescription", + "plot_end": 7.0, + }, + { + "title": "cemergeflag: CE merger treatment", + "flag_key": "cemergeflag", + "flags": { + 0: "0: alternate", + 1: "1: default", + }, + "initial": { + "m1": 8.0, + "m2": 8.0, + "porb": 1500.0, + "kstar1": 1, + "kstar2": 1, + "tphysf": 500.0, + }, + "annotation": "0 permits CE evolution; 1 forces merger for this system", + "plot_end": 45.0, + }, +] + + +def make_binary(initial): + """Build one InitialBinaryTable row from an example dictionary.""" + return InitialBinaryTable.InitialBinaries( + m1=np.array([initial["m1"]]), + m2=np.array([initial["m2"]]), + porb=np.array([initial["porb"]]), + ecc=np.zeros(1), + tphysf=np.ones(1) * initial["tphysf"], + kstar1=np.ones(1) * initial["kstar1"], + kstar2=np.ones(1) * initial["kstar2"], + metallicity=np.ones(1) * 0.014, + ) + + +def evolve_example(example, flag_value): + """Evolve one example binary for one flag value.""" + settings = copy.deepcopy(BSEDict) + settings[example["flag_key"]] = flag_value + + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + message="At least one of your initial binaries is starting in Roche Lobe Overflow:*", + category=UserWarning, + ) + bpp, bcm, initC, _ = Evolve.evolve( + initialbinarytable=make_binary(example["initial"]), + BSEDict=settings, + dtp=0.2, + randomseed=1, + ) + + return bpp, bcm + + +def orbital_period_track(bpp, bcm): + """Combine sampled evolution rows with event rows into one real track.""" + sampled_rows = bcm[["tphys", "porb", "bin_state"]].copy() + sampled_rows["source_order"] = 0 + + event_rows = bpp[["tphys", "porb"]].copy() + event_rows["bin_state"] = np.nan + event_rows["source_order"] = np.arange(len(event_rows)) + 1 + + track = pd.concat([sampled_rows, event_rows], ignore_index=True) + track = track[track["porb"].gt(0)].drop_duplicates( + subset=["tphys", "porb", "bin_state"], keep="first" + ) + track = track.sort_values(["tphys", "source_order"]).reset_index(drop=True) + + time = track["tphys"].to_numpy() + period = track["porb"].to_numpy() + if not np.all(np.diff(time) >= 0): + raise RuntimeError("Orbital-period track is not time sorted.") + if not np.all(np.isfinite(period)) or not np.all(period > 0): + raise RuntimeError("Orbital-period track contains invalid period values.") + + return track + + +def final_status(bcm): + """Return a short final-status string.""" + final_row = bcm.sort_values("tphys").iloc[-1] + if final_row["bin_state"] == 0 and final_row["porb"] > 0: + return "survives" + if final_row["bin_state"] == 1: + return "merges" + return "disrupts" + + +def plotted_track(track, plot_end): + """Return only real COSMIC rows inside the panel time window.""" + if track.empty: + return track + + plot_track = track[track["tphys"].le(plot_end)].copy() + if plot_track.empty: + plot_track = track.iloc[[0]].copy() + + return plot_track + + +def mark_final_outcome( + ax, + track, + status, + color, + marker_size=70, + zorder=5, +): + """Mark the plotted endpoint with a consistent simple symbol.""" + if track.empty: + return + + final_point = track.iloc[-1] + if status == "survives": + ax.scatter( + final_point["tphys"], + final_point["porb"], + marker="o", + s=marker_size * 0.75, + facecolor="white", + edgecolor=color, + linewidths=1.8, + zorder=zorder, + ) + elif status == "merges": + ax.scatter( + final_point["tphys"], + final_point["porb"], + marker="x", + s=marker_size, + linewidths=2.2, + color=color, + zorder=zorder, + ) + else: + ax.scatter( + final_point["tphys"], + final_point["porb"], + marker="^", + s=marker_size, + facecolor="white", + edgecolor=color, + linewidths=2.0, + zorder=zorder, + ) + + +def first_event_time(bpp, event_type): + """Return the first time for an event type, or NaN if absent.""" + rows = bpp[bpp["evol_type"].eq(event_type)] + if rows.empty: + return np.nan + return rows["tphys"].iloc[0] + + +def first_rlo_time(bpp): + """Return the first recorded RLO time, or NaN if absent.""" + rrlo = bpp[["RRLO_1", "RRLO_2"]].max(axis=1) + rows = bpp[rrlo.ge(1.0)] + if rows.empty: + return np.nan + return rows["tphys"].iloc[0] + + +def primary_event(bpp): + """Return the single event label/time most relevant to the panel.""" + ce_time = first_event_time(bpp, 7) + if np.isfinite(ce_time): + return "CE/RLO", ce_time + + rlo_time = first_rlo_time(bpp) + if np.isfinite(rlo_time): + return "CE/RLO", rlo_time + + return "", np.nan + + +def annotate_event(ax, label, event_time): + """Draw one readable event marker.""" + if not label or not np.isfinite(event_time): + return + + ax.axvline(event_time, color="0.35", ls=":", lw=1.4, zorder=4) + + +def line_width(flag_value): + """Use a consistent width for the two CE-flag tracks.""" + return 2.0 + + +def endpoint_marker_size(flag_value): + """Use a consistent endpoint marker size.""" + return 70 + + +fig, axes = plt.subplots(1, 2, figsize=(13, 4.8), constrained_layout=True) +axes = axes.ravel() + +for ax, example in zip(axes, EXAMPLES): + records = [] + reference_record = None + + for flag_value, label in example["flags"].items(): + bpp, bcm = evolve_example(example, flag_value) + track = orbital_period_track(bpp, bcm) + status = final_status(bcm) + record = { + "flag_value": flag_value, + "label": label, + "bpp": bpp, + "bcm": bcm, + "track": track, + "status": status, + } + records.append(record) + + if flag_value == max(example["flags"]): + reference_record = record + + event_label, event_time = primary_event(reference_record["bpp"]) + all_plotted_tracks = [] + plot_end = example["plot_end"] + for plot_index, record in enumerate(records): + track = plotted_track(record["track"], plot_end) + all_plotted_tracks.append(track) + linewidth = line_width(record["flag_value"]) + zorder = 3 + record["flag_value"] + line, = ax.plot( + track["tphys"], + track["porb"], + lw=linewidth, + ls="-", + label=record["label"], + zorder=zorder, + ) + record["color"] = line.get_color() + mark_final_outcome( + ax, + track, + record["status"], + line.get_color(), + marker_size=endpoint_marker_size(record["flag_value"]), + zorder=5, + ) + + annotate_event(ax, event_label, event_time) + combined = pd.concat(all_plotted_tracks, ignore_index=True) + ax.set_yscale("log") + if not combined.empty: + max_time = combined["tphys"].max() + ax.set_xlim(0, max_time * 1.08) + min_porb = combined["porb"].min() + max_porb = combined["porb"].max() + ax.set_ylim(min_porb * 0.65, max_porb * 1.6) + ax.set_title(example["title"], fontsize=11) + ax.set_xlabel("Time [Myr]") + ax.set_ylabel("Orbital period [days]") + ax.grid(True, which="both", alpha=0.25) + flag_legend = ax.legend( + fontsize=7, + loc="best", + title=example["flag_key"], + title_fontsize=8, + ) + ax.text( + 0.02, + 0.03, + example["annotation"], + transform=ax.transAxes, + fontsize=7, + color="0.25", + ) +fig.legend( + handles=[ + Line2D([0], [0], marker="o", color="0.2", markerfacecolor="white", + lw=0, markersize=6, label="survives"), + Line2D([0], [0], marker="x", color="0.2", + lw=0, markersize=7, label="merges"), + Line2D([0], [0], marker="^", color="0.2", markerfacecolor="white", + lw=0, markersize=7, label="disrupts"), + Line2D([0], [0], color="0.35", ls=":", lw=1.4, label="CE/RLO"), + ], + loc="lower center", + ncol=5, + bbox_to_anchor=(0.5, -0.02), + fontsize=7, +) +plt.show()