fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(15,5))
my_cmap = copy.copy(mpl.cm.get_cmap('inferno'))
my_cmap.set_bad((0,0,0))

ax[0].set_xlabel('$x_c$ [$mathrm{kpc}$]')
ax[0].set_ylabel('$y_c$ [$mathrm{kpc}$]')
ax[0].set_aspect('equal')
ax[0].hist2d(x_plot, y_plot, bins=(bin_num, bin_num), norm=mpl.colors.LogNorm(), cmap=my_cmap, range=[[-box_size, box_size], [-box_size, box_size]])

ax[1].set_xlabel('$y_c$ [$mathrm{kpc}$]')
ax[1].set_ylabel('$z_c$ [$mathrm{kpc}$]')
ax[1].set_aspect('equal')
ax[1].hist2d(y_plot, z_plot, bins=(bin_num, bin_num), norm=mpl.colors.LogNorm(), cmap=my_cmap, range=[[-box_size, box_size], [-box_size, box_size]])

ax[2].set_xlabel('$x_c$ [$mathrm{kpc}$]')
ax[2].set_ylabel('$z_c$ [$mathrm{kpc}$]')
ax[2].set_aspect('equal')
ax[2].hist2d(x_plot, z_plot, bins=(bin_num, bin_num), norm=mpl.colors.LogNorm(), cmap=my_cmap, range=[[-box_size, box_size], [-box_size, box_size]])

fig.suptitle('Dark matter density', y=0.99)
plt.tight_layout()
plt.savefig('dm_density.pdf', dpi=300)
plt.show()