Added Fig 9b

Fig9b_gallery_stars.py

+import os
+import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+from matplotlib import patches
+
+import numpy as np
+np.seterr(divide = 'ignore')
+
+from swiftsimio import load
+
+import matplotlib.pylab as pylab
+params = {'legend.fontsize': 'large',
+         'axes.labelsize': 'x-large',
+         'axes.titlesize':'x-large',
+         'figure.titlesize':'x-large',
+         'legend.title_fontsize':'large',
+         'legend.fontsize':'large',
+         'xtick.labelsize':'large',
+         'ytick.labelsize':'large'}
+pylab.rcParams.update(params)
+
+import sys
+sys.path.insert(0, "helpers/")
+
+from plotannotation import add_arrows
+
+from stellar_surface_density_plot import add_stellar_surface_density
+from stellar_surface_density_plot import add_fof_groups
+
+from simulations_details import get_simulation_time_Myr
+from simulations_details import get_simulation_details_all
+
+from select_simulations_for_figures import Fig9_runnames as runnames
+from select_simulations_for_figures import Fig9_filebase_runs as filebase_runs
+from select_simulations_for_figures import Fig9_filebase_reruns as filebase_reruns
+from select_simulations_for_figures import Fig9_filebase_fofruns as filebase_fofruns
+
+from simulation_path import snapshotnumber
+
+
+snapshotnumber_rerun = 0
+############################################################################
+# Set output filename
+############################################################################
+outputname = "Fig7_gallery_stars.png"
+
+############################################################################
+# Some extra plotting settings 
+############################################################################
+
+axis_xsize = 0.40
+axis_ysize = axis_xsize
+axis_dy    = 0.015
+axis_dx_left   = 0.02
+axis_dx        = 0.2
+axis_dy_bottom = 0.02
+axis_overlap   = 0.1
+axis_overlap_x = 0.02
+
+# #################################################################
+# Creates 5 stellar mass density plots 
+# of the 5 runs identified in helpers/Fig7_runs_to_plot.py
+# Labels need to be updated if different runs are selected
+# #################################################################
+
+def make_stars_gallery_plot(outputname):
+    mB_runs, eps_runs, h_min_runs, h_min_ratio_runs, XH_runs, esf_runs, \
+    mB_reruns, eps_reruns, h_min_reruns, h_min_ratio_reruns, XH_reruns, esf_reruns = \
+                                     get_simulation_details_all(filebase_runs, filebase_reruns, \
+                                                           snapshotnumber)
+
+    #get simulation time
+    time_Myr = get_simulation_time_Myr(filebase_runs[0] + '%4.4i'%(snapshotnumber) + '.hdf5')
+
+    fig = plt.figure(figsize = (7.5*0.9,7*0.9))
+    axframe = fig.add_axes([0,0,1,1])
+    for axis in ['top','bottom','left','right']:
+        axframe.spines[axis].set_linewidth(10)
+        axframe.spines[axis].set_color("black")
+        axframe.spines[axis].set_zorder(0)
+
+    ax = []
+
+    ####################################################### 
+    # Setup axes
+    ####################################################### 
+    # 0: top left panel
+    ax.append(fig.add_axes([axis_dx_left, \
+                            1. - axis_dy_bottom - axis_ysize, \
+                            axis_xsize, \
+                            axis_ysize]))
+    # 1: top right panel
+    ax.append(fig.add_axes([1. - axis_dx_left - axis_xsize , \
+                            1. - axis_dy_bottom - axis_ysize, \
+                            axis_xsize, \
+                            axis_ysize]))
+    # 2: middle panel
+    ax.append(fig.add_axes([0.5 - 0.5 * axis_xsize, \
+                            0.5 - 0.5 * axis_ysize, \
+                            axis_xsize, \
+                            axis_ysize]))
+    # 3: bottom left panel
+    ax.append(fig.add_axes([axis_dx_left, \
+                            axis_dy_bottom, \
+                            axis_xsize, \
+                            axis_ysize]))
+    # 4: bottom right panel
+    ax.append(fig.add_axes([1. - axis_dx_left - axis_xsize , \
+                            axis_dy_bottom, \
+                            axis_xsize, \
+                            axis_ysize]))
+
+    for iplot in range(len(filebase_reruns)):
+        snapfilename = filebase_reruns[iplot] + '%4.4i'%(snapshotnumber_rerun) + '.hdf5'
+        mass_map_plot = add_stellar_surface_density(ax[iplot], snapfilename)
+
+        foffilename  = filebase_fofruns[iplot]
+        add_fof_groups(ax[iplot], foffilename)
+   
+    ax[2].text(1.2, 0.5, "Increase m$_{\mathrm{B}}$\n (x 8)", transform = ax[2].transAxes, fontsize = 'x-large', va = 'center', ha = 'left')
+    ax[2].text(0.5, 1.35, "Increase $\epsilon$\n (x 2)", transform = ax[2].transAxes, fontsize = 'x-large', va = 'bottom', ha = 'center')
+    ax[2].text(0.5,-0.35, "Decrease h$_{\mathrm{min}}$\n (x 0.1)", transform = ax[2].transAxes, fontsize = 'x-large', va = 'top', ha = 'center')
+    ax[2].text(-0.2,0.5, "Increase e$_{\mathrm{sf}}$\n (x 3)", transform = ax[2].transAxes, fontsize = 'x-large', va = 'center', ha = 'right')
+
+    add_arrows(fig, ax)
+
+    plt.savefig(outputname, dpi = 150)
+    plt.close()
+    print ("Figure saved: "+outputname)
+    
+    return
+
+make_stars_gallery_plot(outputname)

helpers/stellar_surface_density_plot.py

+import os
+import numpy as np
+np.seterr(divide = 'ignore')
+
+import h5py
+import matplotlib
+import matplotlib.pyplot as plt
+import matplotlib.pyplot as plt
+from matplotlib import gridspec
+from mpl_toolkits.axes_grid1.inset_locator import InsetPosition
+
+from unyt import unyt_array, unyt_quantity
+from unyt import msun, kpc, g, cm, K, pc
+from unyt import boltzmann_constant_cgs as kB
+from unyt import proton_mass_cgs as mH
+
+import swiftsimio as sw
+from swiftsimio import load
+from swiftsimio.visualisation.projection import project_pixel_grid
+from swiftsimio.visualisation.slice import kernel_gamma
+from swiftsimio.visualisation.smoothing_length_generation import (
+    generate_smoothing_lengths,
+)
+# ignore cosmology warnings 
+import warnings
+warnings.filterwarnings('ignore', module='swiftsimio')
+
+from simulations_details import get_simulation_time_Myr
+from simulations_details import get_simulation_details_all
+
+from simulation_plotparameter import MaxAgeYoungStarsMyr
+from simulation_plotparameter import r_img_kpc
+from simulation_plotparameter import slice_thickness_in_kpc
+from simulation_plotparameter import npix
+from simulation_plotparameter import nbins
+from simulation_plotparameter import vmin
+from simulation_plotparameter import vmax
+from simulation_plotparameter import vmax_sigma_line
+from simulation_plotparameter import cmap_stars
+from simulation_plotparameter import colors
+from simulation_plotparameter import sigma_min
+from simulation_plotparameter import eta_res 
+from simulation_plotparameter import mask_central_kpc
+
+# #################################################################
+# Calculates the fraction of pixels in a mass map above a 
+# stellar surface density Sigma for various Sigma
+# #################################################################
+def get_cumulative_pixel_fraction(mass_map):
+
+    surface_densities = np.power(10., mass_map.ravel())
+    surface_densities[surface_densities < sigma_min] = 0.
+    indx = np.where(surface_densities >= sigma_min)[0]
+    histbins_sigma = np.logspace(vmin, vmax_sigma_line, num = nbins, endpoint=True, base=10.0)
+    hist, edges = np.histogram(surface_densities, bins = histbins_sigma, density = False)
+    cumhist = 1. -  np.cumsum(hist) / float(len(indx))
+
+    return np.log10(histbins_sigma[:-1]), cumhist
+
+
+# #################################################################
+# Adds stellar surface density image to subplot
+# #################################################################
+def add_stellar_surface_density(ax, filename):
+    data = load(filename)
+    t_Myr = data.metadata.time.to('Myr')
+    data.metadata.boxsize.convert_to_units("kpc")
+    boxsize = data.metadata.boxsize
+
+    data.stars.birth_times.convert_to_units('Myr')
+    data.stars.coordinates.convert_to_units("kpc")
+    data.stars.masses.convert_to_units("Msun")
+    Mstar = np.sum(data.stars.masses)
+
+    data.metadata.time.convert_to_units('Myr')
+    tmin = np.maximum(data.metadata.time.value - MaxAgeYoungStarsMyr, 0.)
+
+    data.stars.masses = data.stars.masses[data.stars.birth_times >= tmin]
+    data.stars.coordinates = data.stars.coordinates[data.stars.birth_times >= tmin, :]
+
+    total_usermass = np.sum(data.stars.masses)
+
+    data.stars.smoothing_lengths = generate_smoothing_lengths(
+        coordinates=data.stars.coordinates,
+        boxsize=data.metadata.boxsize,
+        kernel_gamma=kernel_gamma,
+        neighbours=11,
+        speedup_fac=1,
+        dimension=3,
+    )
+
+    deltax = unyt_quantity(r_img_kpc, 'kpc')
+    deltay = unyt_quantity(r_img_kpc, 'kpc')
+
+    range_plot = [0.5 * boxsize[0] - deltax, 0.5 * boxsize[0] + deltax , 0.5 * boxsize[1] - deltay, 0.5 * boxsize[1] + deltay]
+    extent = [range_plot[0].value-0.5 * boxsize[0].value, \
+    range_plot[1].value-0.5 * boxsize[0].value, \
+    range_plot[2].value-0.5 * boxsize[1].value, \
+    range_plot[3].value-0.5 * boxsize[1].value ]
+
+    pixel_grid = project_pixel_grid(
+            data.stars,
+            boxsize=data.metadata.boxsize,
+            resolution=int(npix),
+            project="masses",
+             region=range_plot,
+            )
+    x_range = range_plot[1] - range_plot[0]
+    y_range = range_plot[3] - range_plot[2]
+    units = 1.0 / (x_range * y_range)
+    # Unfortunately this is required to prevent us from {over,under}flowing
+    # the units...
+    units.convert_to_units(1.0 / (x_range.units * y_range.units))
+    units *= getattr(data.stars, 'masses').units
+
+    mass_map = unyt_array(pixel_grid, units=units)
+    mass_map.convert_to_units(msun / kpc**2)
+
+    mass_map_plot = np.log10(mass_map.value)
+
+    nonzero = mass_map_plot >= 1e-2
+    min_nonzero = np.min(mass_map_plot[nonzero])
+    mass_map_plot[~nonzero] = min_nonzero
+
+    cblab = 'Star'
+    cmap_loc = matplotlib.cm.get_cmap(cmap_stars)
+    tcolor = 'white'
+    ccolor = 'black'
+
+
+    # mask with circle
+    lx, ly = mass_map_plot.shape
+    X, Y = np.ogrid[0:lx, 0:ly]
+    mask = (X - lx / 2) ** 2 + (Y - ly / 2) ** 2 > lx * ly / 4
+    mass_map_plot[mask] = np.nan
+
+    # make plot 
+    ax.tick_params(labelleft = False, labelbottom = False)
+    ax.set_xticks([])
+    ax.set_yticks([])
+    ax.set_axis_off()
+
+    im = ax.imshow(mass_map_plot, cmap=cmap_stars, extent = extent, vmin = vmin, vmax = vmax, origin = 'lower')
+    circle = plt.Circle((0, 0), 0.99 * r_img_kpc, color=ccolor, fill=False, linewidth = 2)
+    ax.add_artist(circle)
+
+    return mass_map_plot
+
+def add_fof_groups(ax, foffile):
+
+    with h5py.File(foffile, 'r') as f:
+        group_centers_kpc = f["Groups/Centres"][:] 
+        group_masses_msun = f["Groups/Masses"][:]  * f["Groups/Masses"].attrs["Conversion factor to physical CGS (including cosmological corrections)"] / msun
+        boxsize_kpc =  f["Header"].attrs["BoxSize"]
+    
+    group_masses_msun.convert_to_units('1/g')
+    group_masses_msun = group_masses_msun.value
+
+    group_centers_kpc[:,0] = group_centers_kpc[:,0] - 0.5 * boxsize_kpc[0]
+    group_centers_kpc[:,1] = group_centers_kpc[:,1] - 0.5 * boxsize_kpc[1]
+    group_centers_kpc[:,2] = group_centers_kpc[:,2] - 0.5 * boxsize_kpc[2]
+
+    r = np.sqrt( group_centers_kpc[:,1]**2 + group_centers_kpc[:,0]**2 )
+
+    ax.scatter(group_centers_kpc[(r>mask_central_kpc) & (group_masses_msun < 1.e8), 1], \
+               group_centers_kpc[(r>mask_central_kpc) & (group_masses_msun < 1.e8), 0], marker = 'o', edgecolor = 'yellow', facecolor = 'none',
+               s = 110, linewidths = 2.5)
+
+    ax.scatter(group_centers_kpc[(r>mask_central_kpc) & (group_masses_msun > 1.e8), 1], \
+               group_centers_kpc[(r>mask_central_kpc) & (group_masses_msun > 1.e8), 0], marker = 'D', edgecolor = 'yellow', facecolor = 'none', \
+               s = 220, linewidths = 2.5)
+    return
+
+