Source code for deprojected_sersic_models.paper.paper_plots

##################################################################################
# deprojected_sersic_models/paper/paper_plots.py                                 #
#                                                                                #
# Copyright 2018-2022 Sedona Price <sedona.price@gmail.com> / MPE IR/Submm Group #
# Licensed under a 3-clause BSD style license - see LICENSE.rst                  #
##################################################################################

import os
import copy

import dill as pickle

import numpy as np
import scipy.interpolate as scp_interp

from astropy.io import ascii

import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib.ticker import MultipleLocator, FixedLocator, FixedFormatter
import matplotlib.cm as cm
from matplotlib.colors import LinearSegmentedColormap
import matplotlib.patches as mpatches

from deprojected_sersic_models import core, io, table_generation
from deprojected_sersic_models.utils import calcs as util_calcs
from deprojected_sersic_models.paper import plot_calcs
from deprojected_sersic_models.paper import scaling_rel
from deprojected_sersic_models.utils import interp_profiles

__all__ = [ 'make_all_paper_plots', 'list_table1_values' ]

# ---------------------
# Plot settings

mpl.rcParams['text.usetex'] = True

fontsize_leg_sm = 10.
fontsize_leg = 11.
fontsize_labels = 14.
fontsize_labels_sm = 12.
fontsize_ann = 11.
fontsize_title= 14.
fontsize_title_lg = 16.
fontsize_ticks = 12.
fontsize_ticks_sm = 11.

fontsize_ann_lg = 13.
fontsize_ann_latex_sm = 12.
fontsize_ann_latex = 13.

cmap_mass = cm.OrRd
cmap_q = cm.magma_r
cmap_n = cm.viridis_r
cmapg = cm.Greys


cmap_bt = LinearSegmentedColormap.from_list('cmap_bt',
            [(0, 0, 0), (0.9568627450980393, 0.4745098039215686, 0.4745098039215686)],
             N=101)
cmap_D2t_ref = cm.GnBu_r
ref_vals = np.arange(0.1, 0.7, 0.05)
arr_tuples = [(0,0,0,1)]
for rv in ref_vals:
    arr_tuples.append(cmap_D2t_ref(rv))
cmap_D2t = LinearSegmentedColormap.from_list('cmap_D2t', arr_tuples,N=101)


_dir_deprojected_sersic_models = os.getenv('DEPROJECTED_SERSIC_MODELS_DATADIR', None)



# ----------------------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------------

[docs] def list_table1_values(table_path=None, n_arr=[1.,4.], q_arr=[0.4,1.,1.5], rfac_arr=[1.,2.2/1.678], recalculate_all=False): """ Wrapper to list ktot, k3D values for Table 1. Prints values to output. Parameters ---------- table_path: str Path to the directory where the Sersic profile tables are located. Default: system variable `DEPROJECTED_SERSIC_MODELS_DATADIR`, if specified. n_arr: array_like, optional Array of Sersic indices to show q_arr: array_like, optional Array of intrinsic axis ratios to show rfac_arr: array_like, optional Array of radii to show, with R=R_fac_arr*Re Returns ------- """ if table_path is None: table_path = _dir_deprojected_sersic_models print("-----------------------------------------------------") print(" || ktot || k3D ") print("-----------------||----------------||----------------") strout = "n q ||" strrs = " " for j, rf in enumerate(rfac_arr): if j == 0: strrs += "R=" else: strrs += " " if rf == 1.: strrs += "Re" else: strrs += "{:0.1f}Re".format(rf) strrs += " " strout += strrs + "||" + strrs print(strout) k3d_ellip_arr = [] for n in n_arr: print("-----------------------------------------------------") for i, q in enumerate(q_arr): table = io.read_profile_table(n=n, invq=1./q, path=table_path) if recalculate_all: sprof = core.DeprojSersicModel(total_mass=table['total_mass'], Reff=table['Reff'], n=n, q=q) ktot_arr = [] k3d_arr = [] for rfac in rfac_arr: if recalculate_all: ktot_arr.append(sprof.virial_coeff_tot(rfac*table['Reff'])) k3d_arr.append(sprof.virial_coeff_3D(rfac*table['Reff'])) else: if rfac == 1.: # Use table values ktot_arr.append(table['ktot_Reff']) k3d_arr.append(table['k3D_sph_Reff']) else: # Interpolate & calculate: R = rfac*table['Reff'] table_rfac = interp_profiles.interpolate_entire_table(R=R, total_mass=table['total_mass'], Reff=table['Reff'], n=n, invq=1./q, path=table_path) m3D = table_rfac['menc3D_sph'] vc = table_rfac['vcirc'] ktot_arr.append(util_calcs.virial_coeff_tot(R, total_mass=table['total_mass'], vc=vc)) k3d_arr.append(util_calcs.virial_coeff_3D(R, m3D=m3D, vc=vc)) if (q == 0.4): m3D_ellip = table_rfac['menc3D_ellipsoid'] k3d_ellip_arr.append(util_calcs.virial_coeff_3D(R, m3D=m3D_ellip, vc=vc)) if i == 0: strout = 'n={:4.1f} '.format(n) else: strout = ' ' strout += 'q={:0.1f} || '.format(q) for kt in ktot_arr: strout += "{:0.3f} ".format(kt) strout += "|| " for k3d in k3d_arr: strout += "{:0.3f} ".format(k3d) print(strout) print("-----------------------------------------------------") strout = '' for k3d_ellip in k3d_ellip_arr: strout += "{:0.3f} ".format(k3d_ellip) print("Miller+11 comparison points: "+strout) return None
# ---------------------------------------------------------------------------------------------------- # ----------------------------------------------------------------------------------------------------
[docs] def make_all_paper_plots(output_path=None, table_path=None): """ Wrapper function to make all plots for the paper. Saves plots in PDF format. Parameters ---------- output_path: str Path to the directory where the plots will be saved. table_path: str Path to the directory where the Sersic profile tables are located. Default: system variable `DEPROJECTED_SERSIC_MODELS_DATADIR`, if specified. Returns ------- """ if table_path is None: table_path = _dir_deprojected_sersic_models # Figure 1 plot_compare_mencl(output_path=output_path, table_path=table_path, q_arr=[1., 0.4, 0.2]) # Figure 2 plot_rhalf_3D_sersic_potential(output_path=output_path, table_path=table_path, q_arr_plotq = np.array([0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 1., 1.5, 2.]), n_arr_plotn = np.array([0.5, 1., 2., 4., 8.])) # Figure 3 plot_composite_menc_vcirc_profile_invert(output_path=output_path, table_path=table_path, total_mass = 5.e10, Reff = 1., q_arr=[1., 0.4, 0.2], n_arr=[1., 4.]) # Figure 4 plot_virial_coeff(output_path=output_path, table_path=table_path, q_arr=[0.2, 0.4, 0.6, 0.8, 1., 1.5, 2.]) # Figure 5 # Use toy scaling relations / empirical approximations: # z = 2, lmstar=10.5, which gives: # Mstar = 3.2e10 (3.1623e10) # Mbar = 6.6e10 (6.64947e10) / lMbar = 10.82, # Mhalo = 8.9e11 (8.8618527e11) / lMhalo = 11.9475, # halo_conc = 4.2 (4.1931) # Reff_disk = 3.4 (3.400257) # qdisk = 0.25 or invq = 4 plot_example_galaxy_mencl_vcirc(bt_arr=[0., 0.25, 0.5, 0.75, 1.0], output_path=output_path, table_path=table_path, Mbaryon=6.6e10, z=2., Mhalo=8.9e11, halo_conc=4.2, Reff_disk=3.4, n_disk=1., invq_disk=4., Reff_bulge=1., n_bulge=4., invq_bulge=1., logradius=False, rmin=0., rmax=15., rstep=0.01) # Figure 6 plot_fdm_calibration(Mbaryon_arr=[6.6e10], Reff_disk_arr=[3.4], q_disk_arr = [0.01, 0.05, 0.1, 0.2, 0.25, 0.4, 0.6, 0.8, 1., 1.5, 2.], Mhalo_arr=[8.9e11], halo_conc_arr=[4.2], output_path=output_path, table_path=table_path, del_fDM=False) # Figure 7 plot_rhalf3DReD_fdm_calibration_ReB_ReD_rhalf3D(Mbaryon=6.6e10, q_disk_arr = [0.05, 0.1, 0.2, 0.25, 0.4, 0.6, 0.8, 1., 1.5, 2.], ReB_to_ReD_arr=[0.2, 0.5, 1.], n_disk=1., n_bulge=4., invq_bulge=1., Reff_bulge=1., z=2., Mhalo=8.9e11, halo_conc=4.2, output_path=output_path, table_path=table_path, del_fDM=True) # Figure 8 plot_toy_impl_fDM_calibration_z_evol(output_path=output_path, table_path=table_path, n_disk=1., Reff_bulge=1., n_bulge=4., invq_bulge=1., del_fDM=True) # Figure 9 plot_alpha_vs_R(output_path=output_path, table_path=table_path, n_arr=[0.5, 1., 2., 4., 8.], show_literature=True) # Figure 10 plot_AD_sersic_potential_alpha_vs_R(output_path=output_path, table_path=table_path, sigma0_arr = [30., 60., 90.], q_arr=[1., 0.2], n_arr=[0.5, 1., 2., 4.], total_mass = np.power(10.,10.5), Reff=1.) # Figure 11 plot_composite_alpha_bulge_2disk_vs_r(output_path=output_path, table_path=table_path, nDisk=1., nBulge=4., bt_arr=[0.,0.25,0.5,0.75,1.], ReB_to_ReD_arr=[0.2, 0.5, 1.], show_BT=True, show_D2T=False) return None
# ---------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------- # Figure 1 def plot_compare_mencl(fileout=None, output_path=None, table_path=None, q_arr=[1., 0.4, 0.2]): """ Plot the fractional 3D spherical enclosed mass profile versus radius, for a range of Sersic indices n (n=1..8) and intrinsic axis ratios q. Saves plot to PDF. Parameters ---------- output_path: str Path to directory where the output plot will be saved. table_path: str Path to directory containing the Sersic profile tables. q_arr: array_like, optional Range of intrinsic axis ratios to plot. Default: q_arr = [1., 0.4, 0.2] fileout: str, optional Override the default filename and explicitly choose the output filename (must include full path). """ if (output_path is None) & (fileout is None): raise ValueError("Must set 'output_path' if 'fileout' is not set !") if table_path is None: raise ValueError("Must set 'table_path' !") if (fileout is None): # Ensure trailing slash: if output_path[-1] != '/': output_path += '/' fileout = output_path+'mencl_sersic_potential' for q in q_arr: fileout += '_q{:0.2f}'.format(q) fileout += '.pdf' nmin = 1 nmax = 8 nstep = 1 nrange = nmax-nmin n_arr = np.linspace(nmin,nmax, num=np.int(np.round(((nrange)/nstep)))+1) n_arr = np.append(np.array([0.5]), n_arr) color_arr, ls_arr, labels = ([] for _ in range(3)) nextra = -0.5 for n in n_arr: if n == 8.: color_arr.append('#2d0037') else: color_arr.append(cmap_n((n+nextra)/(nrange+nextra))) ls_arr.append('-') if n == 0.5: labels.append(r'$n={:0.1f}$'.format(n)) else: labels.append(r'$n={:0.0f}$'.format(n)) # Load files: ks_dicts = [] for q in q_arr: invq = 1./q ks_dict = {} ks_dict['3D'], ks_dict['3D_ellip'], ks_dict['vcirc'] = ({} for mm in range(3)) for n in n_arr: tab = io.read_profile_table(n=n, invq=invq, path=table_path) if n == 1.: ks_dict['R_arr'] = tab['R'] ks_dict['Reff'] = tab['Reff'] ks_dict['3D']['n={}'.format(n)] = tab['menc3D_sph'] / tab['total_mass'] ks_dict['3D_ellip']['n={}'.format(n)] = tab['menc3D_ellipsoid'] / tab['total_mass'] ks_dict['vcirc']['n={}'.format(n)] = tab['vcirc'] ks_dicts.append(ks_dict) # ++++++++++++++++ # plot 3D: xlabel = r'$\log_{10}(R/R_{\mathrm{e}})$' types, ylabels, titles, vline_loc, vline_lss, vline_labels = ([] for _ in range(6)) for i, q in enumerate(q_arr): types.append('3D') ylabels.append(r'$M_{\mathrm{sph}}(<r=R)/M_{\mathrm{tot}}$') if q < 1.: titles.append(r'$q_0={:0.1f}$'.format(q)) else: titles.append(r'$q_0={:0.0f}$'.format(q)) vline_loc.append([0., np.log10(2.2/1.676)]) vline_lss.append(['--', '-.']) if i == 0: vline_labels.append([r'$R=R_{\mathrm{e}}$', r'$R=1.3 R_{\mathrm{e}}$']) else: vline_labels.append([None, None]) xlim = [-2.0, 2.0] ylims = [[0., 1.], [0., 1.], [0., 1.]] spec_pairs = [[1., 0.], [4., 0.]] lw = 1.25 ###################################### # Setup plot: f = plt.figure() scale = 3.5 n_cols = len(types) fac = 1.155 f.set_size_inches(fac*scale*n_cols,scale) pad_outer = 0.25 gs = gridspec.GridSpec(1, n_cols, wspace=pad_outer) axes = [] for i in range(n_cols): axes.append(plt.subplot(gs[0,i])) for i in range(n_cols): ax = axes[i] ylim = ylims[i] ks_dict = ks_dicts[i] q = q_arr[i] if ks_dict is not None: n_cnt = 0 for j, n in enumerate(n_arr): menc_arr = ks_dict[types[i]]['n={}'.format(n)] if menc_arr is not None: n_cnt += 1 Rarr_plot = np.log10(ks_dict['R_arr']/ks_dict['Reff']) if len(menc_arr) != len(ks_dict['R_arr']): raise ValueError ax.plot(Rarr_plot, menc_arr, ls=ls_arr[j], color=color_arr[j], lw=lw, label=labels[j]) for mm, sp_p in enumerate(spec_pairs): if sp_p[0] == n: wh = np.where(Rarr_plot == sp_p[1])[0] if len(wh) > 0: ax.axhline(y=menc_arr[wh[0]], ls='--', lw=0.8, color=color_arr[j], zorder=-20.) fracx = 0.015 if (q == 1.): delt = 0.015 elif (q == 0.4): delt = -0.07 - 0.04 * (len(spec_pairs)-mm-1) elif (q == 0.2): delt = -0.057 - 0.05 * (len(spec_pairs)-mm-1) xypos = (xlim[1]-fracx*(xlim[1]-xlim[0]), menc_arr[wh[0]]+delt) ax.annotate(r'{:0.1f}'.format(menc_arr[wh[0]]*100)+'\%', xy=xypos, ha='right', color=color_arr[j], fontsize=fontsize_ann) ax.axhline(y=0.5, ls='-.', color='grey', zorder=-20.) delt = 0.015 fracx = 0.015 ax.annotate(r'50\%', xy=(xlim[1]-fracx*(xlim[1]-xlim[0]), 0.5+delt), ha='right', color='dimgrey', fontsize=fontsize_ann) if vline_loc[i] is not None: for vl, vls, vlb in zip(vline_loc[i], vline_lss[i], vline_labels[i]): ax.axvline(x=vl, ls=vls, color='grey', label=vlb, zorder=-20.) ax.set_xlim(xlim) ax.set_ylim(ylim) ax.set_xlabel(xlabel, fontsize=fontsize_labels) ax.set_ylabel(ylabels[i], fontsize=fontsize_labels) ax.tick_params(labelsize=fontsize_ticks) if titles[i] is not None: ax.set_title(titles[i], fontsize=fontsize_title) ax.xaxis.set_minor_locator(MultipleLocator(0.25)) ax.xaxis.set_major_locator(MultipleLocator(1.)) ax.yaxis.set_minor_locator(MultipleLocator(0.05)) ax.yaxis.set_major_locator(MultipleLocator(0.2)) if i == 0: handles, labels_leg = ax.get_legend_handles_labels() neworder = range(n_cnt) handles_arr, labels_arr = ([] for _ in range(2)) for i in neworder: handles_arr.append(handles[i]) labels_arr.append(labels_leg[i]) neworder2 = range(n_cnt, len(handles)) handles_arr2, labels_arr2 = ([] for _ in range(2)) for i in neworder2: handles_arr2.append(handles[i]) labels_arr2.append(labels_leg[i]) frameon = True framealpha = 1. edgecolor = 'none' borderpad = 0.25 fontsize_leg_tmp = fontsize_leg labelspacing=0.15 handletextpad=0.25 legend1 = ax.legend(handles_arr, labels_arr, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, loc='upper left', numpoints=1, scatterpoints=1, frameon=frameon, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp) legend2 = ax.legend(handles_arr2, labels_arr2, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, loc='lower right', numpoints=1, scatterpoints=1, frameon=frameon, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp) ax.add_artist(legend1) ax.add_artist(legend2) plt.savefig(fileout, bbox_inches='tight', dpi=600) plt.close() return None # ---------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------- # Figure 2 def plot_rhalf_3D_sersic_potential(output_path=None, table_path=None, fileout=None, q_arr_plotq = np.array([0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 1., 1.5, 2.]), n_arr_plotn = np.array([0.5, 1., 2., 4., 8.])): """ Plot the ratio of the 3D half mass radius to the 2D projected major axis effective radius (Reff), for a range of Sersic indices n and range of intrinsic axis ratios q. Saves plot to PDF. Parameters ---------- output_path: str Path to directory where the output plot will be saved. table_path: str Path to directory containing the Sersic profile tables. fileout: str, optional Override the default filename and explicitly choose the output filename (must include full path). """ if (output_path is None) & (fileout is None): raise ValueError("Must set 'output_path' if 'fileout' is not set !") if table_path is None: raise ValueError("Must set 'table_path' !") if (fileout is None): # Ensure trailing slash: if output_path[-1] != '/': output_path += '/' fileout = output_path+'rhalf_3D_sersic_potential' fileout += '.pdf' # -------------------------------------------------------------- # Arrays: # PLOTTING COLORED BY q: nmin = 0.5 nmax = 8. nstep = 0.1 nrange = 7. n_arr_plotq = np.linspace(nmin,nmax, num=np.int(np.round((nmax-nmin)/nstep))+1) color_arr_plotq, ls_arr_plotq, labels_plotq = ([] for _ in range(3)) for q in q_arr_plotq: if q <= 1.: color_arr_plotq.append(cmap_q(q)) ls_arr_plotq.append('-') else: color_arr_plotq.append(cmapg(1./q)) ls_arr_plotq.append('--') if q != 1.: labels_plotq.append(r'$q_0={}$'.format(q)) else: labels_plotq.append(r'$q_0={:0.0f}$'.format(q)) # -------------------------------------------------------------- # PLOTTING COLORED BY n: nrange = 8.-1. q_arr_plotn = np.array([0.01, 0.05, 0.1, 0.125, 0.1428, 0.1667, 0.2, 0.25, 0.3, 0.333, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.5, 2.]) color_arr_plotn, ls_arr_plotn, labels_plotn = ([] for _ in range(3)) nextra = -0.5 for n in n_arr_plotn: if n == 8.: color_arr_plotn.append('#2d0037') else: color_arr_plotn.append(cmap_n((n+nextra)/(nrange+nextra))) ls_arr_plotn.append('-') if n < 1: labels_plotn.append(r'$n={}$'.format(n)) else: labels_plotn.append(r'$n={:0.0f}$'.format(n)) # -------------------------------------------------------------- # COMPOSITE STACK: n_arrs = [n_arr_plotq, n_arr_plotn] q_arrs = [q_arr_plotq, q_arr_plotn] color_arrs = [color_arr_plotq, color_arr_plotn] ls_arrs = [ls_arr_plotq, ls_arr_plotn] labels_arrs = [labels_plotq, labels_plotn] plot_types = ['plotq', 'plotn'] # -------------------------------------------------------------- # Load files: ks_dicts_list = [] for i, plot_type in enumerate(plot_types): q_arr = q_arrs[i] n_arr = n_arrs[i] ks_dicts = [] arrs = {'q': q_arr, 'n': n_arr} if plot_type == 'plotq': var_order = ['q', 'n'] elif plot_type == 'plotn': var_order = ['n', 'q'] for v1 in arrs[var_order[0]]: if var_order[0] == 'q': q = v1 elif var_order[0] == 'n': n = v1 ks_dict = {} rhalf3D_arr = np.ones(len(arrs[var_order[1]])) * -99. for j, v2 in enumerate(arrs[var_order[1]]): if var_order[1] == 'q': q = v2 elif var_order[1] == 'n': n = v2 invq = 1./q tab = io.read_profile_table(n=n, invq=invq, path=table_path) rhalf3D_arr[j] = util_calcs.find_rhalf3D_sphere(R=tab['R'], menc3D_sph=tab['menc3D_sph'], total_mass=tab['total_mass']) if plot_type == 'plotq': ks_dict['q'] = q ks_dict['n_arr'] = n_arr elif plot_type == 'plotn': ks_dict['n'] = n ks_dict['q_arr'] = q_arr ks_dict['rhalf3D'] = rhalf3D_arr ks_dicts.append(ks_dict) ### ks_dicts_list.append(ks_dicts) # ++++++++++++++++ # plot: ylim = [0.96, 1.86] xlim_plotq = [0., 8.25] ylim_plotq = ylim xlabel_plotq = r'$n$' xlim_plotn = [-0.05, 2.05] ylim_plotn = ylim xlabel_plotn = r'$q_0$' xlims = [xlim_plotq, xlim_plotn] ylims = [ylim_plotq, ylim_plotn] xlabels = [xlabel_plotq, xlabel_plotn] ylabels = [r'$r_{1/2,\, \mathrm{mass,\,3D}}/R_{\mathrm{e}}$', r'$r_{1/2,\, \mathrm{mass,\,3D}}/R_{\mathrm{e}}$'] titles = [None, None] marker = 'o' marker = None ms = 2 lw = 1.3 ###################################### # Setup plot: f = plt.figure() scale = 3.75 n_cols = len(xlims) n_rows = 1 fac = 1.33 f.set_size_inches(fac*scale*n_cols,scale*n_rows) wspace = 0.6 hspace = 0. gs = gridspec.GridSpec(n_rows, n_cols, wspace=wspace, hspace=hspace) axes = [] for i in range(n_rows): for j in range(n_cols): axes.append(plt.subplot(gs[i,j])) ##### for i, plot_type in enumerate(plot_types): ax = axes[i] if (i == 0): ax.set_zorder(1) ks_dicts = ks_dicts_list[i] q_arr = q_arrs[i] n_arr = n_arrs[i] ls_arr = ls_arrs[i] color_arr = color_arrs[i] labels = labels_arrs[i] xlim = xlims[i] ylim = ylims[i] xlabel = xlabels[i] ylabel = ylabels[i] title = titles[i] if plot_type == 'plotq': for i, q in enumerate(q_arr): ks_dict = ks_dicts[i] wh_fin = np.where(np.isfinite(ks_dict['rhalf3D']))[0] if len(wh_fin) > 0: ax.plot(ks_dict['n_arr'][wh_fin], ks_dict['rhalf3D'][wh_fin], ls=ls_arr[i], color=color_arr[i], lw=lw, label=labels[i], marker=marker, ms=ms) ax.axvline(x=1., ls=':', color='lightgrey', zorder=-20.) ax.axvline(x=4., ls=':', color='lightgrey', zorder=-20.) elif plot_type == 'plotn': for i, n in enumerate(n_arr): ks_dict = ks_dicts[i] wh_fin = np.where(np.isfinite(ks_dict['rhalf3D']))[0] if len(wh_fin) > 0: ax.plot(ks_dict['q_arr'][wh_fin], ks_dict['rhalf3D'][wh_fin], ls=ls_arr[i], color=color_arr[i], lw=lw, label=labels[i], marker=marker, ms=ms) ax.axhline(y=1., ls=":", color='lightgrey', zorder=-20.) ##### if ylim is None: ylim = ax.get_ylim() ax.set_xlim(xlim) ax.set_ylim(ylim) if xlabel is not None: ax.set_xlabel(xlabel, fontsize=fontsize_labels) else: ax.set_xticklabels([]) if ylabel is not None: ax.set_ylabel(ylabel, fontsize=fontsize_labels) else: ax.set_yticklabels([]) ax.tick_params(labelsize=fontsize_ticks) if title is not None: ax.set_title(title, fontsize=fontsize_title) if plot_type == 'plotq': ax.xaxis.set_minor_locator(MultipleLocator(0.2)) ax.xaxis.set_major_locator(MultipleLocator(1.)) elif plot_type == 'plotn': ax.xaxis.set_minor_locator(MultipleLocator(0.05)) ax.xaxis.set_major_locator(MultipleLocator(0.5)) if (ylim[1]-ylim[0]) > 0.6: ax.yaxis.set_minor_locator(MultipleLocator(0.02)) ax.yaxis.set_major_locator(MultipleLocator(0.1)) else: ax.yaxis.set_minor_locator(MultipleLocator(0.01)) ax.yaxis.set_major_locator(MultipleLocator(0.1)) frameon = False #True framealpha = 1. edgecolor = 'none' borderpad = 0.5 fontsize_leg_tmp = fontsize_leg labelspacing=0.15 handletextpad=0.25 if plot_type == 'plotq': loc = 'upper left' bbox_to_anchor = (0.99, 1.02) elif plot_type == 'plotn': loc = 'upper left' bbox_to_anchor = (0., 1.02) legend = ax.legend(labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, loc=loc, bbox_to_anchor=bbox_to_anchor, frameon=frameon, framealpha=framealpha, edgecolor=edgecolor, numpoints=1, scatterpoints=1, fontsize=fontsize_leg_tmp) if fileout is not None: plt.savefig(fileout, bbox_inches='tight', dpi=600) plt.close() else: plt.show() # ---------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------- # Figure 3 def plot_composite_menc_vcirc_profile_invert(q_arr=[1., 0.4, 0.2], n_arr=[1., 4.], total_mass = 5.e10, Reff = 1., output_path=None, table_path=None, fileout=None): """ Plot the enclosed mass and circular velocity profiles versus radius, for a number of direct calculations and inverted assumptions (applying the relation valid for spherically symmetric cases: v^2 = GM/r ). Saves plot to PDF. Parameters ---------- output_path: str Path to directory where the output plot will be saved. table_path: str Path to directory containing the Sersic profile tables. q_arr: array_like, optional Range of intrinsic axis ratios to plot. Default: q_arr = [1., 0.4, 0.2] n_arr: array_like, optional Range of Sersic indices to plot. Default: n_arr = [1., 4.] total_mass: float, optional Total mass of the Sersic component [Msun]. Default: 5e10 Msun Reff: float, optional Effective radius of the Sersic component [kpc]. Default: 1 kpc fileout: str, optional Override the default filename and explicitly choose the output filename (must include full path). """ if (output_path is None) & (fileout is None): raise ValueError("Must set 'output_path' if 'fileout' is not set !") if table_path is None: raise ValueError("Must set 'table_path' !") if (fileout is None): # Ensure trailing slash: if output_path[-1] != '/': output_path += '/' fileout = output_path+'mencl_vcirc_compsite_compare_sersic_potential' for q in q_arr: fileout += '_q{:0.2f}'.format(q) for n in n_arr: fileout += '_n{:1.0f}'.format(n) fileout += '.pdf' # Load files: ks_dicts = [] for n in n_arr: for q in q_arr: # read fits tables, construct ks dict..... invq = 1./q tab_orig = io.read_profile_table(n=n, invq=invq, path=table_path) tab = interp_profiles.interpolate_entire_table(R=tab_orig['R']*Reff, table=tab_orig, total_mass=total_mass, Reff=Reff, n=n, invq=invq) ks_dict = {} ks_dict['R_arr'] = tab['R'] ks_dict['Reff'] = tab['Reff'] ks_dict['total_mass'] = tab['total_mass'] ks_dict['menc2D_ellip'] = util_calcs.total_mass2D_direct(tab['R'], total_mass=tab['total_mass'], q=q, n=n, Reff=tab['Reff'],i=90.) ks_dict['menc3D_sph'] = tab['menc3D_sph'] ks_dict['menc3D_ellip'] = tab['menc3D_ellipsoid'] ks_dict['vcirc'] = tab['vcirc'] ks_dict['v_menc3D_sph'] = util_calcs.vcirc_spherical_symmetry(R=tab['R'], menc=tab['menc3D_sph']) ks_dict['v_menc3D_ellip'] = util_calcs.vcirc_spherical_symmetry(R=tab['R'], menc=tab['menc3D_ellipsoid']) ks_dict['m_invert_vcirc'] = util_calcs.menc_spherical_symmetry(R=tab['R'], vcirc=tab['vcirc']) ks_dict['mass_plot'] = tab['total_mass'] ks_dicts.append(ks_dict) # ++++++++++++++++ # plot: types1 = ['menc2D_ellip', 'menc3D_sph', 'menc3D_ellip', 'm_invert_vcirc'] ls_arr1 = ['-','--','-.',':'] lw_arr1 = [1.1, 1.35, 1.1, 1.35] labels1 = [r'2D, ellipse, Eq. 3', r'3D, sphere, Eq. 8', r'3D, spheroid, Eq. 6', r'$v_{\mathrm{circ}}(R)^2 R/G$, Eq. 5'] color_arr1 = [cmap_q(0.18), cmap_q(1./3.), cmap_q(5./9.), cmap_q(1.)] types2 = ['v_menc3D_sph', 'v_menc3D_ellip', 'vcirc'] ls_arr2 = ['--','-.',':'] lw_arr2 = [1.35, 1.1, 1.35] labels2 = [r'$\sqrt{M_{\mathrm{3D,sph}}(<r=R)G/R}$, Eq. 8', r'$\sqrt{M_{\mathrm{3D,spheroid}}(<m=R)G/R}$, Eq. 6', r'$v_{\mathrm{circ}}(R)$, Eq. 5'] color_arr2 = [cmap_q(1./3.), cmap_q(5./9.), cmap_q(1.)] titles1, vline_loc1, vline_lss1, vline_labels1 = ([] for _ in range(4)) xlims1, ylims1, xlabels1, ylabels1 = ([] for _ in range(4)) titles2, vline_loc2, vline_lss2, vline_labels2 = ([] for _ in range(4)) xlims2, ylims2, xlabels2, ylabels2 = ([] for _ in range(4)) for i, n in enumerate(n_arr): for j, q in enumerate(q_arr): if i == 0: if j == np.int(np.round((len(q_arr)-1)/2.)): titles1.append(r'Fractional mass profile') titles2.append(r'Circular velocity') else: titles1.append(None) titles2.append(None) else: titles1.append(None) titles2.append(None) if i == len(n_arr) - 1: xlabels1.append(r'$\log_{10}(R/R_{\mathrm{e}})$') xlabels2.append(r'$\log_{10}(R/R_{\mathrm{e}})$') else: xlabels1.append(None) xlabels2.append(None) if j == 0: ylabels1.append(r'$M_{\mathrm{enc}}(<R)/M_{\mathrm{total}}$') ylabels2.append(r'$v(R)\ \mathrm{[km\,s^{-1}]}$') else: ylabels1.append(None) ylabels2.append(None) xlims1.append([-2.0, 2.0]) xlims2.append([-2.0, 2.0]) ylims1.append([-0.05, 1.22]) ylims2.append([0., 410.]) vline_loc1.append([0., np.log10(2.2/1.676)]) vline_lss1.append(['--', '-.']) vline_loc2.append([0., np.log10(2.2/1.676)]) vline_lss2.append(['--', '-.']) if ((i == 0) & (j == 0)): vline_labels1.append([r'$R=R_{\mathrm{e}}$', r'$R=1.3 R_{\mathrm{e}}$']) vline_labels2.append([r'$R=R_{\mathrm{e}}$', r'$R=1.3 R_{\mathrm{e}}$']) else: vline_labels1.append([None, None]) vline_labels2.append([None, None]) ###################################### # Setup plot: f = plt.figure() scale = 3.35 n_cols = len(q_arr) n_rows = len(n_arr) fac = 0.98 f.set_size_inches(fac*scale*n_cols,2.*scale*n_rows) pad_outer = 0.175 wspace = 0.025 hspace = wspace gs_outer = gridspec.GridSpec(2, 1, wspace=wspace, hspace=pad_outer) gs1 = gridspec.GridSpecFromSubplotSpec(n_rows,n_cols,subplot_spec=gs_outer[0,0],wspace=wspace, hspace=hspace) gs2 = gridspec.GridSpecFromSubplotSpec(n_rows,n_cols, subplot_spec=gs_outer[1,0], wspace=wspace, hspace=hspace) axes1, axes2 = ([] for _ in range(2)) for i in range(n_rows): for j in range(n_cols): axes1.append(plt.subplot(gs1[i,j])) axes2.append(plt.subplot(gs2[i,j])) for i in range(n_rows): for j in range(n_cols): k = i*n_cols + j n = n_arr[i] q = q_arr[j] ks_dict = ks_dicts[k] ax1 = axes1[k] ax2 = axes2[k] xlim1 = xlims1[k] ylim1= ylims1[k] ylabel1 = ylabels1[k] title1 = titles1[k] xlim2 = xlims2[k] ylim2 = ylims2[k] ylabel2 = ylabels2[k] title2 = titles2[k] if ks_dict is not None: for mm, typ in enumerate(types1): menc_arr = ks_dict[types1[mm]] Rarr_plot = np.log10(ks_dict['R_arr']/ks_dict['Reff']) ax1.plot(Rarr_plot, menc_arr/ks_dict['mass_plot'], ls=ls_arr1[mm], color=color_arr1[mm], lw=lw_arr1[mm], label=labels1[mm], zorder=10.) for mm, typ in enumerate(types2): menc_arr = ks_dict[types2[mm]] Rarr_plot = np.log10(ks_dict['R_arr']/ks_dict['Reff']) ax2.plot(Rarr_plot, menc_arr, ls=ls_arr2[mm], color=color_arr2[mm], lw=lw_arr2[mm], label=labels2[mm]) if ylim1 is None: ylim1 = ax1.get_ylim() if ylim2 is None: ylim2 = ax2.get_ylim() ax1.axhline(y=0.5*ks_dict['total_mass']/ks_dict['mass_plot'], ls='-.', color='darkgrey', zorder=-20.) ax1.axhline(y=ks_dict['total_mass']/ks_dict['mass_plot'], ls=':', color='darkgrey', zorder=-20.) if k == 0: delt = 0.015 fracx = 0.015 yspan = ylim1[1]-ylim1[0] ax1.annotate(r'50\%', xy=(xlim1[1]-fracx*(xlim1[1]-xlim1[0]), (0.5+delt)), va='bottom', ha='right', color='grey', fontsize=fontsize_ann) ax1.annotate(r'100\%', xy=(xlim1[1]-fracx*(xlim1[1]-xlim1[0]), (1.0+delt)), va='bottom', ha='right', color='grey', fontsize=fontsize_ann) if vline_loc1[k] is not None: for vl, vls, vlb in zip(vline_loc1[k], vline_lss1[k], vline_labels1[k]): ax1.axvline(x=vl, ls=vls, color='grey', label=vlb, zorder=-20.) if vline_loc2[k] is not None: for vl, vls, vlb in zip(vline_loc2[k], vline_lss2[k], vline_labels2[k]): ax2.axvline(x=vl, ls=vls, color='grey', label=vlb, zorder=-20.) ###################### # ANNOTATE: ####### posannx = 0.71 posanny = 0.95 dely = 0. bbox_props = dict(boxstyle="square", fc="ghostwhite", ec="dimgrey", lw=1) ann_str = r'$n\ \:'+r'={:0.0f}$'.format(n)+'\n'+r'$q_0={:0.1f}$'.format(q) # Annotate ax1: ax1.annotate(ann_str, xy=(posannx, 1-posanny + dely), xycoords='axes fraction', va='bottom', ha='left', fontsize=fontsize_ann_lg, bbox=bbox_props) # Annotate ax2: if k == 0: posannx = 0.71 posanny = 0.75 ax2.annotate(ann_str, xy=(posannx, posanny), xycoords='axes fraction', va='top', ha='left', fontsize=fontsize_ann_lg, bbox=bbox_props, zorder=10.) ###################### ax1.set_xlim(xlim1) ax1.set_ylim(ylim1) ax2.set_xlim(xlim2) ax2.set_ylim(ylim2) ax1.tick_params(labelsize=fontsize_ticks) ax2.tick_params(labelsize=fontsize_ticks) ax1.xaxis.set_minor_locator(MultipleLocator(0.25)) ax2.xaxis.set_minor_locator(MultipleLocator(0.25)) if (j > 0) & (i == n_rows-1): ax1.xaxis.set_major_locator(FixedLocator([-2., -1., 0., 1., 2.])) ax1.xaxis.set_major_formatter(FixedFormatter(["", "-1", "0", "1", "2"])) ax2.xaxis.set_major_locator(FixedLocator([-2., -1., 0., 1., 2.])) ax2.xaxis.set_major_formatter(FixedFormatter(["", "-1", "0", "1", "2"])) else: ax1.xaxis.set_major_locator(MultipleLocator(1.)) ax2.xaxis.set_major_locator(MultipleLocator(1.)) ax1.yaxis.set_minor_locator(MultipleLocator(0.05)) ax1.yaxis.set_major_locator(MultipleLocator(0.2)) ax2.yaxis.set_minor_locator(MultipleLocator(20.)) ax2.yaxis.set_major_locator(MultipleLocator(100.)) if xlabels1[k] is not None: ax1.set_xlabel(xlabels1[k], fontsize=fontsize_labels) else: ax1.set_xticklabels([]) if ylabels1[k] is not None: ax1.set_ylabel(ylabels1[k], fontsize=fontsize_labels) else: ax1.set_yticklabels([]) if xlabels2[k] is not None: ax2.set_xlabel(xlabels2[k], fontsize=fontsize_labels) else: ax2.set_xticklabels([]) if ylabels2[k] is not None: ax2.set_ylabel(ylabels2[k], fontsize=fontsize_labels) else: ax2.set_yticklabels([]) if titles1[k] is not None: ax1.set_title(titles1[k], fontsize=fontsize_title_lg) if titles2[k] is not None: ax2.set_title(titles2[k], fontsize=fontsize_title_lg) if k == 0: handles1, labels1 = ax1.get_legend_handles_labels() neworder1 = range(len(types1)) handles_arr1 = [] labels_arr1 = [] for ii in neworder1: handles_arr1.append(handles1[ii]) labels_arr1.append(labels1[ii]) neworder12 = range(len(types1), len(handles1)) handles_arr12 = [] labels_arr12 = [] for ii in neworder12: handles_arr12.append(handles1[ii]) labels_arr12.append(labels1[ii]) frameon = True framealpha = 1. edgecolor = 'none' borderpad = 0.25 fontsize_leg_tmp = fontsize_leg_sm labelspacing=0.15 handletextpad=0.25 loc1 = (0.025, 0.25) legend11 = ax1.legend(handles_arr1, labels_arr1, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, handlelength=1.15, loc='upper left', frameon=frameon, numpoints=1, scatterpoints=1, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp) legend12 = ax1.legend(handles_arr12, labels_arr12, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, loc=loc1, frameon=False, numpoints=1, scatterpoints=1, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp) ax1.add_artist(legend11) ax1.add_artist(legend12) ##### handles2, labels2 = ax2.get_legend_handles_labels() neworder2 = range(len(types2)) handles_arr2 = [] labels_arr2 = [] for ii in neworder2: handles_arr2.append(handles2[ii]) labels_arr2.append(labels2[ii]) neworder2 = range(len(types2), len(handles2)) handles_arr22 = [] labels_arr22 = [] for ii in neworder2: handles_arr22.append(handles2[ii]) labels_arr22.append(labels2[ii]) legend21 = ax2.legend(handles_arr2, labels_arr2, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, loc='upper left', frameon=frameon, numpoints=1, scatterpoints=1, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp) legend22 = ax2.legend(handles_arr22, labels_arr22, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, loc='lower right', frameon=False, numpoints=1, scatterpoints=1, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp) ax2.add_artist(legend21) ax2.add_artist(legend22) if fileout is not None: plt.savefig(fileout, bbox_inches='tight', dpi=600) plt.close() else: plt.show() # ---------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------- # Figure 4 def plot_virial_coeff(fileout=None, output_path=None, table_path=None, q_arr=[0.2, 0.4, 0.6, 0.8, 1., 1.5, 2.]): """ Plot the total and 3D virial coefficients at Reff, ktot(Reff) and k3D(Reff), as a function of Sersic index n for a range of intrinsic axis ratios q. Saves plot to PDF. Parameters ---------- output_path: str Path to directory where the output plot will be saved. table_path: str Path to directory containing the Sersic profile tables. q_arr: array_like, optional Range of intrinsic axis ratios to plot. Default: q_arr = [0.2, 0.4, 0.6, 0.8, 1., 1.5, 2.] (mostly oblate, 2 prolate) fileout: str, optional Override the default filename and explicitly choose the output filename (must include full path). """ if (output_path is None) & (fileout is None): raise ValueError("Must set 'output_path' if 'fileout' is not set !") if table_path is None: raise ValueError("Must set 'table_path' !") if (fileout is None): # Ensure trailing slash: if output_path[-1] != '/': output_path += '/' fileout = output_path+'virial_coeff_sersic_potential.pdf' #################### q_arr = np.array(q_arr) color_arr, ls_arr, labels = ([] for _ in range(3)) for q in q_arr: if q <= 1.: color_arr.append(cmap_q(q)) ls_arr.append('-') else: color_arr.append(cmapg(1./q)) ls_arr.append('--') labels.append(r'$q_0={}$'.format(q)) # Load files: # read fits tables, construct ks dict..... invq = 1./q ks_dict = {} ks_dict['tot'], ks_dict['3D'] = ({} for mm in range(2)) nmin = 0.5 nmax = 8. nstep = 0.1 n_arr = np.linspace(nmin,nmax, num=np.int(np.round((nmax-nmin)/nstep))+1) for q in q_arr: invq = 1./q ktot_narr = np.ones(len(n_arr)) * -99. k3D_narr = np.ones(len(n_arr)) * -99. for j, n in enumerate(n_arr): tab = io.read_profile_table(n=n, invq=invq, path=table_path) if (q == 0.2) & (n == 1.): ks_dict['Reff'] = tab['Reff'] ks_dict['tot']['narr'] = n_arr ks_dict['3D']['narr'] = n_arr ktot_narr[j] = util_calcs.virial_coeff_tot(tab['Reff'], total_mass=tab['total_mass'], vc=tab['vcirc_Reff']) k3D_narr[j] = util_calcs.virial_coeff_3D(tab['Reff'], vc=tab['vcirc_Reff'], m3D=tab['menc3D_sph_Reff']) ks_dict['tot']['q={}'.format(q)] = ktot_narr ks_dict['3D']['q={}'.format(q)] = k3D_narr types = ['tot', '3D'] xlabel = r'$n$' ylabels = [r'$k_{\mathrm{tot}}(R=R_{\mathrm{e}})$', r'$k_{\mathrm{3D}}(R=R_{\mathrm{e}})$'] xlim = [0., 8.25] ylims = [[1.75, 5.5], [0.85, 1.15]] titles = [None, None] ann_arr = [ r'$\displaystyle k_{\mathrm{tot}}(R_{\mathrm{e}}) = \frac{M_{\mathrm{tot}} G }{v_{\mathrm{circ}}(R_{\mathrm{e}})^2 R_{\mathrm{e}}}$', r'$\displaystyle k_{\mathrm{3D}}(R_{\mathrm{e}}) = \frac{M_{\mathrm{sph}}(<r=R_{\mathrm{e}}) G }{v_{\mathrm{circ}}(R_{\mathrm{e}})^2 R_{\mathrm{e}}}$' ] ann_arr_pos = ['upperright', 'lowerright'] lw = 1.3 ###################################### # Setup plot: f = plt.figure() scale = 3.85 n_cols = len(types) fac = 1.175 f.set_size_inches(fac*scale*n_cols,scale) pad_outer = 0.25 gs = gridspec.GridSpec(1, n_cols, wspace=pad_outer) axes = [] for i in range(n_cols): axes.append(plt.subplot(gs[0,i])) for i in range(n_cols): ax = axes[i] ylim = ylims[i] for j, q in enumerate(q_arr): k_arr = ks_dict[types[i]]['q={}'.format(q)] n_arr = ks_dict[types[i]]['narr'] if len(k_arr) != len(n_arr): raise ValueError if q > 1: zorder = 5. else: zorder = 10. ax.plot(n_arr, k_arr, ls=ls_arr[j], color=color_arr[j], lw=lw, label=labels[j], zorder=zorder) ax.axvline(x=1., ls=':', color='lightgrey', zorder=-20.) ax.axvline(x=4., ls=':', color='lightgrey', zorder=-20.) ax.set_xlim(xlim) ax.set_ylim(ylim) ax.set_xlabel(xlabel, fontsize=fontsize_labels) ax.set_ylabel(ylabels[i], fontsize=fontsize_labels) ax.tick_params(labelsize=fontsize_ticks) if titles[i] is not None: ax.set_title(titles[i], fontsize=fontsize_title) ax.xaxis.set_minor_locator(MultipleLocator(0.25)) ax.xaxis.set_major_locator(MultipleLocator(1.)) xydelt = 0.04 if ann_arr_pos[i] == 'lowerright': xy = (1.-xydelt, xydelt) va='bottom' ha='right' elif ann_arr_pos[i] == 'upperright': xy = (1.-xydelt, 1.-xydelt) va='top' ha='right' if (i==1): patch = mpatches.FancyBboxPatch((3.5,0.869), 4., 0.025, boxstyle="square,pad=0.", fc="white", ec="dimgrey", lw=0.) ax.add_patch(patch) ax.annotate(ann_arr[i], xy=xy, va=va, ha=ha, fontsize=fontsize_ann, xycoords='axes fraction') if (ylim[1]-ylim[0]) > 3: ax.yaxis.set_minor_locator(MultipleLocator(0.2)) ax.yaxis.set_major_locator(MultipleLocator(1.)) else: ax.yaxis.set_minor_locator(MultipleLocator(0.02)) ax.yaxis.set_major_locator(MultipleLocator(0.1)) if i == 1: frameon = True #False framealpha = 1. edgecolor = 'none' borderpad = 0.25 fontsize_leg_tmp = fontsize_leg labelspacing=0.15 handletextpad=0.25 ncol_leg = 2 # 1 legend = ax.legend(labelspacing=labelspacing, borderpad=borderpad, ncol=ncol_leg, handletextpad=handletextpad, loc='upper right', frameon=frameon, framealpha=framealpha, edgecolor=edgecolor, numpoints=1, scatterpoints=1,fontsize=fontsize_leg_tmp) if fileout is not None: plt.savefig(fileout, bbox_inches='tight', dpi=600) plt.close() else: plt.show() return None # ---------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------- # Figure 5 def plot_example_galaxy_mencl_vcirc(bt_arr=[0., 0.25, 0.5, 0.75, 1.], output_path=None, table_path=None, fileout=None, z=2., Mbaryon=6.6e10, Reff_disk=3.4, n_disk=1., invq_disk=4., Reff_bulge=1., n_bulge=4., invq_bulge=1., Mhalo=8.9e11, halo_conc=4.2, rmin = 0., rmax=15., rstep=0.01, log_rmin=-0.6, log_rmax=2.2, nlogr = 101, logradius=False, ylim_lmenc=[8., 12.], ylim_lmenc_lograd=[6.,12.5], ylim_vcirc=[0., 420.]): """ Plot example enclosed mass and circular velocity profiles for different mass components, over a variety of B/T ratios. Saves plot to PDF. Parameters ---------- output_path: str Path to directory where the output plot will be saved. table_path: str Path to directory containing the Sersic profile tables. bt_arr: array_like, optional Array of B/T ratios to plot, in separate panels. Default: [0., 0.25, 0.5, 0.75] z: float, optional Redshift (to determine NFW halo properties. Default: z=2. Mbaryon: float, optional Total baryon mass [Msun]. Default: 6.6e10 Msun. Reff_disk: float, optional Sersic projected 2D half-light (assumed half-mass) radius of disk component [kpc]. Default: 3.4 kpc n_disk: float, optional Sersic index of disk component. Default: n_disk = 1. invq_disk: float, optional Flattening of disk component. Default: invq_disk = 4. Reff_bulge: float, optional Sersic projected 2D half-light (assumed half-mass) radius of bulge component [kpc]. Default: 1kpc n_bulge: float, optional Sersic index of bulge component. Default: n_bulge = 4. invq_bulge: float, optional Flattening of bulge component. Default: invq_bulge = 1. (spherical) Mhalo: float, optional NFW halo mass within R200 (=M200) [Msun]. Default: 8.9e112 Msun halo_conc: float, optional Concentration of NFW halo. Default: 4.2 logradius: bool, optional Option whether to plot log radius or linear. Default: False (plot linear). rmin: float, optional Minimum radius, if doing linear radius. Ignored if logradius=True. [kpc] rmax: float, optional Maximum radius, if doing linear radius. Ignored if logradius=True. [kpc] rstep: float, optional Radius stepsize, if doing linear radius.Ignored if logradius=True. [kpc] log_rmin: float, optional Log of minimum radius, if doing log radius. Ignored if logradius=False. [log(kpc)] log_rmax: float, optional Log of maximum radius, if doing log radius. Ignored if logradius=False. [log(kpc)] nlogr: int, optional Number of log radius steps, if logradius=True. Ignored if logradius=False. fileout: str, optional Override the default filename and explicitly choose the output filename (must include full path). """ if (output_path is None) & (fileout is None): raise ValueError("Must set 'output_path' if 'fileout' is not set !") if table_path is None: raise ValueError("Must set 'table_path' !") if (fileout is None): # Ensure trailing slash: if output_path[-1] != '/': output_path += '/' fileout = output_path+'plot_example_galaxy_mencl_vcirc' for bt in bt_arr: fileout += '_bt{:0.2f}'.format(bt) if logradius: fileout += '_logradius' fileout += '.pdf' if logradius: r_arr = np.logspace(log_rmin, log_rmax, num=nlogr) Rarr_plot = np.log10(r_arr) else: r_arr = np.arange(rmin, rmax+rstep, rstep) Rarr_plot = r_arr q_disk = 1./invq_disk # ++++++++++++++++ # plot: color_fdm = 'silver' color_vsq = 'dimgrey' types = ['menc3D', 'vcirc', 'fDM'] ls_arr = ['--', '-.', (0, (3, 1, 1, 1, 1, 1)), ':', '-', (0, (10, 4)), (0, (10, 2, 1, 2, 1, 2, 1, 2))] lw_arr = [1.25, 1.25, 1.35, 1.35, 1.5, 1., 1.] labels = ['Disk', 'Bulge', 'Baryons', 'Halo', 'Total', r'$M_{\mathrm{enc,DM}}/M_{\mathrm{enc,tot}}$', r'$v_{\mathrm{circ,DM}}^2/v_{\mathrm{circ,tot}}^2$'] color_arr = ['blue', 'red', 'green', 'purple', 'black', color_fdm, color_vsq] titles, xlims, ylims = ([] for _ in range(3)) xlabels, ylabels, ylims2, ylabels2 = ([] for _ in range(4)) for i, typ in enumerate(types): for j, bt in enumerate(bt_arr): if i == 0: titles.append(r'$B/T={:0.2f}$'.format(bt)) else: titles.append(None) if i == len(types) - 1: if logradius: xlabels.append(r'$\log_{10}(R/\mathrm{[kpc]})$') else: xlabels.append(r'$R\ \mathrm{[kpc]}$') else: xlabels.append(None) if j == 0: if typ == 'menc3D': ylabels.append(r'$\log_{10}(M_{\mathrm{sph}}(<r=R)/M_{\odot})$') elif typ == 'vcirc': ylabels.append(r'$v_{\mathrm{circ}}(R)\ \mathrm{[km\,s^{-1}]}$') elif typ == 'fDM': ylabels.append(r'$f_{\mathrm{DM}}(R)$') else: ylabels.append(None) if j == len(bt_arr) -1: ylabels2.append(r'$\mathrm{frac}$') else: ylabels2.append(None) if logradius: xlims.append([log_rmin, log_rmax]) else: xlims.append([rmin, rmax+0.025*(rmax-rmin)]) if typ == 'menc3D': if logradius: ylims.append(ylim_lmenc_lograd) else: ylims.append(ylim_lmenc) elif typ == 'vcirc': ylims.append(ylim_vcirc) if typ == 'fDM': ylims.append([0.,1.]) ylims2.append([0.,1.]) ###################################### # Setup plot: f = plt.figure() scale = 2.35 n_cols = len(bt_arr) n_rows = len(types) fac = 1.02 f.set_size_inches(fac*scale*n_cols,scale*n_rows) wspace = 0.05 hspace = wspace gs = gridspec.GridSpec(n_rows, n_cols, wspace=wspace, hspace=hspace) axes = [] for i in range(n_rows): for j in range(n_cols): axes.append(plt.subplot(gs[i,j])) for i in range(n_rows): for j in range(n_cols): k = i*n_cols + j typ = types[i] bt = bt_arr[j] ax = axes[k] xlim = xlims[k] ylim = ylims[k] ylabel = ylabels[k] title = titles[k] ############### ## Get mass, velocity components: menc_disk = interp_profiles.interpolate_sersic_profile_menc(R=r_arr, total_mass=(1.-bt)*Mbaryon, Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path) menc_bulge = interp_profiles.interpolate_sersic_profile_menc(R=r_arr, total_mass=(bt)*Mbaryon, Reff=Reff_bulge, n=n_bulge, invq=invq_bulge, path=table_path) nfw = plot_calcs.NFW(z=z, Mvir=Mhalo, conc=halo_conc) menc_halo = nfw.enclosed_mass(r_arr) menc_baryons = menc_disk + menc_bulge menc_tot = menc_baryons + menc_halo fdm_menc = menc_halo/menc_tot vcirc_disk = interp_profiles.interpolate_sersic_profile_VC(R=r_arr, total_mass=(1.-bt)*Mbaryon, Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path) vcirc_bulge = interp_profiles.interpolate_sersic_profile_VC(R=r_arr, total_mass=(bt)*Mbaryon, Reff=Reff_bulge, n=n_bulge, invq=invq_bulge, path=table_path) vcirc_halo = nfw.v_circ(r_arr) vcirc_baryons = np.sqrt(vcirc_disk**2 + vcirc_bulge**2) vcirc_tot = np.sqrt(vcirc_baryons**2 + vcirc_halo**2) fdm_vsq = vcirc_halo**2/vcirc_tot**2 ######### # Get key radii: r3Dhalf_disk = util_calcs.find_rhalf3D_sphere(R=r_arr, menc3D_sph=menc_disk, total_mass=(1.-bt)*Mbaryon) r3Dhalf_bulge = util_calcs.find_rhalf3D_sphere(R=r_arr, menc3D_sph=menc_bulge, total_mass=(bt)*Mbaryon) r3Dhalf_bary = util_calcs.find_rhalf3D_sphere(R=r_arr, menc3D_sph=menc_baryons, total_mass=Mbaryon) # Get fDM at these radii: vcD_Redisk = interp_profiles.interpolate_sersic_profile_VC(R=Reff_disk, total_mass=(1.-bt)*Mbaryon, Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path) vcB_Redisk = interp_profiles.interpolate_sersic_profile_VC(R=Reff_disk, total_mass=(bt)*Mbaryon, Reff=Reff_bulge, n=n_bulge, invq=invq_bulge, path=table_path) vcH_Redisk = nfw.v_circ(Reff_disk) fdm_vsq_Redisk = vcH_Redisk**2/(vcD_Redisk**2 + vcB_Redisk**2 + vcH_Redisk**2) mencD_r3DhalfB = interp_profiles.interpolate_sersic_profile_menc(R=r3Dhalf_bary, total_mass=(1.-bt)*Mbaryon, Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path) mencB_r3DhalfB = interp_profiles.interpolate_sersic_profile_menc(R=r3Dhalf_bary, total_mass=(bt)*Mbaryon, Reff=Reff_bulge, n=n_bulge, invq=invq_bulge, path=table_path) mencH_r3DhalfB = nfw.enclosed_mass(r3Dhalf_bary) fdm_menc_r3Dhalf_bary = mencH_r3DhalfB / (mencD_r3DhalfB + mencB_r3DhalfB + mencH_r3DhalfB) if logradius: vline_loc = [np.log10(Reff_disk), np.log10(r3Dhalf_disk), np.log10(r3Dhalf_bulge), np.log10(r3Dhalf_bary), np.log10(nfw.rvir)] vline_lss = ['-', '--', '-.', (0, (3, 1, 1, 1, 1, 1)), ':'] vline_cols = ['darkgrey', 'blue', 'red', 'green', 'grey'] vline_labels = [r'$R_{\mathrm{e,disk}}$', r'$r_{\mathrm{1/2,3D,disk}}$', r'$r_{\mathrm{1/2,3D,bulge}}$', r'$r_{\mathrm{1/2,3D,baryon}}$', r'$r_{\mathrm{vir}}$'] vline_alphas = [1., 0.5, 0.5, 0.5, 1.] else: vline_loc = [Reff_disk, r3Dhalf_disk, r3Dhalf_bulge, r3Dhalf_bary] vline_lss = ['-', '--', '-.', (0, (3, 1, 1, 1, 1, 1))] vline_cols = ['darkgrey', 'blue', 'red', 'green'] vline_labels = [r'$R_{\mathrm{e,disk}}$', r'$r_{\mathrm{1/2,3D,disk}}$', r'$r_{\mathrm{1/2,3D,bulge}}$', r'$r_{\mathrm{1/2,3D,baryon}}$'] vline_alphas = [1., 0.5, 0.5, 0.5] ######### if typ == 'menc3D': yarrs = [np.log10(menc_disk), np.log10(menc_bulge), np.log10(menc_baryons), np.log10(menc_halo), np.log10(menc_tot)] asymp_vals = [np.log10((1.-bt)*Mbaryon), np.log10(bt*Mbaryon), np.log10(Mbaryon), np.log10(Mhalo), np.NaN] fdm = fdm_menc elif typ == 'vcirc': yarrs = [vcirc_disk, vcirc_bulge, vcirc_baryons, vcirc_halo, vcirc_tot] fdm = fdm_vsq elif typ == 'fDM': yarrs = [] for mm in range(len(yarrs)): ax.plot(Rarr_plot, yarrs[mm], ls=ls_arr[mm], color=color_arr[mm], lw=lw_arr[mm], label=labels[mm]) # Plot fDM: if typ in ['menc3D', 'vcirc']: ax2 = ax.twinx() ax2.set_zorder(-1) ax.patch.set_visible(False) if typ == 'menc3D': mm += 1 ax2.plot(Rarr_plot, fdm_menc, ls=ls_arr[mm], color=color_arr[mm], lw=lw_arr[mm], label=labels[mm], zorder=-10.) elif typ == 'vcirc': mm += 2 ax2.plot(Rarr_plot, fdm_vsq, ls=ls_arr[mm], color=color_arr[mm], lw=lw_arr[mm], label=labels[mm], zorder=-12.) else: ax2 = None if typ == 'fDM': mm = 4 mm += 1 ax.plot(Rarr_plot, fdm_menc, ls=ls_arr[mm], color=color_arr[mm], lw=lw_arr[mm]+.5, label=labels[mm], zorder=-10.) mm += 1 ax.plot(Rarr_plot, fdm_vsq, ls=ls_arr[mm], color=color_arr[mm], lw=lw_arr[mm]+.5, label=labels[mm], zorder=-12.) if ylim is None: ylim = ax.get_ylim() for vl, vls, vlc, vlb, vah in zip(vline_loc, vline_lss, vline_cols, vline_labels, vline_alphas): if ((i == 0) & (j == 1)): pass else: vlb = None ax.axvline(x=vl, ls=vls, color=vlc, label=vlb, zorder=-20., alpha=vah) if typ == 'fDM': # fdm_vsq_Redisk, fdm_menc_r3Dhalf_bary fac = 0.8 ind = 0 ax.axhline(y=fdm_vsq_Redisk, color=vline_cols[ind], ls=vline_lss[ind], alpha=vline_alphas[ind]*fac, zorder=-20) ind = 3 ax.axhline(y=fdm_menc_r3Dhalf_bary, color=vline_cols[ind], ls=vline_lss[ind], alpha=vline_alphas[ind]*fac, zorder=-20) ax.set_xlim(xlim) ax.set_ylim(ylim) ######## if ax2 is not None: ax2.set_ylim(ylims2[k]) ax2.yaxis.set_minor_locator(MultipleLocator(0.05)) ax2.tick_params(axis='y', direction='in', color=color_vsq, which='both') if (j == n_cols-1): if ylabels2[k] is not None: ax2.set_ylabel(ylabels2[k], fontsize=fontsize_labels_sm, color=color_vsq) if (i > 0): ax2.yaxis.set_major_locator(FixedLocator([0., 0.2, 0.4, 0.6, 0.8, 1.])) ax2.yaxis.set_major_formatter(FixedFormatter(["0.0", "0.2", "0.4", "0.6", "0.8", ""])) else: ax2.yaxis.set_major_locator(MultipleLocator(0.2)) ax2.tick_params(axis='y', direction='in', color=color_vsq, labelsize=fontsize_ticks_sm, colors=color_vsq) else: ax2.yaxis.set_major_locator(MultipleLocator(0.2)) ax2.set_yticklabels([]) ######## if logradius: ax.xaxis.set_minor_locator(MultipleLocator(0.1)) if (j > 0) & (i == n_rows-1): ax.xaxis.set_major_locator(FixedLocator([-2., -1., 0., 1., 2.])) ax.xaxis.set_major_formatter(FixedFormatter(["", "-1", "0", "1", "2"])) else: ax.xaxis.set_major_locator(MultipleLocator(1.)) else: ax.xaxis.set_minor_locator(MultipleLocator(1.)) if (j > 0) & (i == n_rows-1): ax.xaxis.set_major_locator(FixedLocator([0., 5., 10., 15.])) ax.xaxis.set_major_formatter(FixedFormatter(["", "5", "10", "15"])) else: ax.xaxis.set_major_locator(MultipleLocator(5.)) if typ == 'menc3D': ax.yaxis.set_minor_locator(MultipleLocator(0.2)) ax.yaxis.set_major_locator(MultipleLocator(1.)) elif typ == 'vcirc': ax.yaxis.set_minor_locator(MultipleLocator(20.)) ax.yaxis.set_major_locator(MultipleLocator(100.)) elif typ == 'fDM': ax.yaxis.set_minor_locator(MultipleLocator(0.05)) if (j==0): ax.yaxis.set_major_locator(FixedLocator([0., 0.2, 0.4, 0.6, 0.8, 1.])) ax.yaxis.set_major_formatter(FixedFormatter(["0.0", "0.2", "0.4", "0.6", "0.8", ""])) else: ax.yaxis.set_major_locator(MultipleLocator(0.2)) ax.tick_params(labelsize=fontsize_ticks) if xlabels[k] is not None: ax.set_xlabel(xlabels[k], fontsize=fontsize_labels) else: ax.set_xticklabels([]) if ylabels[k] is not None: ax.set_ylabel(ylabels[k], fontsize=fontsize_labels) else: ax.set_yticklabels([]) if titles[k] is not None: ax.set_title(titles[k], fontsize=fontsize_title) if (k == 0) | ((i == n_rows -1) & (j == 0)) | ((i == 0) & (j == 1)): handles, labels_leg = ax.get_legend_handles_labels() neworder = range(len(yarrs)) handles_arr = [] labels_arr = [] for ii in neworder: handles_arr.append(handles[ii]) labels_arr.append(labels_leg[ii]) handlesax2, labelsax2 = ax.get_legend_handles_labels() neworder2 = range(len(handlesax2)) handles_arr2 = [] labels_arr2 = [] for ii in neworder2: handles_arr2.append(handlesax2[ii]) labels_arr2.append(labelsax2[ii]) neworder3 = range(len(yarrs), len(handles)) handles_arr3 = [] labels_arr3 = [] for ii in neworder3: handles_arr3.append(handles[ii]) labels_arr3.append(labels_leg[ii]) frameon = True framealpha = 1. edgecolor = 'none' borderpad = 0.25 fontsize_leg_tmp = fontsize_leg labelspacing=0.15 handletextpad=0.25 loc3 = 'lower right' if logradius: loc2 = 'upper left' else: loc2 = 'upper right' bbox2_to_anchor = (1.,1.) if (k == 0): legend1 = ax.legend([handles_arr[-1]], [labels_arr[-1]], labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, loc='upper right', bbox_to_anchor=(0.93,1.02), handlelength=1.55, numpoints=1, scatterpoints=1, frameon=False, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp) legend15 = ax.legend(handles_arr[:-1], labels_arr[:-1], labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, loc='lower right', bbox_to_anchor=(1.04,-0.04), handlelength=1.55, numpoints=1, scatterpoints=1, frameon=False, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp) ax.add_artist(legend1) ax.add_artist(legend15) elif ((i == n_rows -1) & (j == 0)): legend2 = ax.legend(handles_arr2, labels_arr2, labelspacing=labelspacing, borderpad=borderpad, bbox_to_anchor=bbox2_to_anchor, handletextpad=handletextpad, loc=loc2, numpoints=1, scatterpoints=1, handlelength=3.35, frameon=frameon, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp) ax.add_artist(legend2) elif ((i == 0) & (j == 1)): legend3 = ax.legend([handles_arr3[0]], [labels_arr3[0]], labelspacing=0, borderpad=borderpad, handletextpad=handletextpad, loc='upper right', bbox_to_anchor=(0.83, 1.02), handlelength=1.4, numpoints=1, scatterpoints=1, frameon=False, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp) legend35 = ax.legend(handles_arr3[1:], labels_arr3[1:], labelspacing=0, borderpad=borderpad, handletextpad=handletextpad, loc= loc3, bbox_to_anchor=(1.04,-0.05), handlelength=1.4, numpoints=1, scatterpoints=1, frameon=False, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp) # Hack to get alpha more to match for lh, alph in zip(legend3.legendHandles, [vline_alphas[0]]): if alph < 1.: alph *= 0.65 lh.set_alpha(alph) for lh, alph in zip(legend35.legendHandles, vline_alphas[1:]): if alph < 1.: alph *= 0.65 lh.set_alpha(alph) # Change label v alignment for t in legend35.get_texts(): t.set_position((0.,2.)) ax.add_artist(legend3) ax.add_artist(legend35) if fileout is not None: plt.savefig(fileout, bbox_inches='tight', dpi=600) plt.close() else: plt.show() return None # ---------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------- # Figure 6 def plot_fdm_calibration(Mstar_arr=None, Mbaryon_arr=[10**10.5], q_disk_arr = [0.01, 0.05, 0.1, 0.2, 0.25, 0.4, 0.6, 0.8, 1., 1.5, 2.], Mhalo_arr=[1.e12], halo_conc_arr=[4.], Reff_disk_arr = [5.], output_path=None, table_path=None, fileout=None, z=2., n_disk=1., Reff_bulge=1., n_bulge=4., invq_bulge=1., del_fDM=False): """ Plot the calibration ratio between fDM(vcirc,Reff) / fDM(Menc,Reff) Saves plot to PDF. Parameters ---------- output_path: str Path to directory where the output plot will be saved. table_path: str Path to directory containing the Sersic profile tables. z: float, optional Redshift (to determine NFW halo properties). Default: z=2. Mbaryon_arr: array_like, optional Array of total baryon mass to plot [Msun]. Default: [10**10.5] Msun. q_disk_arr: array_like, optional Array of flattening for disk. Default: [0.01, 0.05, 0.1, 0.2, 0.25, 0.4, 0.6, 0.8, 1., 1.5, 2.] Mhalo_arr: array_like, optional Array of NFW halo masses within R200 (=M200) [Msun]. Default: [1.e12] Msun. halo_conc_arr: array_like, optional Array of halo concentrations. Default: [4.] Reff_disk: array_like, optional Sersic projected 2D half-light (assumed half-mass) radius of disk component [kpc]. Default: [5.] kpc n_disk: float, optional Sersic index of disk component. Default: n_disk = 1. Reff_bulge: float, optional Sersic projected 2D half-light (assumed half-mass) radius of bulge component [kpc]. Default: 1kpc n_bulge: Sersic index of bulge component. Default: n_bulge = 4. invq_bulge: Flattening of bulge componen. Default: invq_bulge = 1. (spherical) fileout: str, optional Override the default filename and explicitly choose the output filename (must include full path). """ if (output_path is None) & (fileout is None): raise ValueError("Must set 'output_path' if 'fileout' is not set !") if table_path is None: raise ValueError("Must set 'table_path' !") if (fileout is None): # Ensure trailing slash: if output_path[-1] != '/': output_path += '/' fileout = output_path+'plot_fdm_calibration' if del_fDM: fileout += '_del_fDM' fileout += '.pdf' # ++++++++++++++++ # SETUP: if Mstar_arr is not None: # Override the other settings with "typical" values: Mbaryon_arr, Mhalo_arr, halo_conc_arr, Reff_disk_arr = ([] for _ in range(4)) for Mstar in Mstar_arr: lmstar = np.log10(Mstar) Reff_disk = scaling_rel._mstar_Reff_relation(z=z, lmstar=lmstar, galtype='sf') fgas = scaling_rel._fgas_scaling_relation_MS(z=z, lmstar=lmstar) Mhalo = scaling_rel._smhm_relation(z=z, lmstar=lmstar) halo_conc = scaling_rel._halo_conc_relation(z=z, lmhalo=np.log10(Mhalo)) Mbaryon_arr.append(Mstar/(1.-fgas)) Mhalo_arr.append(Mhalo) halo_conc_arr.append(halo_conc) Reff_disk_arr.append(Reff_disk) # ++++++++++++++++ # plot: btstep = 0.05 bt_arr=np.arange(0.,1.+btstep, btstep) q_disk_arr = np.array(q_disk_arr) invq_disk_arr = 1./q_disk_arr color_arr, ls_arr, labels, lw_arr = ([] for _ in range(4)) for q in q_disk_arr: if q <= 1.: color_arr.append(cmap_q(q)) ls_arr.append('-') else: color_arr.append(cmapg(1./q)) ls_arr.append('--') lw_arr.append(1.25) labels.append(r'$q_{0,\mathrm{disk}}'+r'={}$'.format(q)) lss_bonus = ['-', '--', '-.'] titles, xlims, ylims, xlabels, ylabels = ([] for _ in range(5)) for j, Mbaryon in enumerate(Mbaryon_arr): titles.append(r'$\log_{10}(M_{\mathrm{baryon}}/M_{\odot})'+r'={:0.1f}$'.format(np.log10(Mbaryon))) xlabels.append(r'$B/T$') if j == 0: if del_fDM: ylabels.append(r'$[ (f_{\mathrm{DM}}^{v}-f_{\mathrm{DM}}^{m})/f_{\mathrm{DM}}^{m}](R_{\mathrm{e,disk}})$') else: ylabels.append(r'$f_{\mathrm{DM}}^{v}(R_{\mathrm{e,disk}})/f_{\mathrm{DM}}^{m}(R_{\mathrm{e,disk}})$') else: ylabels.append(None) xlims.append([0., 1.]) if del_fDM: ylims.append([-0.24, 0.07]) else: ylims.append([0.77, 1.045]) ###################################### # Setup plot: f = plt.figure() scale = 4. n_cols = len(Mbaryon_arr) n_rows = 1 fac = 1.02 f.set_size_inches(fac*scale*n_cols,scale*n_rows) wspace = 0.025 hspace = wspace gs = gridspec.GridSpec(n_rows, n_cols, wspace=wspace, hspace=hspace) axes = [] for i in range(n_rows): for j in range(n_cols): axes.append(plt.subplot(gs[i,j])) ###################### # Load bulge: n, invq don't change. tab_bulge = io.read_profile_table(n=n_bulge, invq=invq_bulge, path=table_path) tab_bulge_menc = tab_bulge['menc3D_sph'] tab_bulge_vcirc = tab_bulge['vcirc'] tab_bulge_rad = tab_bulge['R'] tab_bulge_Reff = tab_bulge['Reff'] tab_bulge_mass = tab_bulge['total_mass'] # Clean up values inside rmin: Add the value at r=0: menc=0 if tab_bulge['R'][0] > 0.: tab_bulge_rad = np.append(0., tab_bulge_rad) tab_bulge_menc = np.append(0., tab_bulge_menc) tab_bulge_vcirc = np.append(0., tab_bulge_vcirc) m_interp_bulge = scp_interp.interp1d(tab_bulge_rad, tab_bulge_menc, fill_value=np.NaN, bounds_error=False, kind='cubic') v_interp_bulge = scp_interp.interp1d(tab_bulge_rad, tab_bulge_vcirc, fill_value=np.NaN, bounds_error=False, kind='cubic') ###################### plot_cnt_q = 0 for mm, invq_disk in enumerate(invq_disk_arr): if invq_disk in invq_disk_arr: ## Get mass, velocity components: # FASTER: ###################### tab_disk = io.read_profile_table(n=n_disk, invq=invq_disk, path=table_path) tab_disk_menc = tab_disk['menc3D_sph'] tab_disk_vcirc = tab_disk['vcirc'] tab_disk_rad = tab_disk['R'] tab_disk_Reff = tab_disk['Reff'] tab_disk_mass = tab_disk['total_mass'] # Clean up values inside rmin: Add the value at r=0: menc=0 if tab_disk['R'][0] > 0.: tab_disk_rad = np.append(0., tab_disk_rad) tab_disk_menc = np.append(0., tab_disk_menc) tab_disk_vcirc = np.append(0., tab_disk_vcirc) m_interp_disk = scp_interp.interp1d(tab_disk_rad, tab_disk_menc, fill_value=np.NaN, bounds_error=False, kind='cubic') v_interp_disk = scp_interp.interp1d(tab_disk_rad, tab_disk_vcirc, fill_value=np.NaN, bounds_error=False, kind='cubic') ###################### for i in range(n_rows): for j in range(n_cols): k = i*n_cols + j Mbaryon = Mbaryon_arr[j] Reff_disk = Reff_disk_arr[j] Mhalo = Mhalo_arr[j] halo_conc = halo_conc_arr[j] ax = axes[k] xlim = xlims[k] ylim = ylims[k] ylabel = ylabels[k] title = titles[k] ###################### fdm_menc = np.ones(len(bt_arr))* -99. fdm_vsq = np.ones(len(bt_arr))* -99. for ll, bt in enumerate(bt_arr): menc_disk = (m_interp_disk(Reff_disk / Reff_disk * tab_disk_Reff) * (((1.-bt)*Mbaryon) / tab_disk_mass) ) vcirc_disk = (v_interp_disk(Reff_disk / Reff_disk * tab_disk_Reff) * np.sqrt(((1.-bt)*Mbaryon) / tab_disk_mass) * np.sqrt(tab_disk_Reff / Reff_disk)) menc_bulge = (m_interp_bulge(Reff_disk / Reff_bulge * tab_bulge_Reff) * ((bt*Mbaryon) / tab_bulge_mass) ) vcirc_bulge = (v_interp_bulge(Reff_disk / Reff_bulge * tab_bulge_Reff) * np.sqrt((bt*Mbaryon) / tab_bulge_mass) * np.sqrt(tab_bulge_Reff / Reff_bulge)) nfw = plot_calcs.NFW(z=z, Mvir=Mhalo, conc=halo_conc) menc_halo = nfw.enclosed_mass(Reff_disk) menc_baryons = menc_disk + menc_bulge menc_tot = menc_baryons + menc_halo fdm_menc[ll] = menc_halo/menc_tot vcirc_halo = nfw.v_circ(Reff_disk) vcirc_baryons = np.sqrt(vcirc_disk**2 + vcirc_bulge**2) vcirc_tot = np.sqrt(vcirc_baryons**2 + vcirc_halo**2) fdm_vsq[ll] = vcirc_halo**2/vcirc_tot**2 if (mm == 0) & (bt == 0): menc_baryons = util_calcs.total_mass2D_direct(Reff_disk, total_mass=Mbaryon, Reff=Reff_disk, n=1., q=1., i=0.) vcirc_baryons = plot_calcs.freeman_vcirc(Reff_disk, Mbaryon, Reff_disk) menc_tot = menc_baryons + menc_halo vcirc_tot = np.sqrt(vcirc_baryons**2 + vcirc_halo**2) fdm_menc_tmp = menc_halo/menc_tot fdm_vsq_tmp = vcirc_halo**2/vcirc_tot**2 print("Freeman: fDMv/fDMm = {:0.4f}".format(fdm_vsq_tmp/fdm_menc_tmp)) # # Do another limiting calc: # qprime = 0.001 # sprof = core.DeprojSersicModel(total_mass=Mbaryon, # Reff=Reff_disk, n=n_disk, q=qprime) # # menc_baryons = sprof.enclosed_mass(Reff_disk) # vcirc_baryons = sprof.v_circ(Reff_disk) # menc_tot = menc_baryons + menc_halo # vcirc_tot = np.sqrt(vcirc_baryons**2 + vcirc_halo**2) # fdm_menc_tmp = menc_halo/menc_tot # fdm_vsq_tmp = vcirc_halo**2/vcirc_tot**2 # # print("Oblate q={}: fDMv/fDMm = {:0.4f}".format(qprime, fdm_vsq_tmp/fdm_menc_tmp)) ###################### if invq_disk < 1: zorder = -5. else: zorder = None if del_fDM: ax.plot(bt_arr, (fdm_vsq-fdm_menc)/fdm_menc, ls=ls_arr[mm], color=color_arr[mm], lw=lw_arr[mm], label=labels[mm], zorder=zorder) else: ax.plot(bt_arr, fdm_vsq/fdm_menc, ls=ls_arr[mm], color=color_arr[mm], lw=lw_arr[mm], label=labels[mm], zorder=zorder) if k > 0: if (mm == (len(invq_disk_arr) -1) ): lbl_tmp = r'$\log_{10}(M_{\mathrm{bar}}/M_{\odot})='+r'{:0.1f}$'.format(np.log10(Mbaryon)) else: lbl_tmp = None axes[0].plot(bt_arr, fdm_vsq/fdm_menc, ls=lss_bonus[k], color=color_arr[mm], lw=lw_arr[mm], label=lbl_tmp) plot_cnt_q += 1 else: # Missing invq_disk file pass #### # Go back and do plot formatting: for i in range(n_rows): for j in range(n_cols): k = i*n_cols + j Mbaryon = Mbaryon_arr[j] Reff_disk = Reff_disk_arr[j] Mhalo = Mhalo_arr[j] halo_conc = halo_conc_arr[j] ax = axes[k] xlim = xlims[k] ylim = ylims[k] ylabel = ylabels[k] title = titles[k] if ylim is None: ylim = ax.get_ylim() # Annotate: if k == 0: padx = pady = 0.03 xypos = (padx, 1.-pady) va = 'top' ha = 'left' ann_str_cnst = r'$R_{\mathrm{e,disk}}='+r'{:0.1f}'.format(Reff_disk)+r'\,\mathrm{kpc}$' ann_str_cnst += '\n' ann_str_cnst += r'$n_{\mathrm{disk}}='+r'{:0.1f}$'.format(n_disk) ax.annotate(ann_str_cnst, xy=xypos, xycoords='axes fraction', ha=ha, va=va, color='darkblue', fontsize=fontsize_ann) pady = 0.03 xdelt = 0.35 ann_str_cnst = r'$R_{\mathrm{e,bulge}}='+r'{:0.1f}'.format(Reff_bulge)+r'\,\mathrm{kpc}$' ann_str_cnst += '\n' ann_str_cnst += r'$n_{\mathrm{bulge}}='+r'{:0.1f}$'.format(n_bulge) xypos = (padx + xdelt, 1.-pady) ax.annotate(ann_str_cnst, xy=xypos, xycoords='axes fraction', ha=ha, va=va, color='firebrick', fontsize=fontsize_ann) xdelt = 0.7 xypos = (padx + xdelt, 1.-pady) ann_str_cnst = r'$q_{0,\mathrm{bulge}}='+r'{:0.1f}$'.format(1./invq_bulge) ax.annotate(ann_str_cnst, xy=xypos, xycoords='axes fraction', ha=ha, va=va, color='firebrick', fontsize=fontsize_ann) else: padx = pady = 0.05 xypos = (padx, 1.-pady) va = 'top' ha = 'left' ann_str_cnst = r'$R_{\mathrm{e,disk}}='+r'{:0.1f}'.format(Reff_disk)+r'\,\mathrm{kpc}$' ax.annotate(ann_str_cnst, xy=xypos, xycoords='axes fraction', ha=ha, va=va, color='darkblue', fontsize=fontsize_ann) ############ padx = 0.045 pady = 0.025 xypos = (padx, pady) va = 'bottom' ha = 'left' ann_str = r'$\log_{10}(M_{\mathrm{bar}}/M_{\odot})'+r'={:0.2f}$'.format(np.log10(Mbaryon)) ann_str += '\n' ann_str += r'$\log_{10}(M_{\mathrm{halo}}/M_{\odot})'+r'={:0.2f}$'.format(np.log10(Mhalo)) ann_str += '\n' ann_str += r'$\mathrm{conc}_{\mathrm{halo}}'+r'={:0.1f}$'.format(halo_conc) ax.annotate(ann_str, xy=xypos, xycoords='axes fraction', ha=ha, va=va, color='black', fontsize=fontsize_ann) if del_fDM: ax.axhline(y=0., ls=(0, (5,3)), color='darkgrey', zorder=-20.) else: ax.axhline(y=1., ls=(0, (5,3)), color='darkgrey', zorder=-20.) ax.set_xlim(xlim) ax.set_ylim(ylim) ######## ax.xaxis.set_minor_locator(MultipleLocator(0.05)) if (j > 0) & (i == n_rows-1): bt_loc = [0., 0.2, 0.4, 0.6, 0.8, 1.] bt_loc_str = ["{:0.1f}".format(bt) for bt in bt_loc] bt_loc_str[0] = "" ax.xaxis.set_major_locator(FixedLocator(bt_loc)) ax.xaxis.set_major_formatter(FixedFormatter(bt_loc_str)) else: ax.xaxis.set_major_locator(MultipleLocator(0.2)) ax.yaxis.set_minor_locator(MultipleLocator(0.01)) ax.yaxis.set_major_locator(MultipleLocator(0.05)) ax.tick_params(labelsize=fontsize_ticks) if xlabels[k] is not None: ax.set_xlabel(xlabels[k], fontsize=fontsize_labels) else: ax.set_xticklabels([]) if ylabels[k] is not None: ax.set_ylabel(ylabels[k], fontsize=fontsize_labels) else: ax.set_yticklabels([]) if k == 0: handles, labels_leg = ax.get_legend_handles_labels() neworder = range(plot_cnt_q) handles_arr = [] labels_arr = [] for ii in neworder: handles_arr.append(handles[ii]) labels_arr.append(labels_leg[ii]) neworder2 = range(plot_cnt_q, len(handles)) handles_arr2 = [] labels_arr2 = [] for ii in neworder2: handles_arr2.append(handles[ii]) labels_arr2.append(labels_leg[ii]) frameon = True framealpha = 1. edgecolor = 'none' borderpad = 0.25 fontsize_leg_tmp = fontsize_leg labelspacing=0.15 handletextpad=0.25 legend1 = ax.legend(handles_arr, labels_arr, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, loc='lower right', bbox_to_anchor=(1.01,-0.02), numpoints=1, scatterpoints=1, handlelength=1.65, frameon=False, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp) ax.add_artist(legend1) if len(handles_arr2) > 0: loc2 = (1.05, 0.5) legend2 = ax.legend(handles_arr2, labels_arr2, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, loc=loc2, numpoints=1, scatterpoints=1, frameon=frameon, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp) ax.add_artist(legend2) ax.set_zorder(1) if fileout is not None: plt.savefig(fileout, bbox_inches='tight', dpi=600) plt.close() else: plt.show() return None # ---------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------- # Figure 7 def plot_rhalf3DReD_fdm_calibration_ReB_ReD_rhalf3D(Mbaryon=10**10.5, q_disk_arr = [0.05, 0.1, 0.2, 0.25, 0.4, 0.6, 0.8, 1., 1.5, 2.], ReB_to_ReD_arr=[0.2, 0.5, 1.], Mhalo=1.e12, halo_conc=4., output_path=None, table_path=None, fileout=None, z=2., n_disk=1., n_bulge=4., invq_bulge=1., Reff_bulge=1., rmin = 0., rmax=15., rstep=0.01, btstep = 0.025, del_fDM=True): """ Plot the calibration ratio between fDM(vcirc,Reff) / fDM(Menc,Reff), for a range of ReB/ReD Saves plot to PDF. Parameters ---------- output_path: str Path to directory where the output plot will be saved. table_path: str Path to directory containing the Sersic profile tables. z: float, optional Redshift (to determine NFW halo properties). Default: z=2. Mbaryon: float, optional Total baryon mass to plot [Msun]. Default: 10**10.5 Msun. q_disk_arr: array_like, optional Array of flattening for disk. Default: [0.01, 0.05, 0.1, 0.2, 0.25, 0.4, 0.6, 0.8, 1., 1.5, 2.] Mhalo: float, optional NFW halo mass within R200 (=M200) [Msun]. Default: 1.e12 Msun. halo_conc_arr: flat, optional Halo concentration. Default: 4. n_disk: float, optional Sersic index of disk component. Default: n_disk = 1. Reff_bulge: float, optional Sersic projected 2D half-light (assumed half-mass) radius of bulge component [kpc]. Default: 1kpc n_bulge: Sersic index of bulge component. Default: n_bulge = 4. invq_bulge: Flattening of bulge componen. Default: invq_bulge = 1. (spherical) ReB_to_ReD_arr: array_like, optional Ratio of bulge to disk effective radii. Default: [0.2, 0.5, 1.] fileout: str, optional Override the default filename and explicitly choose the output filename (must include full path). """ if (output_path is None) & (fileout is None): raise ValueError("Must set 'output_path' if 'fileout' is not set !") if table_path is None: raise ValueError("Must set 'table_path' !") if (fileout is None): # Ensure trailing slash: if output_path[-1] != '/': output_path += '/' fileout = output_path+'plot_rhalf3DReD_fdm_calibration_ReB_ReD_rhalf3D' if del_fDM: fileout += '_del_fDM' fileout += '.pdf' # ++++++++++++++++ # plot: bt_arr=np.arange(0.,1.+btstep, btstep) q_disk_arr = np.array(q_disk_arr) invq_disk_arr = 1./q_disk_arr color_arr, ls_arr, labels, lw_arr = ([] for _ in range(4)) for q in q_disk_arr: if q <= 1.: color_arr.append(cmap_q(q)) ls_arr.append('-') else: color_arr.append(cmapg(1./q)) ls_arr.append('--') lw_arr.append(1.25) labels.append(r'$q_{0,\mathrm{disk}}'+r'={}$'.format(q)) lss_bonus = ['-', '--', '-.'] # Redo LS arr: ls_arr, dashes_arr, lw_arr, labels_reratio = ([] for _ in range(4)) dashlen = [5, 10, 20] if len(ReB_to_ReD_arr) > 5: dashlen = np.arange(5,(len(ReB_to_ReD_arr)-1)*5, 5) elif len(ReB_to_ReD_arr) < 5: dashlen = dashlen[:len(ReB_to_ReD_arr)-2] dashlen = dashlen[::-1] # reverse for i, reBtoreD in enumerate(ReB_to_ReD_arr): if i == 0: lw_arr.append(1.3) ls_arr.append('-') dashes_arr.append(None) elif i == len(ReB_to_ReD_arr)-1: lw_arr.append(1.) ls_arr.append(':') dashes_arr.append(None) else: lw_arr.append(1.) ls_arr.append('--') dashes_arr.append((dashlen[i-1], 5)) lblstr = r'$R_{\mathrm{e,bulge}}/R_{\mathrm{e,disk}}=' if reBtoreD % 1. == 0: lblstr += '{:0.0f}$'.format(reBtoreD) elif 10.*reBtoreD % 1. == 0: lblstr += '{:0.1f}$'.format(reBtoreD) else: lblstr += r'{:0.2f}$'.format(reBtoreD) labels_reratio.append(lblstr) xlabels = [r'$B/T$', r'$B/T$'] ylabels = [r'$r_{\mathrm{1/2,3D,baryons}}/R_{\mathrm{e,disk}}$'] if del_fDM: ylabels.append(r'$f_{\mathrm{DM}}^{v}(R_{\mathrm{e,disk}})-f_{\mathrm{DM}}^{m}(r_{1/2,\mathrm{3D,baryons}})$') else: ylabels.append(r'$f_{\mathrm{DM}}^{v}(R_{\mathrm{e,disk}})/f_{\mathrm{DM}}^{m}(r_{1/2,\mathrm{3D,baryons}})$') xlims = [[0., 1.], [0., 1.]] ylims = [[0., 10.]] ylims = [None] if del_fDM: ylims.append([-0.11,0.25]) else: ylims.append(None) ###################################### # Setup plot: f = plt.figure() scale = 3.85 #4. n_cols = 2 n_rows = 1 fac = 1. pad = 0.25 f.set_size_inches(fac*scale*(n_cols+pad*(n_cols-1)),scale*n_rows) wspace = 0.22 hspace = wspace gs = gridspec.GridSpec(n_rows, n_cols, wspace=wspace, hspace=hspace) axes = [] for i in range(n_rows): for j in range(n_cols): axes.append(plt.subplot(gs[i,j])) ###################### # Load bulge: n, invq don't change. tab_bulge = io.read_profile_table(n=n_bulge, invq=invq_bulge, path=table_path) tab_bulge_menc = tab_bulge['menc3D_sph'] tab_bulge_vcirc = tab_bulge['vcirc'] tab_bulge_rad = tab_bulge['R'] tab_bulge_Reff = tab_bulge['Reff'] tab_bulge_mass = tab_bulge['total_mass'] # Clean up values inside rmin: Add the value at r=0: menc=0 if tab_bulge['R'][0] > 0.: tab_bulge_rad = np.append(0., tab_bulge_rad) tab_bulge_menc = np.append(0., tab_bulge_menc) tab_bulge_vcirc = np.append(0., tab_bulge_vcirc) m_interp_bulge = scp_interp.interp1d(tab_bulge_rad, tab_bulge_menc, fill_value=np.NaN, bounds_error=False, kind='cubic') v_interp_bulge = scp_interp.interp1d(tab_bulge_rad, tab_bulge_vcirc, fill_value=np.NaN, bounds_error=False, kind='cubic') r_arr = np.arange(rmin, rmax+rstep, rstep) ###################### plot_cnt_q = 0 for mm, invq_disk in enumerate(invq_disk_arr): if invq_disk in invq_disk_arr: ## Get mass, velocity components: # FASTER: ###################### tab_disk = io.read_profile_table(n=n_disk, invq=invq_disk, path=table_path) tab_disk_menc = tab_disk['menc3D_sph'] tab_disk_vcirc = tab_disk['vcirc'] tab_disk_rad = tab_disk['R'] tab_disk_Reff = tab_disk['Reff'] tab_disk_mass = tab_disk['total_mass'] # Clean up values inside rmin: Add the value at r=0: menc=0 if tab_disk['R'][0] > 0.: tab_disk_rad = np.append(0., tab_disk_rad) tab_disk_menc = np.append(0., tab_disk_menc) tab_disk_vcirc = np.append(0., tab_disk_vcirc) m_interp_disk = scp_interp.interp1d(tab_disk_rad, tab_disk_menc, fill_value=np.NaN, bounds_error=False, kind='cubic') v_interp_disk = scp_interp.interp1d(tab_disk_rad, tab_disk_vcirc, fill_value=np.NaN, bounds_error=False, kind='cubic') ###################### for nn, ReBD in enumerate(ReB_to_ReD_arr): fdm_menc = np.ones(len(bt_arr))* -99. fdm_vsq = np.ones(len(bt_arr))* -99. rhalf3D_arr = np.ones(len(bt_arr))* -99. for ll, bt in enumerate(bt_arr): Reff_disk = Reff_bulge / ReBD menc_bulge_full = interp_profiles.interpolate_sersic_profile_menc(R=r_arr, total_mass=(bt), Reff=Reff_bulge, n=n_bulge, invq=invq_bulge, path=table_path) menc_disk_full = interp_profiles.interpolate_sersic_profile_menc(R=r_arr, total_mass=(1.-bt), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path) rhalf3D = util_calcs.find_rhalf3D_sphere(R=r_arr, menc3D_sph=menc_disk_full+menc_bulge_full, total_mass=1.) rhalf3D_arr[ll] = rhalf3D / Reff_disk menc_disk = (m_interp_disk(rhalf3D / Reff_disk * tab_disk_Reff) * (((1.-bt)*Mbaryon) / tab_disk_mass) ) menc_bulge = (m_interp_bulge(rhalf3D / Reff_bulge * tab_bulge_Reff) * ((bt*Mbaryon) / tab_bulge_mass) ) nfw = plot_calcs.NFW(z=z, Mvir=Mhalo, conc=halo_conc) menc_halo = nfw.enclosed_mass(rhalf3D) menc_baryons = menc_disk + menc_bulge menc_tot = menc_baryons + menc_halo fdm_menc[ll] = menc_halo/menc_tot vcirc_disk = (v_interp_disk(Reff_disk / Reff_disk * tab_disk_Reff) * np.sqrt(((1.-bt)*Mbaryon) / tab_disk_mass) * np.sqrt(tab_disk_Reff / Reff_disk)) vcirc_bulge = (v_interp_bulge(Reff_disk / Reff_bulge * tab_bulge_Reff) * np.sqrt((bt*Mbaryon) / tab_bulge_mass) * np.sqrt(tab_bulge_Reff / Reff_bulge)) vcirc_halo = nfw.v_circ(Reff_disk) vcirc_baryons = np.sqrt(vcirc_disk**2 + vcirc_bulge**2) vcirc_tot = np.sqrt(vcirc_baryons**2 + vcirc_halo**2) fdm_vsq[ll] = vcirc_halo**2/vcirc_tot**2 ###################### for i in range(n_rows): for j in range(n_cols): k = i*n_cols + j ax = axes[k] xlim = xlims[k] ylim = ylims[k] ylabel = ylabels[k] ###################### if invq_disk < 1: zorder = -5. else: zorder = None if k == 0: lbl = None if dashes_arr[nn] is not None: ax.plot(bt_arr, rhalf3D_arr, ls=ls_arr[nn], color=color_arr[mm], lw=lw_arr[nn], label=lbl, dashes=dashes_arr[nn], zorder=zorder) else: ax.plot(bt_arr, rhalf3D_arr, ls=ls_arr[nn], color=color_arr[mm], lw=lw_arr[nn], label=lbl, zorder=zorder) elif k == 1: if nn == 0: lbl = labels[mm] else: lbl = None if del_fDM: delarr =(fdm_vsq-fdm_menc) if dashes_arr[nn] is not None: ax.plot(bt_arr, delarr, ls=ls_arr[nn], color=color_arr[mm], lw=lw_arr[nn], label=lbl, dashes=dashes_arr[nn], zorder=zorder) else: ax.plot(bt_arr, delarr, ls=ls_arr[nn], color=color_arr[mm], lw=lw_arr[nn], label=lbl, zorder=zorder) else: if dashes_arr[nn] is not None: ax.plot(bt_arr, fdm_vsq/fdm_menc, ls=ls_arr[nn], color=color_arr[mm], dashes=dashes_arr[nn], lw=lw_arr[nn], label=lbl, zorder=zorder) else: ax.plot(bt_arr, fdm_vsq/fdm_menc, ls=ls_arr[nn], color=color_arr[mm], lw=lw_arr[nn], label=lbl, zorder=zorder) plot_cnt_q += 1 #### # Go back and do plot formatting: for i in range(n_rows): for j in range(n_cols): k = i*n_cols + j ax = axes[k] xlim = xlims[k] ylim = ylims[k] ylabel = ylabels[k] if ylim is None: ylim = ax.get_ylim() if k == 0: ax.axhline(y=1., ls=(0, (5,3)), color='darkgrey', zorder=-20.) elif k == 1: # Phandom plots for legend: for ii, reBD in enumerate(ReB_to_ReD_arr): lbl = labels_reratio[ii] dashes = dashes_arr[ii] lw = lw_arr[ii] ls = ls_arr[ii] color='black' if dashes is not None: ax.plot(bt_arr, bt_arr*np.NaN, ls=ls, color=color, dashes=dashes, lw=lw, label=lbl) else: ax.plot(bt_arr, bt_arr*np.NaN, ls=ls, color=color, lw=lw, label=lbl) # Annotate: padx = 0.025 pady = 0.07 xypos = (padx, 1-pady-0.125) va = 'top' ha = 'left' ann_str_cnst = r'$n_{\mathrm{disk}}='+r'{:0.1f}$'.format(n_disk) ax.annotate(ann_str_cnst, xy=xypos, xycoords='axes fraction', ha=ha, va=va, color='darkblue', fontsize=fontsize_ann) yoff = 0.185 xypos = (padx, 1-pady-yoff) ann_str_cnst = r'$R_{\mathrm{e,bulge}}='+r'{:0.1f}'.format(Reff_bulge)+r'\,\mathrm{kpc}$' ann_str_cnst += '\n' ann_str_cnst += r'$n_{\mathrm{bulge}}='+r'{:0.1f}$'.format(n_bulge) ann_str_cnst += '\n' ann_str_cnst += r'$q_{0,\mathrm{bulge}}='+r'{:0.1f}$'.format(1./invq_bulge) ax.annotate(ann_str_cnst, xy=xypos, xycoords='axes fraction', ha=ha, va=va, color='firebrick', fontsize=fontsize_ann) ############ if del_fDM: padx = 0.025 pady = 1-0.025 va = 'top' ha = 'left' else: padx = 0.025 pady = 1-0.025 va = 'top' ha = 'left' xypos = (padx, pady) ann_str = r'$\log_{10}(M_{\mathrm{bar}}/M_{\odot})'+r'={:0.2f}$'.format(np.log10(Mbaryon)) ann_str += '\n' ann_str += r'$\log_{10}(M_{\mathrm{halo}}/M_{\odot})'+r'={:0.2f}$'.format(np.log10(Mhalo)) ann_str += '\n' ann_str += r'$\mathrm{conc}_{\mathrm{halo}}'+r'={:0.1f}$'.format(halo_conc) ax.annotate(ann_str, xy=xypos, xycoords='axes fraction', ha=ha, va=va, color='black', fontsize=fontsize_ann) if del_fDM: ax.axhline(y=0., ls=(0, (5,3)), color='lightgrey', zorder=-20.) else: ax.axhline(y=1., ls=(0, (5,3)), color='darkgrey', zorder=-20.) ax.set_xlim(xlim) ax.set_ylim(ylim) ######## ax.xaxis.set_minor_locator(MultipleLocator(0.05)) ax.xaxis.set_major_locator(MultipleLocator(0.2)) if k == 0: ax.yaxis.set_minor_locator(MultipleLocator(0.1)) ax.yaxis.set_major_locator(MultipleLocator(0.5)) elif k == 1: if del_fDM: ax.yaxis.set_minor_locator(MultipleLocator(0.02)) ax.yaxis.set_major_locator(MultipleLocator(0.1)) else: ax.yaxis.set_minor_locator(MultipleLocator(0.2)) ax.yaxis.set_major_locator(MultipleLocator(1.)) ax.tick_params(labelsize=fontsize_ticks) if xlabels[k] is not None: ax.set_xlabel(xlabels[k], fontsize=fontsize_labels) else: ax.set_xticklabels([]) if ylabels[k] is not None: if (k == 1): labelpad = 0. else: labelpad=None ax.set_ylabel(ylabels[k], fontsize=fontsize_labels, labelpad=labelpad) else: ax.set_yticklabels([]) if k == 1: handles, labels_leg = ax.get_legend_handles_labels() neworder = range(plot_cnt_q) handles_arr = [] labels_arr = [] for ii in neworder: handles_arr.append(handles[ii]) labels_arr.append(labels_leg[ii]) neworder2 = range(plot_cnt_q, len(handles)) handles_arr2 = [] labels_arr2 = [] for ii in neworder2: handles_arr2.append(handles[ii]) labels_arr2.append(labels_leg[ii]) frameon = True framealpha = 1. edgecolor = 'none' borderpad = 0.25 fontsize_leg_tmp = fontsize_leg labelspacing=0.15 handletextpad=0.25 loc = 'upper left' bbox_to_anchor = (1.,1.) legend1 = ax.legend(handles_arr, labels_arr, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, loc=loc, bbox_to_anchor=bbox_to_anchor, numpoints=1, scatterpoints=1, frameon=frameon, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp) ax.add_artist(legend1) if len(handles_arr2) > 0: loc2 = 'lower left' bbox_to_anchor2 = (1.,0.) legend2 = ax.legend(handles_arr2, labels_arr2, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, loc=loc2, bbox_to_anchor=bbox_to_anchor2, numpoints=1, scatterpoints=1, frameon=frameon, framealpha=framealpha, edgecolor=edgecolor, #handlelength=1.4, fontsize=fontsize_leg_tmp) ax.add_artist(legend2) ax.set_zorder(1) if fileout is not None: plt.savefig(fileout, bbox_inches='tight', dpi=600) plt.close() else: plt.show() # ---------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------- # Figure 8 def plot_toy_impl_fDM_calibration_z_evol(lmstar_arr=None, output_path=None, table_path=None, table_path_alt=None, fileout=None, n_disk=1., Reff_bulge=1., n_bulge=4., invq_bulge=1., del_fDM=True, save_dict_stack=True, overwrite_dict_stack=False, include_toy_curves=True): """ Plot "typical" change of difference between fDM calibrations fDM(vcirc,Reff) / fDM(Menc,Reff) with redshift for SF, roughly MS galaxies, for different stellar masses. Saves plot to PDF. Parameters ---------- output_path: str Path to directory where the output plot will be saved. table_path: str Path to directory containing the Sersic profile tables. z: float, optional Redshift (to determine NFW halo properties). Default: z=2. lmstar_arr: array_like, optional Array of the log of the total stellar mass [logMsun]. Default: [9.0, 9.5, 10., 10.5, 11.] logMsun. n_disk: float, optional Sersic index of disk component. Default: n_disk = 1. Reff_bulge: float, optional Sersic projected 2D half-light (assumed half-mass) radius of bulge component [kpc]. Default: 1kpc n_bulge: float, optional Sersic index of bulge component. Default: n_bulge = 4. invq_bulge: float, optional Flattening of bulge component. Default: invq_bulge = 1. (spherical) include_toy_curves: bool, optional Include curves showing assumed interpolated / extrapolated parameters curves used to construct the toy MS model. Default: True fileout: str, optional Override the default filename and explicitly choose the output filename (must include full path). """ if (output_path is None) & (fileout is None): raise ValueError("Must set 'output_path' if 'fileout' is not set !") if table_path is None: raise ValueError("Must set 'table_path' !") if (fileout is None): # Ensure trailing slash: if output_path[-1] != '/': output_path += '/' fileout = output_path+'plot_toy_fdm_calibration_z_evol' if del_fDM: fileout += '_del_fDM' fileout += '.pdf' if lmstar_arr is None: lm_step = 0.25 lmstar_arr = np.arange(9.0, 11.+lm_step, lm_step) # Ensure it's an np array, not a list: lmstar_arr = np.array(lmstar_arr) # ++++++++++++++++ # plot: zstep = 0.05 z_arr = np.arange(0., 3.+zstep, zstep) color_arr, ls_arr, labels, lw_arr = ([] for _ in range(4)) for lmstar in lmstar_arr: color_arr.append(cmap_mass( (lmstar-lmstar_arr.min())/(lmstar_arr.max()-lmstar_arr.min()) ) ) if ((lmstar % 1.) == 0): lw_arr.append(2.) ls_arr.append('-') else: lw_arr.append(1.) ls_arr.append('-') if ((lmstar % 1.) == 0): labels.append(r'${:2.0f}$'.format(lmstar)) elif ((lmstar*2. % 1.) == 0): labels.append(r'${:2.1f}$'.format(lmstar)) else: labels.append(r'${:2.2f}$'.format(lmstar)) xlabel = r'$z$' if del_fDM: ylabel = r'$f_{\mathrm{DM}}^{v}(R_{\mathrm{e,disk}})-f_{\mathrm{DM}}^{m}(R_{\mathrm{e,disk}})$' ylim_comp = [-0.04, 0.005] else: ylabel = r'$f_{\mathrm{DM}}^{v}(R_{\mathrm{e,disk}})/f_{\mathrm{DM}}^{m}(R_{\mathrm{e,disk}})$' ylim_comp = [0.9, 1.01] ylabels = [r'$f_{\mathrm{DM}}^{v}(R_{\mathrm{e,disk}})$', ylabel, r'$r_{\mathrm{1/2,3D,baryons}}/R_{\mathrm{e,disk}}$', r'$f_{\mathrm{DM}}^v(R_{\mathrm{e,disk}}) - f_{\mathrm{DM}}^m(r_{\mathrm{1/2,3D,baryons}})$'] xlim = [0., 3.] ylim = [0.9, 1.01] ylims = [[0.,1.], ylim_comp, [0.5,1.1], [-0.05,0.2]] types = ['fdm_vsq', 'fDM_comp', 'rhalf3D_bary_to_Reff_disk', 'fDM_comp_rhalf3D_bary'] ann_arr_n = [2, 1, 0, 0] ann_arr = [ [r'$\displaystyle f_{\mathrm{DM}}^v(R_{\mathrm{e,disk}})$', r'$\displaystyle = \frac{v_{\mathrm{circ,DM}}^2(R_{\mathrm{e,disk}})}{v_{\mathrm{circ,tot}}^2(R_{\mathrm{e,disk}})}$'], r'$\displaystyle f_{\mathrm{DM}}^m(R_{\mathrm{e,disk}}) = \frac{M_{\mathrm{DM,sph}}(<r=R_{\mathrm{e,disk}})}{M_{\mathrm{tot,sph}}(<r=R_{\mathrm{e,disk}})}$', None, None] ann_arr_pos = [['upperleft', 'upperleft'],'lowerright', None, None] if include_toy_curves: keys_toy = ['Reff_disk', 'invq_disk', 'fgas', 'bt', 'lMhalo', 'halo_conc'] ylabels_toy = [] ylims_toy = [] ylocators_toy = [] for keyt in keys_toy: if keyt == 'Reff_disk': ylabels_toy.append(r'$R_{\mathrm{e,disk}}$ [kpc]') ylims_toy.append([0., 12.]) ylocators_toy.append([1.,5.]) elif keyt == 'invq_disk': ylabels_toy.append(r'$1/q_{0,\mathrm{disk}}$') ylims_toy.append([0., 12.]) ylocators_toy.append([1.,5.]) elif keyt == 'fgas': ylabels_toy.append(r'$f_{\mathrm{gas}}$') ylims_toy.append([0., 1.]) ylocators_toy.append([0.05, 0.5]) elif keyt == 'lMbar': ylabels_toy.append(r'$\log_{10}(M_{\mathrm{bar}}/M_{\odot})$') ylims_toy.append([8.8, 11.6]) ylocators_toy.append([0.2,1.]) elif keyt == 'bt': ylabels_toy.append(r'$B/T$') ylims_toy.append([0., 1.]) ylocators_toy.append([0.05, 0.5]) elif keyt == 'lMhalo': ylabels_toy.append(r'$\log_{10}(M_{\mathrm{halo}}/M_{\odot})$') ylims_toy.append([10.8, 13.2]) ylocators_toy.append([0.2,1.]) elif keyt == 'halo_conc': ylabels_toy.append(r'$c_{\mathrm{halo}}$') ylims_toy.append([2., 12.]) ylocators_toy.append([1.,5.]) else: raise ValueError ###################################### # Setup plot: f = plt.figure() scale = 2.5 n_cols = 2 if include_toy_curves: n_rows = 3 if len(keys_toy) == 7: fac = 1.5 elif len(keys_toy) == 6: fac = 1.75 else: n_rows = 1 fac = 1.15 facy = 1.305 pad = 0 f.set_size_inches(fac*scale*(n_cols+(n_cols-1)*pad),facy*scale*(n_rows+(n_rows-1)*pad)) if include_toy_curves: wspace = 0.125 hspace = wspace if len(keys_toy) == 7: height_ratios=[0.4, 2.] elif len(keys_toy) == 6: height_ratios= [2./7., 2.3] else: height_ratios=[len(types)/len(keys_toy), 2.] gs_outer = gridspec.GridSpec(2, 1, wspace=wspace, hspace=hspace, height_ratios=height_ratios) wspace = 0.45 hspace = wspace gs0 = gridspec.GridSpecFromSubplotSpec(1, len(keys_toy),subplot_spec=gs_outer[0,0], wspace=wspace, hspace=hspace ) axes_toy = [] for i in range(1): for j in range(len(keys_toy)): axes_toy.append(plt.subplot(gs0[i,j])) wspace = 0.3 hspace = 0.175 gs = gridspec.GridSpecFromSubplotSpec(n_rows-1, n_cols,subplot_spec=gs_outer[1,0], wspace=wspace, hspace=hspace ) axes = [] for i in range(n_rows-1): for j in range(n_cols): axes.append(plt.subplot(gs[i,j])) else: wspace = 0.25 hspace = wspace gs = gridspec.GridSpec(n_rows, n_cols, wspace=wspace, hspace=hspace) axes = [] for i in range(n_rows): for j in range(n_cols): axes.append(plt.subplot(gs[i,j])) ###################### # Load bulge: n, invq don't change. tab_bulge = io.read_profile_table(n=n_bulge, invq=invq_bulge, path=table_path) tab_bulge_menc = tab_bulge['menc3D_sph'] tab_bulge_vcirc = tab_bulge['vcirc'] tab_bulge_rad = tab_bulge['R'] tab_bulge_Reff = tab_bulge['Reff'] tab_bulge_mass = tab_bulge['total_mass'] # Clean up values inside rmin: Add the value at r=0: menc=0 if tab_bulge['R'][0] > 0.: tab_bulge_rad = np.append(0., tab_bulge_rad) tab_bulge_menc = np.append(0., tab_bulge_menc) tab_bulge_vcirc = np.append(0., tab_bulge_vcirc) m_interp_bulge = scp_interp.interp1d(tab_bulge_rad, tab_bulge_menc, fill_value=np.NaN, bounds_error=False, kind='cubic') v_interp_bulge = scp_interp.interp1d(tab_bulge_rad, tab_bulge_vcirc, fill_value=np.NaN, bounds_error=False, kind='cubic') rmin = 0. rmax=15. rstep=0.1 r_arr = np.arange(rmin, rmax+rstep, rstep) ###################### dict_stack = None f_save = output_path+'toy_impl_fDM_calibration_z' if del_fDM: f_save += '_delfDM' f_save += '.pickle' if save_dict_stack: if (os.path.isfile(f_save)): with open(f_save, 'rb') as f: dict_stack = copy.deepcopy(pickle.load(f)) if overwrite_dict_stack: print("Deleting old file") os.remove(f_save) dict_stack = None if dict_stack is None: dict_stack = [] for mm, lmstar in enumerate(lmstar_arr): print("lmstar={}".format(lmstar)) val_dict = {'z_arr': z_arr, 'lmstar': np.ones(len(z_arr)) * lmstar, 'bt': np.ones(len(z_arr)) * -99., 'Reff_disk': np.ones(len(z_arr)) * -99., 'invq_disk': np.ones(len(z_arr)) * -99., 'invq_near': np.ones(len(z_arr)) * -99., 'fgas': np.ones(len(z_arr)) * -99., 'lMbar': np.ones(len(z_arr)) * -99., 'lMhalo': np.ones(len(z_arr)) * -99., 'Rvir': np.ones(len(z_arr)) * -99., 'halo_conc': np.ones(len(z_arr)) * -99., 'menc_disk': np.ones(len(z_arr)) * -99., 'menc_bulge': np.ones(len(z_arr)) * -99., 'menc_halo': np.ones(len(z_arr)) * -99., 'menc_bar': np.ones(len(z_arr)) * -99., 'menc_tot': np.ones(len(z_arr)) * -99., 'fdm_menc': np.ones(len(z_arr)) * -99., 'vcirc_disk': np.ones(len(z_arr)) * -99., 'vcirc_bulge': np.ones(len(z_arr)) * -99., 'vcirc_halo': np.ones(len(z_arr)) * -99., 'vcirc_bar': np.ones(len(z_arr)) * -99., 'vcirc_tot': np.ones(len(z_arr)) * -99., 'fdm_vsq': np.ones(len(z_arr)) * -99., 'fDM_comp': np.ones(len(z_arr)) * -99., 'rhalf3D_disk': np.ones(len(z_arr)) * -99., 'rhalf3D_bulge': np.ones(len(z_arr)) * -99., 'rhalf3D_bary': np.ones(len(z_arr)) * -99., 'rhalf3D_bary_to_Reff_disk': np.ones(len(z_arr)) * -99., 'fdm_menc_rhalf3D_bary': np.ones(len(z_arr)) * -99., 'fDM_comp_rhalf3D_bary': np.ones(len(z_arr)) * -99. } for ll, z in enumerate(z_arr): print(" z={}".format(z)) Reff_disk = scaling_rel._mstar_Reff_relation(z=z, lmstar=lmstar, galtype='sf') invq_disk = scaling_rel._invq_disk_lmstar_estimate(z=z, lmstar=lmstar) fgas = scaling_rel._fgas_scaling_relation_MS(z=z, lmstar=lmstar) Mstar = np.power(10., lmstar) Mbaryon = Mstar / (1.-fgas) lMbar = np.log10(Mbaryon) bt = scaling_rel._bt_lmstar_relation(z=z, lmstar=lmstar, galtype='sf') Mhalo = scaling_rel._smhm_relation(z=z, lmstar=lmstar) lMhalo = np.log10(Mhalo) halo_conc = scaling_rel._halo_conc_relation(z=z, lmhalo=lMhalo) nfw = plot_calcs.NFW(z=z, Mvir=Mhalo, conc=halo_conc) Rvir = nfw.rvir val_dict['bt'][ll] = bt val_dict['Reff_disk'][ll] = Reff_disk val_dict['invq_disk'][ll] = invq_disk val_dict['fgas'][ll] = fgas val_dict['lMbar'][ll] = lMbar val_dict['lMhalo'][ll] = lMhalo val_dict['Rvir'][ll] = Rvir val_dict['halo_conc'][ll] = halo_conc ###### try: menc_disk = interp_profiles.interpolate_sersic_profile_menc(R=Reff_disk, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path) vcirc_disk = interp_profiles.interpolate_sersic_profile_VC(R=Reff_disk, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path) except: filename_base = 'deproj_sersic_model' fname_table = table_path_alt+filename_base+'_n{:0.1f}_invq{:0.2f}.fits'.format(n_disk, invq_disk) try: menc_disk = interp_profiles.interpolate_sersic_profile_menc(R=Reff_disk, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path_alt, filename=fname_table) vcirc_disk = interp_profiles.interpolate_sersic_profile_VC(R=Reff_disk, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path_alt, filename=fname_table) except: # Write table table_generation.calculate_sersic_profile_table(n=n_disk, invq=invq_disk, output_path=table_path_alt, fileout=fname_table, total_mass=5.e10, Reff=1., overwrite=False) menc_disk = interp_profiles.interpolate_sersic_profile_menc(R=Reff_disk, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path_alt, filename=fname_table) vcirc_disk = interp_profiles.interpolate_sersic_profile_VC(R=Reff_disk, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path_alt, filename=fname_table) menc_bulge = (m_interp_bulge(Reff_disk / Reff_bulge * tab_bulge_Reff) * ((bt*Mbaryon) / tab_bulge_mass) ) vcirc_bulge = (v_interp_bulge(Reff_disk / Reff_bulge * tab_bulge_Reff) * np.sqrt((bt*Mbaryon) / tab_bulge_mass) * np.sqrt(tab_bulge_Reff / Reff_bulge)) nfw = plot_calcs.NFW(z=z, Mvir=Mhalo, conc=halo_conc) menc_halo = nfw.enclosed_mass(Reff_disk) menc_baryons = menc_disk + menc_bulge menc_tot = menc_baryons + menc_halo fdm_menc = menc_halo/menc_tot vcirc_halo = nfw.v_circ(Reff_disk) vcirc_baryons = np.sqrt(vcirc_disk**2 + vcirc_bulge**2) vcirc_tot = np.sqrt(vcirc_baryons**2 + vcirc_halo**2) fdm_vsq = vcirc_halo**2/vcirc_tot**2 if del_fDM: val_dict['fDM_comp'][ll] = (fdm_vsq-fdm_menc) else: val_dict['fDM_comp'][ll] = fdm_vsq/fdm_menc #### val_dict['menc_disk'][ll] = menc_disk val_dict['menc_bulge'][ll] = menc_bulge val_dict['menc_halo'][ll] = menc_halo val_dict['menc_bar'][ll] = menc_baryons val_dict['menc_tot'][ll] = menc_tot val_dict['fdm_menc'][ll] = fdm_menc val_dict['vcirc_disk'][ll] = vcirc_disk val_dict['vcirc_bulge'][ll] = vcirc_bulge val_dict['vcirc_halo'][ll] = vcirc_halo val_dict['vcirc_bar'][ll] = vcirc_baryons val_dict['vcirc_tot'][ll] = vcirc_tot val_dict['fdm_vsq'][ll] = fdm_vsq val_dict['invq_near'][ll] = invq_disk nearest_n, nearest_invq = n_disk, invq_disk try: menc_disk_full = interp_profiles.interpolate_sersic_profile_menc(R=r_arr, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path) except: menc_disk_full = interp_profiles.interpolate_sersic_profile_menc(R=r_arr, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path_alt, filename=fname_table) menc_bulge_full = (m_interp_bulge(r_arr / Reff_bulge * tab_bulge_Reff) * ((bt*Mbaryon) / tab_bulge_mass) ) rhalf3D_disk = util_calcs.find_rhalf3D_sphere(R=r_arr, menc3D_sph=menc_disk_full, total_mass=(1.-bt)*Mbaryon) rhalf3D_bulge = util_calcs.find_rhalf3D_sphere(R=r_arr, menc3D_sph=menc_bulge_full, total_mass=bt*Mbaryon) rhalf3D_bary = util_calcs.find_rhalf3D_sphere(R=r_arr, menc3D_sph=menc_disk_full+menc_bulge_full, total_mass=Mbaryon) try: mD_r3Db = interp_profiles.interpolate_sersic_profile_menc(R=rhalf3D_bary, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path) except: mD_r3Db = interp_profiles.interpolate_sersic_profile_menc(R=rhalf3D_bary, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path_alt, filename=fname_table) mB_r3Db = (m_interp_bulge(rhalf3D_bary / Reff_bulge * tab_bulge_Reff) * ((bt*Mbaryon) / tab_bulge_mass) ) mH_r3Db = nfw.enclosed_mass(rhalf3D_bary) fdm_menc_rhalf3D_bary = mH_r3Db/(mD_r3Db+mB_r3Db+mH_r3Db) val_dict['rhalf3D_disk'][ll] = rhalf3D_disk val_dict['rhalf3D_bulge'][ll] = rhalf3D_bulge val_dict['rhalf3D_bary'][ll] = rhalf3D_bary val_dict['rhalf3D_bary_to_Reff_disk'][ll] = rhalf3D_bary/Reff_disk val_dict['fdm_menc_rhalf3D_bary'][ll] = fdm_menc_rhalf3D_bary val_dict['fDM_comp_rhalf3D_bary'][ll] = (fdm_vsq-fdm_menc_rhalf3D_bary) #### dict_stack.append(val_dict) if save_dict_stack: if not (os.path.isfile(f_save)): with open(f_save, 'wb') as f: pickle.dump(dict_stack, f) ######################## dict_toy = None f_save_toy = output_path+'toy_impl_fDM_calibration_z_MW_M31' if del_fDM: f_save_toy += '_delfDM' f_save_toy += '.pickle' if save_dict_stack: if (os.path.isfile(f_save_toy)): with open(f_save_toy, 'rb') as f: dict_toy = copy.deepcopy(pickle.load(f)) if overwrite_dict_stack: print("Deleting old file") os.remove(f_save_toy) dict_toy = None if dict_toy is None: dict_toy = {} ztoystep = 0.25 z_arr_toy = np.arange(z_arr.min(), z_arr.max()+ztoystep, ztoystep) names = ['MW', 'M31'] ln0_arr = [-2.9, -3.4] n_evol='const' cmf_source='papovich15' for name, ln0 in zip(names, ln0_arr): print("Toy model: {}".format(name)) val_dict = {'z_arr': z_arr_toy, 'lmstar': np.ones(len(z_arr_toy)) * -99, 'bt': np.ones(len(z_arr_toy)) * -99., 'Reff_disk': np.ones(len(z_arr_toy)) * -99., 'invq_disk': np.ones(len(z_arr_toy)) * -99., 'invq_near': np.ones(len(z_arr_toy)) * -99., 'fgas': np.ones(len(z_arr_toy)) * -99., 'lMbar': np.ones(len(z_arr_toy)) * -99., 'lMhalo': np.ones(len(z_arr_toy)) * -99., 'Rvir': np.ones(len(z_arr_toy)) * -99., 'halo_conc': np.ones(len(z_arr_toy)) * -99., 'menc_disk': np.ones(len(z_arr_toy)) * -99., 'menc_bulge': np.ones(len(z_arr_toy)) * -99., 'menc_halo': np.ones(len(z_arr_toy)) * -99., 'menc_bar': np.ones(len(z_arr_toy)) * -99., 'menc_tot': np.ones(len(z_arr_toy)) * -99., 'fdm_menc': np.ones(len(z_arr_toy)) * -99., 'vcirc_disk': np.ones(len(z_arr_toy)) * -99., 'vcirc_bulge': np.ones(len(z_arr_toy)) * -99., 'vcirc_halo': np.ones(len(z_arr_toy)) * -99., 'vcirc_bar': np.ones(len(z_arr_toy)) * -99., 'vcirc_tot': np.ones(len(z_arr_toy)) * -99., 'fdm_vsq': np.ones(len(z_arr_toy)) * -99., 'fDM_comp': np.ones(len(z_arr_toy)) * -99., 'rhalf3D_disk': np.ones(len(z_arr_toy)) * -99., 'rhalf3D_bulge': np.ones(len(z_arr_toy)) * -99., 'rhalf3D_bary': np.ones(len(z_arr_toy)) * -99., 'rhalf3D_bary_to_Reff_disk': np.ones(len(z_arr_toy)) * -99., 'fdm_menc_rhalf3D_bary': np.ones(len(z_arr_toy)) * -99., 'fDM_comp_rhalf3D_bary': np.ones(len(z_arr_toy)) * -99., 'n_evol': n_evol, 'cmf_source': cmf_source, 'name': name, } val_dict['lmstar'] = scaling_rel._mass_progenitor_num_density(ln0, z_arr_toy, n_evol=n_evol, cmf_source=cmf_source) for ll, z in enumerate(z_arr_toy): print(" z={}".format(z)) lmstar = val_dict['lmstar'][ll] Reff_disk = scaling_rel._mstar_Reff_relation(z=z, lmstar=lmstar, galtype='sf') invq_disk = scaling_rel._invq_disk_lmstar_estimate(z=z, lmstar=lmstar) fgas = scaling_rel._fgas_scaling_relation_MS(z=z, lmstar=lmstar) Mstar = np.power(10., lmstar) Mbaryon = Mstar / (1.-fgas) lMbar = np.log10(Mbaryon) bt = scaling_rel._bt_lmstar_relation(z=z, lmstar=lmstar, galtype='sf') Mhalo = scaling_rel._smhm_relation(z=z, lmstar=lmstar) lMhalo = np.log10(Mhalo) halo_conc = scaling_rel._halo_conc_relation(z=z, lmhalo=lMhalo) nfw = plot_calcs.NFW(z=z, Mvir=Mhalo, conc=halo_conc) Rvir = nfw.rvir val_dict['bt'][ll] = bt val_dict['Reff_disk'][ll] = Reff_disk val_dict['invq_disk'][ll] = invq_disk val_dict['fgas'][ll] = fgas val_dict['lMbar'][ll] = lMbar val_dict['lMhalo'][ll] = lMhalo val_dict['Rvir'][ll] = Rvir val_dict['halo_conc'][ll] = halo_conc ###### try: menc_disk = interp_profiles.interpolate_sersic_profile_menc(R=Reff_disk, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path) vcirc_disk = interp_profiles.interpolate_sersic_profile_VC(R=Reff_disk, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path) except: filename_base = 'deproj_sersic_model' fname_table = table_path_alt+filename_base+'_n{:0.1f}_invq{:0.2f}.fits'.format(n_disk, invq_disk) try: menc_disk = interp_profiles.interpolate_sersic_profile_menc(R=Reff_disk, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path_alt, filename=fname_table) vcirc_disk = interp_profiles.interpolate_sersic_profile_VC(R=Reff_disk, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path_alt, filename=fname_table) except: # Write table table_generation.calculate_sersic_profile_table(n=n_disk, invq=invq_disk, output_path=table_path_alt, fileout=fname_table, total_mass=5.e10, Reff=1., overwrite=False) menc_disk = interp_profiles.interpolate_sersic_profile_menc(R=Reff_disk, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path_alt, filename=fname_table) vcirc_disk = interp_profiles.interpolate_sersic_profile_VC(R=Reff_disk, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path_alt, filename=fname_table) menc_bulge = (m_interp_bulge(Reff_disk / Reff_bulge * tab_bulge_Reff) * ((bt*Mbaryon) / tab_bulge_mass) ) vcirc_bulge = (v_interp_bulge(Reff_disk / Reff_bulge * tab_bulge_Reff) * np.sqrt((bt*Mbaryon) / tab_bulge_mass) * np.sqrt(tab_bulge_Reff / Reff_bulge)) nfw = plot_calcs.NFW(z=z, Mvir=Mhalo, conc=halo_conc) menc_halo = nfw.enclosed_mass(Reff_disk) vcirc_halo = nfw.v_circ(Reff_disk) menc_baryons = menc_disk + menc_bulge menc_tot = menc_baryons + menc_halo fdm_menc = menc_halo/menc_tot vcirc_baryons = np.sqrt(vcirc_disk**2 + vcirc_bulge**2) vcirc_tot = np.sqrt(vcirc_baryons**2 + vcirc_halo**2) fdm_vsq = vcirc_halo**2/vcirc_tot**2 if del_fDM: val_dict['fDM_comp'][ll] = (fdm_vsq-fdm_menc) else: val_dict['fDM_comp'][ll] = fdm_vsq/fdm_menc #### val_dict['menc_disk'][ll] = menc_disk val_dict['menc_bulge'][ll] = menc_bulge val_dict['menc_halo'][ll] = menc_halo val_dict['menc_bar'][ll] = menc_baryons val_dict['menc_tot'][ll] = menc_tot val_dict['fdm_menc'][ll] = fdm_menc val_dict['vcirc_disk'][ll] = vcirc_disk val_dict['vcirc_bulge'][ll] = vcirc_bulge val_dict['vcirc_halo'][ll] = vcirc_halo val_dict['vcirc_bar'][ll] = vcirc_baryons val_dict['vcirc_tot'][ll] = vcirc_tot val_dict['fdm_vsq'][ll] = fdm_vsq val_dict['invq_near'][ll] = invq_disk nearest_n, nearest_invq = n_disk, invq_disk try: menc_disk_full = interp_profiles.interpolate_sersic_profile_menc(R=r_arr, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path) except: menc_disk_full = interp_profiles.interpolate_sersic_profile_menc(R=r_arr, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path_alt) menc_bulge_full = (m_interp_bulge(r_arr / Reff_bulge * tab_bulge_Reff) * ((bt*Mbaryon) / tab_bulge_mass) ) rhalf3D_disk = util_calcs.find_rhalf3D_sphere(R=r_arr, menc3D_sph=menc_disk_full, total_mass=(1.-bt)*Mbaryon) rhalf3D_bulge = util_calcs.find_rhalf3D_sphere(R=r_arr, menc3D_sph=menc_bulge_full, total_mass=bt*Mbaryon) rhalf3D_bary = util_calcs.find_rhalf3D_sphere(R=r_arr, menc3D_sph=menc_disk_full+menc_bulge_full, total_mass=Mbaryon) try: mD_r3Db = interp_profiles.interpolate_sersic_profile_menc(R=rhalf3D_bary, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path) except: mD_r3Db = interp_profiles.interpolate_sersic_profile_menc(R=rhalf3D_bary, total_mass=((1.-bt)*Mbaryon), Reff=Reff_disk, n=n_disk, invq=invq_disk, path=table_path_alt) mB_r3Db = (m_interp_bulge(rhalf3D_bary / Reff_bulge * tab_bulge_Reff) * ((bt*Mbaryon) / tab_bulge_mass) ) mH_r3Db = nfw.enclosed_mass(rhalf3D_bary) fdm_menc_rhalf3D_bary = mH_r3Db/(mD_r3Db+mB_r3Db+mH_r3Db) val_dict['rhalf3D_disk'][ll] = rhalf3D_disk val_dict['rhalf3D_bulge'][ll] = rhalf3D_bulge val_dict['rhalf3D_bary'][ll] = rhalf3D_bary val_dict['rhalf3D_bary_to_Reff_disk'][ll] = rhalf3D_bary/Reff_disk val_dict['fdm_menc_rhalf3D_bary'][ll] = fdm_menc_rhalf3D_bary val_dict['fDM_comp_rhalf3D_bary'][ll] = (fdm_vsq-fdm_menc_rhalf3D_bary) #### dict_toy[name] = val_dict if save_dict_stack: if not (os.path.isfile(f_save_toy)): with open(f_save_toy, 'wb') as f: pickle.dump(dict_toy, f) ###################### # FOR TOY PLOTS if include_toy_curves: for i in range(1): for j in range(len(keys_toy)): k = i*n_cols + j ax = axes_toy[k] keyy = keys_toy[j] ylim = ylims_toy[j] ylabel = ylabels_toy[j] plot_cnt_lmstar = 0 for mm, lmstar in enumerate(lmstar_arr): ax.plot(dict_stack[mm]['z_arr'], dict_stack[mm][keyy], ls=ls_arr[mm], color=color_arr[mm], lw=lw_arr[mm], label=labels[mm], zorder=-1.) plot_cnt_lmstar += 1 for name, color, m, s,zord in zip(['MW', 'M31'], ['black', 'grey'], ['*', 's'], [15, 7.5], [-0.5, -0.6]): whshow = np.where(((dict_toy[name]['z_arr'])%0.25 == 0))[0] ax.scatter(dict_toy[name]['z_arr'][whshow], dict_toy[name][keyy][whshow], color=color,facecolor='none', lw=0.5, marker=m, s=s, label=name,zorder=zord) ###################### if ylim is None: ylim = ax.get_ylim() ax.set_xlim(xlim) ax.set_ylim(ylim) ax.xaxis.set_minor_locator(MultipleLocator(0.2)) ax.xaxis.set_major_locator(MultipleLocator(1.)) ax.yaxis.set_minor_locator(MultipleLocator(ylocators_toy[j][0])) ax.yaxis.set_major_locator(MultipleLocator(ylocators_toy[j][1])) if xlabel is not None: ax.set_xlabel(xlabel, fontsize=fontsize_labels_sm-2) else: ax.set_xticklabels([]) if ylabel is not None: ax.set_ylabel(ylabel, fontsize=fontsize_labels_sm-2, labelpad=1) else: ax.set_yticklabels([]) ax.tick_params(labelsize=fontsize_ticks_sm-3) ###################### # FOR JUST THE fDM PLOTS! n_rows = 2 for i in range(n_rows): for j in range(n_cols): k = i*n_cols + j ax = axes[k] keyy = types[k] ylim = ylims[k] ylabel = ylabels[k] plot_cnt_lmstar = 0 for mm, lmstar in enumerate(lmstar_arr): ax.plot(dict_stack[mm]['z_arr'], dict_stack[mm][keyy], ls=ls_arr[mm], color=color_arr[mm], lw=lw_arr[mm], label=labels[mm], zorder=-1.) plot_cnt_lmstar += 1 ##################### for name, color, m, s,zord in zip(['MW', 'M31'], ['black', 'grey'], ['*', 's'], [50, 25], [-0.5, -0.6]): whshow = np.where(((dict_toy[name]['z_arr'])%0.25 == 0))[0] ax.scatter(dict_toy[name]['z_arr'][whshow], dict_toy[name][keyy][whshow], color=color, facecolor='none', marker=m, s=s, lw=1., label=name,zorder=zord) ###################### if ylim is None: ylim = ax.get_ylim() if keyy == 'fDM_comp': if del_fDM: ax.axhline(y=0., ls=(0, (5,3)), color='darkgrey', zorder=-20.) else: ax.axhline(y=1., ls=(0, (5,3)), color='darkgrey', zorder=-20.) elif keyy == 'fDM_comp_rhalf3D_bary': ax.axhline(y=0., ls=(0, (5,3)), color='darkgrey', zorder=-20.) elif keyy == 'rhalf3D_bary_to_Reff_disk': ax.axhline(y=1., ls=(0, (5,3)), color='darkgrey', zorder=-20.) if ann_arr_n[k] > 0: xydelt = 0.04 if (ann_arr_n[k] == 1): if ann_arr_pos[k] == 'lowerright': xy = (1.-xydelt, xydelt) va='bottom' ha='right' elif ann_arr_pos[k] == 'upperright': xy = (1.-xydelt, 1.-xydelt) va='top' ha='right' elif ann_arr_pos[k] == 'upperleft': xy = (xydelt, 1.-xydelt) va='top' ha='left' ax.annotate(ann_arr[k], xy=xy, va=va, ha=ha, fontsize=fontsize_ann-0.5, xycoords='axes fraction') else: for mm in range(ann_arr_n[k]): if ann_arr_pos[k][mm] == 'lowerright': xy = (1.-xydelt, xydelt) va='bottom' ha='right' elif ann_arr_pos[k][mm] == 'upperright': xy = (1.-xydelt, 1.-xydelt) va='top' ha='right' elif ann_arr_pos[k][mm] == 'upperleft': if (mm > 0): xdelt = xydelt * 2. ydelt = xydelt + 0.065 else: xdelt = ydelt = xydelt xy = (xdelt, 1.-ydelt) va='top' ha='left' ax.annotate(ann_arr[k][mm], xy=xy, va=va, ha=ha, fontsize=fontsize_ann-0.5, xycoords='axes fraction') ax.set_xlim(xlim) ax.set_ylim(ylim) ax.xaxis.set_minor_locator(MultipleLocator(0.2)) ax.xaxis.set_major_locator(MultipleLocator(1.)) if keyy == 'fDM_comp': if del_fDM: ax.yaxis.set_minor_locator(MultipleLocator(0.005)) ax.yaxis.set_major_locator(MultipleLocator(0.02)) else: ax.yaxis.set_minor_locator(MultipleLocator(0.01)) ax.yaxis.set_major_locator(MultipleLocator(0.05)) elif keyy == 'fDM_comp_rhalf3D_bary': ax.yaxis.set_minor_locator(MultipleLocator(0.01)) ax.yaxis.set_major_locator(MultipleLocator(0.05)) elif keyy == 'rhalf3D_bary_to_Reff_disk': ax.yaxis.set_minor_locator(MultipleLocator(0.02)) ax.yaxis.set_major_locator(MultipleLocator(0.1)) else: ax.yaxis.set_minor_locator(MultipleLocator(0.05)) ax.yaxis.set_major_locator(MultipleLocator(0.2)) if xlabel is not None: ax.set_xlabel(xlabel, fontsize=fontsize_labels) else: ax.set_xticklabels([]) if ylabel is not None: ax.set_ylabel(ylabel, fontsize=fontsize_labels) else: ax.set_yticklabels([]) ax.tick_params(labelsize=fontsize_ticks) if k == 0: handles, labels_leg = ax.get_legend_handles_labels() neworder = range(plot_cnt_lmstar) handles_arr = [] labels_arr = [] for ii in neworder: handles_arr.append(handles[ii]) labels_arr.append(labels_leg[ii]) neworder2 = range(plot_cnt_lmstar, len(handles)) handles_arr2 = [] labels_arr2 = [] for ii in neworder2: handles_arr2.append(handles[ii]) labels_arr2.append(labels_leg[ii]) frameon = True framealpha = 1. edgecolor = 'none' borderpad = 0.25 fontsize_leg_tmp = fontsize_leg labelspacing=0.15 handletextpad=0.25 loc = 'upper right' leg_title = r'$\log_{10}(M_{\star}/M_{\odot})=$' fontsize_leg_title = fontsize_leg legend1 = ax.legend(handles_arr, labels_arr, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, loc=loc, numpoints=1, scatterpoints=1, frameon=frameon, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp, title=leg_title, title_fontsize=fontsize_leg_title) ax.add_artist(legend1) if len(handles_arr2) > 0: legend2 = ax.legend(handles_arr2, labels_arr2, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, loc='lower left', numpoints=1, scatterpoints=1, frameon=frameon, framealpha=framealpha, edgecolor=edgecolor, fontsize=fontsize_leg_tmp) ax.add_artist(legend2) if fileout is not None: plt.savefig(fileout, bbox_inches='tight', dpi=600) plt.close() else: plt.show() return None # ---------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------- # Figure 9 def plot_alpha_vs_R(fileout=None, output_path=None, table_path=None, n_arr=[0.5, 1., 2., 4., 8.], show_literature=True, print_values=False): """ Plot alpha=-dlnrho_g/dlnr derived for deprojected Sersic distribution, over a range of Sersic index n. Compare to self-gravitating exponential disk case (eg Burkert+10). Saves plot to PDF. Parameters ---------- output_path: str Path to directory where the output plot will be saved. table_path: str Path to directory containing the Sersic profile tables. n_arr: array_like, optional Range of Sersic indices to plot. Default: n_arr = [0.5, 1., 2., 4,, 8.] show_literature: bool, optional Whether to show inferred alpha values from some simulation literature work: Kretschmer et al., 2021, MNRAS, 503, 5238-5253; Dalcanton & Stilp, 2010, ApJ, 721, 547 fileout: str, optional Override the default filename and explicitly choose the output filename (must include full path). """ if (output_path is None) & (fileout is None): raise ValueError("Must set 'output_path' if 'fileout' is not set !") if table_path is None: raise ValueError("Must set 'table_path' !") if (fileout is None): # Ensure trailing slash: if output_path[-1] != '/': output_path += '/' fileout = output_path+'alpha_vs_r_sersic_potential' fileout += '.pdf' #################### n_arr = np.array(n_arr) color_arr, labels = ([] for _ in range(2)) nextra = -0.5 nrange = 7. for n in n_arr: if n == 8.: color_arr.append('#2d0037') elif n < 0.5: color_arr.append('orange') else: color_arr.append(cmap_n((n+nextra)/(nrange+nextra))) if n < 1.: labels.append(r'$n={:0.1f}$'.format(n)) else: labels.append(r'$n={:0.0f}$'.format(n)) q_arr = [1.] # Load files: ks_dict = {} ks_dict['alpha'] = {} invq = 1. q = 1. for j, n in enumerate(n_arr): tab = io.read_profile_table(n=n, invq=invq, path=table_path) if j == 0: ks_dict['Reff'] = tab['Reff'] ks_dict['alpha']['narr'] = n_arr ks_dict['R'] = tab['R'] ks_dict['alpha']['n={}'.format(n)] = -tab['dlnrho_dlnR'] types = ['alpha', 'alpha_by_sg'] xlabel = r'$R/R_{\mathrm{e}}$' ylabels = [r'$\alpha(R)$', r'$\alpha/\alpha_{\mathrm{self-grav}}(R)$'] xlim = [0., 5.] ylims = [[0., 10.], [0., 2.]] titles = [None, None] ann_arr = [ r'$\displaystyle \alpha(R) = -\frac{\mathrm{d}\ln \rho_g}{\mathrm{d}\ln{}R}$', None] ann_arr_pos = ['lowerright', 'lowerright'] lw = 1.3 ls_arr = ['-'] if show_literature: lit_dict = {} del_RtoRe = 0.01 ##### # Eq A.5, Bouche+21 // from Dalcanton & Stilp 2010, Eq 17 (exponential) RtoRe = np.arange(0.,5.+del_RtoRe,del_RtoRe) lit_dict['DS10_ISM_P'] = {'RtoRe': RtoRe, 'alpha': scaling_rel._dalcanton_stilp10_alpha_n(RtoRe, 1.), 'ls': '--', 'lw': 1.25, 'color': 'orange', 'label': r'Dalcanton \& Stilp 2010, Eq. 17, exp.'} ##### # Kretschmer+21, Figure 4 alpha_rho (medians from data from M. Kretschmer) k21_fig4 = scaling_rel._kretschmer21_fig4_alpharho(path=output_path) if print_values: nser_50 = k21_fig4['n_sersic'][0] nser_16 = nser_50 - k21_fig4['n_sersic_err_l68'][0] nser_84 = k21_fig4['n_sersic_err_u68'][0] - nser_50 print("K21: n_sersic: median={:0.2f}, 16,84% range=[{:0.2f},{:0.2f}]".format(nser_50, nser_16, nser_84)) lit_dict['K21_fig4_alpharho'] = { 'RtoRe': k21_fig4['RtoRe'].value, 'alpha': k21_fig4['alpha_rho'].value, 'alpha_err_l68': k21_fig4['alpha_rho_err_l68'].value, 'alpha_err_u68': k21_fig4['alpha_rho_err_u68'].value, 'marker': 'o', 'ms': 25., 'color': 'red', 'label': r'Kretschmer et al. 2021, Fig. 4, $\alpha_{\rho}$'} ###################################### # Setup plot: f = plt.figure() scale = 3.85 n_cols = len(types) f.set_size_inches(1.15*scale*n_cols,scale) pad_outer = 0.2 gs = gridspec.GridSpec(1, n_cols, wspace=pad_outer) axes = [] for i in range(n_cols): axes.append(plt.subplot(gs[0,i])) r_arr = ks_dict['R']/ks_dict['Reff'] for i in range(n_cols): ax = axes[i] ylim = ylims[i] if types[i] == 'alpha': if (i==0): lblSG = 'Self-grav' else: lblSG = None ax.plot(r_arr, 3.36*r_arr, ls='--', color='black', lw=lw, label=lblSG) ls = ls_arr[0] for k, n in enumerate(n_arr): if types[i] == 'alpha': alpha_arr = ks_dict['alpha']['n={}'.format(n)] elif types[i] == 'alpha_by_sg': alpha_arr = ks_dict['alpha']['n={}'.format(n)] / (3.36*r_arr) if len(alpha_arr) != len(r_arr): raise ValueError if (i==0): lbln = labels[k] else: lbln = None ax.plot(r_arr, alpha_arr, ls=ls, color=color_arr[k], lw=lw, label=lbln) if show_literature & (k==len(n_arr)-1): zorder_lit = -10. for key in lit_dict.keys(): if (i==1): lbl = lit_dict[key]['label'] else: lbl = None if 'ls' in lit_dict[key].keys(): if 'frac_uncertainty' in lit_dict[key].keys(): if types[i] == 'alpha': ax.fill_between(lit_dict[key]['RtoRe'], lit_dict[key]['alpha']*(1-lit_dict[key]['frac_uncertainty']), lit_dict[key]['alpha']*(1+lit_dict[key]['frac_uncertainty']), ls=lit_dict[key]['ls'], color=lit_dict[key]['color'], lw=lit_dict[key].get('lw', 1.), alpha=0.2, label=None, zorder=zorder_lit-10) elif types[i] == 'alpha_by_sg': ax.fill_between(lit_dict[key]['RtoRe'], lit_dict[key]['alpha']/(3.36*lit_dict[key]['RtoRe'])*(1-lit_dict[key]['frac_uncertainty']), lit_dict[key]['alpha']/(3.36*lit_dict[key]['RtoRe'])*(1+lit_dict[key]['frac_uncertainty']), ls=lit_dict[key]['ls'], color=lit_dict[key]['color'], lw=lit_dict[key].get('lw', 1.), alpha=0.2, label=None, zorder=zorder_lit-10) if types[i] == 'alpha': ax.plot(lit_dict[key]['RtoRe'], lit_dict[key]['alpha'], ls=lit_dict[key]['ls'], color=lit_dict[key]['color'], lw=lit_dict[key].get('lw', 1.), label=lbl, zorder=zorder_lit) elif types[i] == 'alpha_by_sg': ax.plot(lit_dict[key]['RtoRe'], lit_dict[key]['alpha']/(3.36*lit_dict[key]['RtoRe']), ls=lit_dict[key]['ls'], color=lit_dict[key]['color'], lw=lit_dict[key].get('lw', 1.), label=lbl, zorder=zorder_lit) else: # ------- # Show uncertainties: if 'alpha_err_l68' in lit_dict[key].keys(): if types[i] == 'alpha': ax.errorbar(lit_dict[key]['RtoRe'], lit_dict[key]['alpha'], yerr = [lit_dict[key]['alpha_err_l68'], lit_dict[key]['alpha_err_u68']], marker=None, ls='None', capsize=0., lw=1, color='lightgrey', label=None, zorder=zorder_lit+1) elif types[i] == 'alpha_by_sg': ax.errorbar(lit_dict[key]['RtoRe'], lit_dict[key]['alpha']/(3.36*lit_dict[key]['RtoRe']), yerr = [lit_dict[key]['alpha_err_l68']/(3.36*lit_dict[key]['RtoRe']), lit_dict[key]['alpha_err_u68']/(3.36*lit_dict[key]['RtoRe'])], marker=None, ls='None', capsize=0., lw=1, color='lightgrey', label=None, zorder=zorder_lit+1) # ------- if types[i] == 'alpha': ax.scatter(lit_dict[key]['RtoRe'], lit_dict[key]['alpha'], marker=lit_dict[key]['marker'], color=lit_dict[key]['color'], s=lit_dict[key].get('ms', 30.), label=lbl, zorder=zorder_lit+1) elif types[i] == 'alpha_by_sg': ax.scatter(lit_dict[key]['RtoRe'], lit_dict[key]['alpha']/(3.36*lit_dict[key]['RtoRe']), marker=lit_dict[key]['marker'], color=lit_dict[key]['color'], s=lit_dict[key].get('ms', 30.), label=lbl, zorder=zorder_lit+1) if types[i] == 'alpha_by_sg': ax.axhline(y=1., ls='--', color='black', zorder=-20.) ax.axvline(x=1., ls=':', color='lightgrey', zorder=-20.) ax.set_xlim(xlim) ax.set_ylim(ylim) ax.set_xlabel(xlabel, fontsize=fontsize_labels) ax.set_ylabel(ylabels[i], fontsize=fontsize_labels) ax.tick_params(labelsize=fontsize_ticks) if titles[i] is not None: ax.set_title(titles[i], fontsize=fontsize_title) ax.xaxis.set_minor_locator(MultipleLocator(0.2)) ax.xaxis.set_major_locator(MultipleLocator(1.0)) xydelt = 0.04 if ann_arr_pos[i] == 'lowerright': xy = (1.-xydelt, xydelt) va='bottom' ha='right' elif ann_arr_pos[i] == 'upperright': xy = (1.-xydelt, 1.-xydelt) va='top' ha='right' ax.annotate(ann_arr[i], xy=xy, va=va, ha=ha, fontsize=fontsize_ann_latex_sm, xycoords='axes fraction') if (ylim[1]-ylim[0]) > 9: ax.yaxis.set_minor_locator(MultipleLocator(0.5)) ax.yaxis.set_major_locator(MultipleLocator(2.)) elif (ylim[1]-ylim[0]) > 3: ax.yaxis.set_minor_locator(MultipleLocator(0.2)) ax.yaxis.set_major_locator(MultipleLocator(1.)) elif (ylim[1]-ylim[0]) > 1: ax.yaxis.set_minor_locator(MultipleLocator(0.1)) ax.yaxis.set_major_locator(MultipleLocator(0.5)) else: ax.yaxis.set_minor_locator(MultipleLocator(0.02)) ax.yaxis.set_major_locator(MultipleLocator(0.1)) if i == 0: frameon = False framealpha = 1. borderpad = 0.75 fontsize_leg_tmp = fontsize_leg labelspacing= 0.2 handletextpad= 0.5 fancybox = False edgecolor='None' facecolor = 'white' loc = 'upper left' bbox_to_anchor = (0., 1.) legend = ax.legend(labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, frameon=frameon, loc=loc, bbox_to_anchor=bbox_to_anchor, numpoints=1, scatterpoints=1,fontsize=fontsize_leg_tmp, fancybox=fancybox,edgecolor=edgecolor, facecolor=facecolor, framealpha=framealpha) patch = mpatches.FancyBboxPatch((0.8,5.95), 0.4, 3.55, boxstyle="square,pad=0.", fc="white", lw=0.) ax.add_patch(patch) if show_literature & (i==1): frameon = False framealpha = 1. borderpad = 0.75 fontsize_leg_tmp = fontsize_leg labelspacing= 0.2 handletextpad= 0.5 fancybox = False edgecolor='None' facecolor = 'white' loc = 'upper left' bbox_to_anchor = (0.075, 1.) legend = ax.legend(labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad,frameon=frameon, loc=loc, bbox_to_anchor=bbox_to_anchor, numpoints=1, scatterpoints=1,fontsize=fontsize_leg_tmp, fancybox=fancybox,edgecolor=edgecolor, facecolor=facecolor, framealpha=framealpha) patch = mpatches.FancyBboxPatch((0.8,1.625), 4.0, 0.275, boxstyle="square,pad=0.", fc="white", lw=0., zorder=5.) ax.add_patch(patch) if fileout is not None: plt.savefig(fileout, bbox_inches='tight', dpi=600) plt.close() else: plt.show() return None # ---------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------- # Figure 10 def plot_AD_sersic_potential_alpha_vs_R(fileout=None, output_path=None, table_path=None, sigma0_arr = [30., 60., 90.], q_arr=[1., 0.2], n_arr=[0.5, 1., 2., 4.], total_mass = np.power(10.,10.5), Reff=1., show_sigmar_toy=False, show_sigmar_toy_nosig0=False): """ Plot pressure support using alpha=-dlnrho_g/dlnr derived for deprojected Sersic distributions, over a range of Sersic index n and intrinsic axis ratios q. Compare to self-gravitating exponential disk case (eg Burkert+10). Saves plot to PDF. Parameters ---------- output_path: str Path to directory where the output plot will be saved. table_path: str Path to directory containing the Sersic profile tables. sigma0_arr: array_like, optional Intrinsic constant velocity dispersions to plot (km/s). Default: [30., 60., 90.] q_arr: array_like, optional Range of intrinsic axis ratios to plot. Default: q_arr = [1., 0.2] n_arr: array_like, optional Range of Sersic indices to plot. Default: n_arr = [0.5, 1., 2., 4.] total_mass: float, optional Total mass of the Sersic component (Msun). Default: 10^10.5 Msun. Reff: float, optional Effective radius of the Sersic component (kpc). Default: 1 kpc fileout: str, optional Override the default filename and explicitly choose the output filename (must include full path). """ if (output_path is None) & (fileout is None): raise ValueError("Must set 'output_path' if 'fileout' is not set !") if table_path is None: raise ValueError("Must set 'table_path' !") if (fileout is None): # Ensure trailing slash: if output_path[-1] != '/': output_path += '/' fileout = output_path+'AD_sersic_potential_alpha_vs_r' if show_sigmar_toy: fileout += '_sigmar_toy' if show_sigmar_toy_nosig0: fileout += '_sigmar_toy_nosig0' fileout += '.pdf' #################### n_arr = np.array(n_arr) q_arr = np.array(q_arr) labels, color_arr = ([] for _ in range(2)) color_arr_basic = ['tab:purple', 'tab:cyan', 'tab:orange'] for i, sig0 in enumerate(sigma0_arr): color_arr.append(color_arr_basic[i]) labels.append(r'$\sigma_0='+'{:0.0f}'.format(sig0)+r'\ \mathrm{km\,s^{-1}}$') xlabel = r'$R/R_{\mathrm{e}}$' ylabel = r'$v(R)\ \mathrm{[km\,s^{-1}]}$' xlim = [0., 5.] ylim = [0.,340.] # Load files: delR = 0.001 Rarr = np.arange(xlim[0], xlim[1]+delR, delR)*Reff xlim = [0., 3.] ks_dict = {} for q in q_arr: invq = 1./q ks_dict_q = {} for j, n in enumerate(n_arr): ks_dict_q['n={}'.format(n)] = {} tab = io.read_profile_table(n=n, invq=invq, path=table_path) if (q == q_arr.min()) & (n == n_arr.min()): ks_dict['Reff'] = Reff ks_dict['narr'] = n_arr ks_dict['qarr'] = q_arr ks_dict['R'] = Rarr ks_dict['total_mass'] = total_mass ks_dict_q['n={}'.format(n)]['alpha'] = -1. * interp_profiles.interpolate_sersic_profile_dlnrho_dlnR_nearest(R=Rarr, Reff=Reff, n=n, invq=invq, path=table_path) ks_dict_q['n={}'.format(n)]['vcirc'] = interp_profiles.interpolate_sersic_profile_VC_nearest(R=Rarr, total_mass=total_mass, Reff=Reff, n=n, invq=invq, path=table_path) ks_dict['q={}'.format(q)] = ks_dict_q titles, ann_arr, ann_arr_pos = ([] for _ in range(3)) for n in n_arr: if n < 1.: nann = '{:0.1f}'.format(n) else: nann = '{:0.0f}'.format(n) titles.append(r'$n='+"{}$".format(nann)) for i, q in enumerate(q_arr): if q != 1.: qstr = "{:0.1f}".format(q) else: qstr = "{:0.0f}".format(q) ann_arr.append(r'$q_0='+'{}$'.format(qstr)) ann_arr_pos.append('upperright') lw = 1.3 ###################################### # Setup plot: f = plt.figure() scale = 2.9 n_cols = len(n_arr) n_rows = len(q_arr) f.set_size_inches(1.05*scale*n_cols,scale*n_rows) pad_outer = 0.05 gs = gridspec.GridSpec(n_rows, n_cols, wspace=pad_outer,hspace=pad_outer) axes = [] for j in range(n_rows): for i in range(n_cols): axes.append(plt.subplot(gs[j,i])) r_arr = ks_dict['R']/ks_dict['Reff'] for j, q in enumerate(q_arr): for i, n in enumerate(n_arr): mm = j*n_cols+i ax = axes[mm] for k, sig0 in enumerate(sigma0_arr): vcirc = ks_dict['q={}'.format(q)]['n={}'.format(n)]['vcirc'].copy() alpha = ks_dict['q={}'.format(q)]['n={}'.format(n)]['alpha'].copy() if (mm == 1): lbl_sig0 = labels[k] else: lbl_sig0 = None if (mm == 0): lbl_circ = r'$v_{\mathrm{circ}}$' lbl_alpha = r'$v_{\mathrm{rot}}$, $\alpha(n)$' lbl_SG = r'$v_{\mathrm{rot}}$, Self-grav' lbl_alpha_sigmar = r'$v_{\mathrm{rot}}$, $\alpha(n)+\alpha_{\sigma(r)}$' lbl_alpha_sigmar_nosig0 = r'$v_{\mathrm{rot}}$, $\alpha(n)+\alpha_{\sigma(r)}$, $\sigma_0=0$' else: lbl_circ = None lbl_alpha = None lbl_SG = None lbl_alpha_sigmar = None lbl_alpha_sigmar_nosig0 = None if k == 0: ax.plot(r_arr, vcirc, ls='-', color='black', lw=lw, label=lbl_circ, zorder=10.) if lbl_alpha is not None: ax.plot(r_arr, r_arr*np.NaN, ls='--', color='black', lw=lw, label=lbl_alpha) if show_sigmar_toy: ax.plot(r_arr, r_arr*np.NaN, ls='-.', color='black', lw=lw, label=lbl_alpha_sigmar) if show_sigmar_toy_nosig0: ax.plot(r_arr, r_arr*np.NaN, ls=(0, (5, 2, 1, 2, 1, 2, 1, 2)), color='black', lw=lw, label=lbl_alpha_sigmar_nosig0) ax.plot(r_arr, r_arr*np.NaN, ls=':', color='black', lw=lw, label=lbl_SG) if lbl_sig0 is not None: ax.plot(r_arr, r_arr*np.NaN, ls='-', color=color_arr[k], lw=lw, label=lbl_sig0) ax.plot(r_arr, np.sqrt(vcirc**2-alpha*(sig0**2)), ls='--', color=color_arr[k], lw=lw) if show_sigmar_toy: sigr = _sigr_toy(r_arr, 2.*sig0, sig0, 0.5*Reff) alphasigr = _alpha_sigr_toy(r_arr, 2.*sig0, sig0, 0.5*Reff) ax.plot(r_arr, np.sqrt(vcirc**2-(alpha+alphasigr)*(sigr**2)), ls='-.', color=color_arr[k], lw=lw) if show_sigmar_toy_nosig0: sigr_nosig0 = _sigr_toy(r_arr, np.sqrt(5.)*sig0, 0., 0.5*Reff) alphasigr_nosig0 = _alpha_sigr_toy(r_arr, np.sqrt(5.)*sig0, 0., 0.5*Reff) print(sigr[0], sigr_nosig0[0]) ax.plot(r_arr, np.sqrt(vcirc**2-(alpha+alphasigr_nosig0)*(sigr_nosig0**2)), ls=(0, (5, 2, 1, 2, 1, 2, 1, 2)), color=color_arr[k], lw=lw) ax.plot(r_arr, np.sqrt(vcirc**2-3.36*r_arr*(sig0**2)), ls=':', color=color_arr[k], lw=lw) ax.set_xlim(xlim) ax.set_ylim(ylim) ax.tick_params(labelsize=fontsize_ticks) if (titles[i] is not None) & (j==0): ax.set_title(titles[i], fontsize=fontsize_title) ax.xaxis.set_minor_locator(MultipleLocator(0.2)) ax.xaxis.set_major_locator(MultipleLocator(1.)) xydelt = 0.04 arpos = ann_arr_pos[j] if arpos == 'lowerright': xy = (1.-xydelt, xydelt) va='bottom' ha='right' elif arpos == 'upperright': xy = (1.-xydelt, 1.-xydelt) va='top' ha='right' elif arpos == 'lowerleft': xy = (xydelt, xydelt) va='bottom' ha='left' elif arpos == 'upperleft': xy = (xydelt, 1.-xydelt) va='top' ha='left' ax.annotate(ann_arr[j], xy=xy, va=va, ha=ha, fontsize=fontsize_ann_latex, xycoords='axes fraction') if (ylim[1]-ylim[0]) > 300: ax.yaxis.set_minor_locator(MultipleLocator(20.)) ax.yaxis.set_major_locator(MultipleLocator(100.)) elif (ylim[1]-ylim[0]) > 50: ax.yaxis.set_minor_locator(MultipleLocator(10.)) ax.yaxis.set_major_locator(MultipleLocator(50.)) elif (ylim[1]-ylim[0]) > 3: ax.yaxis.set_minor_locator(MultipleLocator(0.2)) ax.yaxis.set_major_locator(MultipleLocator(1.)) elif (ylim[1]-ylim[0]) > 1: ax.yaxis.set_minor_locator(MultipleLocator(0.1)) ax.yaxis.set_major_locator(MultipleLocator(0.5)) else: ax.yaxis.set_minor_locator(MultipleLocator(0.02)) ax.yaxis.set_major_locator(MultipleLocator(0.1)) if (j == n_rows-1): ax.set_xlabel(xlabel, fontsize=fontsize_labels) else: ax.set_xticklabels([]) if (i == 0): ax.set_ylabel(ylabel, fontsize=fontsize_labels) else: ax.set_yticklabels([]) if (mm <= 1): frameon = False borderpad = 0.1 fontsize_leg_tmp = fontsize_leg labelspacing=0.15 handletextpad= 0.5 fancybox = False edgecolor='None' loc='upper left' bbox_to_anchor = (0.,1.) legend = ax.legend(labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, loc=loc, frameon=frameon, bbox_to_anchor=bbox_to_anchor, numpoints=1, scatterpoints=1,fontsize=fontsize_leg_tmp, fancybox=fancybox,edgecolor=edgecolor) if fileout is not None: plt.savefig(fileout, bbox_inches='tight', dpi=600) plt.close() else: plt.show() return None # ---------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------- # Figure 11 def plot_composite_alpha_bulge_2disk_vs_r(fileout=None, output_path=None, table_path=None, nDisk=1., nBulge=4., bt_arr=[0.,0.25,0.5,0.75,1.], ReB_to_ReD_arr=[0.2, 0.5, 1.], nDisk2=1.2, D2t_arr=[0.,0.25,0.5,0.75,1.], ReD2_to_ReD_arr=[0.5, 1., 1.5], show_BT=True, show_D2T=True): """ Plot composite alpha for bulge+disk over a range of B/T ratios and Rebulge/Redisk, and a 2-disk system over a range of D2/T ratios and ReD2/ReD. Saves plot to PDF. Parameters ---------- output_path: str Path to directory where the output plot will be saved. table_path: str Path to directory containing the Sersic profile tables. bt_arr: array_like, optional Range of B/T ratios to plot. Default: bt_arr=[0.,0.25,0.5,0.75,1.] ReB_to_ReD_arr: array_like, optional Range of Re,bulge/Re,disk ratios to plot. Default: ReB_to_ReD_arr=[0.2, 0.5, 1.] nDisk: float, optional Sersic index of disk. Default: nDisk = 1. nBulge: float, optional Sersic index of bulge. Default: nBulge = 4. nDisk2: float, optional Sersic index of second disk. Defautl: nDisk2 = 1.2 D2t_arr: array_like, optional Range of D2/T ratios to plot. Default: 2t_arr=[0.,0.25,0.5,0.75,1.] ReD2_to_ReD_arr: array_like, optional Range of Re,bulge/Re,disk ratios to plot. Default: ReD2_to_ReD_arr=[0.5, 1., 1.5] fileout: str, optional Override the default filename and explicitly choose the output filename (must include full path). """ if (output_path is None) & (fileout is None): raise ValueError("Must set 'output_path' if 'fileout' is not set !") if table_path is None: raise ValueError("Must set 'table_path' !") if (fileout is None): # Ensure trailing slash: if output_path[-1] != '/': output_path += '/' if show_BT & show_D2T: fileout = output_path+'alpha_composite_vs_r_BT_RebulgeRedisk_D2T_ReD2ReD.pdf' elif (not show_BT) & show_D2T: fileout = output_path+'alpha_composite_vs_r_D2T_Redisk2Redisk.pdf' elif show_BT & (not show_D2T): fileout = output_path+'alpha_composite_vs_r_BT_RebulgeRedisk.pdf' #################### bt_arr = np.array(bt_arr) ReB_to_ReD_arr = np.array(ReB_to_ReD_arr) D2t_arr = np.array(D2t_arr) ReD2_to_ReD_arr = np.array(ReD2_to_ReD_arr) color_arr_bt, labels_bt = ([] for _ in range(2)) color_arr_D2t, labels_D2t = ([] for _ in range(2)) btextra = 0. btrange = 1. for bt in bt_arr: color_arr_bt.append(cmap_bt((bt+btextra)/(btrange+btextra))) if bt % 1. == 0: labels_bt.append(r'$B/T={:0.0f}$'.format(bt)) elif 10.*bt % 1. == 0: labels_bt.append(r'$B/T={:0.1f}$'.format(bt)) else: labels_bt.append(r'$B/T={:0.2f}$'.format(bt)) D2textra = 0. D2trange = 1. for D2t in D2t_arr: color_arr_D2t.append(cmap_D2t((D2t+D2textra)/(D2trange))) if D2t % 1. == 0: labels_D2t.append(r'$D_2/T='+r'{:0.0f}$'.format(D2t)) elif 10.*D2t % 1. == 0: labels_D2t.append(r'$D_2/T='+r'{:0.1f}$'.format(D2t)) else: labels_D2t.append(r'$D_2/T='+r'{:0.2f}$'.format(D2t)) ls_arr_bt, dashes_arr_bt, lw_arr_bt, labels_reratio_bt = ([] for _ in range(4)) dashlen = [5, 10, 20] if len(ReB_to_ReD_arr) > 5: dashlen = np.arange(5,(len(ReB_to_ReD_arr)-1)*5, 5) elif len(ReB_to_ReD_arr) < 5: dashlen = dashlen[:len(ReB_to_ReD_arr)-2] dashlen = dashlen[::-1] # reverse for i, reBtoreD in enumerate(ReB_to_ReD_arr): if i == 0: lw_arr_bt.append(1.3) ls_arr_bt.append('-') dashes_arr_bt.append(None) elif i == len(ReB_to_ReD_arr)-1: lw_arr_bt.append(1.) ls_arr_bt.append(':') dashes_arr_bt.append(None) else: lw_arr_bt.append(1.) ls_arr_bt.append('--') dashes_arr_bt.append((dashlen[i-1], 5)) lblstr = r'$R_{\mathrm{e,bulge}}/R_{\mathrm{e,disk}}=' if reBtoreD % 1. == 0: lblstr += '{:0.0f}$'.format(reBtoreD) elif 10.*reBtoreD % 1. == 0: lblstr += '{:0.1f}$'.format(reBtoreD) else: lblstr += r'{:0.2f}$'.format(reBtoreD) labels_reratio_bt.append(lblstr) ls_arr_D2t, dashes_arr_D2t, lw_arr_D2t, labels_reratio_D2t = ([] for _ in range(4)) dashlen = [5, 10, 20] if len(ReD2_to_ReD_arr) > 5: dashlen = np.arange(5,(len(ReD2_to_ReD_arr)-1)*5, 5) elif len(ReD2_to_ReD_arr) < 5: dashlen = dashlen[:len(ReD2_to_ReD_arr)-2] dashlen = dashlen[::-1] # reverse for i, reD2toreD in enumerate(ReD2_to_ReD_arr): if (reD2toreD == 1.): lw_arr_D2t.append(1.3) ls_arr_D2t.append('-') dashes_arr_D2t.append(None) elif (i == len(ReD2_to_ReD_arr)-1): lw_arr_D2t.append(1.) ls_arr_D2t.append(':') dashes_arr_D2t.append(None) else: lw_arr_D2t.append(1.) ls_arr_D2t.append('--') dashes_arr_D2t.append((dashlen[i-1], 5)) lblstr = r'$R_{\mathrm{e,disk_2}}/R_{\mathrm{e,disk}}=' if reD2toreD % 1. == 0: lblstr += '{:0.0f}$'.format(reD2toreD) elif 10.*reD2toreD % 1. == 0: lblstr += '{:0.1f}$'.format(reD2toreD) else: lblstr += r'{:0.2f}$'.format(reD2toreD) labels_reratio_D2t.append(lblstr) q_arr = [1.] delR = 0.01 xlim = [0., 5.] RtoRed_arr = np.arange(xlim[0], xlim[1]+delR, delR) Redisk = 1. # placeholder r_arr = RtoRed_arr * Redisk # Get values: # read fits tables, construct ks dict..... ks_dict = {} invq_disk = invq_bulge = invq_disk2 = 1. total_mass = 1.e11 for k, bt in enumerate(bt_arr): ks_dict['bt={}'.format(bt)] = {} for j, reBtoreD in enumerate(ReB_to_ReD_arr): alphan = -1. * interp_profiles.interpolate_sersic_profile_dlnrho_dlnR_two_component_nearest(R=r_arr, mass_comp1=total_mass*(1.-bt), mass_comp2=total_mass*bt, Reff_comp1=Redisk, n_comp1=nDisk, invq_comp1=invq_disk, Reff_comp2=reBtoreD*Redisk, n_comp2=nBulge, invq_comp2=invq_bulge, path=table_path) ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)] = {} ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)]['alpha'] = alphan ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)]['color'] = color_arr_bt[k] ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)]['ls'] = ls_arr_bt[j] ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)]['lw'] = lw_arr_bt[j] ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)]['dashes'] = dashes_arr_bt[j] ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)]['label_bt'] = labels_bt[k] ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)]['label_reratio'] = labels_reratio_bt[j] for k, D2t in enumerate(D2t_arr): ks_dict['D2t={}'.format(D2t)] = {} for j, reD2toreD in enumerate(ReD2_to_ReD_arr): alphan = -1. * interp_profiles.interpolate_sersic_profile_dlnrho_dlnR_two_component_nearest(R=r_arr, mass_comp1=total_mass*(1.-D2t), mass_comp2=total_mass*D2t, Reff_comp1=Redisk, n_comp1=nDisk, invq_comp1=invq_disk, Reff_comp2=reD2toreD*Redisk, n_comp2=nDisk2, invq_comp2=invq_disk2, path=table_path) ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)] = {} ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)]['alpha'] = alphan ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)]['color'] = color_arr_D2t[k] ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)]['ls'] = ls_arr_D2t[j] ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)]['lw'] = lw_arr_D2t[j] ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)]['dashes'] = dashes_arr_D2t[j] ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)]['label_D2t'] = labels_D2t[k] ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)]['label_reratio'] = labels_reratio_D2t[j] xlabel = r'$R/R_{\mathrm{e,disk}}$' ylabels = [r'$\alpha_{\mathrm{tot,\ disk+bulge}}(R)$', r'$\alpha_{\mathrm{tot,\ disk+disk_2}}(R)$'] ylims = [[0., 10.], [0., 10.]] titles = [None] lw = 1.3 ###################################### # Setup plot: f = plt.figure() scale = 4. n_cols = 0 if show_BT: n_cols += 1 if show_D2T: n_cols += 1 pad_outer = 0.2 f.set_size_inches(scale*(n_cols + pad_outer*(n_cols-1)),scale) pad_outer = 0.2 gs = gridspec.GridSpec(1, n_cols, wspace=pad_outer) axes = [] for i in range(n_cols): axes.append(plt.subplot(gs[0,i])) i = -1 ############ # BT if show_BT: i += 1 ax = axes[i] ylim = ylims[i] for j, reBtoreD in enumerate(ReB_to_ReD_arr): for k, bt in enumerate(bt_arr): alpha_arr = ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)]['alpha'] color = ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)]['color'] ls = ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)]['ls'] lw = ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)]['lw'] dashes = ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)]['dashes'] if len(alpha_arr) != len(r_arr): raise ValueError if (j==0): if (ls=='-'): lbl = ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)]['label_bt'] else: lbl = None else: lbl = None if dashes is not None: ax.plot(r_arr/Redisk, alpha_arr, ls=ls, color=color, dashes=dashes, lw=lw, label=lbl, zorder=-j) ax.axvline(x=reBtoreD, ls=ls, lw=1., dashes=dashes, color='lightgrey', zorder=-20.) else: ax.plot(r_arr/Redisk, alpha_arr, ls=ls, color=color, lw=lw, label=lbl, zorder=-j) ax.axvline(x=reBtoreD, ls=ls, lw=1., color='lightgrey', zorder=-20.) # Phantom plot: if (j==0) & (ls != '-'): lbl = ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)]['label_bt'] ax.plot(r_arr*np.NaN, r_arr*np.NaN, ls='-', color=color, lw=lw, label=lbl, zorder=-j) # Add self-grav comparison: nextra = 1 if (j==0) & (k==len(bt_arr)-1): lblSG = 'Self-grav' ax.plot(r_arr/Redisk, 3.36*(r_arr/Redisk), ls='-.', color='dimgrey', lw=1., label=lblSG, zorder=-15.) # Phandom plots for legend: if (k==len(bt_arr)-1): indbt_color = 0 color = ks_dict['bt={}'.format(bt_arr[indbt_color])]['reBtoreD={}'.format(reBtoreD)]['color'] lbl = ks_dict['bt={}'.format(bt)]['reBtoreD={}'.format(reBtoreD)]['label_reratio'] if dashes is not None: ax.plot(r_arr/Redisk, r_arr*np.NaN, ls=ls, color=color, dashes=dashes, lw=lw, label=lbl) else: ax.plot(r_arr/Redisk, r_arr*np.NaN, ls=ls, color=color, lw=lw, label=lbl) ax.axvline(x=1., ls=':', color='lightgrey', zorder=-20.) ax.set_xlim(xlim) ax.set_ylim(ylim) ax.set_xlabel(xlabel, fontsize=fontsize_labels) ax.set_ylabel(ylabels[i], fontsize=fontsize_labels) ax.tick_params(labelsize=fontsize_ticks) ax.xaxis.set_minor_locator(MultipleLocator(0.2)) ax.xaxis.set_major_locator(MultipleLocator(1.0)) if (ylim[1]-ylim[0]) > 9: ax.yaxis.set_minor_locator(MultipleLocator(0.5)) ax.yaxis.set_major_locator(MultipleLocator(2.)) elif (ylim[1]-ylim[0]) > 3: ax.yaxis.set_minor_locator(MultipleLocator(0.2)) ax.yaxis.set_major_locator(MultipleLocator(1.)) elif (ylim[1]-ylim[0]) > 1: ax.yaxis.set_minor_locator(MultipleLocator(0.1)) ax.yaxis.set_major_locator(MultipleLocator(0.5)) else: ax.yaxis.set_minor_locator(MultipleLocator(0.02)) ax.yaxis.set_major_locator(MultipleLocator(0.1)) ######################################## # Annotate Sersic indices: ann_str = r'$n_{\mathrm{disk}} = '+r'{:0.0f}$'.format(nDisk) ann_str2 = r'$n_{\mathrm{bulge}} = '+r'{:0.0f}$'.format(nBulge) xdelt = 0.075 ypos = 0.2275 ydelt = 0.055 xy = (1.-xdelt, ypos+ydelt) xy2 = (1.-xdelt, ypos) ax.annotate(ann_str, xy=xy, va='bottom', ha='right', fontsize=fontsize_ann, xycoords='axes fraction', color=cmap_bt(0.+btextra)) ax.annotate(ann_str2, xy=xy2, va='bottom', ha='right', fontsize=fontsize_ann, xycoords='axes fraction', color=cmap_bt(0.9+btextra)) ######################################## frameon = False framealpha = 1. borderpad = 0.2 fontsize_leg_tmp = fontsize_leg labelspacing= 0.2 handletextpad= 0.5 fancybox = False edgecolor='None' facecolor = 'white' # Split into two legends: len_ax1 = len(bt_arr) handles, labels_leg = ax.get_legend_handles_labels() neworder = range(len_ax1) handles_arr, labels_arr = ([] for _ in range(2)) for ii in neworder: handles_arr.append(handles[ii]) labels_arr.append(labels_leg[ii]) # COMPARISON SAMPLE: incl Self-grav neworder3 = range(len_ax1, len_ax1+nextra) handles_arr3, labels_arr3 = ([] for _ in range(2)) for ii in neworder3: handles_arr3.append(handles[ii]) labels_arr3.append(labels_leg[ii]) neworder2 = range(len_ax1+nextra, len(handles)) handles_arr2, labels_arr2 = ([] for _ in range(2)) for ii in neworder2: handles_arr2.append(handles[ii]) labels_arr2.append(labels_leg[ii]) loc1 = 'upper left' bbox_to_anchor1 = (0., 1.) legend1 = ax.legend(handles_arr, labels_arr, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, frameon=frameon, loc=loc1, bbox_to_anchor=bbox_to_anchor1, numpoints=1, scatterpoints=1,fontsize=fontsize_leg_tmp, fancybox=fancybox,edgecolor=edgecolor, facecolor=facecolor, framealpha=framealpha) patch = mpatches.FancyBboxPatch((0.1,6.45), 1.1, 3.55, boxstyle="square,pad=0.", fc="white", lw=0.) ax.add_patch(patch) loc2 = 'lower right' bbox_to_anchor2 = (1., -0.025) handlelength = 4.1 borderpad = 0.5 legend2 = ax.legend(handles_arr2, labels_arr2, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, frameon=frameon, handlelength=handlelength, loc=loc2, bbox_to_anchor=bbox_to_anchor2, numpoints=1, scatterpoints=1,fontsize=fontsize_leg_tmp, fancybox=fancybox,edgecolor=edgecolor, facecolor=facecolor, framealpha=framealpha) loc3 = 'upper right' bbox_to_anchor3 = (1., 1.) legend3 = ax.legend(handles_arr3, labels_arr3, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, frameon=frameon, loc=loc3, bbox_to_anchor=bbox_to_anchor3, handlelength=2.4, numpoints=1, scatterpoints=1,fontsize=fontsize_leg_tmp, fancybox=fancybox,edgecolor=edgecolor, facecolor=facecolor, framealpha=framealpha) ax.add_artist(legend1) ax.add_artist(legend2) ax.add_artist(legend3) ############ # D2T if show_D2T: i += 1 ax = axes[i] ylim = ylims[i] for j, reD2toreD in enumerate(ReD2_to_ReD_arr): for k, D2t in enumerate(D2t_arr): alpha_arr = ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)]['alpha'] color = ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)]['color'] ls = ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)]['ls'] lw = ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)]['lw'] dashes = ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)]['dashes'] if len(alpha_arr) != len(r_arr): raise ValueError if (j==0): if (ls=='-'): lbl = ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)]['label_D2t'] else: lbl = None else: lbl = None if dashes is not None: ax.plot(r_arr/Redisk, alpha_arr, ls=ls, color=color, dashes=dashes, lw=lw, label=lbl, zorder=-j) ax.axvline(x=reD2toreD, ls=ls, lw=1., dashes=dashes, color='lightgrey', zorder=-20.) else: ax.plot(r_arr/Redisk, alpha_arr, ls=ls, color=color, lw=lw, label=lbl, zorder=-j) ax.axvline(x=reD2toreD, ls=ls, lw=1., color='lightgrey', zorder=-20.) # Phantom plot: if (j==0) & (ls != '-'): lbl = ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)]['label_D2t'] ax.plot(r_arr*np.NaN, r_arr*np.NaN, ls='-', color=color, lw=lw, label=lbl, zorder=-j) # Add self-grav comparison: nextra = 1 if (j==0) & (k==len(D2t_arr)-1): lblSG = 'Self-grav' ax.plot(r_arr/Redisk, 3.36*(r_arr/Redisk), ls='-.', color='dimgrey', lw=1., label=lblSG, zorder=-15.) # Phandom plots for legend: if (k==len(D2t_arr)-1): indd2t_color = 0 color = ks_dict['D2t={}'.format(D2t_arr[indd2t_color])]['reD2toreD={}'.format(reD2toreD)]['color'] lbl = ks_dict['D2t={}'.format(D2t)]['reD2toreD={}'.format(reD2toreD)]['label_reratio'] if dashes is not None: ax.plot(r_arr/Redisk, r_arr*np.NaN, ls=ls, color=color, dashes=dashes, lw=lw, label=lbl) else: ax.plot(r_arr/Redisk, r_arr*np.NaN, ls=ls, color=color, lw=lw, label=lbl) ax.set_xlim(xlim) ax.set_ylim(ylim) ax.set_xlabel(xlabel, fontsize=fontsize_labels) ax.set_ylabel(ylabels[i], fontsize=fontsize_labels) ax.tick_params(labelsize=fontsize_ticks) ax.xaxis.set_minor_locator(MultipleLocator(0.2)) ax.xaxis.set_major_locator(MultipleLocator(1.0)) if (ylim[1]-ylim[0]) > 9: ax.yaxis.set_minor_locator(MultipleLocator(0.5)) ax.yaxis.set_major_locator(MultipleLocator(2.)) elif (ylim[1]-ylim[0]) > 3: ax.yaxis.set_minor_locator(MultipleLocator(0.2)) ax.yaxis.set_major_locator(MultipleLocator(1.)) elif (ylim[1]-ylim[0]) > 1: ax.yaxis.set_minor_locator(MultipleLocator(0.1)) ax.yaxis.set_major_locator(MultipleLocator(0.5)) else: ax.yaxis.set_minor_locator(MultipleLocator(0.02)) ax.yaxis.set_major_locator(MultipleLocator(0.1)) ######################################## # Annotate Sersic indices: ann_str = r'$n_D = {:0.0f}$'.format(nDisk) ann_str2 = r'$n_{D_2} = '+r'{:0.1f}$'.format(nDisk2) xdelt = 0.21 ypos = 0.2 ydelt = 0.05 xy = (1.-xdelt, ypos+ydelt) xy2 = (1.-xdelt, ypos) ax.annotate(ann_str, xy=xy, va='bottom', ha='left', fontsize=fontsize_ann, xycoords='axes fraction', color=cmap_D2t(0.+D2textra)) ax.annotate(ann_str2, xy=xy2, va='bottom', ha='left', fontsize=fontsize_ann, xycoords='axes fraction', color=cmap_D2t(0.9+D2textra)) ######################################## frameon = True framealpha = 1. borderpad = 0.75 fontsize_leg_tmp = fontsize_leg labelspacing= 0.2 handletextpad= 0.5 fancybox = False edgecolor='None' facecolor = 'white' # Split into two legends: len_ax1 = len(D2t_arr) handles, labels_leg = ax.get_legend_handles_labels() neworder = range(len_ax1) handles_arr, labels_arr = ([] for _ in range(2)) for ii in neworder: handles_arr.append(handles[ii]) labels_arr.append(labels_leg[ii]) # COMPARISON SAMPLE: incl Self-grav neworder3 = range(len_ax1, len_ax1+nextra) handles_arr3, labels_arr3 = ([] for _ in range(2)) for ii in neworder3: handles_arr3.append(handles[ii]) labels_arr3.append(labels_leg[ii]) neworder2 = range(len_ax1+nextra, len(handles)) handles_arr2, labels_arr2 = ([] for _ in range(2)) for ii in neworder2: handles_arr2.append(handles[ii]) labels_arr2.append(labels_leg[ii]) loc1 = 'upper left' bbox_to_anchor1 = (0., 1.) legend1 = ax.legend(handles_arr, labels_arr, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, frameon=frameon, loc=loc1, bbox_to_anchor=bbox_to_anchor1, numpoints=1, scatterpoints=1,fontsize=fontsize_leg_tmp, fancybox=fancybox,edgecolor=edgecolor, facecolor=facecolor, framealpha=framealpha) loc2 = 'lower right' bbox_to_anchor2 = (1., -0.025) handlelength = 5. frameon = False legend2 = ax.legend(handles_arr2, labels_arr2, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, frameon=frameon, handlelength=handlelength, loc=loc2, bbox_to_anchor=bbox_to_anchor2, numpoints=1, scatterpoints=1,fontsize=fontsize_leg_tmp, fancybox=fancybox,edgecolor=edgecolor, facecolor=facecolor, framealpha=framealpha) loc3 = 'upper right' bbox_to_anchor3 = (1., 1.) frameon = False legend3 = ax.legend(handles_arr3, labels_arr3, labelspacing=labelspacing, borderpad=borderpad, handletextpad=handletextpad, frameon=frameon, loc=loc3, bbox_to_anchor=bbox_to_anchor3, numpoints=1, scatterpoints=1,fontsize=fontsize_leg_tmp, fancybox=fancybox,edgecolor=edgecolor, facecolor=facecolor, framealpha=framealpha) ax.add_artist(legend1) ax.add_artist(legend2) ax.add_artist(legend3) if fileout is not None: plt.savefig(fileout, bbox_inches='tight', dpi=600) plt.close() else: plt.show() return None if __name__ == "__main__": # From the command line, call the wrapper to make all plots. # Input args: output_path, table_path output_path = sys.argv[1] table_path = sys.argv[2] make_all_paper_plots(output_path=output_path, table_path=table_path)