"""Demonstration of how to use cross-matching to combine two catalogues. The example here is how to obtain the stellar formation time in a snipshot, for which this information is not written out. This is essentially a more complex version of the star_density.py example, which you may want to have a look at first. This script also demonstrates some basic approaches for working with snipshots in general (viz. finding the right one, finding galaxy positions). """ # Import required packages import numpy as np import hydrangea as hy import matplotlib.pyplot as plt from pdb import set_trace # Set script parameters sim_index = 0 # Which simulation do we want? target_z = 0.31 # Target redshift, we'll find closest snepshot gal_id = 7679 # Which galaxy do we want? imsize = 30 # (Half-)size of analysis region, in kpc nbins = 30 # Number of bins for plotting ptype = 4 # Look at stars here plotloc = 'galaxy_age.png' # Where to save the output plot? # First, let's find the target snepshot. Load the full list of snepshot # times, including snipshots (including the 25-Myr-filler ones): snep_zred = hy.snep_times(time_type='zred', snep_list='default_long') # Find the one closest to our target redshift ind_snep = np.argmin(np.abs(snep_zred - target_z)) print(f"Closest snepshot to target z ({target_z:.3f}) is {ind_snep}, " f"at z={snep_zred[ind_snep]}") # Translate this into snepshot type and number for the particular simulation sim = hy.Simulation(index=sim_index) root_index, aexp, source_type, source_num = hy.get_snepshot_indices( sim.run_dir, 'default_long', index=ind_snep) print(f"Target snepshot identified as {source_type} {source_num}.") # Next step: find the location of the galaxy in this snepshot # NOTE: these are stored in code units in the file. We must explicitly # convert it to have the inputs to ReadRegion all in pMpc. dataset = f'Snepshot_{root_index:04d}/Coordinates' gal_position = hy.hdf5.read_data(sim.galaxy_path_loc, dataset, read_index=gal_id) aexp_factor = hy.hdf5.read_attribute(sim.galaxy_path_loc, dataset, 'aexp-factor') h_factor = hy.hdf5.read_attribute(sim.galaxy_path_loc, dataset, 'h-factor') gal_position *= (aexp_factor * h_factor) print(f"\nGalaxy position is ({gal_position[0]:.2f}, {gal_position[1]:.2f}, " f"{gal_position[2]:.2f}) pMpc.") # Now we can set up a ReadRegion object to load the stars close to the # position of the galaxy... print("\nSetting up reader:") if source_type == 'snap': snepshot_file = sim.get_snap_file(source_num) elif source_type == 'snip': snepshot_file = sim.get_snipshot_file(source_num) else: print(f"Unexpected source type '{source_type}'!") set_trace() parts = hy.ReadRegion(snepshot_file, # (First) file to read ptype, # Particle type to read gal_position, # Centre of region imsize/1e3, # Size of region [pMpc] shape='cube') # Shape of region print("") # To get the stellar age, we need to link each particle to a subsequent # snapshot. ref_snap_index = hy.get_next_snapshot(parts.aexp) ref_snap_file = sim.get_snap_file(ref_snap_index) print(f"Next snapshot is {ref_snap_index}") # For simplicity, we will only load particles around the target galaxy. # In principle, we could also set up a SplitFile reader to load the entire # star catalogue in the reference snapshot. gal_position_ref = hy.hdf5.read_data( sim.gps_loc, 'Centre', read_index=gal_id)[ref_snap_index, :] print("\nSetting up reader:") parts_ref = hy.ReadRegion(ref_snap_file, ptype, gal_position_ref, imsize/10) # Now set up a 'gate' to cross-match particles between the target snepshot # and the reference snapshot: print("\nNote how and now read 'ParticleIDs':") gate = hy.crossref.Gate(parts.ParticleIDs, parts_ref.ParticleIDs) print("") # We can now access the indices of particles in the (selected region of the) # reference snapshot corresponding to indices in the target snipshot indices_in_ref_all, subind_matched = gate.in_int() snip_formation_aexp = parts_ref.StellarFormationTime[indices_in_ref_all] # Make sure we could match all particles -- otherwise we loaded too small # a region in the ref snapshot! if len(subind_matched) != len(indices_in_ref_all): print("Could not match all stars, must load wider region in ref snap...") set_trace() # Aside: if we only wanted the (ref-snapshot) indices for e.g. the first 100 # particles in the target snepshot read-region, we could get these directly # by specifying them in the call to in_int(): indices_in_ref_first = gate.in_int(np.arange(100)) # Translate formation aexp to age at the time of the target snepshot # (aexp_to_time gives the cosmic time at star formation) snip_stars_age = (parts.time - hy.aexp_to_time(snip_formation_aexp)) # Get particle coordinates relative to subhalo centre (in kpc) part_relpos = (parts.Coordinates - gal_position[None, :]) * 1e3 # Make a simple 2D histogram map of mass-weighted stellar age # Note that we here don't convert the (temporary) mass_map to proper # surface density units (c.f. find_particles_near_subhalo.py example), # because we only need it to normalise the age map mass_map = np.histogram2d(part_relpos[:, 1], part_relpos[:, 0], bins=nbins, range=[[-imsize, imsize], [-imsize, imsize]], weights=parts.Mass)[0] print("\nMay get 'division-by-zero' warning now, which we can ignore") age_map = np.histogram2d(part_relpos[:, 1], part_relpos[:, 0], bins=nbins, range=[[-imsize, imsize], [-imsize, imsize]], weights=parts.Mass*snip_stars_age)[0] / mass_map # Plot the age map. Age is not defined for pixels without any stars, # which will have NaN value in age_map (left white in image) fig = plt.figure(figsize=(4/0.8, 4)) ax1 = fig.add_axes([0.15, 0.15, 0.65, 0.8]) im = plt.imshow(age_map, extent=[-imsize, imsize, -imsize, imsize], aspect='equal', interpolation='nearest', origin='lower', alpha=1.0, cmap=plt.cm.inferno, vmin=0, vmax=10) ax1.set_xlim((-imsize, imsize)) ax1.set_ylim((-imsize, imsize)) ax1.set_xlabel(r'$\Delta x$ [kpc]') ax1.set_ylabel(r'$\Delta y$ [kpc]') # Add a colour bar on the right side ax2 = fig.add_axes([0.81, 0.15, 0.05, 0.8]) ax2.set_xticks([]) ax2.set_yticks([]) cbar = plt.colorbar(im, cax=ax2, orientation='vertical') fig.text(0.95, 0.55, r'Stellar age [Gyr]', rotation=90.0, va='center', ha='left') # Save the figure plt.show plt.savefig(plotloc, transparent=False)