""" General-purpose program to make image of a galaxy. -- Started 21-Mar-2019 -- Expanded 14-May-2019 """ import sys import socket hostname = socket.gethostname() if hostname == 'quasar.strw.leidenuniv.nl': sys.path.insert(0, '/net/quasar/data3/Hydrangea/OLD_EXAMPLE_CODE/' 'TOOLS_QUASAR/sim-utils/') else: sys.path.insert(0, '/net/quasar/data3/Hydrangea/OLD_EXAMPLE_CODE/sim-utils/') import sim_tools as st import yb_utils as yb import hydrangea_tools as ht import image_routines as ir from pdb import set_trace from astropy.io import ascii import numpy as np import os import matplotlib matplotlib.use('pdf') import matplotlib.pyplot as plt plt.style.use('dark_background') matplotlib.rcParams['font.family'] = 'serif' matplotlib.rcParams['font.serif'][0] = 'palatino' # Set up the simulation to image rootdir = '/net/quasar/data3/Hydrangea/' simtype = 'HYDRO' # ---------- Basic settings: simulation, galaxy, size, type ------------- isim = 0 # Simulation index (CE-xx) galID = 126 # Galaxy to center on (may be None) galSH = None # If galID is None, *subhalo* index to center on plot_snap = 29 # Snapshot index to plot ptype = 4 # Particle type to plot imsize = 0.1 # x/y size in Mpc (+/- vs. central galaxy) zsize = 0.1 # z size in Mpc (+/- vs. mid-plane); same as x/y if None # ---------- Projection settings --------------------------------------- projectionPlane = 'xy' # Image projection plane (may be None) camDir = [0, 0, -1] # Camera direction vector (may be left at [0,0,0]) camAngle = [0, 0, 0] # Camera angles (may left at [0,0,0]) align = None#'MIT' # Experimental mode to align frame to AM vector alignDir = 'face' # 'edge' or 'face' show = None#'MIT' # (only possible with Cantor as subhalo catalogue) # ------------- Image settings ---------------------------------------- imtype = 'gri' # 'gri' or 'mass', only for ptype=4 scale_hsml_with_age = False mag_range = [28, 17] lum_source = 'direct' gas_tmax = None # Only consider gas below a given T (can be None for all) fixedSmoothingLength = 0 # Leave at 0 to compute adaptive smoothing lengths desNGB = 32 # Number of neighbours for smoothing calculation numPix = 2500#2400 # Image sidelength in pixels inch = 6 # Image size in inches # ------------ Plot one galaxy's particles only ---------------------- ref_galaxy = None#4#3272#48 # Galaxy whose particles to consider exclusively ref_shi = None # Subhalo whose particles to consider exclusively (unless # ref_galaxy != None, or it is None) ref_snap = plot_snap # Snapshot in which particle membership is tested sh_catalogue = 'Subfind' # Subhalo catalogue to use ('Subfind' or 'Cantor') # ------------- Galaxy location indicators ---------------------------- mark_galaxies = False # Highlight locations of galaxies? mass_crit = 'Msub' # Selection criterion for which galaxies to show # ('Mstar', 'MstarPeak', 'Msub', or 'MsubPeak') mass_threshold = 10.5 # Minimum mass for galaxies to show plot_dead_galaxies = False # True: also show dead galaxies (in red). # NB: only effective if mass_crit == '[...]Peak' label_galaxies = True # True: also print galaxy IDs in image # (otherwise only a circle is drawn) labelSize = 4.0 # Specify size of galaxy labels print_galaxy_id = False # Print ID of central galaxy in top-right corner # ------------- Output settings -------------------------------------- save_maps = False # Store images as HDF5 as well as picture if ref_galaxy is None and ref_shi is None: refCode = '' elif ref_galaxy == galID or (ref_shi == galSH and ref_shi is not None): refCode = 'x' elif ref_galaxy is not None: refCode = 'x{:d}' .format(ref_galaxy) else: refCode = 'xS{:d}' .format(ref_shi) if len(refCode) and sh_catalogue == 'Subfind': refCode = refCode + 'SF' if galID is None: galCode = 'S-{:d}' .format(galSH) else: galCode = 'G-{:d}' .format(galID) sizeStr = '{:.2f}' .format(imsize) sizeStr = sizeStr.replace('.', 'p') if align is None: alCode = '' else: alCode = '_' + align + alignDir[0].upper() plotloc = ('/data/Hydrangea/' 'CE-{:d}_{:s}{:s}{:s}_PT-{:d}_{:s}_MR_{:s}_' .format(isim, galCode, refCode, alCode, ptype, sizeStr, imtype)) # ------------------------------------------------------------------------ # -------- Should not have to adjust anything below this line ------------ # ------------------------------------------------------------------------ rundir = rootdir + 'CE-' + str(isim) + '/' + simtype + '/' hldir = rundir + '/highlev/' fgtloc = hldir + 'FullGalaxyTables.hdf5' posloc = hldir + 'GalaxyPositionsSnap.hdf5' spiderloc = hldir + 'SpiderwebTables.hdf5' if sh_catalogue == 'Cantor': cantorloc = hldir + 'Cantor/GalaxyTables.hdf5' cantorloc_plot = hldir + 'Cantor/Cantor_{:03d}.hdf5' .format(plot_snap) shi = yb.read_hdf5(cantorloc, 'SubhaloIndex') sh2gal_plt = yb.read_hdf5(cantorloc_plot, 'Subhalo/Galaxy') else: subdir_plt = st.form_files(rundir, plot_snap) shi = yb.read_hdf5(fgtloc, 'SHI') sh2gal_plt = yb.read_hdf5( spiderloc, 'Subhalo/Snapshot_' + str(plot_snap).zfill(3) + '/Galaxy') if galID is None: galID = sh2gal_plt[galSH] pos_gal_all = yb.read_hdf5(posloc, 'Centre')[:, plot_snap, :] galpos = np.copy(pos_gal_all[galID, :]) #galpos += np.array((0.19, 0.175, -0.2)) # Set up reading region if zsize is None: zsize = imsize readRad = max([imsize, zsize])*np.sqrt(3) snaplistloc = rundir + "/sneplists/allsnaps.dat" snaplist = ascii.read(snaplistloc) aexpSnap = np.array(snaplist['aexp']) aexp_factor = aexpSnap[plot_snap] conv_astro_pos = aexpSnap[plot_snap]/0.6777 readRad_sim = readRad/conv_astro_pos galpos_sim = galpos/conv_astro_pos snapdir = st.form_files(rundir, plot_snap, 'snap') if ref_galaxy is None: if ref_shi is not None: ref_galaxy = yb.read_hdf5( spiderloc, 'Subhalo/Snapshot_' + str(ref_snap).zfill(3) + '/Galaxy')[ref_shi] if ref_galaxy is not None: ref_shi = shi[ref_galaxy, ref_snap] if sh_catalogue == 'Cantor': cantorloc_ref = hldir + 'Cantor/Cantor_{:03d}.hdf5' .format(ref_snap) cantoridloc_ref = (hldir + 'Cantor/Cantor_{:03d}_IDs.hdf5' .format(ref_snap)) ref_ids_all = yb.read_hdf5(cantoridloc_ref, 'IDs') ref_off_all = yb.read_hdf5(cantorloc_ref, 'Subhalo/OffsetType') #ref_len_all = yb.read_hdf5( # cantorloc, 'Snapshot_' + str(ref_snap).zfill(3) # + '/Subhalo/LengthType') ref_off = ref_off_all[ref_shi, ptype] ref_end = ref_off_all[ref_shi, ptype+1]# + ref_len_all[ref_shi] ref_ids = ref_ids_all[ref_off : ref_end] else: subdir = st.form_files(rundir, ref_snap) ref_ids_all = st.eagleread(subdir, 'IDs/ParticleID', astro=False) ref_off_all = st.eagleread(subdir, 'Subhalo/SubOffset', astro=False) ref_len_all = st.eagleread(subdir, 'Subhalo/SubLength', astro=False) ref_off = ref_off_all[ref_shi] ref_len = ref_len_all[ref_shi] ref_ids = ref_ids_all[ref_off : ref_off+ref_len] xBase = None if (align is not None or show is not None): if sh_catalogue != 'Cantor': cantorloc = hldir + 'Cantor/GalaxyTables.hdf5' cantorloc_snap = hldir + 'Cantor/Cantor_{:03d}.hdf5' .format(plot_snap) shi_align = yb.read_hdf5(cantorloc, 'SubhaloIndex') else: shi_align = shi gal_shi = shi_align[galID, plot_snap] extra_id = yb.read_hdf5( cantorloc_snap, '/Subhalo/Extra/ExtraIDs')[gal_shi] if extra_id < 0: edge_on = False show = False else: if align == 'AngMom' or show == 'AngMom': #am = yb.read_hdf5(cantorloc, 'Snapshot_' + str(plot_snap).zfill(3) # + '/Subhalo/Extra/Stars/AngularMomentum')[ # extra_id, 0, :] am = yb.read_hdf5(cantorloc_snap, '/Subhalo/AngularMomentum_Stars')[gal_shi, :] if align == 'AngMom': camDir = am if alignDir == 'edge': edge_on = True else: edge_on = False projectionPlane = None rho = 0 else: edge_on = False rho = 0 if align == 'MIT' or show == 'MIT': ax = yb.read_hdf5(cantorloc_snap, '/Subhalo/Extra/Stars/Axes')[extra_id, 0, :, :] axRat = yb.read_hdf5(cantorloc_snap, '/Subhalo/Extra/Stars/AxisRatios')[ extra_id, 0, :] if align == 'MIT': if alignDir == 'edge': xBase = ax[2, :].astype(np.float32) yBase = ax[0, :].astype(np.float32) zBase = ax[1, :].astype(np.float32) else: xBase = ax[2, :].astype(np.float32) yBase = ax[1, :].astype(np.float32) zBase = ax[0, :].astype(np.float32) """ camDir = ax[:, 1] n_vec = np.array((-camDir[0], -camDir[1], 1-camDir[2])) cz = -(n_vec[0]*camDir[0] + n_vec[1]*camDir[1])/camDir[2] c_vec = np.array((n_vec[0], n_vec[1], cz)) c_vec = c_vec / np.linalg.norm(c_vec) cosrho = np.dot(c_vec, ax[:, 0]) rho = np.arccos(cosrho)*180/np.pi """ edge_on = False projectionPlane = None rho = 0 else: edge_on = False rho = 0 else: edge_on = False else: edge_on = False rho = 0 #edge_on = False #projectionPlane = 'xy' readReg = ht.ReadRegion(snapdir, ptype, [*galpos_sim, readRad_sim]) pos = readReg.read_data("Coordinates", astro = True) vel = readReg.read_data("Velocity", astro = True) if ptype != 1: if ptype == 0 and imtype == 'hi': hiloc = ('/virgo/scratch/ybahe/HYDRANGEA/ANALYSIS/10r200/CE-{:d}/' .format(isim) + simtype + '/highlev/' 'HI_particle_masses_18Oct19.hdf5') if not os.path.exists(hiloc): print("Could not find HI file.") set_trace() mass = readReg.read_data( "H_Neutral", astro = False, filename = hiloc, PTName = '', singleFile = True) mass = np.clip(mass, 0, None) elif ptype == 0 and imtype == 'sfr': mass = readReg.read_data("StarFormationRate", astro = True) else: mass = readReg.read_data("Mass", astro = True) else: mass = np.zeros(pos.shape[0])+st.m_dm(rundir, astro = True, isnap=29) if ptype == 4 and imtype == 'gri': if lum_source == 'andrea_old': if plot_snap != 29: set_trace() mag_file = ht.clone_dir(hldir) + '/StellarMagnitudes_029.hdf5' load_mag = True lsf = True lptn = None elif lum_source == 'andrea_new': mag_file = st.form_files(rundir, plot_snap, 'lum') load_mag = True lsf = False lptn = 'restFrame' elif lum_source == 'direct': load_mag = False m_init = readReg.read_data('InitialMass', astro = True) z_star = readReg.read_data('SmoothedMetallicity') sft = readReg.read_data('StellarFormationTime', astro = False) if load_mag: mag_g = readReg.read_data( 'g-Magnitude', astro = False, filename = mag_file, singleFile = lsf, PTName = lptn) mag_r = readReg.read_data( 'r-Magnitude', astro = False, filename = mag_file, singleFile = lsf, PTName = lptn) mag_i = readReg.read_data( 'i-Magnitude', astro = False, filename = mag_file, singleFile = lsf, PTName = lptn) lum_g = 10.0**(-0.4*mag_g) lum_r = 10.0**(-0.4*mag_r) lum_i = 10.0**(-0.4*mag_i) else: age_star = st.aexp_to_age(sft, lbt_aexp = aexp_factor) lum_g = st.stellar_luminosity(m_init*1e10, z_star, age_star*1e9, 'g') lum_r = st.stellar_luminosity(m_init*1e10, z_star, age_star*1e9, 'r') lum_i = st.stellar_luminosity(m_init*1e10, z_star, age_star*1e9, 'i') if ptype == 0: hsml = readReg.read_data("SmoothingLength", astro = True) if imtype not in ['hi', 'sfr']: temp = readReg.read_data("Temperature", astro = True) hsml = readReg.read_data("SmoothingLength", astro = True) sfr = readReg.read_data("StarFormationRate", astro = True) temp[np.nonzero(sfr > 0)] = 1e4 if ref_galaxy is not None: ids = readReg.read_data("ParticleIDs", astro = False) gate = st.Gate(ids, ref_ids) ref_index = gate.in2() ind_in_ref = np.nonzero(ref_index >= 0)[0] vel = vel[ind_in_ref, :] pos = pos[ind_in_ref, :] mass = mass[ind_in_ref] if ptype == 0: hsml = hsml[ind_in_ref] if imtype not in ['hi', 'sfr']: temp = temp[ind_in_ref] sfr = sfr[ind_in_ref] if ptype == 4 and imtype == 'gri': sft = sft[ind_in_ref] lum_g = lum_g[ind_in_ref] lum_r = lum_r[ind_in_ref] lum_i = lum_i[ind_in_ref] # Select only those particles that are actually within target sphere deltaPos = pos - galpos[None, :] rad = np.linalg.norm(deltaPos, axis = 1) ind_sphere = np.nonzero((rad < readRad*np.sqrt(3)))[0]# & (np.abs(deltaPos[:, 2]) < zsize))[0] #ind_sphere = np.arange(len(mass)) pos = pos[ind_sphere, :] mass = mass[ind_sphere] if ptype == 4 and imtype == 'gri': lum_g = lum_g[ind_sphere] lum_r = lum_r[ind_sphere] lum_i = lum_i[ind_sphere] sft = sft[ind_sphere] if ptype == 0: hsml = hsml[ind_sphere] if imtype not in ['hi', 'sfr']: if gas_tmax is not None: ind_tsel = np.nonzero(temp[ind_sphere] < gas_tmax)[0] quant = temp[ind_sphere[ind_tsel]] mass = mass[ind_tsel] pos = pos[ind_tsel, :] else: quant = temp[ind_sphere] else: quant = mass else: quant = mass if ptype != 0: if fixedSmoothingLength > 0: hsml = np.zeros(len(mass), dtype = np.float32) + fixedSmoothingLength else: hsml = None if xBase is None: xBase = np.zeros(3, dtype = np.float32) yBase = np.copy(xBase) zBase = np.copy(xBase) # Generate actual image if ptype == 4 and imtype == 'gri': image_weight_all, image_quant, hsml_true = ir.make_sph_image_new_3d( pos, mass, quant, hsml, DesNgb=desNGB, imsize=numPix, zpix = 1, boxsize = imsize, CamPos = galpos, CamDir = camDir, ProjectionPlane = projectionPlane, verbose=True, CamAngle = [0,0,rho], rollMode = 0, edge_on = edge_on, treeAllocFac = 10, xBase=xBase, yBase=yBase, zBase=zBase, return_hsml = True) if scale_hsml_with_age: hsml = (hsml_true/1e-3)**((1.7-sft)) * 1e-3 #ind_young = np.nonzero(sft > 0.9)[0] #hsml[ind_young] = (hsml_true[ind_young]*1e3)**(1.0) * 1e3 hsml = np.clip(hsml, None, np.clip(0.1*(1-sft)**(0.5), 0.003, None)) #hsml[ind_young] = np.clip(hsml[ind_young], None, 0.005) else: hsml = hsml_true image_weight_all_g, image_quant = ir.make_sph_image_new_3d( pos, lum_g, lum_g, hsml, DesNgb=desNGB, imsize=numPix, zpix = 1, boxsize = imsize, CamPos = galpos, CamDir = camDir, ProjectionPlane = projectionPlane, verbose=True, CamAngle = [0,0,rho], rollMode = 0, edge_on = edge_on, treeAllocFac = 10, xBase=xBase, yBase=yBase, zBase=zBase, return_hsml = False) image_weight_all_r, image_quant = ir.make_sph_image_new_3d( pos, lum_r, lum_r, hsml, DesNgb=desNGB, imsize=numPix, zpix = 1, boxsize = imsize, CamPos = galpos, CamDir = camDir, ProjectionPlane = projectionPlane, verbose=True, CamAngle = [0,0,rho], rollMode = 0, edge_on = edge_on, treeAllocFac = 10, xBase=xBase, yBase=yBase, zBase=zBase, return_hsml = False) image_weight_all_i, image_quant = ir.make_sph_image_new_3d( pos, lum_i, lum_i, hsml, DesNgb=desNGB, imsize=numPix, zpix = 1, boxsize = imsize, CamPos = galpos, CamDir = camDir, ProjectionPlane = projectionPlane, verbose=True, CamAngle = [0,0,rho], rollMode = 0, edge_on = edge_on, treeAllocFac = 10, xBase=xBase, yBase=yBase, zBase=zBase, return_hsml = False) map_maas_g = -5/2*np.log10( image_weight_all_g[:, :, 1]+1e-5)+5*np.log10(180*3600/np.pi) + 25 map_maas_r = -5/2*np.log10( image_weight_all_r[:, :, 1]+1e-5)+5*np.log10(180*3600/np.pi) + 25 map_maas_i = -5/2*np.log10( image_weight_all_i[:, :, 1]+1e-5)+5*np.log10(180*3600/np.pi) + 25 else: image_weight_all, image_quant = ir.make_sph_image_new_3d( pos, mass, quant, hsml, DesNgb=desNGB, imsize=numPix, zpix = 1, boxsize = imsize, CamPos = galpos, CamDir = camDir, ProjectionPlane = projectionPlane, verbose=True, CamAngle = [0,0,rho], rollMode = 0, edge_on = edge_on, treeAllocFac = 10, xBase=xBase, yBase=yBase, zBase=zBase, zrange = [-zsize, zsize]) if save_maps: maploc = plotloc + str(plot_snap).zfill(4) + '.hdf5' yb.write_hdf5(image_weight_all, maploc, 'Sigma', new = True) if ptype == 4 and imtype == 'gri': yb.write_hdf5(map_maas_g, maploc, 'Sigma_g') yb.write_hdf5(map_maas_r, maploc, 'Sigma_r') yb.write_hdf5(map_maas_i, maploc, 'Sigma_i') if ptype == 0: yb.write_hdf5(image_quant, maploc, 'Temperature') yb.write_hdf5(np.array((-imsize, imsize, -imsize, imsize)), maploc, 'Extent') print("Obtained image, plotting...") # Plot image... fig = plt.figure(figsize = (inch, inch)) if ptype == 4 and imtype == 'gri': vmin = -mag_range[0] + np.array([-0.5, -0.25, 0.0]) vmax = -mag_range[1] + np.array([-0.5, -0.25, 0.0]) clmap_rgb = np.zeros((numPix, numPix, 3)) clmap_rgb[:, :, 2] = np.clip( ((-map_maas_g)-vmin[0])/((vmax[0]-vmin[0])), 0, 1) clmap_rgb[:, :, 1] = np.clip( ((-map_maas_r)-vmin[1])/((vmax[1]-vmin[1])), 0, 1) clmap_rgb[:, :, 0] = np.clip( ((-map_maas_i)-vmin[2])/((vmax[2]-vmin[2])), 0, 1) im = plt.imshow(clmap_rgb, extent = [-imsize, imsize, -imsize, imsize], aspect = 'equal', interpolation = 'nearest', origin = 'lower', alpha = 1.0) else: sigma = np.log10(image_weight_all[:, :, 1]+1e-5)-2 # in M_sun / pc^2 if ptype == 0 and imtype == 'temp': temp = np.log10(image_quant[:, :, 1]) #clmap_rgb = ir.make_double_image(sigma, temp, percSigma = [-2, 3], absSigma = True) clmap_rgb = ir.make_double_image(sigma, temp, percSigma = [0.01, 99.99], absSigma = False, rangeQuant = [4.0, 8.0], cmap = None) im = plt.imshow(clmap_rgb, extent = [-imsize, imsize, -imsize, imsize], aspect = 'equal', interpolation = 'nearest', origin = 'lower', alpha = 1.0) elif (pos.shape[0] > 32): if ptype == 1: sigRange_use = np.percentile(sigma[np.nonzero(sigma >= 1e-4)], [0.01, 99.9]) vmin, vmax = sigRange_use[0], sigRange_use[1]#0, 4.5 #vmin, vmax = np.min(sigma), np.max(sigma) cmap = plt.cm.Greys_r else: if ptype == 4: vmin, vmax = -1, 3.5 cmap = plt.cm.bone elif ptype == 0 and imtype == 'hi': #sigRange_use = np.percentile(sigma[np.nonzero(sigma >= 1e-4)], [0.01, 99.9]) vmin, vmax = -1, 1#sigRange_use[0], sigRange_use[1] cmap = plt.cm.bone elif ptype == 0 and imtype == 'sfr': vmin, vmax = -2, 3 cmap = plt.cm.viridis else: vmin, vmax = -3, 1 cmap = plt.cm.inferno im = plt.imshow(sigma, cmap = cmap, origin = 'lower', extent = [-imsize, imsize, -imsize, imsize], vmin = vmin, vmax = vmax, interpolation = 'none') else: plt.scatter(deltaPos[:, 0], deltaPos[:, 1], color = 'white') if show == 'MIT': ax_min = ax[0, :]#0, :] ax_int = ax[1, :] ax_maj = ax[2, :] print("|min|={:.2f}, |int|={:.2f}, |maj|={:.2f}" .format(np.linalg.norm(ax_min), np.linalg.norm(ax_int), np.linalg.norm(ax_maj))) print("ab={:.2f}, ac={:.2f}, bc={:.2f}" .format(np.dot(ax_min, ax_int), np.dot(ax_min, ax_maj), np.dot(ax_int, ax_maj))) plt.plot((0, np.sum(ax_min*xBase)*0.01), (0, np.sum(ax_min*yBase)*0.01), color = 'red') plt.plot((0, np.sum(ax_int*xBase)*0.01), (0, np.sum(ax_int*yBase)*0.01), color = 'green') plt.plot((0, np.sum(ax_maj*xBase)*0.01), (0, np.sum(ax_maj*yBase)*0.01), color = 'blue') phivec = np.arange(0, np.pi*2.01, 0.001) plt.plot(0.01*np.cos(phivec), 0.01*np.sin(phivec), color = 'grey', linestyle = ':') plt.text(-0.045/0.05*imsize, -0.04/0.05*imsize, 'a/c = {:.3f}\nb/c = {:.3f}' .format(*axRat), va = 'top', ha = 'left', color = 'white') elif show == 'AngMom': aam = np.linalg.norm(am) plt.plot((0, np.sum(am*xBase)*0.01/aam), (0, np.sum(am*yBase)*0.01/aam), color = 'red') phivec = np.arange(0, np.pi*2.01, 0.001) plt.plot(0.01*np.cos(phivec), 0.01*np.sin(phivec), color = 'grey', linestyle = ':') if mark_galaxies: if mass_crit == 'MstarPeak': mpeak = yb.read_hdf5(fgtloc, 'Full/Mstar') elif mass_crit == 'MsubPeak': mpeak = yb.read_hdf5(fgtloc, 'Full/Msub') elif mass_crit == 'Mstar': mpeak = yb.read_hdf5(fgtloc, 'Mstar')[:, plot_snap] elif mass_crit == 'Msub': mpeak = yb.read_hdf5(fgtloc, 'Msub')[:, plot_snap] else: print("I do not understand mass_crit = '" + mass_crit + "'") set_trace() spiderloc = rundir + '/highlev/SpiderwebTables.hdf5' lastsnap = yb.read_hdf5(spiderloc, 'LastSnap') firstsnap = yb.read_hdf5(spiderloc, 'FirstSnap') satFlag = yb.read_hdf5(fgtloc, 'SatFlag')[:, plot_snap] relPosGal = pos_gal_all - galpos[None, :] xpos = np.sum(relPosGal * xBase[None, :], axis = 1) ypos = np.sum(relPosGal * yBase[None, :], axis = 1) zpos = np.sum(relPosGal * zBase[None, :], axis = 1) ind_in_field = np.nonzero((abs(xpos) < imsize) & (abs(ypos) < imsize) & (abs(zpos) < imsize) & (mpeak >= mass_threshold) & (aexpSnap[firstsnap] <= aexpSnap[plot_snap]))[0] sorter_lastsnap = np.argsort(lastsnap[ind_in_field]) for iigal, igal in enumerate(ind_in_field[sorter_lastsnap]): if aexpSnap[lastsnap[igal]] < aexpSnap[plot_snap]: if plot_dead_galaxies: plotcol = 'red' else: continue else: plotcol = 'limegreen' if satFlag[igal]: marker = 'o' else: marker = 'D' plt.scatter(xpos[igal], ypos[igal], 5, edgecolor = plotcol, facecolor = 'none', alpha = 0.5, marker = marker, linewidth = 0.2) if label_galaxies: plt.text(xpos[igal]+imsize/100, ypos[igal]+imsize/100, str(igal), color = plotcol, va = 'bottom', ha = 'left', alpha = 0.5, fontsize = labelSize) # Some embellishments on the image plt.text(-0.045/0.05*imsize, 0.045/0.05*imsize, 'z = {:.3f}' .format(1/aexp_factor - 1), va = 'center', ha = 'left', color = 'white') if print_galaxy_id: plt.text(0.045/0.05*imsize, 0.045/0.05*imsize, 'galID = {:d}' .format(galID), va = 'center', ha = 'right', color = 'white') ax = plt.gca() #ax.set_xlabel(r'$\Delta x$ [pMpc]') #ax.set_ylabel(r'$\Delta y$ [pMpc]') """Plot length bar in bottom-left corner of specified axis.""" lengths = [10, 5, 2, 1, 0.5, 0.2, 0.1, 0.05, 0.02, 0.01, 0.005, 0.002, 0.001] for ilen in lengths: frac_len = ilen / imsize if frac_len > 0.5: continue if ilen >= 1: label = '{:.0f} pMpc' .format(ilen) else: label = '{:.0f} pkpc' .format(ilen*1000) lower_pt = -imsize*0.045/0.05 upper_pt = lower_pt + ilen cen_pt = lower_pt + ilen/2 plt.plot((lower_pt, upper_pt), (-0.045/0.05*imsize, -0.045/0.05*imsize), color = 'white', linewidth = 2) plt.text(cen_pt, -0.0435/0.05*imsize, label, va = 'bottom', ha = 'center', color = 'white', fontsize = 9) break ax.set_xlim((-imsize, imsize)) ax.set_ylim((-imsize, imsize)) #plt.subplots_adjust(left = 0.13, right = 0.95, bottom = 0.12, top = 0.92) #plt.subplots_adjust(left = 0.1, right = 0.95, bottom = 0.1, top = 0.95) plt.subplots_adjust(left = 0.0, right = 1.0, bottom = 0.0, top = 1.0) plt.savefig(plotloc + str(plot_snap).zfill(4) + '.png', dpi = 1*numPix/inch) plt.close() print("Image located in {:s}" .format(plotloc + '{:04d}.png' .format(plot_snap))) print("Done!")