diff --git a/examples/GEAR/AgoraDisk/README b/examples/GEAR/AgoraDisk/README
new file mode 100644
index 0000000000000000000000000000000000000000..c51d39ac476250eed6d2032053cc037e73eb188e
--- /dev/null
+++ b/examples/GEAR/AgoraDisk/README
@@ -0,0 +1,6 @@
+An analysis script for this test case can be found in the AGORA project repository:
+https://bitbucket.org/mornkr/agora-analysis-script/
+
+A modified version of the script able to read in SWIFT data is available upon request.
+
+(Note to the SWIFT authors: The modified script is stored with the other reference solutions.)
diff --git a/examples/GEAR/AgoraDisk/plotSolution.py b/examples/GEAR/AgoraDisk/plotSolution.py
deleted file mode 100644
index 0e3f4d14f7a7f291e10c3f331816c95f25aeacc5..0000000000000000000000000000000000000000
--- a/examples/GEAR/AgoraDisk/plotSolution.py
+++ /dev/null
@@ -1,6065 +0,0 @@
-#!/usr/bin/env python3
-#######################################################################
-#
-# UNIFIED ANALYSIS SCRIPT FOR DISK SIMULATION FOR THE AGORA PROJECT
-#
-# FOR SCRIPT HISTORY SEE VERSION CONTROL CHANGELOG
-#
-# Note: This script works with yt-3.5.1.
-#
-#######################################################################
-# This script is a copy of the AGORA project (https://bitbucket.org/mornkr/agora-analysis-script/)
-# modified in order to take into account SWIFT
-import matplotlib
-
-matplotlib.use("Agg")
-import sys
-import math
-import copy
-import csv
-import numpy as np
-import matplotlib.colorbar as cb
-import matplotlib.lines as ln
-import yt.utilities.physical_constants as constants
-import matplotlib.pyplot as plt
-import matplotlib.gridspec as gridspec
-from yt.mods import *
-from yt.units.yt_array import YTQuantity
-from mpl_toolkits.axes_grid1 import AxesGrid
-from matplotlib.offsetbox import AnchoredText
-from matplotlib.ticker import MaxNLocator
-from yt.fields.particle_fields import add_volume_weighted_smoothed_field
-from yt.fields.particle_fields import add_nearest_neighbor_field
-from yt.analysis_modules.star_analysis.api import StarFormationRate
-from yt.analysis_modules.halo_analysis.api import *
-from yt.data_objects.particle_filters import add_particle_filter
-from scipy.stats import kde
-from subprocess import call
-
-# mylog.setLevel(1)
-
-from yt.utilities.math_utils import get_cyl_r, get_cyl_theta, get_cyl_z
-
-import sys
-
-draw_density_map = 1 # 1 # 0/1 = OFF/ON
-draw_temperature_map = 1 # 1 # 0/1 = OFF/ON
-draw_cellsize_map = (
- 2
-) # 2 # 0/1/2/3 = OFF/ON/ON now with resolution map where particle_size is defined as [M/rho]^(1/3) for SPH/ON with both cell_size and resolution maps
-draw_elevation_map = 1 # 1 # 0/1 = OFF/ON
-draw_metal_map = (
- 0
-) # 0 # 0/1/2 = OFF/ON/ON with total metal mass in the simulation annotated (this will be inaccurate for SPH)
-draw_zvel_map = 0 # 0 # 0/1 = OFF/ON
-draw_star_map = 1 # 1 # 0/1 = OFF/ON
-draw_star_clump_stats = (
- 1
-) # 1 # 0/1/2/3 = OFF/ON/ON with additional star_map with annotated clumps/ON with addtional star_map + extra clump mass function with GIZMO-ps2
-draw_SFR_map = (
- 2
-) ###2 # 0/1/2 = OFF/ON/ON with local K-S plot using patches
-draw_PDF = (
- 2
-) ###2 # 0/1/2/3 = OFF/ON/ON with constant pressure/entropy lines/ON with additional annotations such as 1D profile from a specific code (CHANGA)
-draw_pos_vel_PDF = (
- 4
-) ##4 # 0/1/2/3/4 = OFF/ON/ON with 1D profile/ON also with 1D dispersion profile/ON also with separate 1D vertical dispersion profiles
-draw_star_pos_vel_PDF = (
- 4
-) ##4 # 0/1/2/3/4 = OFF/ON/ON with 1D profile/ON also with 1D dispersion profile/ON also with separate 1D vertical dispersion profiles
-draw_rad_height_PDF = (
- 2
-) ##2 # 0/1/2/3 = OFF/ON/ON with 1D profile/ON with analytic ftn subtracted
-draw_metal_PDF = 1 ##1 # 0/1 = OFF/ON
-draw_density_DF = (
- 2
-) # 2 # 0/1/2 = OFF/ON/ON with difference plot between 1st and 2nd datasets (when 2, dataset_num should be set to 2)
-draw_radius_DF = 1 # 1 # 0/1 = OFF/ON
-draw_star_radius_DF = (
- 2
-) # 2 # 0/1/2 = OFF/ON/ON with SFR profile and K-S plot (when 2, this automatically turns on draw_radius_DF)
-draw_height_DF = 1 # 1 # 0/1 = OFF/ON
-draw_SFR = 1 # 1 # 0/1/2 = OFF/ON/ON with extra line with GIZMO-ps2
-draw_cut_through = 0 # 0 # 0/1 = OFF/ON
-add_nametag = 1 # 1 # 0/1 = OFF/ON
-add_mean_fractional_dispersion = 1 # 1 # 0/1 = OFF/ON
-
-dataset_num = (
- 2
-) # 1/2 = 1st dataset(Grackle+noSF)/2nd dataset(Grackle+SF+ThermalFbck) for AGORA Paper 4
-yt_version_pre_3_2_3 = (
- 0
-) # 0/1 = NO/YES to "Is the yt version being used pre yt-dev-3.2.3?"
-times = [0, 500] # in Myr
-figure_width = 30 # in kpc
-n_ref = 64 # 8; for SPH codes
-over_refine_factor = 1 # 2; for SPH codes
-aperture_size_SFR_map = (
- 750
-) # in pc = Used if draw_SFR_map = 1, 750 matches Bigiel et al. (2008)
-young_star_cutoff_SFR_map = 20 # in Myr = Used if draw_SFR_map = 1
-young_star_cutoff_star_radius_DF = 20 # in Myr = Used if draw_star_radius_DF = 1
-mean_dispersion_radius_range = [
- 2,
- 10,
-] # in kpc = Used if add_mean_fractional_dispersion = 1, for draw_pos_vel_PDF, draw_star_pos_vel_PDF, draw_rad_height_PDF, draw_radius_DF, draw_star_radius_DF
-mean_dispersion_height_range = [
- 0,
- 0.6,
-] # in kpc = Used if add_mean_fractional_dispersion = 1, for draw_height_DF
-mean_dispersion_time_range = [
- 50,
- 500,
-] # in Myr = Used if add_mean_fractional_dispersion = 1, for draw_SFR
-mean_dispersion_density_range = [
- 1e-25,
- 1e-22,
-] # in g/cm3 = Used if add_mean_fractional_dispersion = 1, for draw_density_DF
-disk_normal_vector = [0.0, 0.0, 1.0]
-
-gadget_default_unit_base = {
- "UnitLength_in_cm": 3.08567758e21,
- "UnitMass_in_g": 1.98848e43,
- "UnitVelocity_in_cm_per_s": 100000,
-}
-color_names = [
- "red",
- "magenta",
- "orange",
- "gold",
- "green",
- "cyan",
- "blue",
- "blueviolet",
- "black",
-]
-linestyle_names = ["-"]
-marker_names = ["s", "o", "p", "v", "^", "<", ">", "h", "*"]
-
-# file_location = ['/lustre/ki/pfs/mornkr/080112_CHaRGe/pfs-hyades/AGORA-DISK-repository-for-use/Grackle+noSF/',
-# '/lustre/ki/pfs/mornkr/080112_CHaRGe/pfs-hyades/AGORA-DISK-repository-for-use/Grackle+SF+ThermalFbck/']
-# file_location = ['/global/project/projectdirs/agora/AGORA-DISK-repository-for-use/Grackle+noSF/',
-# '/global/project/projectdirs/agora/AGORA-DISK-repository-for-use/Grackle+SF+ThermalFbck/']
-# codes = ['ART-I', 'ART-II', 'ENZO', 'RAMSES', 'CHANGA', 'GASOLINE', 'GADGET-3', 'GEAR', 'GIZMO']
-# filenames = [[[file_location[0]+'ART-I/IC/AGORA_Galaxy_LOW.d', file_location[0]+'ART-I/t0.5Gyr/10MpcBox_csf512_02350.d'],
-# [file_location[0]+'ART-II/noSF_def_2p/OUT/AGORA_LOW_000000.art', file_location[0]+'ART-II/noSF_def_2p/OUT/AGORA_LOW_000098.art'],
-# [file_location[0]+'ENZO/DD0000/DD0000', file_location[0]+'ENZO/DD0100/DD0100'],
-# [file_location[0]+'RAMSES/output_00001/info_00001.txt', file_location[0]+'RAMSES/output_00068/info_00068.txt'],
-# [file_location[0]+'CHANGA/disklow/disklow.000000', file_location[0]+'CHANGA/disklow/disklow.000500'],
-# [file_location[0]+'GASOLINE/LOW_dataset1.00001', file_location[0]+'GASOLINE/LOW_dataset1.00335'],
-# [file_location[0]+'GADGET-3/AGORA_ISO_LOW_DRY/snap_iso_dry_000.hdf5', file_location[0]+'GADGET-3/AGORA_ISO_LOW_DRY/snap_iso_dry_010.hdf5'],
-# [file_location[0]+'GEAR/snapshot_0000', file_location[0]+'GEAR/snapshot_0500'],
-# [file_location[0]+'GIZMO/snapshot_temp_000', file_location[0]+'GIZMO/snapshot_temp_100']],
-# [[file_location[1]+'ART-I/IC/AGORA_Galaxy_LOW.d', file_location[1]+'ART-I/t0.5Gyr/10MpcBox_csf512_02350.d'],
-# [file_location[1]+'ART-II/SF_FBth_def_2p/OUT/AGORA_LOW_000000.art', file_location[1]+'ART-II/SF_FBth_def_2p/OUT/AGORA_LOW_000311.art'],
-# [file_location[1]+'ENZO/DD0000/DD0000', file_location[1]+'ENZO/DD0050/DD0050'],
-# [file_location[1]+'RAMSES/output_00001/info_00001.txt', file_location[1]+'RAMSES/output_00216/info_00216.txt'],
-# [file_location[1]+'CHANGA/disklow/disklow.000000', file_location[1]+'CHANGA/disklow/disklow.000500'],
-# [file_location[1]+'GASOLINE/LOW_dataset1.00001', file_location[1]+'GASOLINE/LOW_dataset2.00335'],
-# [file_location[1]+'GADGET-3/AGORA_ISO_LOW_SF_SNII_Thermal_Chevalier_SFT10/snap_iso_sf_000.hdf5', file_location[1]+'GADGET-3/AGORA_ISO_LOW_SF_SNII_Thermal_Chevalier_SFT10/snap_iso_sf_010.hdf5'],
-# [file_location[1]+'GEAR/snapshot_0000', file_location[1]+'GEAR/snapshot_0500'],
-# [file_location[1]+'GIZMO/snapshot_temp_000', file_location[1]+'GIZMO/snapshot_temp_100']]]
-codes = ["SWIFT", "GEAR"]
-filenames = [
- [
- ["./agora_disk_IC.hdf5", "./agora_disk_500Myr.hdf5"], # Sim-noSFF
- ["./snapshot_0000.hdf5", "./snapshot_0050.hdf5"],
- ],
- [
- [
- "./agora_disk_IC.hdf5",
- "./agora_disk_500Myr.hdf5",
- ], # Sim-SFF (with star formation and feedback)
- ["./snapshot_0000.hdf5", "./snapshot_0050.hdf5"],
- ],
-] # I did not check the order, they can be switched
-
-# codes = ["SWIFT"]
-# filenames = [[["./agora_disk_0000.hdf5", "./agora_disk_0050.hdf5"]],
-# [["./agora_disk_0000.hdf5", "./agora_disk_0050.hdf5"]]]
-# codes = ['ART-I']
-# filenames = [[[file_location[0]+'ART-I/IC/AGORA_Galaxy_LOW.d', file_location[0]+'ART-I/t0.5Gyr/10MpcBox_csf512_02350.d']],
-# [[file_location[1]+'ART-I/IC/AGORA_Galaxy_LOW.d', file_location[1]+'ART-I/t0.5Gyr/10MpcBox_csf512_02350.d']]]
-# codes = ['ART-II']
-# filenames = [[[file_location[0]+'ART-II/noSF_def_2p/OUT/AGORA_LOW_000000.art', file_location[0]+'ART-II/noSF_def_2p/OUT/AGORA_LOW_000098.art']],
-# [[file_location[1]+'ART-II/SF_FBth_def_2p/OUT/AGORA_LOW_000000.art', file_location[1]+'ART-II/SF_FBth_def_2p/OUT/AGORA_LOW_000311.art']]]
-# codes = ['ENZO']
-# filenames = [[[file_location[0]+'ENZO/DD0000/DD0000', file_location[0]+'ENZO/DD0100/DD0100']],
-# [[file_location[1]+'ENZO/DD0000/DD0000', file_location[1]+'ENZO/DD0050/DD0050']]]
-# codes = ['RAMSES']
-# filenames = [[[file_location[0]+'RAMSES/output_00001/info_00001.txt', file_location[0]+'RAMSES/output_00068/info_00068.txt']],
-# [[file_location[1]+'RAMSES/output_00001/info_00001.txt', file_location[1]+'RAMSES/output_00216/info_00216.txt']]]
-# codes = ['CHANGA']
-# filenames = [[[file_location[0]+'CHANGA/disklow/disklow.000000', file_location[0]+'CHANGA/disklow/disklow.000500']],
-# [[file_location[1]+'CHANGA/disklow/disklow.000000', file_location[1]+'CHANGA/disklow/disklow.000500']]]
-# codes = ['GASOLINE']
-# filenames = [[[file_location[0]+'GASOLINE/LOW_dataset1.00001', file_location[0]+'GASOLINE/LOW_dataset1.00335']],
-# [[file_location[1]+'GASOLINE/LOW_dataset1.00001', file_location[1]+'GASOLINE/LOW_dataset2.00335']]]
-# codes = ['GADGET-3']
-# filenames = [[[file_location[0]+'GADGET-3/AGORA_ISO_LOW_DRY/snap_iso_dry_000.hdf5', file_location[0]+'GADGET-3/AGORA_ISO_LOW_DRY/snap_iso_dry_010.hdf5']],
-# [[file_location[1]+'GADGET-3/AGORA_ISO_LOW_SF_SNII_Thermal_Chevalier_SFT10/snap_iso_sf_000.hdf5', file_location[1]+'GADGET-3/AGORA_ISO_LOW_SF_SNII_Thermal_Chevalier_SFT10/snap_iso_sf_010.hdf5']]]
-# codes = ['GEAR']
-# filenames = [[['snapshot_0000.hdf5', 'snapshot_0500.hdf5']],
-# [['snapshot_0000.hdf5', 'snapshot_0500']]]
-# codes = ['GIZMO']
-# filenames = [[[file_location[0]+'GIZMO/snapshot_temp_000', file_location[0]+'GIZMO/snapshot_temp_100']],
-# [[file_location[1]+'GIZMO/snapshot_temp_000', file_location[1]+'GIZMO/snapshot_temp_100']]]
-
-
-def load_dataset(dataset_num, time, code, codes, filenames_entry):
- if (
- codes[code] == "ART-I"
- ): # ART frontend doesn't find the accompanying files, so we specify them; see http://yt-project.org/docs/dev/examining/loading_data.html#art-data
- dirnames = filenames_entry[code][time][
- : filenames_entry[code][time].rfind("/") + 1
- ]
- if time == 0:
- timestamp = ""
- else:
- timestamp = filenames_entry[code][time][
- filenames_entry[code][time]
- .rfind("_") : filenames_entry[code][time]
- .rfind(".")
- ]
- pf = load(
- filenames_entry[code][time],
- file_particle_header=dirnames + "PMcrd" + timestamp + ".DAT",
- file_particle_data=dirnames + "PMcrs0" + timestamp + ".DAT",
- file_particle_stars=dirnames + "stars" + timestamp + ".dat",
- )
- elif (
- codes[code] == "CHANGA" or codes[code] == "GASOLINE"
- ): # For TIPSY frontend, always make sure to place your parameter file in the same directory as your datasets
- pf = load(
- filenames_entry[code][time],
- n_ref=n_ref,
- over_refine_factor=over_refine_factor,
- )
- elif (
- codes[code] == "GADGET-3"
- ): # For GADGET-3 2nd dataset, there somehow exist particles very far from the center; so we use [-2000, 2000] for a bounding_box
- pf = load(
- filenames_entry[code][time],
- unit_base=gadget_default_unit_base,
- bounding_box=[[-2000.0, 2000.0], [-2000.0, 2000.0], [-2000.0, 2000.0]],
- n_ref=n_ref,
- over_refine_factor=over_refine_factor,
- )
- elif codes[code] == "SWIFT":
- pf = load(
- filenames_entry[code][time],
- unit_base=gadget_default_unit_base,
- bounding_box=[[-2000.0, 2000.0], [-2000.0, 2000.0], [-2000.0, 2000.0]],
- n_ref=n_ref,
- over_refine_factor=over_refine_factor,
- )
- elif codes[code] == "GEAR":
- pf = load(
- filenames_entry[code][time],
- unit_base=gadget_default_unit_base,
- bounding_box=[[-2000.0, 2000.0], [-2000.0, 2000.0], [-2000.0, 2000.0]],
- n_ref=n_ref,
- over_refine_factor=over_refine_factor,
- )
- # from yt.frontends.gadget.definitions import gadget_header_specs
- # from yt.frontends.gadget.definitions import gadget_ptype_specs
- # from yt.frontends.gadget.definitions import gadget_field_specs
- # if with_cooling:
- # gadget_header_specs["chemistry"] = (('h1', 4, 'c'),('h2', 4, 'c'),('empty', 256, 'c'),)
- # header_spec = "default+chemistry"
- # else:
- # header_spec = "default"
- # gear_ptype_specs = ("Gas", "Stars", "Halo", "Disk", "Bulge", "Bndry")
- # if dataset_num == 1:
- # pf = GadgetDataset(filenames_entry[code][time], unit_base = gadget_default_unit_base, bounding_box=[[-1000.0, 1000.0], [-1000.0, 1000.0], [-1000.0, 1000.0]], header_spec=header_spec,
- # ptype_spec=gear_ptype_specs, n_ref=n_ref, over_refine_factor=over_refine_factor)
- # elif dataset_num == 2:
- # # For GEAR 2nd dataset (Grackle+SF+ThermalFbck), "Metals" is acutally 10-species field; check how metallicity field is added below.
- # if with_cooling:
- # agora_gear = ( "Coordinates",
- # "Velocities",
- # "ParticleIDs",
- # "Mass",
- # ("InternalEnergy", "Gas"),
- # ("Density", "Gas"),
- # ("SmoothingLength", "Gas"),
- # ("Metals", "Gas"),
- # ("StellarFormationTime", "Stars"),
- # ("StellarInitMass", "Stars"),
- # ("StellarIDs", "Stars"),
- # ("StellarDensity", "Stars"),
- # ("StellarSmoothingLength", "Stars"),
- # ("StellarMetals", "Stars"),
- # ("Opt1", "Stars"),
- # ("Opt2", "Stars"),
- # )
- # else:
- # agora_gear = ( "Coordinates",
- # "Velocities",
- # "ParticleIDs",
- # "Mass",
- # ("InternalEnergy", "Gas"),
- # ("Density", "Gas"),
- # ("SmoothingLength", "Gas"),
- # )
- # gadget_field_specs["agora_gear"] = agora_gear
- # pf = GadgetDataset(filenames_entry[code][time], unit_base = gadget_default_unit_base, bounding_box=[[-1000.0, 1000.0], [-1000.0, 1000.0], [-1000.0, 1000.0]], header_spec=header_spec,
- # ptype_spec=gear_ptype_specs, field_spec="agora_gear", n_ref=n_ref, over_refine_factor=over_refine_factor)
- elif codes[code] == "GIZMO":
- from yt.frontends.gadget.definitions import gadget_field_specs
-
- if dataset_num == 1:
- agora_gizmo = (
- "Coordinates",
- "Velocities",
- "ParticleIDs",
- "Mass",
- ("Temperature", "Gas"),
- ("Density", "Gas"),
- ("Electron_Number_Density", "Gas"),
- ("HI_NumberDensity", "Gas"),
- ("SmoothingLength", "Gas"),
- )
- elif dataset_num == 2:
- agora_gizmo = (
- "Coordinates",
- "Velocities",
- "ParticleIDs",
- "Mass",
- ("Temperature", "Gas"),
- ("Density", "Gas"),
- ("ElectronAbundance", "Gas"),
- ("NeutralHydrogenAbundance", "Gas"),
- ("SmoothingLength", "Gas"),
- ("StarFormationRate", "Gas"),
- ("StellarFormationTime", "Stars"),
- ("Metallicity", ("Gas", "Stars")),
- ("DelayTime", "Stars"),
- ("StellarInitMass", "Stars"),
- )
- gadget_field_specs["agora_gizmo"] = agora_gizmo
- pf = load(
- filenames_entry[code][time],
- unit_base=gadget_default_unit_base,
- bounding_box=[[-1000.0, 1000.0], [-1000.0, 1000.0], [-1000.0, 1000.0]],
- field_spec="agora_gizmo",
- n_ref=n_ref,
- over_refine_factor=over_refine_factor,
- )
- else:
- pf = load(filenames_entry[code][time])
- # Bug with some fields, need to redefine them in order to get correct units
-
- def _radius(field, data):
- """The cylindrical radius component of the particle positions
-
- Relative to the coordinate system defined by the *normal* vector
- and *center* field parameters.
- """
- normal = data.get_field_parameter("normal")
- pos = data["relative_particle_position"].T
- return data.ds.arr(get_cyl_r(pos, normal), "cm")
-
- pf.add_field(
- ("PartType0", "radius"),
- sampling_type="particle",
- function=_radius,
- units="cm",
- validators=[ValidateParameter("normal"), ValidateParameter("center")],
- )
- pf.add_field(
- ("PartType4", "radius"),
- sampling_type="particle",
- function=_radius,
- units="cm",
- validators=[ValidateParameter("normal"), ValidateParameter("center")],
- )
- pf.add_field(
- ("PartType1", "radius"),
- sampling_type="particle",
- function=_radius,
- units="cm",
- validators=[ValidateParameter("normal"), ValidateParameter("center")],
- )
- pf.add_field(
- ("all", "radius"),
- sampling_type="particle",
- function=_radius,
- units="cm",
- validators=[ValidateParameter("normal"), ValidateParameter("center")],
- )
-
- def _height(field, data):
- """The cylindrical radius component of the particle positions
-
- Relative to the coordinate system defined by the *normal* vector
- and *center* field parameters.
- """
- normal = data.get_field_parameter("normal")
- pos = data["relative_particle_position"].T
- return data.ds.arr(get_cyl_z(pos, normal), "cm")
-
- pf.add_field(
- ("PartType0", "height"),
- sampling_type="particle",
- function=_height,
- units="cm",
- validators=[ValidateParameter("normal"), ValidateParameter("center")],
- )
- pf.add_field(
- ("PartType1", "height"),
- sampling_type="particle",
- function=_height,
- units="cm",
- validators=[ValidateParameter("normal"), ValidateParameter("center")],
- )
- pf.add_field(
- ("PartType4", "height"),
- sampling_type="particle",
- function=_height,
- units="cm",
- validators=[ValidateParameter("normal"), ValidateParameter("center")],
- )
- pf.add_field(
- ("all", "height"),
- sampling_type="particle",
- function=_height,
- units="cm",
- validators=[ValidateParameter("normal"), ValidateParameter("center")],
- )
-
- return pf
-
-
-fig_density_map = []
-fig_temperature_map = []
-fig_cellsize_map = []
-fig_cellsize_map_2 = []
-fig_elevation_map = []
-fig_metal_map = []
-fig_zvel_map = []
-fig_star_map = []
-fig_star_map_2 = []
-fig_star_map_3 = []
-fig_degr_sfr_map = []
-fig_degr_density_map = []
-fig_PDF = []
-fig_pos_vel_PDF = []
-fig_star_pos_vel_PDF = []
-fig_rad_height_PDF = []
-fig_metal_PDF = []
-grid_density_map = []
-grid_temperature_map = []
-grid_cellsize_map = []
-grid_cellsize_map_2 = []
-grid_elevation_map = []
-grid_metal_map = []
-grid_zvel_map = []
-grid_star_map = []
-grid_star_map_2 = []
-grid_star_map_3 = []
-grid_degr_sfr_map = []
-grid_degr_density_map = []
-grid_PDF = []
-grid_pos_vel_PDF = []
-grid_star_pos_vel_PDF = []
-grid_rad_height_PDF = []
-grid_metal_PDF = []
-star_clump_masses_hop = []
-star_clump_masses_fof = []
-star_clump_masses_hop_ref = []
-star_clump_masses_fof_ref = []
-pos_vel_xs = []
-pos_vel_profiles = []
-pos_disp_xs = []
-pos_disp_profiles = []
-pos_disp_vert_xs = []
-pos_disp_vert_profiles = []
-star_pos_vel_xs = []
-star_pos_vel_profiles = []
-star_pos_disp_xs = []
-star_pos_disp_profiles = []
-star_pos_disp_vert_xs = []
-star_pos_disp_vert_profiles = []
-rad_height_xs = []
-rad_height_profiles = []
-density_DF_xs = []
-density_DF_profiles = []
-density_DF_1st_xs = []
-density_DF_1st_profiles = []
-radius_DF_xs = []
-radius_DF_profiles = []
-surface_density = []
-star_radius_DF_xs = []
-star_radius_DF_profiles = []
-star_surface_density = []
-sfr_radius_DF_xs = []
-sfr_radius_DF_profiles = []
-sfr_surface_density = []
-height_DF_xs = []
-height_DF_profiles = []
-height_surface_density = []
-sfr_ts = []
-sfr_cum_masses = []
-sfr_sfrs = []
-surf_dens_SFR_map = []
-sfr_surf_dens_SFR_map = []
-cut_through_zs = []
-cut_through_zvalues = []
-cut_through_xs = []
-cut_through_xvalues = []
-
-for time in range(len(times)):
- if draw_density_map == 1:
- fig_density_map += [plt.figure(figsize=(100, 20))]
- grid_density_map += [
- AxesGrid(
- fig_density_map[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(2, len(codes)),
- axes_pad=0.02,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="2%",
- cbar_pad=0.02,
- )
- ]
- if draw_temperature_map == 1:
- fig_temperature_map += [plt.figure(figsize=(100, 20))]
- grid_temperature_map += [
- AxesGrid(
- fig_temperature_map[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(2, len(codes)),
- axes_pad=0.02,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="2%",
- cbar_pad=0.02,
- )
- ]
- if draw_cellsize_map >= 1:
- fig_cellsize_map += [plt.figure(figsize=(100, 20))]
- grid_cellsize_map += [
- AxesGrid(
- fig_cellsize_map[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(2, len(codes)),
- axes_pad=0.02,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="2%",
- cbar_pad=0.02,
- )
- ]
- fig_cellsize_map_2 += [plt.figure(figsize=(100, 20))]
- grid_cellsize_map_2 += [
- AxesGrid(
- fig_cellsize_map_2[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(2, len(codes)),
- axes_pad=0.02,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="2%",
- cbar_pad=0.02,
- )
- ]
- if draw_elevation_map == 1:
- fig_elevation_map += [plt.figure(figsize=(100, 20))]
- grid_elevation_map += [
- AxesGrid(
- fig_elevation_map[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(1, len(codes)),
- axes_pad=0.02,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="6%",
- cbar_pad=0.02,
- )
- ]
- if draw_metal_map >= 1:
- fig_metal_map += [plt.figure(figsize=(100, 20))]
- grid_metal_map += [
- AxesGrid(
- fig_metal_map[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(2, len(codes)),
- axes_pad=0.02,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="2%",
- cbar_pad=0.02,
- )
- ]
- if draw_zvel_map == 1:
- fig_zvel_map += [plt.figure(figsize=(100, 20))]
- grid_zvel_map += [
- AxesGrid(
- fig_zvel_map[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(2, len(codes)),
- axes_pad=0.02,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="2%",
- cbar_pad=0.02,
- )
- ]
- if draw_star_map == 1:
- fig_star_map += [plt.figure(figsize=(100, 20))]
- grid_star_map += [
- AxesGrid(
- fig_star_map[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(2, len(codes)),
- axes_pad=0.02,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="2%",
- cbar_pad=0.02,
- )
- ]
- if draw_star_clump_stats >= 1:
- star_clump_masses_hop.append([])
- star_clump_masses_fof.append([])
- fig_star_map_2 += [plt.figure(figsize=(100, 20))]
- grid_star_map_2 += [
- AxesGrid(
- fig_star_map_2[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(2, len(codes)),
- axes_pad=0.02,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="2%",
- cbar_pad=0.02,
- )
- ]
- fig_star_map_3 += [plt.figure(figsize=(100, 20))]
- grid_star_map_3 += [
- AxesGrid(
- fig_star_map_3[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(2, len(codes)),
- axes_pad=0.02,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="2%",
- cbar_pad=0.02,
- )
- ]
- if draw_SFR_map >= 1:
- fig_degr_density_map += [plt.figure(figsize=(100, 20))]
- grid_degr_density_map += [
- AxesGrid(
- fig_degr_density_map[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(1, len(codes)),
- axes_pad=0.02,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="6%",
- cbar_pad=0.02,
- )
- ]
- fig_degr_sfr_map += [plt.figure(figsize=(100, 20))]
- grid_degr_sfr_map += [
- AxesGrid(
- fig_degr_sfr_map[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(1, len(codes)),
- axes_pad=0.02,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="6%",
- cbar_pad=0.02,
- )
- ]
- surf_dens_SFR_map.append([])
- sfr_surf_dens_SFR_map.append([])
- if draw_SFR_map >= 2 or draw_star_radius_DF >= 2:
-
- def import_text(filename, separator):
- for line in csv.reader(
- open(filename), delimiter=separator, skipinitialspace=True
- ):
- if line:
- yield line
-
- if draw_PDF >= 1:
- fig_PDF += [plt.figure(figsize=(50, 80))]
- grid_PDF += [
- AxesGrid(
- fig_PDF[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(2, int(math.ceil(len(codes) / 2.0))),
- axes_pad=0.05,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="2%",
- cbar_pad=0.05,
- aspect=False,
- )
- ]
- if draw_pos_vel_PDF >= 1:
- fig_pos_vel_PDF += [plt.figure(figsize=(50, 80))]
- grid_pos_vel_PDF += [
- AxesGrid(
- fig_pos_vel_PDF[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(3, int(math.ceil(len(codes) / 3.0))),
- axes_pad=0.05,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="2%",
- cbar_pad=0.05,
- aspect=False,
- )
- ]
- pos_vel_xs.append([])
- pos_vel_profiles.append([])
- pos_disp_xs.append([])
- pos_disp_profiles.append([])
- pos_disp_vert_xs.append([])
- pos_disp_vert_profiles.append([])
- if draw_star_pos_vel_PDF >= 1:
- fig_star_pos_vel_PDF += [plt.figure(figsize=(50, 80))]
- grid_star_pos_vel_PDF += [
- AxesGrid(
- fig_star_pos_vel_PDF[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(3, int(math.ceil(len(codes) / 3.0))),
- axes_pad=0.05,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="2%",
- cbar_pad=0.05,
- aspect=False,
- )
- ]
- star_pos_vel_xs.append([])
- star_pos_vel_profiles.append([])
- star_pos_disp_xs.append([])
- star_pos_disp_profiles.append([])
- star_pos_disp_vert_xs.append([])
- star_pos_disp_vert_profiles.append([])
- if draw_rad_height_PDF >= 1:
- fig_rad_height_PDF += [plt.figure(figsize=(50, 80))]
- grid_rad_height_PDF += [
- AxesGrid(
- fig_rad_height_PDF[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(3, int(math.ceil(len(codes) / 3.0))),
- axes_pad=0.05,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="2%",
- cbar_pad=0.05,
- aspect=False,
- )
- ]
- rad_height_xs.append([])
- rad_height_profiles.append([])
- if draw_metal_PDF == 1:
- fig_metal_PDF += [plt.figure(figsize=(50, 80))]
- grid_metal_PDF += [
- AxesGrid(
- fig_metal_PDF[time],
- (0.01, 0.01, 0.99, 0.99),
- nrows_ncols=(2, int(math.ceil(len(codes) / 2.0))),
- axes_pad=0.05,
- add_all=True,
- share_all=True,
- label_mode="1",
- cbar_mode="single",
- cbar_location="right",
- cbar_size="2%",
- cbar_pad=0.05,
- aspect=False,
- )
- ]
- if draw_density_DF >= 1:
- density_DF_xs.append([])
- density_DF_profiles.append([])
- density_DF_1st_xs.append([])
- density_DF_1st_profiles.append([])
- if draw_density_DF == 2 and dataset_num == 1:
- print("This won't work; resetting draw_density_DF to 1...")
- draw_density_DF == 1
- if draw_star_radius_DF >= 1:
- star_radius_DF_xs.append([])
- star_radius_DF_profiles.append([])
- star_surface_density.append([])
- sfr_radius_DF_xs.append([])
- sfr_radius_DF_profiles.append([])
- sfr_surface_density.append([])
- if draw_star_radius_DF == 2:
- draw_radius_DF = 1
- if draw_radius_DF == 1:
- radius_DF_xs.append([])
- radius_DF_profiles.append([])
- surface_density.append([])
- if draw_height_DF == 1:
- height_DF_xs.append([])
- height_DF_profiles.append([])
- height_surface_density.append([])
- if draw_SFR >= 1:
- sfr_ts.append([])
- sfr_cum_masses.append([])
- sfr_sfrs.append([])
- if draw_cut_through == 1:
- cut_through_zs.append([])
- cut_through_zvalues.append([])
- cut_through_xs.append([])
- cut_through_xvalues.append([])
-
-for time in range(len(times)):
-
- for code in range(len(codes)):
-
- if (dataset_num != 1 and dataset_num != 2) or (
- filenames[dataset_num - 1][code] == []
- ):
- continue
-
- ####################################
- # PRE-ANALYSIS STEPS #
- ####################################
-
- # LOAD DATASETS
- pf = load_dataset(dataset_num, time, code, codes, filenames[dataset_num - 1])
-
- # PARTICLE FILED NAMES FOR SPH CODES, AND STELLAR PARTICLE FILTERS
- PartType_Gas_to_use = "Gas"
- PartType_Star_to_use = "Stars"
- PartType_StarBeforeFiltered_to_use = "Stars"
- MassType_to_use = "Mass"
- MetallicityType_to_use = "Metallicity"
- FormationTimeType_to_use = (
- "StellarFormationTime"
- ) # for GADGET/GEAR/GIZMO, this field has to be added in frontends/sph/fields.py, in which only "FormationTime" can be recognized
- SmoothingLengthType_to_use = "SmoothingLength"
-
- if codes[code] == "CHANGA" or codes[code] == "GASOLINE":
- MetallicityType_to_use = "Metals"
- PartType_Star_to_use = "NewStars"
- PartType_StarBeforeFiltered_to_use = "Stars"
- FormationTimeType_to_use = "FormationTime"
- SmoothingLengthType_to_use = "smoothing_length"
-
- def NewStars(
- pfilter, data
- ): # see http://yt-project.org/docs/dev/analyzing/filtering.html#filtering-particle-fields
- return data[(pfilter.filtered_type, FormationTimeType_to_use)] > 0
-
- add_particle_filter(
- PartType_Star_to_use,
- function=NewStars,
- filtered_type=PartType_StarBeforeFiltered_to_use,
- requires=[FormationTimeType_to_use],
- )
- pf.add_particle_filter(PartType_Star_to_use)
- pf.periodicity = (
- True,
- True,
- True,
- ) # this is needed especially when bPeriodic = 0 in GASOLINE, to avoid RuntimeError in geometry/selection_routines.pyx:855
- elif codes[code] == "ART-I":
- PartType_StarBeforeFiltered_to_use = "stars"
- FormationTimeType_to_use = "particle_creation_time"
-
- def Stars(pfilter, data):
- return data[(pfilter.filtered_type, FormationTimeType_to_use)] > 0
-
- # return ((data[(pfilter.filtered_type, FormationTimeType_to_use)] > 0) & (data[(pfilter.filtered_type, "particle_index")] >= 212500)) # without the fix metioned above, the above line won't work because all "stars"="specie1" have the same wrong particle_creation_time of 6.851 Gyr; so you will have to use this quick and dirty patch
- add_particle_filter(
- PartType_Star_to_use,
- function=Stars,
- filtered_type=PartType_StarBeforeFiltered_to_use,
- requires=[FormationTimeType_to_use, "particle_index"],
- )
- pf.add_particle_filter(PartType_Star_to_use)
- elif codes[code] == "ART-II":
- PartType_StarBeforeFiltered_to_use = "STAR"
- if (
- time != 0
- ): # BIRTH_TIME field exists in ART-II but as a dimensionless quantity for some reason in frontends/artio/fields.py; so we create StellarFormationTime field
-
- def _FormationTime(field, data):
- return pf.arr(data["STAR", "BIRTH_TIME"].d, "code_time")
-
- pf.add_field(
- ("STAR", FormationTimeType_to_use),
- function=_FormationTime,
- particle_type=True,
- take_log=False,
- units="code_time",
- )
-
- def Stars(pfilter, data):
- return data[(pfilter.filtered_type, FormationTimeType_to_use)] > 0
-
- add_particle_filter(
- PartType_Star_to_use,
- function=Stars,
- filtered_type=PartType_StarBeforeFiltered_to_use,
- requires=[FormationTimeType_to_use],
- )
- pf.add_particle_filter(PartType_Star_to_use)
- elif codes[code] == "ENZO":
- PartType_StarBeforeFiltered_to_use = "all"
- FormationTimeType_to_use = "creation_time"
-
- def Stars(pfilter, data):
- return (data[(pfilter.filtered_type, "particle_type")] == 2) & (
- data[(pfilter.filtered_type, FormationTimeType_to_use)] > 0
- )
-
- add_particle_filter(
- PartType_Star_to_use,
- function=Stars,
- filtered_type=PartType_StarBeforeFiltered_to_use,
- requires=["particle_type", FormationTimeType_to_use],
- )
- pf.add_particle_filter(PartType_Star_to_use)
- elif codes[code] == "GADGET-3":
- PartType_Gas_to_use = "PartType0"
- PartType_Star_to_use = "PartType4"
- PartType_StarBeforeFiltered_to_use = "PartType4"
- MassType_to_use = "Masses"
- elif codes[code] == "SWIFT":
- PartType_Gas_to_use = "PartType0"
- PartType_Star_to_use = "PartType4"
- PartType_StarBeforeFiltered_to_use = "PartType4"
- if time != 0:
-
- def _FormationTime(field, data):
- return pf.arr(data["PartType4", "BirthTimes"].d, "code_time")
-
- pf.add_field(
- ("PartType4", FormationTimeType_to_use),
- function=_FormationTime,
- particle_type=True,
- take_log=False,
- units="code_time",
- )
-
- MassType_to_use = "Masses"
- elif codes[code] == "GEAR":
- PartType_Gas_to_use = "PartType0"
- PartType_Star_to_use = "PartType1"
- PartType_StarBeforeFiltered_to_use = "PartType1"
- MassType_to_use = "Masses"
- if time != 0:
-
- def _FormationTime(field, data):
- return pf.arr(data["PartType1", "StarFormationTime"].d, "code_time")
-
- pf.add_field(
- ("PartType1", FormationTimeType_to_use),
- function=_FormationTime,
- particle_type=True,
- take_log=False,
- units="code_time",
- )
- elif codes[code] == "RAMSES":
- PartType_StarBeforeFiltered_to_use = "all"
- pf.current_time = pf.arr(
- pf.parameters["time"], "code_time"
- ) # reset pf.current_time because it is incorrectly set up in frontends/ramses/data_structure.py, and I don't wish to mess with units there
- FormationTimeType_to_use = (
- "particle_age"
- ) # particle_age field actually means particle creation time, at least for this particular dataset, so the new field below is not needed
- # if time != 0: # Only particle_age field exists in RAMSES (only for new stars + IC stars), so we create StellarFormationTime field
- # def _FormationTime(field, data):
- # return pf.current_time - data["all", "particle_age"].in_units("s")
- # pf.add_field(("all", FormationTimeType_to_use), function=_FormationTime, particle_type=True, take_log=False, units="code_time")
- def Stars(pfilter, data):
- return data[(pfilter.filtered_type, "particle_age")] > 0
-
- add_particle_filter(
- PartType_Star_to_use,
- function=Stars,
- filtered_type=PartType_StarBeforeFiltered_to_use,
- requires=["particle_age"],
- )
- pf.add_particle_filter(PartType_Star_to_use)
-
- # AXIS SWAP FOR PLOT COLLECTION
- pf.coordinates.x_axis[1] = 0
- pf.coordinates.y_axis[1] = 2
- pf.coordinates.x_axis["y"] = 0
- pf.coordinates.y_axis["y"] = 2
-
- # ADDITIONAL FIELDS I: GLOBALLY USED FIELDS
- def _density_squared(field, data):
- return data[("gas", "density")] ** 2
-
- pf.add_field(
- ("gas", "density_squared"), function=_density_squared, units="g**2/cm**6"
- )
-
- # ADDITIONAL FIELDS II: FOR CELL SIZE AND RESOLUTION MAPS
- def _CellSizepc(field, data):
- return (data[("index", "cell_volume")].in_units("pc**3")) ** (1 / 3.0)
-
- pf.add_field(
- ("index", "cell_size"),
- function=_CellSizepc,
- units="pc",
- display_name="Resolution $\Delta$ x",
- take_log=True,
- )
-
- def _Inv2CellVolumeCode(field, data):
- return data[("index", "cell_volume")] ** -2
-
- pf.add_field(
- ("index", "cell_volume_inv2"),
- function=_Inv2CellVolumeCode,
- units="code_length**(-6)",
- display_name="Inv2CellVolumeCode",
- take_log=True,
- )
- if draw_cellsize_map >= 2:
- if (
- codes[code] == "CHANGA"
- or codes[code] == "GEAR"
- or codes[code] == "GADGET-3"
- or codes[code] == "GASOLINE"
- or codes[code] == "GIZMO"
- or codes[code] == "SWIFT"
- ):
-
- def _ParticleSizepc(field, data):
- return (
- data[(PartType_Gas_to_use, MassType_to_use)]
- / data[(PartType_Gas_to_use, "Density")]
- ) ** (1.0 / 3.0)
-
- pf.add_field(
- (PartType_Gas_to_use, "particle_size"),
- function=_ParticleSizepc,
- units="pc",
- display_name="Resolution $\Delta$ x",
- particle_type=True,
- take_log=True,
- )
-
- def _Inv2ParticleVolumepc(field, data):
- return (
- data[(PartType_Gas_to_use, MassType_to_use)]
- / data[(PartType_Gas_to_use, "Density")]
- ) ** (-2.0)
-
- pf.add_field(
- (PartType_Gas_to_use, "particle_volume_inv2"),
- function=_Inv2ParticleVolumepc,
- units="pc**(-6)",
- display_name="Inv2ParticleVolumepc",
- particle_type=True,
- take_log=True,
- )
- # Also creating smoothed field following an example in yt-project.org/docs/dev/cookbook/calculating_information.html; use hardcoded num_neighbors as in frontends/gadget/fields.py
- fn = add_volume_weighted_smoothed_field(
- PartType_Gas_to_use,
- "Coordinates",
- MassType_to_use,
- SmoothingLengthType_to_use,
- "Density",
- "particle_size",
- pf.field_info,
- nneighbors=64,
- )
- fn = add_volume_weighted_smoothed_field(
- PartType_Gas_to_use,
- "Coordinates",
- MassType_to_use,
- SmoothingLengthType_to_use,
- "Density",
- "particle_volume_inv2",
- pf.field_info,
- nneighbors=64,
- )
-
- # Alias doesn't work -- e.g. pf.field_info.alias(("gas", "particle_size"), fn[0]) -- check alias below; so I simply add ("gas", "particle_size")
- def _ParticleSizepc_2(field, data):
- return data[
- "deposit", PartType_Gas_to_use + "_smoothed_" + "particle_size"
- ]
-
- pf.add_field(
- ("gas", "particle_size"),
- function=_ParticleSizepc_2,
- units="pc",
- force_override=True,
- display_name="Resolution $\Delta$ x",
- particle_type=False,
- take_log=True,
- )
-
- def _Inv2ParticleVolumepc_2(field, data):
- return data[
- "deposit",
- PartType_Gas_to_use + "_smoothed_" + "particle_volume_inv2",
- ]
-
- pf.add_field(
- ("gas", "particle_volume_inv2"),
- function=_Inv2ParticleVolumepc_2,
- units="pc**(-6)",
- force_override=True,
- display_name="Inv2ParticleVolumepc",
- particle_type=False,
- take_log=True,
- )
-
- # ADDITIONAL FIELDS III: TEMPERATURE
- if (
- codes[code] == "GEAR"
- or codes[code] == "GADGET-3"
- or codes[code] == "RAMSES"
- or codes[code] == "SWIFT"
- ):
- # From grackle/src/python/utilities/convenience.py: Calculate a tabulated approximation to mean molecular weight (valid for data that used Grackle 2.0 or below)
- def calc_mu_table_local(temperature):
- tt = np.array(
- [
- 1.0e01,
- 1.0e02,
- 1.0e03,
- 1.0e04,
- 1.3e04,
- 2.1e04,
- 3.4e04,
- 6.3e04,
- 1.0e05,
- 1.0e09,
- ]
- )
- mt = np.array(
- [
- 1.18701555,
- 1.15484424,
- 1.09603514,
- 0.9981496,
- 0.96346395,
- 0.65175895,
- 0.6142901,
- 0.6056833,
- 0.5897776,
- 0.58822635,
- ]
- )
- logttt = np.log(temperature)
- logmu = np.interp(
- logttt, np.log(tt), np.log(mt)
- ) # linear interpolation in log-log space
- return np.exp(logmu)
-
- temperature_values = []
- mu_values = []
- T_over_mu_values = []
- current_temperature = 1e1
- final_temperature = 1e7
- dlogT = 0.1
- while current_temperature < final_temperature:
- temperature_values.append(current_temperature)
- current_mu = calc_mu_table_local(current_temperature)
- mu_values.append(current_mu)
- T_over_mu_values.append(current_temperature / current_mu)
- current_temperature = np.exp(np.log(current_temperature) + dlogT)
-
- def convert_T_over_mu_to_T(T_over_mu):
- logT_over_mu = np.log(T_over_mu)
- logT = np.interp(
- logT_over_mu, np.log(T_over_mu_values), np.log(temperature_values)
- ) # linear interpolation in log-log space
- return np.exp(logT)
-
- if (
- codes[code] == "GEAR"
- or codes[code] == "GADGET-3"
- or codes[code] == "SWIFT"
- ):
-
- def _Temperature_3(field, data):
- gamma = 5.0 / 3.0
- T_over_mu = (
- (
- data[PartType_Gas_to_use, "InternalEnergy"]
- * (gamma - 1)
- * constants.mass_hydrogen_cgs
- / constants.boltzmann_constant_cgs
- )
- .in_units("K")
- .d
- ) # T/mu
- return YTArray(convert_T_over_mu_to_T(T_over_mu), "K") # now T
-
- pf.add_field(
- (PartType_Gas_to_use, "Temperature"),
- function=_Temperature_3,
- particle_type=True,
- force_override=True,
- units="K",
- )
- elif codes[code] == "RAMSES":
- # The pressure field includes the artificial pressure support term, so one needs to be careful (compare with the exsiting frontends/ramses/fields.py)
- def _temperature_3(field, data):
- T_J = 1800.0 # in K
- n_J = 8.0 # in H/cc
- gamma_0 = 2.0
- x_H = 0.76
- mH = (
- 1.66e-24
- ) # from pymses/utils/constants/__init__.py (vs. in yt, mass_hydrogen_cgs = 1.007947*amu_cgs = 1.007947*1.660538921e-24 = 1.6737352e-24)
- kB = (
- 1.3806504e-16
- ) # from pymses/utils/constants/__init__.py (vs. in yt, boltzmann_constant_cgs = 1.3806488e-16)
- if time != 0:
- T_over_mu = data["gas", "pressure"].d / data[
- "gas", "density"
- ].d * mH / kB - T_J * (
- data["gas", "density"].d * x_H / mH / n_J
- ) ** (
- gamma_0 - 1.0
- ) # T/mu = T2 in Ramses
- else:
- T_over_mu = (
- data["gas", "pressure"].d
- / data["gas", "density"].d
- * mH
- / kB
- ) # IC: no pressure support
- return YTArray(convert_T_over_mu_to_T(T_over_mu), "K") # now T
-
- pf.add_field(
- ("gas", "temperature"),
- function=_temperature_3,
- force_override=True,
- units="K",
- )
-
- # ADDITIONAL FIELDS IV: METALLICITY (IN MASS FRACTION, NOT IN ZSUN)
- if draw_metal_map >= 1 or draw_metal_PDF == 1:
- if (
- codes[code] == "ART-I"
- ): # metallicity field in ART-I has a different meaning (see frontends/art/fields.py), and metallicity field in ART-II is missing
-
- def _metallicity_2(field, data):
- return (
- (
- data["gas", "metal_ii_density"]
- + data["gas", "metal_ia_density"]
- )
- / data["gas", "density"]
- ) # metal_ia_density needs to be added to account for initial 0.5 Zsun, even though we don't have SNe Ia
-
- pf.add_field(
- ("gas", "metallicity"),
- function=_metallicity_2,
- force_override=True,
- display_name="Metallicity",
- take_log=True,
- units="",
- )
- elif codes[code] == "ART-II":
-
- def _metallicity_2(field, data):
- return data["gas", "metal_ii_density"] / data["gas", "density"]
-
- pf.add_field(
- ("gas", "metallicity"),
- function=_metallicity_2,
- force_override=True,
- display_name="Metallicity",
- take_log=True,
- units="",
- )
- elif (
- codes[code] == "ENZO"
- ): # metallicity field in ENZO is in Zsun, so we create a new field
-
- def _metallicity_2(field, data):
- return data["gas", "metal_density"] / data["gas", "density"]
-
- pf.add_field(
- ("gas", "metallicity"),
- function=_metallicity_2,
- force_override=True,
- display_name="Metallicity",
- take_log=True,
- units="",
- )
- elif (
- codes[code] == "GEAR"
- ): # "Metals" in GEAR is 10-species field ([:,9] is the total metal fraction), so requires a change in _vector_fields in frontends/gadget/io.py: added ("Metals", 10)
-
- def _metallicity_2(field, data):
- if len(data[PartType_Gas_to_use, "Metals"].shape) == 1:
- return data[PartType_Gas_to_use, "Metals"]
- else:
- return data[PartType_Gas_to_use, "Metals"][:, 9].in_units(
- ""
- ) # in_units("") turned out to be crucial!; otherwise code_metallicity will be used and it will mess things up
-
- # We are creating ("Gas", "Metallicity") here, different from ("Gas", "metallicity") which is auto-generated by yt but doesn't work properly
- pf.add_field(
- (PartType_Gas_to_use, MetallicityType_to_use),
- function=_metallicity_2,
- display_name="Metallicity",
- particle_type=True,
- take_log=True,
- units="",
- )
- # Also creating smoothed field following an example in yt-project.org/docs/dev/cookbook/calculating_information.html; use hardcoded num_neighbors as in frontends/gadget/fields.py
- fn = add_volume_weighted_smoothed_field(
- PartType_Gas_to_use,
- "Coordinates",
- MassType_to_use,
- "SmoothingLength",
- "Density",
- MetallicityType_to_use,
- pf.field_info,
- nneighbors=64,
- )
- # Alias doesn't work -- e.g. pf.field_info.alias(("gas", "metallicity"), fn[0]) -- probably because pf=GadgetDataset()?, not load()?; so I add and replace existing ("gas", "metallicity")
- def _metallicity_3(field, data):
- return data[
- "deposit",
- PartType_Gas_to_use + "_smoothed_" + MetallicityType_to_use,
- ]
-
- pf.add_field(
- ("gas", "metallicity"),
- function=_metallicity_3,
- force_override=True,
- display_name="Metallicity",
- particle_type=False,
- take_log=True,
- units="",
- )
- elif (
- codes[code] == "SWIFT"
- ): # "Metals" in SWIFT is 10-species field ([:,9] is the total metal fraction), so requires a change in _vector_fields in frontends/gadget/io.py: added ("Metals", 10)
-
- def _metallicity_2(field, data):
- if (
- len(
- data[
- PartType_Gas_to_use, "SmoothedMetalMassFractions"
- ].shape
- )
- == 1
- ):
- return data[PartType_Gas_to_use, "SmoothedMetalMassFractions"]
- else:
- return data[PartType_Gas_to_use, "SmoothedMetalMassFractions"][
- :, 9
- ].in_units(
- ""
- ) # in_units("") turned out to be crucial!; otherwise code_metallicity will be used and it will mess things up
-
- # We are creating ("Gas", "Metallicity") here, different from ("Gas", "metallicity") which is auto-generated by yt but doesn't work properly
- pf.add_field(
- (PartType_Gas_to_use, MetallicityType_to_use),
- function=_metallicity_2,
- display_name="Metallicity",
- particle_type=True,
- take_log=True,
- units="",
- )
- # Also creating smoothed field following an example in yt-project.org/docs/dev/cookbook/calculating_information.html; use hardcoded num_neighbors as in frontends/gadget/fields.py
- fn = add_volume_weighted_smoothed_field(
- PartType_Gas_to_use,
- "Coordinates",
- MassType_to_use,
- "SmoothingLength",
- "Density",
- MetallicityType_to_use,
- pf.field_info,
- nneighbors=64,
- )
- # Alias doesn't work -- e.g. pf.field_info.alias(("gas", "metallicity"), fn[0]) -- probably because pf=GadgetDataset()?, not load()?; so I add and replace existing ("gas", "metallicity")
- def _metallicity_3(field, data):
- return data[
- "deposit",
- PartType_Gas_to_use + "_smoothed_" + MetallicityType_to_use,
- ]
-
- pf.add_field(
- ("gas", "metallicity"),
- function=_metallicity_3,
- force_override=True,
- display_name="Metallicity",
- particle_type=False,
- take_log=True,
- units="",
- )
- if draw_metal_map >= 2:
-
- def _metal_mass(field, data):
- return data["gas", "cell_mass"] * data["gas", "metallicity"]
-
- pf.add_field(
- ("gas", "metal_mass"),
- function=_metal_mass,
- force_override=True,
- display_name="Metal Mass",
- take_log=True,
- units="Msun",
- )
-
- # ADDITIONAL FIELDS V: FAKE PARTICLE FIELDS, CYLINDRICAL COORDINATES, etc.
- def rho_agora_disk(r, z):
- r_d = YTArray(3.432, "kpc")
- z_d = 0.1 * r_d
- M_d = YTArray(4.297e10, "Msun")
- f_gas = 0.2
- r_0 = (f_gas * M_d / (4.0 * np.pi * r_d ** 2 * z_d)).in_units("g/cm**3")
- return r_0 * numpy.exp(-r / r_d) * numpy.exp(-z / z_d)
-
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
-
- def _cylindrical_z_abs(field, data):
- return numpy.abs(data[("index", "cylindrical_z")])
-
- pf.add_field(
- ("index", "cylindrical_z_abs"),
- function=_cylindrical_z_abs,
- take_log=False,
- particle_type=False,
- units="cm",
- )
-
- def _density_minus_analytic(field, data):
- return data[("gas", "density")] - rho_agora_disk(
- data[("index", "cylindrical_r")],
- data[("index", "cylindrical_z_abs")],
- )
-
- pf.add_field(
- ("gas", "density_minus_analytic"),
- function=_density_minus_analytic,
- take_log=False,
- particle_type=False,
- display_name="density residual abs",
- units="g/cm**3",
- )
- else:
-
- def _particle_position_cylindrical_z_abs(field, data):
- return numpy.abs(data[(PartType_Gas_to_use, "height")])
-
- pf.add_field(
- (PartType_Gas_to_use, "particle_position_cylindrical_z_abs"),
- function=_particle_position_cylindrical_z_abs,
- take_log=False,
- particle_type=True,
- units="cm",
- )
- # particle_type=False doesn't make sense, but is critical for PhasePlot/ProfilePlot to work for particle fields
- # requires a change in data_objects/data_container.py: remove raise YTFieldTypeNotFound(ftype)
- if draw_metal_map >= 1 or draw_metal_PDF == 1:
-
- def _Metallicity_2(field, data):
- return data[(PartType_Gas_to_use, MetallicityType_to_use)]
-
- pf.add_field(
- (PartType_Gas_to_use, "Metallicity_2"),
- function=_Metallicity_2,
- take_log=True,
- particle_type=True,
- display_name="Metallicity",
- units="",
- )
-
- def _Density_2_minus_analytic(field, data):
- return data[(PartType_Gas_to_use, "density")] - rho_agora_disk(
- data[(PartType_Gas_to_use, "radius")],
- data[(PartType_Gas_to_use, "particle_position_cylindrical_z_abs")],
- )
-
- pf.add_field(
- (PartType_Gas_to_use, "Density_2_minus_analytic"),
- function=_Density_2_minus_analytic,
- take_log=False,
- particle_type=False,
- display_name="Density Residual",
- units="g/cm**3",
- )
-
- ####################################
- # MAIN ANALYSIS ROUTINES #
- ####################################
-
- # FIND CENTER AND PROJ_REGION
- v, cen = pf.h.find_max(
- ("gas", "density")
- ) # find the center to keep the galaxy at the center of all the images; here we assume that the gas disk is no bigger than 30 kpc in radius
- sp = pf.sphere(cen, (30.0, "kpc"))
- cen2 = sp.quantities.center_of_mass(use_gas=True, use_particles=False).in_units(
- "kpc"
- )
- sp2 = pf.sphere(cen2, (1.0, "kpc"))
- cen3 = sp2.quantities.max_location(("gas", "density"))
- center = pf.arr(
- [cen3[1].d, cen3[2].d, cen3[3].d], "cm"
- ) # naive usage such as YTArray([cen3[1], cen3[2], cen3[3]]) doesn't work somehow for ART-II data
- if yt_version_pre_3_2_3 == 1:
- center = pf.arr(
- [cen3[2].d, cen3[3].d, cen3[4].d], "code_length"
- ) # for yt-3.2.3 or before
- if codes[code] == "GASOLINE" and dataset_num == 2 and time == 1:
- # center = pf.arr([2.7912903206399826e+20, 1.5205303149849894e+21, 1.5398968883245956e+21], 'cm') # a temporary hack for GASOLINE (center of the most massive stellar clump)
- center = pf.arr([0.09045956, 0.49277032, 0.49904659], "code_length")
- # if codes[code] == "GADGET-3" and dataset_num == 2 and time == 1:
- # #center = pf.arr([6.184520935812296e+21, 4.972678132728082e+21, 6.559067311284074e+21], 'cm') # a temporary hack for GADGET-3 (center of the most massive stellar clump)
- # center = pf.arr([ 2.00426673674, 1.61153523084, 2.12564895041], 'code_length')
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
- proj_region = pf.box(
- center - YTArray([figure_width, figure_width, figure_width], "kpc"),
- center + YTArray([figure_width, figure_width, figure_width], "kpc"),
- ) # projected images made using a (2*figure_width)^3 box for AMR codes
- else:
- proj_region = pf.all_data()
-
- # DENSITY MAPS
- if draw_density_map == 1:
- my_cmap = copy.copy(matplotlib.cm.get_cmap("algae"))
- my_cmap.set_bad(my_cmap(0)) # cmap range [0, 1)
- my_cmap.set_under(my_cmap(0))
- for ax in range(1, 3):
- p = ProjectionPlot(
- pf,
- ax,
- ("gas", "density"),
- center=center,
- data_source=proj_region,
- width=(figure_width, "kpc"),
- weight_field=None,
- fontsize=9,
- )
- p.set_zlim(("gas", "density"), 1e-5, 1e-1)
- p.set_cmap(("gas", "density"), my_cmap)
- plot = p.plots[("gas", "density")]
-
- plot.figure = fig_density_map[time]
- plot.axes = grid_density_map[time][(ax - 1) * len(codes) + code].axes
- if code == 0:
- plot.cax = grid_density_map[time].cbar_axes[0]
- p._setup_plots()
-
- if add_nametag == 1:
- at = AnchoredText(
- "%s" % codes[code], loc=2, prop=dict(size=6), frameon=True
- )
- grid_density_map[time][code].axes.add_artist(at)
-
- # TEMPERATURE MAPS
- if draw_temperature_map == 1:
- my_cmap = copy.copy(matplotlib.cm.get_cmap("algae"))
- my_cmap.set_bad(
- my_cmap(0.6)
- ) # (log(1e4) - log(1e1)) / (log(1e6) - log(1e1)) = 0.6
- my_cmap.set_under(my_cmap(0))
- for ax in range(1, 3):
- p2 = ProjectionPlot(
- pf,
- ax,
- ("gas", "temperature"),
- center=center,
- data_source=proj_region,
- width=(figure_width, "kpc"),
- weight_field=("gas", "density_squared"),
- fontsize=9,
- )
- p2.set_zlim(("gas", "temperature"), 1e1, 1e6)
- p2.set_cmap(("gas", "temperature"), my_cmap)
- plot2 = p2.plots[("gas", "temperature")]
-
- plot2.figure = fig_temperature_map[time]
- plot2.axes = grid_temperature_map[time][
- (ax - 1) * len(codes) + code
- ].axes
- if code == 0:
- plot2.cax = grid_temperature_map[time].cbar_axes[0]
- p2._setup_plots()
-
- if add_nametag == 1:
- at = AnchoredText(
- "%s" % codes[code], loc=2, prop=dict(size=6), frameon=True
- )
- grid_temperature_map[time][code].axes.add_artist(at)
-
- # CELL-SIZE MAPS
- if draw_cellsize_map >= 1:
- my_cmap = copy.copy(matplotlib.cm.get_cmap("algae"))
- my_cmap.set_bad(my_cmap(0.99999)) # 1 doesn't work!
- my_cmap.set_under(my_cmap(0))
- if draw_cellsize_map == 1 or draw_cellsize_map == 3:
- for ax in range(1, 3):
- p25 = ProjectionPlot(
- pf,
- ax,
- ("index", "cell_size"),
- center=center,
- data_source=proj_region,
- width=(figure_width, "kpc"),
- weight_field=("index", "cell_volume_inv2"),
- fontsize=9,
- )
- p25.set_zlim(("index", "cell_size"), 10, 1e3)
- p25.set_cmap(("index", "cell_size"), my_cmap)
- plot25 = p25.plots[("index", "cell_size")]
-
- plot25.figure = fig_cellsize_map[time]
- plot25.axes = grid_cellsize_map[time][
- (ax - 1) * len(codes) + code
- ].axes
- if code == 0:
- plot25.cax = grid_cellsize_map[time].cbar_axes[0]
- p25._setup_plots()
-
- # Create another map with a different resolution definition if requested
- if draw_cellsize_map == 2 or draw_cellsize_map == 3:
- for ax in range(1, 3):
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
- p251 = ProjectionPlot(
- pf,
- ax,
- ("index", "cell_size"),
- center=center,
- data_source=proj_region,
- width=(figure_width, "kpc"),
- weight_field=("index", "cell_volume_inv2"),
- fontsize=9,
- )
- p251.set_zlim(("index", "cell_size"), 10, 1e3)
- p251.set_cmap(("index", "cell_size"), my_cmap)
- plot251 = p251.plots[("index", "cell_size")]
- else:
- p251 = ProjectionPlot(
- pf,
- ax,
- ("gas", "particle_size"),
- center=center,
- data_source=proj_region,
- width=(figure_width, "kpc"),
- weight_field=("gas", "particle_volume_inv2"),
- fontsize=9,
- )
- p251.set_zlim(("gas", "particle_size"), 10, 1e3)
- p251.set_cmap(("gas", "particle_size"), my_cmap)
- plot251 = p251.plots[("gas", "particle_size")]
-
- plot251.figure = fig_cellsize_map_2[time]
- plot251.axes = grid_cellsize_map_2[time][
- (ax - 1) * len(codes) + code
- ].axes
- if code == 0:
- plot251.cax = grid_cellsize_map_2[time].cbar_axes[0]
- p251._setup_plots()
-
- if add_nametag == 1:
- if draw_cellsize_map == 1 or draw_cellsize_map == 3:
- at = AnchoredText(
- "%s" % codes[code], loc=2, prop=dict(size=6), frameon=True
- )
- grid_cellsize_map[time][code].axes.add_artist(at)
- if draw_cellsize_map == 2 or draw_cellsize_map == 3:
- at = AnchoredText(
- "%s" % codes[code], loc=2, prop=dict(size=6), frameon=True
- )
- grid_cellsize_map_2[time][code].axes.add_artist(at)
-
- # ELEVATION MAPS
- if draw_elevation_map == 1:
-
- def _CellzElevationpc(field, data):
- return (data[("index", "z")] - center[2]).in_units("pc")
-
- pf.add_field(
- ("index", "z_elevation"),
- function=_CellzElevationpc,
- units="kpc",
- display_name="$z$ Elevation",
- take_log=False,
- )
- my_cmap = copy.copy(matplotlib.cm.get_cmap("algae"))
- my_cmap.set_bad(my_cmap(0.5))
- my_cmap.set_under(my_cmap(0))
- for ax in range(2, 3):
- p255 = ProjectionPlot(
- pf,
- ax,
- ("index", "z_elevation"),
- center=center,
- data_source=proj_region,
- width=(figure_width, "kpc"),
- weight_field=("gas", "density"),
- fontsize=9,
- )
- p255.set_zlim(("index", "z_elevation"), -1, 1)
- p255.set_cmap(("index", "z_elevation"), my_cmap)
- plot255 = p255.plots[("index", "z_elevation")]
-
- plot255.figure = fig_elevation_map[time]
- plot255.axes = grid_elevation_map[time][
- (ax - 2) * len(codes) + code
- ].axes
- if code == 0:
- plot255.cax = grid_elevation_map[time].cbar_axes[0]
- p255._setup_plots()
-
- if add_nametag == 1:
- at = AnchoredText(
- "%s" % codes[code], loc=2, prop=dict(size=6), frameon=True
- )
- grid_elevation_map[time][code].axes.add_artist(at)
-
- # Z-VELOCITY MAPS
- if draw_zvel_map == 1:
- my_cmap = copy.copy(matplotlib.cm.get_cmap("algae"))
- my_cmap.set_bad(my_cmap(0.5))
- my_cmap.set_under(my_cmap(0))
- for ax in range(1, 3):
- p26 = ProjectionPlot(
- pf,
- ax,
- ("gas", "velocity_z"),
- center=center,
- data_source=proj_region,
- width=(figure_width, "kpc"),
- weight_field=("gas", "density_squared"),
- fontsize=9,
- )
- p26.set_zlim(("gas", "velocity_z"), -1e7, 1e7)
- p26.set_cmap(("gas", "velocity_z"), my_cmap)
- plot26 = p26.plots[("gas", "velocity_z")]
-
- plot26.figure = fig_zvel_map[time]
- plot26.axes = grid_zvel_map[time][(ax - 1) * len(codes) + code].axes
- if code == 0:
- plot26.cax = grid_zvel_map[time].cbar_axes[0]
- p26._setup_plots()
-
- if add_nametag == 1:
- at = AnchoredText(
- "%s" % codes[code], loc=2, prop=dict(size=6), frameon=True
- )
- grid_zvel_map[time][code].axes.add_artist(at)
-
- # METAL MAPS
- if draw_metal_map >= 1:
- my_cmap = copy.copy(matplotlib.cm.get_cmap("algae"))
- my_cmap.set_bad(my_cmap(0.347)) # (0.02041 - 0.01) / (0.04 - 0.01) = 0.347
- my_cmap.set_under(my_cmap(0))
- for ax in range(1, 3):
- pf.field_info[("gas", "metallicity")].take_log = False
- p26 = ProjectionPlot(
- pf,
- ax,
- ("gas", "metallicity"),
- center=center,
- data_source=proj_region,
- width=(figure_width, "kpc"),
- weight_field=("gas", "density_squared"),
- fontsize=9,
- )
- p26.set_zlim(("gas", "metallicity"), 0.01, 0.04)
- p26.set_cmap(("gas", "metallicity"), my_cmap)
- plot26 = p26.plots[("gas", "metallicity")]
-
- plot26.figure = fig_metal_map[time]
- plot26.axes = grid_metal_map[time][(ax - 1) * len(codes) + code].axes
- if code == 0:
- plot26.cax = grid_metal_map[time].cbar_axes[0]
- p26._setup_plots()
-
- if add_nametag == 1:
- if draw_metal_map == 2:
- sp = pf.sphere(center, (1, "Mpc"))
- total_metal_mass = sp.quantities.total_quantity(
- [("gas", "metal_mass")]
- )
- at = AnchoredText(
- "%s - %e" % (codes[code], total_metal_mass),
- loc=2,
- prop=dict(size=6),
- frameon=True,
- )
- else:
- at = AnchoredText(
- "%s" % codes[code], loc=2, prop=dict(size=6), frameon=True
- )
- grid_metal_map[time][code].axes.add_artist(at)
-
- # STELLAR MAPS
- if draw_star_map == 1 and time != 0:
- for ax in range(1, 3):
- p27 = ParticleProjectionPlot(
- pf,
- ax,
- (PartType_Star_to_use, "particle_mass"),
- center=center,
- data_source=proj_region,
- width=(figure_width, "kpc"),
- weight_field=None,
- fontsize=9,
- )
- p27.set_unit((PartType_Star_to_use, "particle_mass"), "Msun")
- p27.set_zlim((PartType_Star_to_use, "particle_mass"), 1e4, 1e7)
- p27.set_cmap((PartType_Star_to_use, "particle_mass"), "algae")
- p27.set_buff_size(400) # default is 800
- p27.set_colorbar_label(
- (PartType_Star_to_use, "particle_mass"),
- "Stellar Mass Per Pixel ($\mathrm{M}_{\odot}$)",
- )
- plot27 = p27.plots[(PartType_Star_to_use, "particle_mass")]
- # p27 = ParticleProjectionPlot(pf, ax, (PartType_Star_to_use, "particle_velocity_cylindrical_theta"), center = center, data_source=proj_region, width = (figure_width, 'kpc'), weight_field = (PartType_Star_to_use, "particle_mass"), fontsize=9)
- # p27.set_unit((PartType_Star_to_use, "particle_velocity_cylindrical_theta"), 'km/s')
- # p27.set_buff_size(400)
- # p27.set_colorbar_label((PartType_Star_to_use, "particle_velocity_cylindrical_theta"), "Rotational Velocity (km/s)")
- # plot27 = p27.plots[(PartType_Star_to_use, "particle_velocity_cylindrical_theta")]
-
- plot27.figure = fig_star_map[time]
- plot27.axes = grid_star_map[time][(ax - 1) * len(codes) + code].axes
- if code == 0:
- plot27.cax = grid_star_map[time].cbar_axes[0]
- p27._setup_plots()
-
- if add_nametag == 1:
- at = AnchoredText(
- "%s" % codes[code], loc=2, prop=dict(size=6), frameon=True
- )
- grid_star_map[time][code].axes.add_artist(at)
-
- # STELLAR CLUMP STATISTICS AND ANNOTATED STELLAR MAPS
- if draw_star_clump_stats >= 1 and time != 0:
- # For make yt's HaloCatalog to work with non-cosmological dataset, a fix needed to be applied to analysis_modules/halo_finding/halo_objects.py: self.period = ds.arr([3.254, 3.254, 3.254], 'Mpc')
- pf.hubble_constant = 0.71
- pf.omega_lambda = 0.73
- pf.omega_matter = 0.27
- pf.omega_curvature = 0.0
- pf.current_redshift = (
- 0.0
- ) # another trick to make HaloCatalog work especially for ART-I dataset
- if (
- os.path.exists(
- "./halo_catalogs/hop_%s_%05d/hop_%s_%05d.0.h5"
- % (codes[code], times[time], codes[code], times[time])
- )
- == False
- ):
- hc = HaloCatalog(
- data_ds=pf,
- finder_method="hop",
- output_dir="./halo_catalogs/hop_%s_%05d"
- % (codes[code], times[time]),
- finder_kwargs={
- "threshold": 2e8,
- "dm_only": False,
- "ptype": PartType_Star_to_use,
- },
- )
- # finder_kwargs={'threshold': 2e5, 'dm_only': False, 'ptype': "all"})
- hc.add_filter(
- "quantity_value", "particle_mass", ">", 2.6e6, "Msun"
- ) # exclude halos with less than 30 particles
- hc.add_filter(
- "quantity_value", "particle_mass", "<", 2.6e8, "Msun"
- ) # exclude the most massive halo (threshold 1e8.4 is hand-picked, so one needs to be careful!)
- hc.create()
- halo_ds = load(
- "./halo_catalogs/hop_%s_%05d/hop_%s_%05d.0.h5"
- % (codes[code], times[time], codes[code], times[time])
- )
- hc = HaloCatalog(
- halos_ds=halo_ds,
- output_dir="./halo_catalogs/hop_%s_%05d" % (codes[code], times[time]),
- )
- hc.load()
-
- halo_ad = hc.halos_ds.all_data()
- star_clump_masses_hop[time].append(
- np.log10(halo_ad["particle_mass"][:].in_units("Msun"))
- )
-
- if (
- os.path.exists(
- "./halo_catalogs/fof_%s_%05d/fof_%s_%05d.0.h5"
- % (codes[code], times[time], codes[code], times[time])
- )
- == False
- ):
- hc2 = HaloCatalog(
- data_ds=pf,
- finder_method="fof",
- output_dir="./halo_catalogs/fof_%s_%05d"
- % (codes[code], times[time]),
- finder_kwargs={
- "link": 0.0025,
- "dm_only": False,
- "ptype": PartType_Star_to_use,
- },
- )
- # finder_kwargs={'link': 0.02, 'dm_only': False, 'ptype': "all"})
- hc2.add_filter(
- "quantity_value", "particle_mass", ">", 2.6e6, "Msun"
- ) # exclude halos with less than 30 particles
- hc2.add_filter(
- "quantity_value", "particle_mass", "<", 2.6e8, "Msun"
- ) # exclude the most massive halo (threshold 1e8.4 is hand-picked, so one needs to be careful!)
- hc2.create()
- halo_ds2 = load(
- "./halo_catalogs/fof_%s_%05d/fof_%s_%05d.0.h5"
- % (codes[code], times[time], codes[code], times[time])
- )
- hc2 = HaloCatalog(
- halos_ds=halo_ds2,
- output_dir="./halo_catalogs/fof_%s_%05d" % (codes[code], times[time]),
- )
- hc2.load()
-
- halo_ad2 = hc2.halos_ds.all_data()
- star_clump_masses_fof[time].append(
- np.log10(halo_ad2["particle_mass"][:].in_units("Msun"))
- )
- # most_massive = halo_ad['particle_mass'].max().in_units('Msun')
- # cen4 = [halo_ad['particle_position_x'][halo_ad['particle_mass'] == most_massive][0].in_units('code_length').d,
- # halo_ad['particle_position_y'][halo_ad['particle_mass'] == most_massive][0].in_units('code_length').d,
- # halo_ad['particle_position_z'][halo_ad['particle_mass'] == most_massive][0].in_units('code_length').d]
-
- # Add additional star_map with annotated clumps if requested
- if draw_star_clump_stats >= 2:
- for ax in range(1, 3):
- p271 = ParticleProjectionPlot(
- pf,
- ax,
- (PartType_Star_to_use, "particle_mass"),
- center=center,
- data_source=proj_region,
- width=(figure_width, "kpc"),
- weight_field=None,
- fontsize=9,
- )
- p271.set_unit((PartType_Star_to_use, "particle_mass"), "Msun")
- p271.set_zlim((PartType_Star_to_use, "particle_mass"), 1e4, 1e7)
- p271.set_cmap((PartType_Star_to_use, "particle_mass"), "algae")
- p271.set_buff_size(400) # default is 800
- p271.set_colorbar_label(
- (PartType_Star_to_use, "particle_mass"),
- "Stellar Mass Per Pixel ($\mathrm{M}_{\odot}$)",
- )
- plot271 = p271.plots[(PartType_Star_to_use, "particle_mass")]
- if ax == 2:
- p271.annotate_halos(
- hc,
- factor=1.0,
- circle_args={
- "linewidth": 0.7,
- "alpha": 0.8,
- "facecolor": "none",
- "edgecolor": "k",
- },
- )
-
- plot271.figure = fig_star_map_2[time]
- plot271.axes = grid_star_map_2[time][
- (ax - 1) * len(codes) + code
- ].axes
- if code == 0:
- plot271.cax = grid_star_map_2[time].cbar_axes[0]
- p271._setup_plots()
-
- if add_nametag == 1:
- at = AnchoredText(
- "%s" % codes[code], loc=2, prop=dict(size=6), frameon=True
- )
- grid_star_map_2[time][code].axes.add_artist(at)
-
- for ax in range(1, 3):
- p272 = ParticleProjectionPlot(
- pf,
- ax,
- (PartType_Star_to_use, "particle_mass"),
- center=center,
- data_source=proj_region,
- width=(figure_width, "kpc"),
- weight_field=None,
- fontsize=9,
- )
- p272.set_unit((PartType_Star_to_use, "particle_mass"), "Msun")
- p272.set_zlim((PartType_Star_to_use, "particle_mass"), 1e4, 1e7)
- p272.set_cmap((PartType_Star_to_use, "particle_mass"), "algae")
- p272.set_buff_size(400) # default is 800
- p272.set_colorbar_label(
- (PartType_Star_to_use, "particle_mass"),
- "Stellar Mass Per Pixel ($\mathrm{M}_{\odot}$)",
- )
- plot272 = p272.plots[(PartType_Star_to_use, "particle_mass")]
- # p272 = ParticleProjectionPlot(pf, ax, ("all", "particle_mass"), center = center, data_source=proj_region, width = (figure_width, 'kpc'), weight_field = None, fontsize=9)
- # p272.set_unit(("all", "particle_mass"), 'Msun')
- # p272.set_zlim(("all", "particle_mass"), 1e4, 1e9)
- # p272.set_cmap(("all", "particle_mass"), 'algae')
- # p272.set_buff_size(400) # default is 800
- # p272.set_colorbar_label(("all", "particle_mass"), "Mass Per Pixel ($\mathrm{M}_{\odot}$)")
- # plot272 = p272.plots[("all", "particle_mass")]
- if ax == 2:
- p272.annotate_halos(
- hc2,
- factor=1.0,
- circle_args={
- "linewidth": 0.7,
- "alpha": 0.8,
- "facecolor": "none",
- "edgecolor": "k",
- },
- ) # , annotate_field='particle_mass')
-
- plot272.figure = fig_star_map_3[time]
- plot272.axes = grid_star_map_3[time][
- (ax - 1) * len(codes) + code
- ].axes
- if code == 0:
- plot272.cax = grid_star_map_3[time].cbar_axes[0]
- p272._setup_plots()
-
- if add_nametag == 1:
- at = AnchoredText(
- "%s" % codes[code], loc=2, prop=dict(size=6), frameon=True
- )
- grid_star_map_3[time][code].axes.add_artist(at)
-
- # SFR MAPS
- if draw_SFR_map >= 1 and time != 0:
- sp = pf.sphere(center, (1.0 * figure_width, "kpc"))
- xy_bins = int(float(figure_width) / (float(aperture_size_SFR_map / 1000.0)))
- my_cmap = copy.copy(matplotlib.cm.get_cmap("algae"))
- my_cmap.set_bad(my_cmap(0))
- my_cmap.set_under(my_cmap(0))
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
-
- def _position_x_from_center(field, data):
- return (data[("index", "x")] - center[0]).in_units("kpc")
-
- pf.add_field(
- ("index", "position_x_from_center"),
- function=_position_x_from_center,
- take_log=False,
- particle_type=False,
- units="kpc",
- )
-
- def _position_y_from_center(field, data):
- return (data[("index", "y")] - center[1]).in_units("kpc")
-
- pf.add_field(
- ("index", "position_y_from_center"),
- function=_position_y_from_center,
- take_log=False,
- particle_type=False,
- units="kpc",
- )
-
- def _cell_mass_surface_density_SFR_map(field, data):
- trans = np.zeros(data[("gas", "cell_mass")].shape)
- ind = np.where(data[("gas", "cell_mass")] > 0)
- trans[ind] = (
- data[("gas", "cell_mass")][ind].in_units("Msun").d
- / (float(aperture_size_SFR_map)) ** 2
- )
- return data.ds.arr(trans, "Msun/pc**2").in_base(
- data.ds.unit_system.name
- )
-
- pf.add_field(
- ("gas", "cell_mass_surface_density_SFR_map"),
- function=_cell_mass_surface_density_SFR_map,
- display_name="$\mathrm{\Sigma}_\mathrm{gas}$",
- units="Msun/pc**2",
- particle_type=False,
- take_log=True,
- )
-
- p28 = PhasePlot(
- sp,
- ("index", "position_x_from_center"),
- ("index", "position_y_from_center"),
- ("gas", "cell_mass_surface_density_SFR_map"),
- weight_field=None,
- fontsize=9,
- x_bins=xy_bins,
- y_bins=xy_bins,
- )
- p28.set_zlim(("gas", "cell_mass_surface_density_SFR_map"), 1e0, 1e3)
- p28.set_cmap(("gas", "cell_mass_surface_density_SFR_map"), my_cmap)
- plot28 = p28.plots[("gas", "cell_mass_surface_density_SFR_map")]
- else:
-
- def _particle_position_x_from_center(field, data):
- return (
- data[(PartType_Gas_to_use, "particle_position_x")] - center[0]
- ).in_units("kpc")
-
- pf.add_field(
- (PartType_Gas_to_use, "particle_position_x_from_center"),
- function=_particle_position_x_from_center,
- take_log=False,
- particle_type=True,
- units="kpc",
- )
-
- def _particle_position_y_from_center(field, data):
- return (
- data[(PartType_Gas_to_use, "particle_position_y")] - center[0]
- ).in_units("kpc")
-
- pf.add_field(
- (PartType_Gas_to_use, "particle_position_y_from_center"),
- function=_particle_position_y_from_center,
- take_log=False,
- particle_type=True,
- units="kpc",
- )
-
- def _particle_mass_surface_density_SFR_map(field, data):
- trans = np.zeros(data[(PartType_Gas_to_use, "particle_mass")].shape)
- ind = np.where(data[(PartType_Gas_to_use, "particle_mass")] > 0)
- trans[ind] = (
- data[(PartType_Gas_to_use, "particle_mass")][ind]
- .in_units("Msun")
- .d
- / (float(aperture_size_SFR_map)) ** 2
- )
- return data.ds.arr(trans, "Msun/pc**2").in_base(
- data.ds.unit_system.name
- )
-
- pf.add_field(
- (PartType_Gas_to_use, "particle_mass_surface_density_SFR_map"),
- function=_particle_mass_surface_density_SFR_map,
- display_name="$\mathrm{\Sigma}_\mathrm{gas}$",
- units="Msun/pc**2",
- particle_type=True,
- take_log=True,
- )
-
- p28 = ParticlePhasePlot(
- sp,
- (PartType_Gas_to_use, "particle_position_x_from_center"),
- (PartType_Gas_to_use, "particle_position_y_from_center"),
- (PartType_Gas_to_use, "particle_mass_surface_density_SFR_map"),
- weight_field=None,
- deposition="cic",
- fontsize=9,
- x_bins=xy_bins,
- y_bins=xy_bins,
- )
- p28.set_zlim(
- (PartType_Gas_to_use, "particle_mass_surface_density_SFR_map"),
- 1e0,
- 1e3,
- )
- p28.set_cmap(
- (PartType_Gas_to_use, "particle_mass_surface_density_SFR_map"),
- my_cmap,
- )
- plot28 = p28.plots[
- (PartType_Gas_to_use, "particle_mass_surface_density_SFR_map")
- ]
-
- p28.set_xlabel("x (kpc)")
- p28.set_ylabel("y (kpc)")
- p28.set_xlim(-15, 15)
- p28.set_ylim(-15, 15)
- plot28.figure = fig_degr_density_map[time]
- plot28.axes = grid_degr_density_map[time][code].axes
- if code == 0:
- plot28.cax = grid_degr_density_map[time].cbar_axes[0]
- p28._setup_plots()
-
- def _particle_position_x_from_center_2(field, data):
- return (
- data[(PartType_StarBeforeFiltered_to_use, "particle_position_x")]
- - center[0]
- ).in_units("kpc")
-
- pf.add_field(
- (PartType_StarBeforeFiltered_to_use, "particle_position_x_from_center"),
- function=_particle_position_x_from_center_2,
- take_log=False,
- particle_type=True,
- units="kpc",
- )
- pf.add_particle_filter(
- PartType_Star_to_use
- ) # This is needed for a filtered particle type PartType_Star_to_use to work, because we have just created new particle fields.
-
- def _particle_position_y_from_center_2(field, data):
- return (
- data[(PartType_StarBeforeFiltered_to_use, "particle_position_y")]
- - center[0]
- ).in_units("kpc")
-
- pf.add_field(
- (PartType_StarBeforeFiltered_to_use, "particle_position_y_from_center"),
- function=_particle_position_y_from_center_2,
- take_log=False,
- particle_type=True,
- units="kpc",
- )
- pf.add_particle_filter(PartType_Star_to_use)
-
- def _particle_mass_young_stars_SFR_surface_density_SFR_map(field, data):
- trans = np.zeros(
- data[(PartType_StarBeforeFiltered_to_use, "particle_mass")].shape
- )
- ind = np.where(
- data[
- (PartType_StarBeforeFiltered_to_use, FormationTimeType_to_use)
- ].in_units("Myr")
- > (pf.current_time.in_units("Myr").d - young_star_cutoff_SFR_map)
- ) # mass for young stars only
- trans[ind] = (
- data[(PartType_StarBeforeFiltered_to_use, "particle_mass")][ind]
- .in_units("Msun")
- .d
- / (young_star_cutoff_SFR_map * 1e6)
- / (float(aperture_size_SFR_map) / 1000.0) ** 2
- )
- return data.ds.arr(trans, "Msun/yr/kpc**2").in_base(
- data.ds.unit_system.name
- )
-
- pf.add_field(
- (
- PartType_StarBeforeFiltered_to_use,
- "particle_mass_young_stars_SFR_surface_density_SFR_map",
- ),
- function=_particle_mass_young_stars_SFR_surface_density_SFR_map,
- display_name="$\mathrm{\Sigma}_\mathrm{SFR}$",
- units="Msun/yr/kpc**2",
- particle_type=True,
- take_log=True,
- )
- pf.add_particle_filter(PartType_Star_to_use)
-
- p281 = ParticlePhasePlot(
- sp,
- (PartType_Star_to_use, "particle_position_x_from_center"),
- (PartType_Star_to_use, "particle_position_y_from_center"),
- (
- PartType_Star_to_use,
- "particle_mass_young_stars_SFR_surface_density_SFR_map",
- ),
- deposition="cic",
- weight_field=None,
- fontsize=9,
- x_bins=xy_bins,
- y_bins=xy_bins,
- )
- p281.set_zlim(
- (
- PartType_Star_to_use,
- "particle_mass_young_stars_SFR_surface_density_SFR_map",
- ),
- 3e-4,
- 3e-1,
- )
- p281.set_cmap(
- (
- PartType_Star_to_use,
- "particle_mass_young_stars_SFR_surface_density_SFR_map",
- ),
- my_cmap,
- )
- plot281 = p281.plots[
- (
- PartType_Star_to_use,
- "particle_mass_young_stars_SFR_surface_density_SFR_map",
- )
- ]
-
- p281.set_xlabel("x (kpc)")
- p281.set_ylabel("y (kpc)")
- p281.set_xlim(-15, 15)
- p281.set_ylim(-15, 15)
- plot281.figure = fig_degr_sfr_map[time]
- plot281.axes = grid_degr_sfr_map[time][code].axes
- if code == 0:
- plot281.cax = grid_degr_sfr_map[time].cbar_axes[0]
- p281._setup_plots()
-
- if add_nametag == 1:
- at = AnchoredText(
- "%s" % codes[code], loc=2, prop=dict(size=6), frameon=True
- )
- grid_degr_density_map[time][code].axes.add_artist(at)
- at2 = AnchoredText(
- "%s" % codes[code], loc=2, prop=dict(size=6), frameon=True
- )
- grid_degr_sfr_map[time][code].axes.add_artist(at2)
-
- # Record degraded maps for the later use in a local K-S plot
- if draw_SFR_map == 2 and time != 0:
- fn = (
- p28.profile.save_as_dataset()
- ) # see http://yt-project.org/docs/dev/reference/api/generated/yt.data_objects.profiles.ProfileND.save_as_dataset.html
- phaseplot_ds = load(fn)
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
- phaseplot_surf_dens = phaseplot_ds.data[
- ("data", "cell_mass_surface_density_SFR_map")
- ]
- else:
- phaseplot_surf_dens = phaseplot_ds.data[
- ("data", "particle_mass_surface_density_SFR_map")
- ]
- fn_1 = p281.profile.save_as_dataset()
- phaseplot_ds_1 = load(fn_1)
- phaseplot_sfr_surf_dens = phaseplot_ds_1.data[
- ("data", "particle_mass_young_stars_SFR_surface_density_SFR_map")
- ]
- xy_shape = phaseplot_surf_dens.shape[0]
- surf_dens_SFR_map[time].append(
- phaseplot_surf_dens.reshape(1, xy_shape ** 2)[0]
- )
- sfr_surf_dens_SFR_map[time].append(
- phaseplot_sfr_surf_dens.reshape(1, xy_shape ** 2)[0]
- )
- call(
- ["rm", "%s_%s.h5" % (str(pf), "ParticleProfile")]
- ) # Remove unwanted .h5 files saved
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
- call(["rm", "%s_%s.h5" % (str(pf), "Profile2D")])
-
- # DENSITY-TEMPERATURE PDF
- if draw_PDF >= 1:
- sp = pf.sphere(center, (0.5 * figure_width, "kpc"))
- my_cmap2 = copy.copy(matplotlib.cm.get_cmap("algae"))
- my_cmap2.set_under(my_cmap2(0))
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
- p3 = PhasePlot(
- sp,
- ("gas", "density"),
- ("gas", "temperature"),
- ("gas", "cell_mass"),
- weight_field=None,
- fontsize=12,
- x_bins=300,
- y_bins=300,
- )
- p3.set_unit("cell_mass", "Msun")
- p3.set_zlim(("gas", "cell_mass"), 1e3, 1e8)
- p3.set_cmap(("gas", "cell_mass"), my_cmap2)
- p3.set_colorbar_label(
- ("gas", "cell_mass"), "Mass ($\mathrm{M}_{\odot}$)"
- )
- plot3 = p3.plots[("gas", "cell_mass")]
- else:
- # Because ParticlePhasePlot doesn't yet work for a log-log PDF for some reason, I will do the following trick.
- p3 = PhasePlot(
- sp,
- (PartType_Gas_to_use, "density"),
- (PartType_Gas_to_use, "temperature"),
- (PartType_Gas_to_use, MassType_to_use),
- weight_field=None,
- fontsize=12,
- x_bins=300,
- y_bins=300,
- )
- p3.set_unit(MassType_to_use, "Msun")
- p3.set_zlim((PartType_Gas_to_use, MassType_to_use), 1e3, 1e8)
- p3.set_cmap((PartType_Gas_to_use, MassType_to_use), my_cmap2)
- plot3 = p3.plots[(PartType_Gas_to_use, MassType_to_use)]
-
- p3.set_unit("density", "g/cm**3")
- p3.set_xlim(1e-29, 1e-21)
- p3.set_unit("temperature", "K")
- p3.set_ylim(10, 1e7)
- p3.set_log("density", True)
- p3.set_log("temperature", True)
-
- plot3.figure = fig_PDF[time]
- plot3.axes = grid_PDF[time][code].axes
- if code == 0:
- plot3.cax = grid_PDF[time].cbar_axes[0]
- p3._setup_plots()
-
- # Add constant pressure and constant entropy lines to guide the eyes
- if draw_PDF >= 2:
- dummy_density = np.logspace(-29, -21, num=4)
- for constant_i in range(1, 100):
- constant = 1e-100 * 100 ** constant_i
- line = ln.Line2D(
- dummy_density,
- constant / dummy_density,
- linestyle=":",
- linewidth=0.6,
- color="k",
- alpha=0.7,
- ) # constant pressure P ~ n*T
- grid_PDF[time][code].axes.add_line(line)
- line2 = ln.Line2D(
- dummy_density,
- constant * dummy_density ** (2.0 / 3.0),
- linestyle="-.",
- linewidth=0.6,
- color="k",
- alpha=0.7,
- ) # constant entropy S ~ n**(-2/3)*T
- grid_PDF[time][code].axes.add_line(line2)
-
- if add_nametag == 1:
- at = AnchoredText(
- "%s" % codes[code], loc=3, prop=dict(size=10), frameon=True
- )
- grid_PDF[time][code].axes.add_artist(at)
-
- # POSITION-VELOCITY PDF FOR GAS
- if draw_pos_vel_PDF >= 1:
- sp = pf.sphere(center, (0.5 * figure_width, "kpc"))
- sp.set_field_parameter("normal", disk_normal_vector)
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
- sp_dense = sp.cut_region(
- ["obj['gas', 'density'].in_units('g/cm**3') > 1.e-25"]
- ) # For AMR codes, consider only the cells that are dense enough
- pf.field_info[("index", "cylindrical_r")].take_log = False
- pf.field_info[
- ("gas", "cylindrical_tangential_velocity")
- ].take_log = False
- p4 = PhasePlot(
- sp_dense,
- ("index", "cylindrical_r"),
- ("gas", "cylindrical_tangential_velocity"),
- ("gas", "cell_mass"),
- weight_field=None,
- fontsize=12,
- x_bins=300,
- y_bins=300,
- )
- p4.set_unit("cylindrical_r", "kpc")
- p4.set_unit("cylindrical_tangential_velocity", "km/s")
- p4.set_unit("cell_mass", "Msun")
- p4.set_zlim(("gas", "cell_mass"), 1e3, 1e8)
- p4.set_cmap(("gas", "cell_mass"), "algae")
- p4.set_colorbar_label(
- ("gas", "cell_mass"), "Mass ($\mathrm{M}_{\odot}$)"
- )
- plot4 = p4.plots[("gas", "cell_mass")]
- else:
- pf.field_info[(PartType_Gas_to_use, "radius")].take_log = False
- pf.field_info[
- (PartType_Gas_to_use, "particle_velocity_cylindrical_theta")
- ].take_log = False
- p4 = ParticlePhasePlot(
- sp,
- (PartType_Gas_to_use, "radius"),
- (PartType_Gas_to_use, "particle_velocity_cylindrical_theta"),
- (PartType_Gas_to_use, MassType_to_use),
- deposition="ngp",
- weight_field=None,
- fontsize=12,
- x_bins=300,
- y_bins=300,
- )
- p4.set_unit("radius", "kpc")
- p4.set_unit("particle_velocity_cylindrical_theta", "km/s")
- p4.set_unit(MassType_to_use, "Msun")
- p4.set_zlim((PartType_Gas_to_use, MassType_to_use), 1e3, 1e8)
- p4.set_cmap((PartType_Gas_to_use, MassType_to_use), "algae")
- p4.set_colorbar_label(
- (PartType_Gas_to_use, MassType_to_use),
- "Mass ($\mathrm{M}_{\odot}$)",
- )
- plot4 = p4.plots[(PartType_Gas_to_use, MassType_to_use)]
-
- p4.set_xlabel("Cylindrical Radius (kpc)")
- p4.set_ylabel("Rotational Velocity (km/s)")
- p4.set_xlim(0, 14)
- p4.set_ylim(-100, 350)
-
- plot4.figure = fig_pos_vel_PDF[time]
- plot4.axes = grid_pos_vel_PDF[time][code].axes
- if code == 0:
- plot4.cax = grid_pos_vel_PDF[time].cbar_axes[0]
- p4._setup_plots()
-
- # Add 1D profile line if requested
- if draw_pos_vel_PDF >= 2 and time != 0:
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
- p5 = ProfilePlot(
- sp_dense,
- ("index", "cylindrical_r"),
- ("gas", "cylindrical_tangential_velocity"),
- weight_field=("gas", "cell_mass"),
- n_bins=50,
- x_log=False,
- )
- p5.set_log("cylindrical_tangential_velocity", False)
- p5.set_log("cylindrical_r", False)
- p5.set_unit("cylindrical_r", "kpc")
- p5.set_xlim(1e-3, 14)
- p5.set_ylim("cylindrical_tangential_velocity", -100, 350)
- line = ln.Line2D(
- p5.profiles[0].x.in_units("kpc"),
- p5.profiles[0]["cylindrical_tangential_velocity"].in_units(
- "km/s"
- ),
- linestyle="-",
- linewidth=2,
- color="k",
- alpha=0.7,
- )
- pos_vel_xs[time].append(p5.profiles[0].x.in_units("kpc").d)
- pos_vel_profiles[time].append(
- p5.profiles[0]["cylindrical_tangential_velocity"]
- .in_units("km/s")
- .d
- )
- else:
- p5 = ProfilePlot(
- sp,
- (PartType_Gas_to_use, "radius"),
- (PartType_Gas_to_use, "particle_velocity_cylindrical_theta"),
- weight_field=(PartType_Gas_to_use, MassType_to_use),
- n_bins=50,
- x_log=False,
- )
- p5.set_log("particle_velocity_cylindrical_theta", False)
- p5.set_log("radius", False)
- p5.set_unit("radius", "kpc")
- p5.set_xlim(1e-3, 14)
- p5.set_ylim("particle_velocity_cylindrical_theta", -100, 350)
- line = ln.Line2D(
- p5.profiles[0].x.in_units("kpc"),
- p5.profiles[0]["particle_velocity_cylindrical_theta"].in_units(
- "km/s"
- ),
- linestyle="-",
- linewidth=2,
- color="k",
- alpha=0.7,
- )
- pos_vel_xs[time].append(p5.profiles[0].x.in_units("kpc").d)
- pos_vel_profiles[time].append(
- p5.profiles[0]["particle_velocity_cylindrical_theta"]
- .in_units("km/s")
- .d
- )
- grid_pos_vel_PDF[time][code].axes.add_line(line)
-
- # Add velocity dispersion profile if requested
- if draw_pos_vel_PDF >= 3 and time != 0:
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
- # To calculate velocity dispersion -- residual velocity components other than the rotational velocity found above -- we will need a mass-weighted average of [v_i - v_rot_local]**2
- # Here we assumed that the disk is on x-y plane; i.e. the variable disk_normal_vector above needs to be [0.0, 0.0, 1.0]
- def _local_rotational_velocity_x(field, data):
- trans = np.zeros(data[("gas", "velocity_x")].shape)
- dr = 0.5 * (
- pos_vel_xs[time][code][1] - pos_vel_xs[time][code][0]
- )
- for radius, v_rot_local in zip(
- pos_vel_xs[time][code], pos_vel_profiles[time][code]
- ):
- ind = np.where(
- (
- data[("index", "cylindrical_r")].in_units("kpc")
- >= (radius - dr)
- )
- & (
- data[("index", "cylindrical_r")].in_units("kpc")
- < (radius + dr)
- )
- )
- trans[ind] = (
- -np.sin(data["index", "cylindrical_theta"][ind])
- * v_rot_local
- * 1e5
- ) # in cm/s
- return data.ds.arr(trans, "cm/s").in_base(
- data.ds.unit_system.name
- )
-
- pf.add_field(
- ("gas", "local_rotational_velocity_x"),
- function=_local_rotational_velocity_x,
- take_log=False,
- particle_type=False,
- units="cm/s",
- )
-
- def _local_rotational_velocity_y(field, data):
- trans = np.zeros(data[("gas", "velocity_y")].shape)
- dr = 0.5 * (
- pos_vel_xs[time][code][1] - pos_vel_xs[time][code][0]
- )
- for radius, v_rot_local in zip(
- pos_vel_xs[time][code], pos_vel_profiles[time][code]
- ):
- ind = np.where(
- (
- data[("index", "cylindrical_r")].in_units("kpc")
- >= (radius - dr)
- )
- & (
- data[("index", "cylindrical_r")].in_units("kpc")
- < (radius + dr)
- )
- )
- trans[ind] = (
- np.cos(data["index", "cylindrical_theta"][ind])
- * v_rot_local
- * 1e5
- )
- return data.ds.arr(trans, "cm/s").in_base(
- data.ds.unit_system.name
- )
-
- pf.add_field(
- ("gas", "local_rotational_velocity_y"),
- function=_local_rotational_velocity_y,
- take_log=False,
- particle_type=False,
- units="cm/s",
- )
-
- def _velocity_minus_local_rotational_velocity_squared(field, data):
- return (
- (
- data[("gas", "velocity_x")]
- - data[("gas", "local_rotational_velocity_x")]
- )
- ** 2
- + (
- data[("gas", "velocity_y")]
- - data[("gas", "local_rotational_velocity_y")]
- )
- ** 2
- + (data[("gas", "velocity_z")]) ** 2
- )
-
- pf.add_field(
- ("gas", "velocity_minus_local_rotational_velocity_squared"),
- function=_velocity_minus_local_rotational_velocity_squared,
- take_log=False,
- particle_type=False,
- units="cm**2/s**2",
- )
- # slc = SlicePlot(pf, 'z', ("gas", "local_rotational_velocity_y"), center = center, width = (figure_width, 'kpc'))
- # slc = SlicePlot(pf, 'z', ("gas", "velocity_minus_local_rotational_velocity_squared"), center = center, width = (figure_width, 'kpc')) # this should give zeros everywhere for IC at t=0
- # slc.save()
-
- p55 = ProfilePlot(
- sp_dense,
- ("index", "cylindrical_r"),
- ("gas", "velocity_minus_local_rotational_velocity_squared"),
- weight_field=("gas", "cell_mass"),
- n_bins=50,
- x_log=False,
- )
- p55.set_log("cylindrical_r", False)
- p55.set_unit("cylindrical_r", "kpc")
- p55.set_xlim(1e-3, 14)
- pos_disp_xs[time].append(p55.profiles[0].x.in_units("kpc").d)
- pos_disp_profiles[time].append(
- np.sqrt(
- p55.profiles[0][
- "velocity_minus_local_rotational_velocity_squared"
- ]
- )
- .in_units("km/s")
- .d
- )
- else:
-
- def _particle_local_rotational_velocity_x(field, data):
- trans = np.zeros(
- data[(PartType_Gas_to_use, "particle_velocity_x")].shape
- )
- dr = 0.5 * (
- pos_vel_xs[time][code][1] - pos_vel_xs[time][code][0]
- )
- for radius, v_rot_local in zip(
- pos_vel_xs[time][code], pos_vel_profiles[time][code]
- ):
- ind = np.where(
- (
- data[(PartType_Gas_to_use, "radius")].in_units(
- "kpc"
- )
- >= (radius - dr)
- )
- & (
- data[(PartType_Gas_to_use, "radius")].in_units(
- "kpc"
- )
- < (radius + dr)
- )
- )
- trans[ind] = (
- -np.sin(
- data[
- (
- PartType_Gas_to_use,
- "particle_position_cylindrical_theta",
- )
- ][ind]
- )
- * v_rot_local
- * 1e5
- )
- return data.ds.arr(trans, "cm/s").in_base(
- data.ds.unit_system.name
- )
-
- pf.add_field(
- (PartType_Gas_to_use, "particle_local_rotational_velocity_x"),
- function=_particle_local_rotational_velocity_x,
- take_log=False,
- particle_type=True,
- units="cm/s",
- )
-
- def _particle_local_rotational_velocity_y(field, data):
- trans = np.zeros(
- data[(PartType_Gas_to_use, "particle_velocity_y")].shape
- )
- dr = 0.5 * (
- pos_vel_xs[time][code][1] - pos_vel_xs[time][code][0]
- )
- for radius, v_rot_local in zip(
- pos_vel_xs[time][code], pos_vel_profiles[time][code]
- ):
- ind = np.where(
- (
- data[(PartType_Gas_to_use, "radius")].in_units(
- "kpc"
- )
- >= (radius - dr)
- )
- & (
- data[(PartType_Gas_to_use, "radius")].in_units(
- "kpc"
- )
- < (radius + dr)
- )
- )
- trans[ind] = (
- np.cos(
- data[
- (
- PartType_Gas_to_use,
- "particle_position_cylindrical_theta",
- )
- ][ind]
- )
- * v_rot_local
- * 1e5
- )
- return data.ds.arr(trans, "cm/s").in_base(
- data.ds.unit_system.name
- )
-
- pf.add_field(
- (PartType_Gas_to_use, "particle_local_rotational_velocity_y"),
- function=_particle_local_rotational_velocity_y,
- take_log=False,
- particle_type=True,
- units="cm/s",
- )
-
- def _particle_velocity_minus_local_rotational_velocity_squared(
- field, data
- ):
- return (
- (
- data[(PartType_Gas_to_use, "particle_velocity_x")]
- - data[
- (
- PartType_Gas_to_use,
- "particle_local_rotational_velocity_x",
- )
- ]
- )
- ** 2
- + (
- data[(PartType_Gas_to_use, "particle_velocity_y")]
- - data[
- (
- PartType_Gas_to_use,
- "particle_local_rotational_velocity_y",
- )
- ]
- )
- ** 2
- + (data[(PartType_Gas_to_use, "particle_velocity_z")]) ** 2
- )
-
- pf.add_field(
- (
- PartType_Gas_to_use,
- "particle_velocity_minus_local_rotational_velocity_squared",
- ),
- function=_particle_velocity_minus_local_rotational_velocity_squared,
- take_log=False,
- particle_type=True,
- units="cm**2/s**2",
- )
-
- p55 = ProfilePlot(
- sp,
- (PartType_Gas_to_use, "radius"),
- (
- PartType_Gas_to_use,
- "particle_velocity_minus_local_rotational_velocity_squared",
- ),
- weight_field=(PartType_Gas_to_use, MassType_to_use),
- n_bins=50,
- x_log=False,
- )
- p55.set_log("radius", False)
- p55.set_unit("radius", "kpc")
- p55.set_xlim(1e-3, 14)
- pos_disp_xs[time].append(p55.profiles[0].x.in_units("kpc").d)
- pos_disp_profiles[time].append(
- np.sqrt(
- p55.profiles[0][
- "particle_velocity_minus_local_rotational_velocity_squared"
- ]
- )
- .in_units("km/s")
- .d
- )
-
- # Add vertical velocity dispersion profile if requested
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
-
- def _velocity_z_squared(field, data):
- return (data[("gas", "velocity_z")]) ** 2
-
- pf.add_field(
- ("gas", "velocity_z_squared"),
- function=_velocity_z_squared,
- take_log=False,
- particle_type=False,
- units="cm**2/s**2",
- )
-
- p551 = ProfilePlot(
- sp_dense,
- ("index", "cylindrical_r"),
- ("gas", "velocity_z_squared"),
- weight_field=("gas", "cell_mass"),
- n_bins=50,
- x_log=False,
- )
- p551.set_log("cylindrical_r", False)
- p551.set_unit("cylindrical_r", "kpc")
- p551.set_xlim(1e-3, 14)
- pos_disp_vert_xs[time].append(p551.profiles[0].x.in_units("kpc").d)
- pos_disp_vert_profiles[time].append(
- np.sqrt(p551.profiles[0]["velocity_z_squared"])
- .in_units("km/s")
- .d
- )
- else:
-
- def _particle_velocity_z_squared(field, data):
- return (data[(PartType_Gas_to_use, "particle_velocity_z")]) ** 2
-
- pf.add_field(
- (PartType_Gas_to_use, "particle_velocity_z_squared"),
- function=_particle_velocity_z_squared,
- take_log=False,
- particle_type=True,
- units="cm**2/s**2",
- )
-
- p551 = ProfilePlot(
- sp,
- (PartType_Gas_to_use, "radius"),
- (PartType_Gas_to_use, "particle_velocity_z_squared"),
- weight_field=(PartType_Gas_to_use, MassType_to_use),
- n_bins=50,
- x_log=False,
- )
- p551.set_log("radius", False)
- p551.set_unit("radius", "kpc")
- p551.set_xlim(1e-3, 14)
- pos_disp_vert_xs[time].append(p551.profiles[0].x.in_units("kpc").d)
- pos_disp_vert_profiles[time].append(
- np.sqrt(p551.profiles[0]["particle_velocity_z_squared"])
- .in_units("km/s")
- .d
- )
-
- if add_nametag == 1:
- at = AnchoredText(
- "%s" % codes[code], loc=4, prop=dict(size=10), frameon=True
- )
- grid_pos_vel_PDF[time][code].axes.add_artist(at)
-
- # POSITION-VELOCITY PDF FOR NEW STARS
- if draw_star_pos_vel_PDF >= 1 and time != 0:
- sp = pf.sphere(center, (0.5 * figure_width, "kpc"))
- sp.set_field_parameter("normal", disk_normal_vector)
- pf.field_info[(PartType_Star_to_use, "radius")].take_log = False
- pf.field_info[
- (PartType_Star_to_use, "particle_velocity_cylindrical_theta")
- ].take_log = False
- pf.field_info[
- (PartType_Star_to_use, "particle_mass")
- ].output_units = (
- "code_mass"
- ) # this turned out to be crucial!; otherwise wrong output_unit 'g' is assumed in ParticlePhasePlot->create_profile in visualizaiton/particle_plots.py for ART-I/ENZO/RAMSES
- p41 = ParticlePhasePlot(
- sp,
- (PartType_Star_to_use, "radius"),
- (PartType_Star_to_use, "particle_velocity_cylindrical_theta"),
- (PartType_Star_to_use, "particle_mass"),
- deposition="ngp",
- weight_field=None,
- fontsize=12,
- x_bins=300,
- y_bins=300,
- )
- p41.set_unit("radius", "kpc")
- p41.set_unit("particle_velocity_cylindrical_theta", "km/s")
- p41.set_unit(
- "particle_mass", "Msun"
- ) # requires a change in set_unit in visualization/profile_plotter.py: remove self.plots[field].zmin, self.plots[field].zmax = (None, None)
- # p41.set_unit((PartType_Star_to_use, "particle_mass"), 'Msun') # Neither this nor above works without such change
- p41.set_zlim((PartType_Star_to_use, "particle_mass"), 1e3, 1e7)
- p41.set_cmap((PartType_Star_to_use, "particle_mass"), "algae")
-
- p41.set_colorbar_label(
- (PartType_Star_to_use, "particle_mass"),
- "Newly Formed Stellar Mass ($\mathrm{M}_{\odot}$)",
- )
- plot41 = p41.plots[(PartType_Star_to_use, "particle_mass")]
-
- p41.set_xlabel("Cylindrical Radius (kpc)")
- p41.set_ylabel("Rotational Velocity (km/s)")
- p41.set_xlim(0, 14)
- p41.set_ylim(-100, 350)
-
- plot41.figure = fig_star_pos_vel_PDF[time]
- plot41.axes = grid_star_pos_vel_PDF[time][code].axes
- if code == 0:
- plot41.cax = grid_star_pos_vel_PDF[time].cbar_axes[0]
- p41._setup_plots()
-
- # Add 1D profile line if requested
- if draw_star_pos_vel_PDF >= 2 and time != 0:
- p51 = ProfilePlot(
- sp,
- (PartType_Star_to_use, "radius"),
- (PartType_Star_to_use, "particle_velocity_cylindrical_theta"),
- weight_field=(PartType_Star_to_use, "particle_mass"),
- n_bins=50,
- x_log=False,
- )
- p51.set_log(
- (PartType_Star_to_use, "particle_velocity_cylindrical_theta"), False
- )
- p51.set_log((PartType_Star_to_use, "radius"), False)
- p51.set_unit("radius", "kpc")
- p51.set_xlim(1e-3, 14)
- p51.set_ylim("particle_velocity_cylindrical_theta", -100, 350)
- line = ln.Line2D(
- p51.profiles[0].x.in_units("kpc"),
- p51.profiles[0]["particle_velocity_cylindrical_theta"].in_units(
- "km/s"
- ),
- linestyle="-",
- linewidth=2,
- color="k",
- alpha=0.7,
- )
- star_pos_vel_xs[time].append(p51.profiles[0].x.in_units("kpc").d)
- star_pos_vel_profiles[time].append(
- p51.profiles[0]["particle_velocity_cylindrical_theta"]
- .in_units("km/s")
- .d
- )
- grid_star_pos_vel_PDF[time][code].axes.add_line(line)
-
- # Add velocity dispersion profile if requested
- if draw_star_pos_vel_PDF >= 3 and time != 0:
-
- def _particle_local_rotational_velocity_x(field, data):
- trans = np.zeros(
- data[
- (PartType_StarBeforeFiltered_to_use, "particle_velocity_x")
- ].shape
- )
- dr = 0.5 * (
- star_pos_vel_xs[time][code][1] - star_pos_vel_xs[time][code][0]
- )
- for radius, v_rot_local in zip(
- star_pos_vel_xs[time][code], star_pos_vel_profiles[time][code]
- ):
- ind = np.where(
- (
- data[
- (PartType_StarBeforeFiltered_to_use, "radius")
- ].in_units("kpc")
- >= (radius - dr)
- )
- & (
- data[
- (PartType_StarBeforeFiltered_to_use, "radius")
- ].in_units("kpc")
- < (radius + dr)
- )
- )
- trans[ind] = (
- -np.sin(
- data[
- (
- PartType_StarBeforeFiltered_to_use,
- "particle_position_cylindrical_theta",
- )
- ][ind]
- )
- * v_rot_local
- * 1e5
- )
- return data.ds.arr(trans, "cm/s").in_base(data.ds.unit_system.name)
-
- pf.add_field(
- (
- PartType_StarBeforeFiltered_to_use,
- "particle_local_rotational_velocity_x",
- ),
- function=_particle_local_rotational_velocity_x,
- take_log=False,
- particle_type=True,
- units="cm/s",
- )
-
- def _particle_local_rotational_velocity_y(field, data):
- trans = np.zeros(
- data[
- (PartType_StarBeforeFiltered_to_use, "particle_velocity_y")
- ].shape
- )
- dr = 0.5 * (
- star_pos_vel_xs[time][code][1] - star_pos_vel_xs[time][code][0]
- )
- for radius, v_rot_local in zip(
- star_pos_vel_xs[time][code], star_pos_vel_profiles[time][code]
- ):
- ind = np.where(
- (
- data[
- (PartType_StarBeforeFiltered_to_use, "radius")
- ].in_units("kpc")
- >= (radius - dr)
- )
- & (
- data[
- (PartType_StarBeforeFiltered_to_use, "radius")
- ].in_units("kpc")
- < (radius + dr)
- )
- )
- trans[ind] = (
- np.cos(
- data[
- (
- PartType_StarBeforeFiltered_to_use,
- "particle_position_cylindrical_theta",
- )
- ][ind]
- )
- * v_rot_local
- * 1e5
- )
- return data.ds.arr(trans, "cm/s").in_base(data.ds.unit_system.name)
-
- pf.add_field(
- (
- PartType_StarBeforeFiltered_to_use,
- "particle_local_rotational_velocity_y",
- ),
- function=_particle_local_rotational_velocity_y,
- take_log=False,
- particle_type=True,
- units="cm/s",
- )
-
- def _particle_velocity_minus_local_rotational_velocity_squared(
- field, data
- ):
- return (
- (
- data[
- (
- PartType_StarBeforeFiltered_to_use,
- "particle_velocity_x",
- )
- ]
- - data[
- (
- PartType_StarBeforeFiltered_to_use,
- "particle_local_rotational_velocity_x",
- )
- ]
- )
- ** 2
- + (
- data[
- (
- PartType_StarBeforeFiltered_to_use,
- "particle_velocity_y",
- )
- ]
- - data[
- (
- PartType_StarBeforeFiltered_to_use,
- "particle_local_rotational_velocity_y",
- )
- ]
- )
- ** 2
- + (
- data[
- (
- PartType_StarBeforeFiltered_to_use,
- "particle_velocity_z",
- )
- ]
- )
- ** 2
- )
-
- pf.add_field(
- (
- PartType_StarBeforeFiltered_to_use,
- "particle_velocity_minus_local_rotational_velocity_squared",
- ),
- function=_particle_velocity_minus_local_rotational_velocity_squared,
- take_log=False,
- particle_type=True,
- units="cm**2/s**2",
- )
- pf.add_particle_filter(
- PartType_Star_to_use
- ) # This is needed for a filtered particle type PartType_Star_to_use to work, because we have just created new particle fields.
-
- p515 = ProfilePlot(
- sp,
- (PartType_Star_to_use, "radius"),
- (
- PartType_Star_to_use,
- "particle_velocity_minus_local_rotational_velocity_squared",
- ),
- weight_field=(PartType_Star_to_use, "particle_mass"),
- n_bins=50,
- x_log=False,
- )
- p515.set_log((PartType_Star_to_use, "radius"), False)
- p515.set_unit("radius", "kpc")
- p515.set_xlim(1e-3, 14)
- star_pos_disp_xs[time].append(p515.profiles[0].x.in_units("kpc").d)
- star_pos_disp_profiles[time].append(
- np.sqrt(
- p515.profiles[0][
- "particle_velocity_minus_local_rotational_velocity_squared"
- ]
- )
- .in_units("km/s")
- .d
- )
-
- # Add vertical velocity dispersion profile if requested
- if draw_star_pos_vel_PDF >= 4 and time != 0:
-
- def _particle_velocity_z_squared(field, data):
- return (
- data[
- (PartType_StarBeforeFiltered_to_use, "particle_velocity_z")
- ]
- ) ** 2
-
- pf.add_field(
- (PartType_StarBeforeFiltered_to_use, "particle_velocity_z_squared"),
- function=_particle_velocity_z_squared,
- take_log=False,
- particle_type=True,
- units="cm**2/s**2",
- )
- pf.add_particle_filter(
- PartType_Star_to_use
- ) # This is needed for a filtered particle type PartType_Star_to_use to work, because we have just created new particle fields.
-
- p5151 = ProfilePlot(
- sp,
- (PartType_Star_to_use, "radius"),
- (PartType_Star_to_use, "particle_velocity_z_squared"),
- weight_field=(PartType_Star_to_use, "particle_mass"),
- n_bins=50,
- x_log=False,
- )
- p5151.set_log((PartType_Star_to_use, "radius"), False)
- p5151.set_unit("radius", "kpc")
- p5151.set_xlim(1e-3, 14)
- star_pos_disp_vert_xs[time].append(
- p5151.profiles[0].x.in_units("kpc").d
- )
- star_pos_disp_vert_profiles[time].append(
- np.sqrt(p5151.profiles[0]["particle_velocity_z_squared"])
- .in_units("km/s")
- .d
- )
-
- if add_nametag == 1:
- at = AnchoredText(
- "%s" % codes[code], loc=4, prop=dict(size=10), frameon=True
- )
- grid_star_pos_vel_PDF[time][code].axes.add_artist(at)
-
- # RADIUS-HEIGHT PDF
- if draw_rad_height_PDF >= 1:
- sp = pf.sphere(center, (0.5 * figure_width, "kpc"))
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
- if draw_rad_height_PDF == 1 or draw_rad_height_PDF == 2:
- p55 = PhasePlot(
- sp,
- ("index", "cylindrical_r"),
- ("index", "cylindrical_z_abs"),
- ("gas", "density"),
- weight_field=("gas", "cell_mass"),
- fontsize=12,
- x_bins=200,
- y_bins=200,
- )
- p55.set_zlim(("gas", "density"), 1e-26, 1e-21)
- p55.set_cmap(("gas", "density"), "algae")
- p55.set_log("cylindrical_r", False)
- p55.set_log("cylindrical_z_abs", False)
- p55.set_unit("cylindrical_r", "kpc")
- p55.set_unit("cylindrical_z_abs", "kpc")
- plot55 = p55.plots[("gas", "density")]
- elif draw_rad_height_PDF == 3:
- p55 = PhasePlot(
- sp,
- ("index", "cylindrical_r"),
- ("index", "cylindrical_z_abs"),
- ("gas", "density_minus_analytic"),
- weight_field=("gas", "cell_mass"),
- fontsize=12,
- x_bins=200,
- y_bins=200,
- )
- p55.set_zlim(("gas", "density_minus_analytic"), -1e-24, 1e-24)
- p55.set_cmap(("gas", "density_minus_analytic"), "algae")
- p55.set_log("cylindrical_r", False)
- p55.set_log("cylindrical_z_abs", False)
- p55.set_unit("cylindrical_r", "kpc")
- p55.set_unit("cylindrical_z_abs", "kpc")
- plot55 = p55.plots[("gas", "density_minus_analytic")]
- else:
- # Because ParticlePhasePlot doesn't work for these newly created fields for some reason, I will do the following trick.
- if draw_rad_height_PDF == 1 or draw_rad_height_PDF == 2:
- p55 = PhasePlot(
- sp,
- (PartType_Gas_to_use, "radius"),
- (PartType_Gas_to_use, "particle_position_cylindrical_z_abs"),
- (PartType_Gas_to_use, "density"),
- weight_field=(PartType_Gas_to_use, MassType_to_use),
- fontsize=12,
- x_bins=200,
- y_bins=200,
- )
- p55.set_zlim((PartType_Gas_to_use, "density"), 1e-26, 1e-21)
- p55.set_cmap((PartType_Gas_to_use, "density"), "algae")
- p55.set_log("radius", False)
- p55.set_log("particle_position_cylindrical_z_abs", False)
- p55.set_unit("radius", "kpc")
- p55.set_unit("particle_position_cylindrical_z_abs", "kpc")
- plot55 = p55.plots[(PartType_Gas_to_use, "density")]
- elif draw_rad_height_PDF == 3:
- p55 = PhasePlot(
- sp,
- (PartType_Gas_to_use, "radius"),
- (PartType_Gas_to_use, "particle_position_cylindrical_z_abs"),
- (PartType_Gas_to_use, "Density_2_minus_analytic"),
- weight_field=(PartType_Gas_to_use, MassType_to_use),
- fontsize=12,
- x_bins=200,
- y_bins=200,
- )
- p55.set_zlim(
- (PartType_Gas_to_use, "Density_2_minus_analytic"), -1e-24, 1e-24
- )
- p55.set_cmap(
- (PartType_Gas_to_use, "Density_2_minus_analytic"), "algae"
- )
- p55.set_log("radius", False)
- p55.set_log("particle_position_cylindrical_z_abs", False)
- p55.set_unit("radius", "kpc")
- p55.set_unit("particle_position_cylindrical_z_abs", "kpc")
- plot55 = p55.plots[
- (PartType_Gas_to_use, "Density_2_minus_analytic")
- ]
-
- p55.set_xlabel("Cylindrical Radius (kpc)")
- p55.set_ylabel("Vertical Height (kpc)")
- p55.set_xlim(0, 14)
- p55.set_ylim(0, 1.4)
-
- plot55.figure = fig_rad_height_PDF[time]
- plot55.axes = grid_rad_height_PDF[time][code].axes
- if code == 0:
- plot55.cax = grid_rad_height_PDF[time].cbar_axes[0]
- p55._setup_plots()
-
- # Add 1D profile line if requested
- if draw_rad_height_PDF == 2:
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
- p56 = ProfilePlot(
- sp,
- ("index", "cylindrical_r"),
- ("index", "cylindrical_z_abs"),
- weight_field=("gas", "cell_mass"),
- n_bins=50,
- x_log=False,
- )
- p56.set_log("cylindrical_r", False)
- p56.set_log("cylindrical_z_abs", False)
- p56.set_unit("cylindrical_r", "kpc")
- p56.set_unit("cylindrical_z_abs", "kpc")
- p56.set_xlim(0, 14)
- p56.set_ylim("cylindrical_z_abs", 0, 1.4)
- line = ln.Line2D(
- p56.profiles[0].x.in_units("kpc"),
- p56.profiles[0]["cylindrical_z_abs"].in_units("kpc"),
- linestyle="-",
- linewidth=2,
- color="k",
- alpha=0.7,
- )
- rad_height_xs[time].append(p56.profiles[0].x.in_units("kpc").d)
- rad_height_profiles[time].append(
- p56.profiles[0]["cylindrical_z_abs"].in_units("kpc").d
- )
- else:
- p56 = ProfilePlot(
- sp,
- (PartType_Gas_to_use, "radius"),
- (PartType_Gas_to_use, "particle_position_cylindrical_z_abs"),
- weight_field=(PartType_Gas_to_use, MassType_to_use),
- n_bins=50,
- x_log=False,
- )
- p56.set_log("radius", False)
- p56.set_log("particle_position_cylindrical_z_abs", False)
- p56.set_unit("radius", "kpc")
- p56.set_unit("particle_position_cylindrical_z_abs", "kpc")
- p56.set_xlim(0, 14)
- p56.set_ylim("particle_position_cylindrical_z_abs", 0, 1.4)
- line = ln.Line2D(
- p56.profiles[0].x.in_units("kpc"),
- p56.profiles[0]["particle_position_cylindrical_z_abs"].in_units(
- "kpc"
- ),
- linestyle="-",
- linewidth=2,
- color="k",
- alpha=0.7,
- )
- rad_height_xs[time].append(p56.profiles[0].x.in_units("kpc").d)
- rad_height_profiles[time].append(
- p56.profiles[0]["particle_position_cylindrical_z_abs"]
- .in_units("kpc")
- .d
- )
- grid_rad_height_PDF[time][code].axes.add_line(line)
-
- if add_nametag == 1:
- at = AnchoredText(
- "%s" % codes[code], loc=1, prop=dict(size=10), frameon=True
- )
- grid_rad_height_PDF[time][code].axes.add_artist(at)
-
- # DENSITY-TEMPERATURE-METALLICITY PDF
- if draw_metal_PDF == 1:
- sp = pf.sphere(center, (0.5 * figure_width, "kpc"))
- my_cmap2 = copy.copy(matplotlib.cm.get_cmap("algae"))
- my_cmap2.set_under("w")
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
- pf.field_info[("gas", "metallicity")].take_log = False
- p3 = PhasePlot(
- sp,
- ("gas", "density"),
- ("gas", "temperature"),
- ("gas", "metallicity"),
- weight_field=("gas", "cell_mass"),
- fontsize=12,
- x_bins=300,
- y_bins=300,
- )
- p3.set_zlim(("gas", "metallicity"), 0.01, 0.04)
- p3.set_cmap(("gas", "metallicity"), my_cmap2)
- p3.set_colorbar_label(
- ("gas", "metallicity"),
- "Metallicity (Mass-weighted average of mass fraction)",
- )
- plot3 = p3.plots[("gas", "metallicity")]
- else:
- # Because ParticlePhasePlot doesn't yet work for a log-log PDF for some reason, I will do the following trick.
- pf.field_info[(PartType_Gas_to_use, "Metallicity_2")].take_log = False
- p3 = ParticlePhasePlot(
- sp,
- (PartType_Gas_to_use, "density"),
- (PartType_Gas_to_use, "temperature"),
- (PartType_Gas_to_use, "Metallicity_2"),
- weight_field=(PartType_Gas_to_use, MassType_to_use),
- fontsize=12,
- x_bins=300,
- y_bins=300,
- )
- p3.set_zlim((PartType_Gas_to_use, "Metallicity_2"), 0.01, 0.04)
- p3.set_cmap((PartType_Gas_to_use, "Metallicity_2"), my_cmap2)
- plot3 = p3.plots[(PartType_Gas_to_use, "Metallicity_2")]
-
- p3.set_unit("density", "g/cm**3")
- p3.set_xlim(1e-29, 1e-21)
- p3.set_unit("temperature", "K")
- p3.set_ylim(10, 1e7)
- p3.set_log("density", True)
- p3.set_log("temperature", True)
-
- plot3.figure = fig_metal_PDF[time]
- plot3.axes = grid_metal_PDF[time][code].axes
- if code == 0:
- plot3.cax = grid_metal_PDF[time].cbar_axes[0]
- p3._setup_plots()
-
- if add_nametag == 1:
- at = AnchoredText(
- "%s" % codes[code], loc=3, prop=dict(size=10), frameon=True
- )
- grid_metal_PDF[time][code].axes.add_artist(at)
-
- # DENSITY DF (DISTRIBUTION FUNCTION)
- if draw_density_DF >= 1:
- sp = pf.sphere(center, (0.5 * figure_width, "kpc"))
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
- p6 = ProfilePlot(
- sp,
- ("gas", "density"),
- ("gas", "cell_mass"),
- weight_field=None,
- n_bins=50,
- x_log=True,
- accumulation=False,
- )
- # p6 = ProfilePlot(sp, ("gas", "density"), ("gas", "cell_mass"), weight_field=None, n_bins=50, x_log=True, accumulation=True)
- p6.set_log("cell_mass", True)
- p6.set_xlim(1e-29, 1e-21)
- density_DF_xs[time].append(p6.profiles[0].x.in_units("g/cm**3").d)
- density_DF_profiles[time].append(
- p6.profiles[0]["cell_mass"].in_units("Msun").d
- )
- else:
- # Because ParticleProfilePlot doesn't exist, I will do the following trick.
- p6 = ProfilePlot(
- sp,
- (PartType_Gas_to_use, "density"),
- (PartType_Gas_to_use, MassType_to_use),
- weight_field=None,
- n_bins=50,
- x_log=True,
- accumulation=False,
- )
- # p6 = ProfilePlot(sp, (PartType_Gas_to_use, "density"), (PartType_Gas_to_use, MassType_to_use), weight_field=None, n_bins=50, x_log=True, accumulation=True)
- p6.set_log((PartType_Gas_to_use, MassType_to_use), True)
- density_DF_xs[time].append(p6.profiles[0].x.in_units("g/cm**3").d)
- density_DF_profiles[time].append(
- p6.profiles[0][MassType_to_use].in_units("Msun").d
- )
-
- # Add difference plot between 1st and 2nd datasets, if requested
- if draw_density_DF == 2 and time != 0:
- if dataset_num == 2:
- pf_1st = load_dataset(
- 1, time, code, codes, filenames[0]
- ) # load 1st datasets
- v, cen = pf_1st.h.find_max(("gas", "density"))
- sp = pf_1st.sphere(cen, (30.0, "kpc"))
- cen2 = sp.quantities.center_of_mass(
- use_gas=True, use_particles=False
- ).in_units("kpc")
- sp2 = pf_1st.sphere(cen2, (1.0, "kpc"))
- cen3 = sp2.quantities.max_location(("gas", "density"))
- center_1st = pf_1st.arr(
- [cen3[1].d, cen3[2].d, cen3[3].d], "code_length"
- )
- if yt_version_pre_3_2_3 == 1:
- center_1st = pf_1st.arr(
- [cen3[2].d, cen3[3].d, cen3[4].d], "code_length"
- ) # for yt-3.2.3 or before
- sp_1st = pf_1st.sphere(center_1st, (0.5 * figure_width, "kpc"))
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
- p61 = ProfilePlot(
- sp_1st,
- ("gas", "density"),
- ("gas", "cell_mass"),
- weight_field=None,
- n_bins=50,
- x_log=True,
- accumulation=False,
- )
- p61.set_log("cell_mass", True)
- p61.set_xlim(1e-29, 1e-21)
- density_DF_1st_xs[time].append(
- p61.profiles[0].x.in_units("g/cm**3").d
- )
- density_DF_1st_profiles[time].append(
- p61.profiles[0]["cell_mass"].in_units("Msun").d
- )
- else:
-
- def _Density_2(field, data):
- return data[(PartType_Gas_to_use, "Density")].in_units(
- "g/cm**3"
- )
-
- pf_1st.add_field(
- (PartType_Gas_to_use, "density"),
- function=_Density_2,
- take_log=True,
- particle_type=False,
- display_name="Density",
- units="g/cm**3",
- )
-
- def _Mass_2(field, data):
- return data[
- (PartType_Gas_to_use, MassType_to_use)
- ].in_units("Msun")
-
- pf_1st.add_field(
- (PartType_Gas_to_use, MassType_to_use),
- function=_Mass_2,
- take_log=True,
- particle_type=False,
- display_name="Mass",
- units="Msun",
- )
- p61 = ProfilePlot(
- sp_1st,
- (PartType_Gas_to_use, "density"),
- (PartType_Gas_to_use, MassType_to_use),
- weight_field=None,
- n_bins=50,
- x_log=True,
- accumulation=False,
- )
- p61.set_log((PartType_Gas_to_use, MassType_to_use), True)
- p61.set_xlim(1e-29, 1e-21)
- density_DF_1st_xs[time].append(
- p61.profiles[0].x.in_units("g/cm**3").d
- )
- density_DF_1st_profiles[time].append(
- p61.profiles[0][MassType_to_use].in_units("Msun").d
- )
- else:
- print("This won't work; consider setting dataset_num to 2...")
- continue
-
- # CYLINDRICAL RADIUS DF + RADIALLY-BINNED GAS SURFACE DENSITY
- if draw_radius_DF == 1:
- sp = pf.sphere(center, (0.5 * figure_width, "kpc"))
- sp.set_field_parameter("normal", disk_normal_vector)
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
- p7 = ProfilePlot(
- sp,
- ("index", "cylindrical_r"),
- ("gas", "cell_mass"),
- weight_field=None,
- n_bins=50,
- x_log=False,
- accumulation=False,
- )
- p7.set_log("cell_mass", True)
- p7.set_log("cylindrical_r", False)
- p7.set_unit("cylindrical_r", "kpc")
- p7.set_xlim(1e-3, 15)
- radius_DF_xs[time].append(p7.profiles[0].x.in_units("kpc").d)
- radius_DF_profiles[time].append(
- p7.profiles[0]["cell_mass"].in_units("Msun").d
- )
- else:
-
- p7 = ProfilePlot(
- sp,
- (PartType_Gas_to_use, "radius"),
- (PartType_Gas_to_use, MassType_to_use),
- weight_field=None,
- n_bins=50,
- x_log=False,
- accumulation=False,
- )
- p7.set_log(MassType_to_use, True)
- p7.set_log("radius", False)
- p7.set_unit("radius", "kpc")
- p7.set_xlim(1e-3, 15)
- radius_DF_xs[time].append(p7.profiles[0].x.in_units("kpc").d)
- radius_DF_profiles[time].append(
- p7.profiles[0][MassType_to_use].in_units("Msun").d
- )
-
- # CYLINDRICAL RADIUS DF + RADIALLY-BINNED SURFACE DENSITY FOR NEW STARS
- if draw_star_radius_DF >= 1 and time != 0:
- sp = pf.sphere(center, (0.5 * figure_width, "kpc"))
- sp.set_field_parameter("normal", disk_normal_vector)
- pf.field_info[(PartType_Star_to_use, "particle_mass")].take_log = True
- pf.field_info[
- (PartType_Star_to_use, "particle_mass")
- ].output_units = (
- "code_mass"
- ) # this turned out to be crucial!; check output_units above
- pf.field_info[(PartType_Star_to_use, "radius")].take_log = False
- p71 = ProfilePlot(
- sp,
- (PartType_Star_to_use, "radius"),
- (PartType_Star_to_use, "particle_mass"),
- weight_field=None,
- n_bins=50,
- x_log=False,
- accumulation=False,
- )
- p71.set_unit("radius", "kpc")
- p71.set_xlim(1e-3, 15)
- star_radius_DF_xs[time].append(p71.profiles[0].x.in_units("kpc").d)
- star_radius_DF_profiles[time].append(
- p71.profiles[0]["particle_mass"].in_units("Msun").d
- )
-
- # Add RADIALLY-BINNED SFR SURFACE DENSITY PROFILE if requested (SFR estimated using stars younger than 20 Myrs old)
- if draw_star_radius_DF == 2 and time != 0:
-
- def _particle_mass_young_stars_star_radius_DF(field, data):
- trans = np.zeros(
- data[
- (PartType_StarBeforeFiltered_to_use, "particle_mass")
- ].shape
- )
- ind = np.where(
- data[
- (
- PartType_StarBeforeFiltered_to_use,
- FormationTimeType_to_use,
- )
- ].in_units("Myr")
- > (
- pf.current_time.in_units("Myr").d
- - young_star_cutoff_star_radius_DF
- )
- )
- trans[ind] = data[
- (PartType_StarBeforeFiltered_to_use, "particle_mass")
- ][ind].in_units("code_mass")
- return data.ds.arr(trans, "code_mass").in_base(
- data.ds.unit_system.name
- )
-
- pf.add_field(
- (
- PartType_StarBeforeFiltered_to_use,
- "particle_mass_young_stars_star_radius_DF",
- ),
- function=_particle_mass_young_stars_star_radius_DF,
- display_name="Young Stellar Mass",
- units="code_mass",
- particle_type=True,
- take_log=True,
- )
- pf.add_particle_filter(PartType_Star_to_use)
-
- pf.field_info[
- (PartType_Star_to_use, "particle_mass_young_stars_star_radius_DF")
- ].take_log = True
- pf.field_info[
- (PartType_Star_to_use, "particle_mass_young_stars_star_radius_DF")
- ].output_units = (
- "code_mass"
- ) # this turned out to be crucial!; check output_units above
- pf.field_info[(PartType_Star_to_use, "radius")].take_log = False
- p72 = ProfilePlot(
- sp,
- (PartType_Star_to_use, "radius"),
- (PartType_Star_to_use, "particle_mass_young_stars_star_radius_DF"),
- weight_field=None,
- n_bins=50,
- x_log=False,
- accumulation=False,
- )
- p72.set_unit("radius", "kpc")
- p72.set_xlim(1e-3, 15)
- sfr_radius_DF_xs[time].append(p72.profiles[0].x.in_units("kpc").d)
- sfr_radius_DF_profiles[time].append(
- p72.profiles[0]["particle_mass_young_stars_star_radius_DF"]
- .in_units("Msun")
- .d
- / young_star_cutoff_star_radius_DF
- / 1e6
- ) # in Msun/yr
-
- # VERTICAL HEIGHT DF + VERTICALLY-BINNED GAS SURFACE DENSITY
- if draw_height_DF == 1:
- sp = pf.sphere(center, (0.5 * figure_width, "kpc"))
- sp.set_field_parameter("normal", disk_normal_vector)
- if (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- ):
- p8 = ProfilePlot(
- sp,
- ("index", "cylindrical_z_abs"),
- ("gas", "cell_mass"),
- weight_field=None,
- n_bins=10,
- x_log=False,
- accumulation=False,
- )
- p8.set_log("cell_mass", True)
- p8.set_log("cylindrical_z_abs", False)
- p8.set_unit("cylindrical_z_abs", "kpc")
- p8.set_xlim(1e-3, 1.4)
- height_DF_xs[time].append(p8.profiles[0].x.in_units("kpc").d)
- height_DF_profiles[time].append(
- p8.profiles[0]["cell_mass"].in_units("Msun").d
- )
- else:
- p8 = ProfilePlot(
- sp,
- (PartType_Gas_to_use, "particle_position_cylindrical_z_abs"),
- (PartType_Gas_to_use, MassType_to_use),
- weight_field=None,
- n_bins=10,
- x_log=False,
- accumulation=False,
- )
- p8.set_log(MassType_to_use, True)
- p8.set_log("particle_position_cylindrical_z_abs", False)
- p8.set_unit("particle_position_cylindrical_z_abs", "kpc")
- p8.set_xlim(1e-3, 1.4)
- height_DF_xs[time].append(p8.profiles[0].x.in_units("kpc").d)
- height_DF_profiles[time].append(
- p8.profiles[0][MassType_to_use].in_units("Msun").d
- )
-
- # STAR FORMATION RATE + CUMULATIVE STELLAR MASS GROWTH IN TIME
- if draw_SFR >= 1 and time != 0:
- from yt.units.dimensions import (
- length,
- ) # Below are tricks to make StarFormationRate() work, particularly with "volume" argument, as it currently works only with comoving datasets
-
- pf.unit_registry.add(
- "pccm", pf.unit_registry.lut["pc"][0], length, "\\rm{pc}/(1+z)"
- )
- pf.hubble_constant = 0.71
- pf.omega_lambda = 0.73
- pf.omega_matter = 0.27
- pf.omega_curvature = 0.0
-
- sp = pf.sphere(center, (0.5 * figure_width, "kpc"))
- draw_SFR_mass = sp[(PartType_Star_to_use, "particle_mass")].in_units("Msun")
- draw_SFR_ct = sp[(PartType_Star_to_use, FormationTimeType_to_use)].in_units(
- "Myr"
- )
- sfr = StarFormationRate(
- pf,
- star_mass=draw_SFR_mass,
- star_creation_time=draw_SFR_ct,
- volume=sp.volume(),
- bins=25,
- ) # 25 bins hardcoded; see: http://yt-project.org/docs/dev/analyzing/analysis_modules/star_analysis.html
-
- sfr_ts[time].append(sfr.time.in_units("Myr")) # in Myr
- sfr_cum_masses[time].append(sfr.Msol_cumulative) # in Msun
- sfr_sfrs[time].append(sfr.Msol_yr) # in Msun/yr
-
- # DENSITY ALONG THE ORTHO-RAY OBJECT CUTTING THROUGH THE CENTER
- if draw_cut_through == 1:
- ray = pf.ortho_ray(
- 2,
- (center[0].in_units("code_length"), center[1].in_units("code_length")),
- ) # see: http://yt-project.org/doc/visualizing/manual_plotting.html#line-plots
- srt = np.argsort(ray["z"])
- cut_through_zs[time].append(
- np.array(ray["z"][srt].in_units("kpc").d - center[2].in_units("kpc").d)
- )
- cut_through_zvalues[time].append(
- np.array(ray[("gas", "density")][srt].in_units("g/cm**3").d)
- )
-
- ray = pf.ortho_ray(
- 0,
- (center[1].in_units("code_length"), center[2].in_units("code_length")),
- )
- srt = np.argsort(ray["x"])
- cut_through_xs[time].append(
- np.array(ray["x"][srt].in_units("kpc").d - center[2].in_units("kpc").d)
- )
- cut_through_xvalues[time].append(
- np.array(ray[("gas", "density")][srt].in_units("g/cm**3").d)
- )
-
- ####################################
- # POST-ANALYSIS STEPS #
- ####################################
-
- # SAVE FIGURES
- if draw_density_map == 1:
- fig_density_map[time].savefig(
- "Sigma_%dMyr" % times[time], bbox_inches="tight", pad_inches=0.03, dpi=300
- )
- if draw_temperature_map == 1:
- fig_temperature_map[time].savefig(
- "Temp_%dMyr" % times[time], bbox_inches="tight", pad_inches=0.03, dpi=300
- )
- if draw_cellsize_map == 1 or draw_cellsize_map == 3:
- fig_cellsize_map[time].savefig(
- "Cell_%dMyr" % times[time], bbox_inches="tight", pad_inches=0.03, dpi=300
- )
- if draw_cellsize_map == 2 or draw_cellsize_map == 3:
- fig_cellsize_map_2[time].savefig(
- "Resolution_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- if draw_elevation_map == 1:
- fig_elevation_map[time].savefig(
- "Elevation_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- if draw_metal_map >= 1:
- fig_metal_map[time].savefig(
- "Metal_%dMyr" % times[time], bbox_inches="tight", pad_inches=0.03, dpi=300
- )
- if draw_zvel_map == 1:
- fig_zvel_map[time].savefig(
- "z-vel_%dMyr" % times[time], bbox_inches="tight", pad_inches=0.03, dpi=300
- )
- if draw_star_map == 1 and time != 0:
- fig_star_map[time].savefig(
- "Star_%dMyr" % times[time], bbox_inches="tight", pad_inches=0.03, dpi=300
- )
- if draw_star_clump_stats >= 1 and time != 0:
- # if draw_star_clump_stats >= 2:
- # fig_star_map_2[time].savefig("Star_with_clumps_HOP_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
- # fig_star_map_3[time].savefig("Star_with_clumps_FOF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
- if draw_star_clump_stats == 3:
- pf_clump_ref = load_dataset(
- 2,
- 1,
- 0,
- ["GIZMO"],
- [
- [
- "dummy",
- "/lustre/ki/pfs/mornkr/080112_CHaRGe/pfs-hyades/AGORA-DISK-repository/Grackle+SF+ThermalFbck/GIZMO/AGORA_disk_second_ps1.5/snapshot_temp_100_last",
- ]
- ],
- )
- PartType_Star_to_use = "Stars"
- if (
- os.path.exists(
- "./halo_catalogs/hop_%s_%05d_high_floor/hop_%s_%05d_high_floor.0.h5"
- % ("GIZMO", 500, "GIZMO", 500)
- )
- == False
- ):
- hc_clump_ref = HaloCatalog(
- data_ds=pf_clump_ref,
- finder_method="hop",
- output_dir="./halo_catalogs/hop_%s_%05d_high_floor"
- % ("GIZMO", 500),
- finder_kwargs={
- "threshold": 2e8,
- "dm_only": False,
- "ptype": PartType_Star_to_use,
- },
- )
- hc_clump_ref.add_filter(
- "quantity_value", "particle_mass", ">", 2.6e6, "Msun"
- )
- hc_clump_ref.add_filter(
- "quantity_value", "particle_mass", "<", 2.6e8, "Msun"
- )
- hc_clump_ref.create()
- halo_ds_clump_ref = load(
- "./halo_catalogs/hop_%s_%05d_high_floor/hop_%s_%05d_high_floor.0.h5"
- % ("GIZMO", 500, "GIZMO", 500)
- )
- hc_clump_ref = HaloCatalog(
- halos_ds=halo_ds_clump_ref,
- output_dir="./halo_catalogs/hop_%s_%05d_high_floor" % ("GIZMO", 500),
- )
- hc_clump_ref.load()
-
- halo_ad_clump_ref = hc_clump_ref.halos_ds.all_data()
- star_clump_masses_hop_ref.append(
- np.log10(halo_ad_clump_ref["particle_mass"][:].in_units("Msun"))
- )
-
- if (
- os.path.exists(
- "./halo_catalogs/fof_%s_%05d_high_floor/fof_%s_%05d_high_floor.0.h5"
- % ("GIZMO", 500, "GIZMO", 500)
- )
- == False
- ):
- hc2_clump_ref = HaloCatalog(
- data_ds=pf_clump_ref,
- finder_method="fof",
- output_dir="./halo_catalogs/fof_%s_%05d_high_floor"
- % ("GIZMO", 500),
- finder_kwargs={
- "link": 0.0025,
- "dm_only": False,
- "ptype": PartType_Star_to_use,
- },
- )
- hc2_clump_ref.add_filter(
- "quantity_value", "particle_mass", ">", 2.6e6, "Msun"
- )
- hc2_clump_ref.add_filter(
- "quantity_value", "particle_mass", "<", 2.6e8, "Msun"
- )
- hc2_clump_ref.create()
- halo_ds2_clump_ref = load(
- "./halo_catalogs/fof_%s_%05d_high_floor/fof_%s_%05d_high_floor.0.h5"
- % ("GIZMO", 500, "GIZMO", 500)
- )
- hc2_clump_ref = HaloCatalog(
- halos_ds=halo_ds2_clump_ref,
- output_dir="./halo_catalogs/fof_%s_%05d_high_floor" % ("GIZMO", 500),
- )
- hc2_clump_ref.load()
-
- halo_ad2_clump_ref = hc2_clump_ref.halos_ds.all_data()
- star_clump_masses_fof_ref.append(
- np.log10(halo_ad2_clump_ref["particle_mass"][:].in_units("Msun"))
- )
- plt.clf()
- fig = plt.figure(figsize=(8, 4))
- gridspec.GridSpec(1, 2)
- for star_clump_stats_i in range(1, 3, 1):
- codes_plotted = []
- # plt.subplot(1,2,star_clump_stats_i, aspect=0.25)
- plt.subplot2grid((1, 2), (0, star_clump_stats_i - 1))
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- if (
- star_clump_stats_i == 1
- and (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- )
- ) or (
- star_clump_stats_i == 2
- and (
- codes[code] == "CHANGA"
- or codes[code] == "GASOLINE"
- or codes[code] == "GADGET-3"
- or codes[code] == "GEAR"
- or codes[code] == "GIZMO"
- or codes[code] == "SWIFT"
- )
- ):
- hist = np.histogram(
- star_clump_masses_hop[time][code], bins=10, range=(6.0, 8.5)
- )
- dbin = 0.5 * (hist[1][1] - hist[1][0])
- lines = plt.plot(
- hist[1][:-1] + dbin,
- hist[0],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- codes_plotted.append(codes[code])
- plt.xlim(6, 8.5)
- plt.ylim(-0.1, 15)
- plt.grid(True)
- plt.xlabel(
- "$\mathrm{log[Newly\ Formed\ Stellar\ Clump\ Mass\ (M_{\odot})]}$"
- )
- if star_clump_stats_i == 1:
- plt.ylabel("$\mathrm{Stellar\ Clump\ Counts, \ \ N_{clump}(M)}$")
- plt.legend(codes_plotted, loc=1, frameon=True, ncol=1, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="xx-small")
- # plt.savefig("star_clump_stats_HOP_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
- plt.clf()
- # Reiterate for cumulative plots
- fig = plt.figure(figsize=(8, 4))
- gridspec.GridSpec(1, 2)
- for star_clump_stats_i in range(1, 3, 1):
- codes_plotted = []
- # plt.subplot(1,2,star_clump_stats_i, aspect=0.15)
- plt.subplot2grid((1, 2), (0, star_clump_stats_i - 1))
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- if (
- star_clump_stats_i == 1
- and (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- )
- ) or (
- star_clump_stats_i == 2
- and (
- codes[code] == "CHANGA"
- or codes[code] == "GASOLINE"
- or codes[code] == "GADGET-3"
- or codes[code] == "GEAR"
- or codes[code] == "GIZMO"
- or codes[code] == "SWIFT"
- )
- ):
- hist = np.histogram(
- star_clump_masses_hop[time][code], bins=10, range=(6.0, 8.5)
- )
- dbin = 0.5 * (hist[1][1] - hist[1][0])
- lines = plt.plot(
- hist[1][:-1] + dbin,
- np.cumsum(hist[0][::-1])[::-1],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- codes_plotted.append(codes[code])
- if draw_star_clump_stats == 3 and star_clump_stats_i == 2:
- hist_clump_ref = np.histogram(
- star_clump_masses_hop_ref, bins=10, range=(6.0, 8.5)
- )
- dbin = 0.5 * (hist[1][1] - hist[1][0])
- lines = plt.plot(
- hist_clump_ref[1][:-1] + dbin,
- np.cumsum(hist_clump_ref[0][::-1])[::-1],
- color="k",
- linestyle="--",
- marker="*",
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- codes_plotted.append("GIZMO-ps2")
- plt.xlim(6, 8.5)
- # plt.ylim(-0.1, 30)
- plt.ylim(0.9, 50)
- plt.semilogy()
- plt.grid(True)
- plt.xlabel(
- "$\mathrm{log[Newly\ Formed\ Stellar\ Clump\ Mass\ (M_{\odot})]}$"
- )
- if star_clump_stats_i == 1:
- plt.ylabel("$\mathrm{Cumulative\ Clump\ Counts, \ \ N_{clump}(>M)}}$")
- plt.legend(
- codes_plotted,
- loc=1,
- frameon=True,
- ncol=1,
- fancybox=True,
- labelspacing=0.15,
- borderpad=0.3,
- )
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="xx-small")
- # plt.savefig("star_clump_stats_HOP_cumulative_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
- plt.clf()
-
- fig = plt.figure(figsize=(8, 4))
- gridspec.GridSpec(1, 2)
- for star_clump_stats_i in range(1, 3, 1):
- codes_plotted = []
- # plt.subplot(1,2,star_clump_stats_i, aspect=0.25)
- plt.subplot2grid((1, 2), (0, star_clump_stats_i - 1))
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- if (
- star_clump_stats_i == 1
- and (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- )
- ) or (
- star_clump_stats_i == 2
- and (
- codes[code] == "CHANGA"
- or codes[code] == "GASOLINE"
- or codes[code] == "GADGET-3"
- or codes[code] == "GEAR"
- or codes[code] == "GIZMO"
- or codes[code] == "SWIFT"
- )
- ):
- hist = np.histogram(
- star_clump_masses_fof[time][code], bins=10, range=(6.0, 8.5)
- )
- dbin = 0.5 * (hist[1][1] - hist[1][0])
- lines = plt.plot(
- hist[1][:-1] + dbin,
- hist[0],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- codes_plotted.append(codes[code])
- plt.xlim(6, 8.5)
- plt.ylim(-0.1, 15)
- plt.grid(True)
- plt.xlabel(
- "$\mathrm{log[Newly\ Formed\ Stellar\ Clump\ Mass\ (M_{\odot})]}$"
- )
- if star_clump_stats_i == 1:
- plt.ylabel("$\mathrm{Stellar\ Clump\ Counts, \ \ N_{clump}(M)}$")
- plt.legend(codes_plotted, loc=1, frameon=True, ncol=1, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="xx-small")
- # plt.savefig("star_clump_stats_FOF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
- plt.clf()
- # Reiterate for cumulative plots
- fig = plt.figure(figsize=(8, 4))
- gridspec.GridSpec(1, 2)
- for star_clump_stats_i in range(1, 3, 1):
- codes_plotted = []
- # plt.subplot(1,2,star_clump_stats_i, aspect=0.15)
- plt.subplot2grid((1, 2), (0, star_clump_stats_i - 1))
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- if (
- star_clump_stats_i == 1
- and (
- codes[code] == "ART-I"
- or codes[code] == "ART-II"
- or codes[code] == "ENZO"
- or codes[code] == "RAMSES"
- )
- ) or (
- star_clump_stats_i == 2
- and (
- codes[code] == "CHANGA"
- or codes[code] == "GASOLINE"
- or codes[code] == "GADGET-3"
- or codes[code] == "GEAR"
- or codes[code] == "GIZMO"
- or codes[code] == "SWIFT"
- )
- ):
- hist = np.histogram(
- star_clump_masses_fof[time][code], bins=10, range=(6.0, 8.5)
- )
- dbin = 0.5 * (hist[1][1] - hist[1][0])
- lines = plt.plot(
- hist[1][:-1] + dbin,
- np.cumsum(hist[0][::-1])[::-1],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- codes_plotted.append(codes[code])
- if draw_star_clump_stats == 3 and star_clump_stats_i == 2:
- hist_clump_ref = np.histogram(
- star_clump_masses_fof_ref, bins=10, range=(6.0, 8.5)
- )
- dbin = 0.5 * (hist[1][1] - hist[1][0])
- lines = plt.plot(
- hist_clump_ref[1][:-1] + dbin,
- np.cumsum(hist_clump_ref[0][::-1])[::-1],
- color="k",
- linestyle="--",
- marker="*",
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- codes_plotted.append("GIZMO-ps2")
- plt.xlim(6, 8.5)
- # plt.ylim(-0.1, 30)
- plt.ylim(0.9, 50)
- plt.semilogy()
- plt.grid(True)
- plt.xlabel(
- "$\mathrm{log[Newly\ Formed\ Stellar\ Clump\ Mass\ (M_{\odot})]}$"
- )
- if star_clump_stats_i == 1:
- plt.ylabel("$\mathrm{Cumulative\ Clump\ Counts, \ \ N_{clump}(>M)}}$")
- plt.legend(
- codes_plotted,
- loc=1,
- frameon=True,
- ncol=1,
- fancybox=True,
- labelspacing=0.15,
- borderpad=0.3,
- )
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="xx-small")
- plt.savefig(
- "star_clump_stats_FOF_cumulative_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
- if draw_SFR_map >= 1 and time != 0:
- fig_degr_density_map[time].savefig(
- "degraded_Sigma_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- fig_degr_sfr_map[time].savefig(
- "degraded_SFR_map_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- if draw_SFR_map == 2 and time != 0:
- # If requested, draw the local K-S plot using the degraded maps created above
- plt.clf()
- # plt.subplot(111)
- fig = plt.figure(figsize=(8, 7))
- Bigiel_surface_density = []
- Bigiel_sfr_surface_density = []
- for code in range(len(codes)):
- # Remove bins where SFR surface density is zero
- KS_x = np.array(surf_dens_SFR_map[time][code])
- KS_x[np.where(np.array(sfr_surf_dens_SFR_map[time][code]) < 1e-10)] = 0
- KS_x = np.log10(KS_x)
- KS_y = np.log10(np.array(sfr_surf_dens_SFR_map[time][code]))
- ind = np.logical_or(np.isinf(KS_x), np.isinf(KS_y))
- ind2 = np.logical_or(np.isnan(KS_x), np.isnan(KS_y))
- ind = np.logical_or(ind, ind2)
- KS_x = KS_x[~ind]
- KS_y = KS_y[~ind]
- if codes[code] == "SWIFT":
- plt.scatter(
- KS_x,
- KS_y,
- color="k",
- edgecolor="k",
- s=20,
- linewidth=0.7,
- marker=marker_names[code],
- alpha=0.1,
- )
- # Draw a 80% contour rather than scattering all the datapoints; see http://stackoverflow.com/questions/19390320/scatterplot-contours-in-matplotlib
- Gaussian_density_estimation_nbins = 20
- kernel = kde.gaussian_kde(np.vstack([KS_x, KS_y]))
- xi, yi = np.mgrid[
- KS_x.min() : KS_x.max() : Gaussian_density_estimation_nbins * 1j,
- KS_y.min() : KS_y.max() : Gaussian_density_estimation_nbins * 1j,
- ]
- zi = np.reshape(
- kernel(np.vstack([xi.flatten(), yi.flatten()])), xi.shape
- )
- contours = plt.contour(
- xi,
- yi,
- zi,
- np.array([0.2]),
- colors=color_names[code],
- linewidths=1.2,
- alpha=1.0,
- ) # 80% percentile contour
- # contours = plt.contour(xi, yi, zi, np.array([0.3173]), colors=color_names[code], linewidths=1.2, alpha=1.0) # 68.27% percentile contour (1-sigma)
- contours.collections[0].set_label(
- codes[code]
- ) # setting names for legend
- plt.clabel(
- contours, fmt=codes[code], inline=True, fontsize=7
- ) # setting labels
- plt.xlim(0, 3)
- plt.ylim(-4, 1)
- plt.grid(True)
- plt.xlabel("$\mathrm{log[Gas\ Surface\ Density\ (M_{\odot}/pc^2)]}$")
- plt.ylabel(
- "$\mathrm{log[Star\ Formation\ Rate\ Surface\ Density\ (M_{\odot}/yr/kpc^2)]}$"
- )
- plt.legend(
- loc=2, frameon=True, ncol=2, fancybox=True
- ) # note difference from others since we want the legends to list contours, not scattered points
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
- t = np.arange(-2, 5, 0.01)
- # for line in import_text("./Bigieletal2008_Fig8_Contour.txt", " "):
- # Bigiel_surface_density.append(float(line[0]) + np.log10(1.36)) # for the factor 1.36, see Section 2.3.1 of Bigiel et al. 2008
- # Bigiel_sfr_surface_density.append(float(line[1]))
- # plt.fill(Bigiel_surface_density, Bigiel_sfr_surface_density, fill=True, color='b', alpha = 0.1, hatch='\\') # contour by Bigiel et al. 2008 (Fig 8), or Feldmann et al. 2012 (Fig 1)
- plt.plot(
- t, 1.37 * t - 3.78, "k--", linewidth=2, alpha=0.7
- ) # observational fit by Kennicutt et al. 2007
- # plt.axhline(y = np.log10(8.593e4/young_star_cutoff_SFR_map/1e6/(float(aperture_size_SFR_map)/1000.)**2),
- # color='k', linestyle ='-.', linewidth=1, alpha=0.7) # SFR surface density cutoff due to limited star particle mass resolution
- plt.savefig(
- "K-S_local_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
- if draw_PDF >= 1:
- if draw_PDF == 3 and time != 0:
- # If requested, find a 1D profile line from a specific snapshot as a reference (in this case ENZO or CHANGA noSF run at 500 Myr), then we add it to all panels
- # pf_profile_ref = load_dataset(1, 1, 0, ['ENZO'], [['dummy', file_location[0]+'ENZO/DD0100/DD0100']])
- pf_profile_ref = load_dataset(
- 1,
- 1,
- 0,
- ["CHANGA"],
- [["dummy", file_location[0] + "CHANGA/disklow/disklow.000500"]],
- )
-
- def _Density_2(field, data):
- return data[(PartType_Gas_to_use, "Density")].in_units("g/cm**3")
-
- pf_profile_ref.add_field(
- (PartType_Gas_to_use, "density"),
- function=_Density_2,
- take_log=True,
- particle_type=False,
- display_name="Density",
- units="g/cm**3",
- )
-
- def _Temperature_2(field, data):
- return data[(PartType_Gas_to_use, "Temperature")].in_units("K")
-
- pf_profile_ref.add_field(
- (PartType_Gas_to_use, "temperature"),
- function=_Temperature_2,
- take_log=True,
- particle_type=False,
- display_name="Temperature",
- units="K",
- )
-
- def _Mass_2(field, data):
- return data[(PartType_Gas_to_use, MassType_to_use)].in_units("Msun")
-
- pf_profile_ref.add_field(
- (PartType_Gas_to_use, MassType_to_use),
- function=_Mass_2,
- take_log=True,
- particle_type=False,
- display_name="Mass",
- units="Msun",
- )
-
- v, cen = pf_profile_ref.h.find_max(("gas", "density"))
- sp = pf_profile_ref.sphere(cen, (30.0, "kpc"))
- cen2 = sp.quantities.center_of_mass(
- use_gas=True, use_particles=False
- ).in_units("kpc")
- sp2 = pf_profile_ref.sphere(cen2, (1.0, "kpc"))
- cen3 = sp2.quantities.max_location(("gas", "density"))
- center_profile_ref = pf_profile_ref.arr(
- [cen3[1].d, cen3[2].d, cen3[3].d], "code_length"
- )
- if yt_version_pre_3_2_3 == 1:
- center_profile_ref = pf_profile_ref.arr(
- [cen3[2].d, cen3[3].d, cen3[4].d], "code_length"
- )
- sp_profile_ref = pf_profile_ref.sphere(
- center_profile_ref, (0.5 * figure_width, "kpc")
- )
-
- # p35 = ProfilePlot(sp_profile_ref, ("gas", "density"), ("gas", "temperature"), weight_field=("gas", "cell_mass"), n_bins=30)
- # p35.set_xlim(1e-29, 1e-21)
- # p35.set_ylim("temperature", 10, 1e7)
- p35 = ProfilePlot(
- sp_profile_ref,
- (PartType_Gas_to_use, "density"),
- (PartType_Gas_to_use, "temperature"),
- weight_field=(PartType_Gas_to_use, MassType_to_use),
- n_bins=30,
- )
- p35.set_xlim(1e-26, 1e-21)
- p35.set_ylim("temperature", 10, 1e7)
- for code in range(len(codes)):
- # line_profile_ref = ln.Line2D(p35.profiles[0].x.in_units('g/cm**3'), p35.profiles[0]["temperature"].in_units('K'), linestyle="--", linewidth=2, color='k', alpha=0.7)
- line_profile_ref = ln.Line2D(
- p35.profiles[0].x.in_units("g/cm**3"),
- p35.profiles[0]["temperature"].in_units("K"),
- linestyle="--",
- linewidth=2,
- color="k",
- alpha=0.7,
- )
- grid_PDF[time][code].axes.add_line(line_profile_ref)
- fig_PDF[time].savefig(
- "PDF_%dMyr" % times[time], bbox_inches="tight", pad_inches=0.03, dpi=300
- )
- if draw_pos_vel_PDF >= 1:
- # fig_pos_vel_PDF[time].savefig("pos_vel_PDF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
- if draw_pos_vel_PDF >= 2 and time != 0:
- plt.clf()
- # plt.subplot(111, aspect=0.04)
- fig = plt.figure(figsize=(8, 8))
- gridspec.GridSpec(4, 1)
- plt.subplot2grid((4, 1), (0, 0), rowspan=3)
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- lines = plt.plot(
- pos_vel_xs[time][code],
- pos_vel_profiles[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(0, 270)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel("$\mathrm{Rotational\ Velocity\ (km/s)}$", fontsize="large")
- plt.legend(codes, loc=4, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
-
- plt.subplot2grid((4, 1), (3, 0), rowspan=1)
- ave_profiles = np.mean(np.array(pos_vel_profiles[time]), axis=0)
- for code in range(len(codes)):
- lines = plt.plot(
- pos_vel_xs[time][code],
- (pos_vel_profiles[time][code] - ave_profiles) / ave_profiles,
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(-0.5, 0.5)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Residual,}\/(\mathrm{v} - \mathrm{\bar v})/\mathrm{\bar v}$",
- fontsize="large",
- )
- if add_mean_fractional_dispersion == 1:
- mean_fractional_dispersion = (
- np.std(np.array(pos_vel_profiles[time]), axis=0) / ave_profiles
- )[
- (mean_dispersion_radius_range[0] < pos_vel_xs[time][0])
- & (pos_vel_xs[time][0] < mean_dispersion_radius_range[1])
- ].mean()
- mean_fractional_dispersion_in_dex = np.log10(
- 1.0 + mean_fractional_dispersion
- )
- plt.annotate(
- "Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) "
- % (
- mean_dispersion_radius_range[0],
- mean_dispersion_radius_range[1],
- mean_fractional_dispersion,
- mean_fractional_dispersion_in_dex,
- ),
- xy=(0.35, 0.08),
- size=10,
- xycoords="axes fraction",
- )
- # plt.text(12, 0.3, "mean dispersion = %.3f or %.3f dex " % (mean_fractional_dispersion, mean_fractional_dispersion_in_dex), ha="center", va="center", size=8, bbox=dict(boxstyle="round", fc="w", ec="0.5", alpha=0.9))
- plt.savefig(
- "pos_vel_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
- if draw_pos_vel_PDF >= 3 and time != 0:
- plt.clf()
- # plt.subplot(111, aspect=0.064)
- fig = plt.figure(figsize=(8, 10))
- gridspec.GridSpec(5, 1)
- plt.subplot2grid((5, 1), (0, 0), rowspan=3)
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- lines = plt.plot(
- pos_disp_xs[time][code],
- pos_disp_profiles[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(0, 170)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel("$\mathrm{Velocity\ Dispersion\ (km/s)}$", fontsize="large")
- plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
-
- plt.subplot2grid((5, 1), (3, 0), rowspan=1)
- ave_profiles = np.mean(np.array(pos_disp_profiles[time]), axis=0)
- for code in range(len(codes)):
- lines = plt.plot(
- pos_disp_xs[time][code],
- (pos_disp_profiles[time][code] - ave_profiles) / ave_profiles,
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(-1.0, 1.0)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Residual,}\/(\mathrm{\sigma} - \mathrm{\bar \sigma})/\mathrm{\bar \sigma}$",
- fontsize="large",
- )
- if add_mean_fractional_dispersion == 1:
- mean_fractional_dispersion = (
- np.std(np.array(pos_disp_profiles[time]), axis=0) / ave_profiles
- )[
- (mean_dispersion_radius_range[0] < pos_disp_xs[time][0])
- & (pos_disp_xs[time][0] < mean_dispersion_radius_range[1])
- ].mean()
- mean_fractional_dispersion_in_dex = np.log10(
- 1.0 + mean_fractional_dispersion
- )
- plt.annotate(
- "Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) "
- % (
- mean_dispersion_radius_range[0],
- mean_dispersion_radius_range[1],
- mean_fractional_dispersion,
- mean_fractional_dispersion_in_dex,
- ),
- xy=(0.35, 0.08),
- size=10,
- xycoords="axes fraction",
- )
- plt.savefig(
- "pos_disp_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.subplot2grid((5, 1), (4, 0), rowspan=1)
- for code in range(len(codes)):
- lines = plt.plot(
- pos_disp_xs[time][code],
- pos_disp_vert_profiles[time][code] / pos_disp_profiles[time][code],
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(0.0, 1.0)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Ratio,}\/\mathrm{\sigma}_{\perp}/\mathrm{\sigma}$",
- fontsize="large",
- )
- plt.savefig(
- "pos_disp_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
- if draw_pos_vel_PDF == 4 and time != 0:
- plt.clf()
- # plt.subplot(111, aspect=0.064)
- fig = plt.figure(figsize=(8, 8))
- gridspec.GridSpec(4, 1)
- plt.subplot2grid((4, 1), (0, 0), rowspan=3)
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- lines = plt.plot(
- pos_disp_vert_xs[time][code],
- pos_disp_vert_profiles[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(0, 170)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- "$\mathrm{Vertical\ Velocity\ Dispersion\ (km/s)}$", fontsize="large"
- )
- plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
-
- plt.subplot2grid((4, 1), (3, 0), rowspan=1)
- ave_profiles = np.mean(np.array(pos_disp_vert_profiles[time]), axis=0)
- for code in range(len(codes)):
- lines = plt.plot(
- pos_disp_vert_xs[time][code],
- (pos_disp_vert_profiles[time][code] - ave_profiles) / ave_profiles,
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(-1.0, 1.0)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Residual,}\/(\mathrm{\sigma} - \mathrm{\bar \sigma})/\mathrm{\bar \sigma}$",
- fontsize="large",
- )
- if add_mean_fractional_dispersion == 1:
- mean_fractional_dispersion = (
- np.std(np.array(pos_disp_vert_profiles[time]), axis=0)
- / ave_profiles
- )[
- (mean_dispersion_radius_range[0] < pos_disp_vert_xs[time][0])
- & (pos_disp_vert_xs[time][0] < mean_dispersion_radius_range[1])
- ].mean()
- mean_fractional_dispersion_in_dex = np.log10(
- 1.0 + mean_fractional_dispersion
- )
- plt.annotate(
- "Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) "
- % (
- mean_dispersion_radius_range[0],
- mean_dispersion_radius_range[1],
- mean_fractional_dispersion,
- mean_fractional_dispersion_in_dex,
- ),
- xy=(0.35, 0.08),
- size=10,
- xycoords="axes fraction",
- )
- # plt.savefig("pos_disp_vert_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
- plt.clf()
- if draw_star_pos_vel_PDF >= 1 and time != 0:
- # fig_star_pos_vel_PDF[time].savefig("star_pos_vel_PDF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
- if draw_star_pos_vel_PDF >= 2 and time != 0:
- plt.clf()
- # plt.subplot(111, aspect=0.04)
- fig = plt.figure(figsize=(8, 8))
- gridspec.GridSpec(4, 1)
- plt.subplot2grid((4, 1), (0, 0), rowspan=3)
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- lines = plt.plot(
- star_pos_vel_xs[time][code],
- star_pos_vel_profiles[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(0, 270)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel("$\mathrm{Rotational\ Velocity\ (km/s)}$", fontsize="large")
- plt.legend(codes, loc=4, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
-
- plt.subplot2grid((4, 1), (3, 0), rowspan=1)
- ave_profiles = np.mean(np.array(star_pos_vel_profiles[time]), axis=0)
- for code in range(len(codes)):
- lines = plt.plot(
- star_pos_vel_xs[time][code],
- (star_pos_vel_profiles[time][code] - ave_profiles) / ave_profiles,
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(-0.5, 0.5)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Residual,}\/(\mathrm{v} - \mathrm{\bar v})/\mathrm{\bar v}$",
- fontsize="large",
- )
- if add_mean_fractional_dispersion == 1:
- mean_fractional_dispersion = (
- np.std(np.array(star_pos_vel_profiles[time]), axis=0) / ave_profiles
- )[
- (mean_dispersion_radius_range[0] < star_pos_vel_xs[time][0])
- & (star_pos_vel_xs[time][0] < mean_dispersion_radius_range[1])
- ].mean()
- mean_fractional_dispersion_in_dex = np.log10(
- 1.0 + mean_fractional_dispersion
- )
- plt.annotate(
- "Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) "
- % (
- mean_dispersion_radius_range[0],
- mean_dispersion_radius_range[1],
- mean_fractional_dispersion,
- mean_fractional_dispersion_in_dex,
- ),
- xy=(0.35, 0.08),
- size=10,
- xycoords="axes fraction",
- )
- plt.savefig(
- "star_pos_vel_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
- if draw_star_pos_vel_PDF >= 3 and time != 0:
- plt.clf()
- # plt.subplot(111, aspect=0.064)
- fig = plt.figure(figsize=(8, 10))
- gridspec.GridSpec(5, 1)
- plt.subplot2grid((5, 1), (0, 0), rowspan=3)
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- lines = plt.plot(
- star_pos_disp_xs[time][code],
- star_pos_disp_profiles[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(0, 170)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel("$\mathrm{Velocity\ Dispersion\ (km/s)}$", fontsize="large")
- plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
-
- plt.subplot2grid((5, 1), (3, 0), rowspan=1)
- ave_profiles = np.mean(np.array(star_pos_disp_profiles[time]), axis=0)
- for code in range(len(codes)):
- lines = plt.plot(
- star_pos_disp_xs[time][code],
- (star_pos_disp_profiles[time][code] - ave_profiles) / ave_profiles,
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(-1.0, 1.0)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Residual,}\/(\mathrm{\sigma} - \mathrm{\bar \sigma})/\mathrm{\bar \sigma}$",
- fontsize="large",
- )
- if add_mean_fractional_dispersion == 1:
- mean_fractional_dispersion = (
- np.std(np.array(star_pos_disp_profiles[time]), axis=0)
- / ave_profiles
- )[
- (mean_dispersion_radius_range[0] < star_pos_disp_xs[time][0])
- & (star_pos_disp_xs[time][0] < mean_dispersion_radius_range[1])
- ].mean()
- mean_fractional_dispersion_in_dex = np.log10(
- 1.0 + mean_fractional_dispersion
- )
- plt.annotate(
- "Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) "
- % (
- mean_dispersion_radius_range[0],
- mean_dispersion_radius_range[1],
- mean_fractional_dispersion,
- mean_fractional_dispersion_in_dex,
- ),
- xy=(0.35, 0.08),
- size=10,
- xycoords="axes fraction",
- )
- plt.savefig(
- "star_pos_disp_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.subplot2grid((5, 1), (4, 0), rowspan=1)
- for code in range(len(codes)):
- lines = plt.plot(
- star_pos_disp_xs[time][code],
- star_pos_disp_vert_profiles[time][code]
- / star_pos_disp_profiles[time][code],
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(0.0, 1.0)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Ratio,}\/\mathrm{\sigma}_{\perp}/\mathrm{\sigma}$",
- fontsize="large",
- )
- plt.savefig(
- "star_pos_disp_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
- if draw_star_pos_vel_PDF == 4 and time != 0:
- plt.clf()
- # plt.subplot(111, aspect=0.064)
- fig = plt.figure(figsize=(8, 8))
- gridspec.GridSpec(4, 1)
- plt.subplot2grid((4, 1), (0, 0), rowspan=3)
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- lines = plt.plot(
- star_pos_disp_vert_xs[time][code],
- star_pos_disp_vert_profiles[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(0, 170)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- "$\mathrm{Vertical\ Velocity\ Dispersion\ (km/s)}$", fontsize="large"
- )
- plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
-
- plt.subplot2grid((4, 1), (3, 0), rowspan=1)
- ave_profiles = np.mean(np.array(star_pos_disp_vert_profiles[time]), axis=0)
- for code in range(len(codes)):
- lines = plt.plot(
- star_pos_disp_vert_xs[time][code],
- (star_pos_disp_vert_profiles[time][code] - ave_profiles)
- / ave_profiles,
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(-1.0, 1.0)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Residual,}\/(\mathrm{\sigma} - \mathrm{\bar \sigma})/\mathrm{\bar \sigma}$",
- fontsize="large",
- )
- if add_mean_fractional_dispersion == 1:
- mean_fractional_dispersion = (
- np.std(np.array(star_pos_disp_vert_profiles[time]), axis=0)
- / ave_profiles
- )[
- (mean_dispersion_radius_range[0] < star_pos_disp_vert_xs[time][0])
- & (star_pos_disp_vert_xs[time][0] < mean_dispersion_radius_range[1])
- ].mean()
- mean_fractional_dispersion_in_dex = np.log10(
- 1.0 + mean_fractional_dispersion
- )
- plt.annotate(
- "Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) "
- % (
- mean_dispersion_radius_range[0],
- mean_dispersion_radius_range[1],
- mean_fractional_dispersion,
- mean_fractional_dispersion_in_dex,
- ),
- xy=(0.35, 0.08),
- size=10,
- xycoords="axes fraction",
- )
- # plt.savefig("star_pos_disp_vert_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
- plt.clf()
- if draw_rad_height_PDF >= 1:
- # fig_rad_height_PDF[time].savefig("rad_height_PDF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
- if draw_rad_height_PDF == 2 and time != 0:
- plt.clf()
- # plt.subplot(111, aspect=18)
- fig = plt.figure(figsize=(8, 8))
- gridspec.GridSpec(4, 1)
- plt.subplot2grid((4, 1), (0, 0), rowspan=3)
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- lines = plt.plot(
- rad_height_xs[time][code],
- rad_height_profiles[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(0, 0.45)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel("$\mathrm{Average\ Vertical\ Height\ (kpc)}$", fontsize="large")
- plt.legend(codes, loc=2, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
-
- plt.subplot2grid((4, 1), (3, 0), rowspan=1)
- ave_profiles = np.mean(np.array(rad_height_profiles[time]), axis=0)
- for code in range(len(codes)):
- lines = plt.plot(
- rad_height_xs[time][code],
- (rad_height_profiles[time][code] - ave_profiles) / ave_profiles,
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(-1.0, 1.0)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Residual,}\/(\mathrm{h} - \mathrm{\bar h})/\mathrm{\bar h}$",
- fontsize="large",
- )
- if add_mean_fractional_dispersion == 1:
- mean_fractional_dispersion = (
- np.std(np.array(rad_height_profiles[time]), axis=0) / ave_profiles
- )[
- (mean_dispersion_radius_range[0] < rad_height_xs[time][0])
- & (rad_height_xs[time][0] < mean_dispersion_radius_range[1])
- ].mean()
- mean_fractional_dispersion_in_dex = np.log10(
- 1.0 + mean_fractional_dispersion
- )
- plt.annotate(
- "Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) "
- % (
- mean_dispersion_radius_range[0],
- mean_dispersion_radius_range[1],
- mean_fractional_dispersion,
- mean_fractional_dispersion_in_dex,
- ),
- xy=(0.35, 0.08),
- size=10,
- xycoords="axes fraction",
- )
- plt.savefig(
- "rad_height_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
- if draw_metal_PDF == 1:
- fig_metal_PDF[time].savefig(
- "metal_PDF_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- if draw_density_DF >= 1:
- # plt.clf()
- # plt.subplot(111, aspect=1)
- fig = plt.figure(figsize=(8, 8))
- gridspec.GridSpec(4, 1)
- plt.subplot2grid((4, 1), (0, 0), rowspan=3)
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- lines = plt.plot(
- density_DF_xs[time][code],
- density_DF_profiles[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.semilogx()
- plt.semilogy()
- plt.xlim(1e-28, 1e-21)
- plt.ylim(1e4, 1e9) # accumulation=False
- # plt.ylim(1e5, 2e10) # accumulation=True
- if time != 0 and dataset_num == 2:
- plt.axvline(
- x=1.67e-24 * 10, color="k", linestyle="--", linewidth=2, alpha=0.7
- ) # SF threshold density
- plt.grid(True)
- plt.xlabel("$\mathrm{Density\ (g/cm^3)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Mass,}\/\mathrm{d}M\mathrm{/dlog}\/\mathrm{\rho}\/\mathrm{(M_{\odot})}$",
- fontsize="large",
- )
- plt.legend(codes, loc=4, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
-
- plt.subplot2grid((4, 1), (3, 0), rowspan=1)
- ave_profiles = np.mean(np.array(density_DF_profiles[time]), axis=0)
- for code in range(len(codes)):
- lines = plt.plot(
- density_DF_xs[time][code],
- (density_DF_profiles[time][code] - ave_profiles) / ave_profiles,
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.semilogx()
- plt.xlim(1e-28, 1e-21)
- plt.ylim(-1.0, 3.0)
- plt.grid(True)
- plt.locator_params(axis="y", nbins=4)
- plt.xlabel("$\mathrm{Density\ (g/cm^3)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Residual,}\/(\mathrm{m} - \mathrm{\bar m})/\mathrm{\bar m}$",
- fontsize="large",
- )
- if add_mean_fractional_dispersion == 1:
- mean_fractional_dispersion = (
- np.std(np.array(density_DF_profiles[time]), axis=0) / ave_profiles
- )[
- (mean_dispersion_density_range[0] < density_DF_xs[time][0])
- & (density_DF_xs[time][0] < mean_dispersion_density_range[1])
- ].mean()
- mean_fractional_dispersion_in_dex = np.log10(
- 1.0 + mean_fractional_dispersion
- )
- plt.annotate(
- "Mean fractional dispersion for %.1e < rho < %.1e = %.3f (%.3f dex) "
- % (
- mean_dispersion_density_range[0],
- mean_dispersion_density_range[1],
- mean_fractional_dispersion,
- mean_fractional_dispersion_in_dex,
- ),
- xy=(0.27, 0.84),
- size=10,
- xycoords="axes fraction",
- )
- plt.savefig(
- "density_DF_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
- if draw_density_DF == 2 and time != 0 and dataset_num == 2:
- fig = plt.figure(figsize=(8, 8))
- gridspec.GridSpec(16, 1)
- plt.subplot2grid((16, 1), (0, 0), rowspan=9)
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- lines = plt.plot(
- density_DF_xs[time][code],
- (
- density_DF_profiles[time][code]
- - density_DF_1st_profiles[time][code]
- )
- / 1e8,
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xscale("log")
- plt.xlim(1e-28, 1e-21)
- plt.ylim(-7, 3)
- plt.axvline(
- x=1.67e-24 * 10, color="k", linestyle="--", linewidth=2, alpha=0.7
- ) # SF threshold density
- plt.grid(True)
- plt.xlabel("$\mathrm{Density\ (g/cm^3)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Mass\ Change,}\/\mathrm{\Delta}(\mathrm{d}M\mathrm{/dlog}\/\mathrm{\rho})\/\mathrm{(10^8 M_{\odot})}$",
- fontsize="large",
- )
- plt.legend(codes, loc=3, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
-
- plt.subplot2grid((16, 1), (9, 0), rowspan=7)
- for code in range(len(codes)):
- lines = plt.plot(
- density_DF_xs[time][code],
- (
- density_DF_profiles[time][code]
- - density_DF_1st_profiles[time][code]
- )
- / 1e8,
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlabel("$\mathrm{Density\ (g/cm^3)}$", fontsize="large")
- plt.ylabel(
- r"$\mathtt{symlog}[\/\mathrm{\Delta}(\mathrm{d}M\mathrm{/dlog}\/\mathrm{\rho})\/\mathrm{(10^8 M_{\odot})}\/]$",
- fontsize="large",
- )
- plt.xscale("log")
- plt.yscale("symlog", linthreshy=0.01)
- plt.xlim(1e-28, 1e-21)
- plt.ylim(-10, 10)
- plt.axvline(
- x=1.67e-24 * 10, color="k", linestyle="--", linewidth=2, alpha=0.7
- ) # SF threshold density
- plt.grid(True)
- plt.savefig(
- "density_DF_change_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
- if draw_radius_DF == 1:
- plt.clf()
- # plt.subplot(111, aspect=1)
- fig = plt.figure(figsize=(8, 6))
- for code in range(len(codes)):
- lines = plt.plot(
- radius_DF_xs[time][code],
- np.add.accumulate(radius_DF_profiles[time][code]),
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.semilogy()
- plt.xlim(0, 14)
- plt.ylim(1e7, 2e10)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel("$\mathrm{Mass\ (M_{\odot})}$", fontsize="large")
- plt.legend(codes, loc=4, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
- # plt.savefig("radius_DF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
-
- plt.clf()
- # plt.subplot(111, aspect=1)
- fig = plt.figure(figsize=(8, 8))
- gridspec.GridSpec(4, 1)
- plt.subplot2grid((4, 1), (0, 0), rowspan=3)
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- temp = []
- dr = 0.5 * (
- radius_DF_xs[time][code][1] - radius_DF_xs[time][code][0]
- ) # Here we assume that ProfilePlot was made with linearly binned radius_DF_xs
- for radius in range(len(radius_DF_profiles[time][code])):
- surface_area = np.pi * (
- ((radius_DF_xs[time][code][radius] + dr) * 1e3) ** 2
- - ((radius_DF_xs[time][code][radius] - dr) * 1e3) ** 2
- )
- temp.append(radius_DF_profiles[time][code][radius] / surface_area)
- surface_density[time].append(temp)
- lines = plt.plot(
- radius_DF_xs[time][code],
- surface_density[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.semilogy()
- plt.xlim(0, 14)
- plt.ylim(1e-1, 2e3)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel("$\mathrm{Surface\ Density\ (M_{\odot}/pc^2)}$", fontsize="large")
- plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
-
- plt.subplot2grid((4, 1), (3, 0), rowspan=1)
- ave_profiles = np.mean(np.array(surface_density[time]), axis=0)
- for code in range(len(codes)):
- lines = plt.plot(
- radius_DF_xs[time][code],
- (surface_density[time][code] - ave_profiles) / ave_profiles,
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(-1.0, 1.0)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Residual,}\/(\mathrm{\Sigma} - \mathrm{\bar \Sigma})/\mathrm{\bar \Sigma}$",
- fontsize="large",
- )
- if add_mean_fractional_dispersion == 1:
- mean_fractional_dispersion = (
- np.std(np.array(surface_density[time]), axis=0) / ave_profiles
- )[
- (mean_dispersion_radius_range[0] < radius_DF_xs[time][0])
- & (radius_DF_xs[time][0] < mean_dispersion_radius_range[1])
- ].mean()
- mean_fractional_dispersion_in_dex = np.log10(
- 1.0 + mean_fractional_dispersion
- )
- plt.annotate(
- "Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) "
- % (
- mean_dispersion_radius_range[0],
- mean_dispersion_radius_range[1],
- mean_fractional_dispersion,
- mean_fractional_dispersion_in_dex,
- ),
- xy=(0.35, 0.08),
- size=10,
- xycoords="axes fraction",
- )
- plt.savefig(
- "gas_surface_density_radial_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
- if draw_star_radius_DF >= 1 and time != 0:
- plt.clf()
- # plt.subplot(111, aspect=1)
- fig = plt.figure(figsize=(8, 6))
- for code in range(len(codes)):
- lines = plt.plot(
- star_radius_DF_xs[time][code],
- np.add.accumulate(star_radius_DF_profiles[time][code]),
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.semilogy()
- plt.xlim(0, 14)
- plt.ylim(1e7, 2e10)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- "$\mathrm{Newly\ Formed\ Stellar\ Mass\ (M_{\odot})}$", fontsize="large"
- )
- plt.legend(codes, loc=4, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
- # plt.savefig("star_radius_DF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
-
- plt.clf()
- # plt.subplot(111, aspect=1)
- fig = plt.figure(figsize=(8, 8))
- gridspec.GridSpec(4, 1)
- plt.subplot2grid((4, 1), (0, 0), rowspan=3)
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- temp = []
- dr = 0.5 * (
- star_radius_DF_xs[time][code][1] - star_radius_DF_xs[time][code][0]
- )
- for radius in range(len(star_radius_DF_profiles[time][code])):
- surface_area = np.pi * (
- ((star_radius_DF_xs[time][code][radius] + dr) * 1e3) ** 2
- - ((star_radius_DF_xs[time][code][radius] - dr) * 1e3) ** 2
- )
- temp.append(star_radius_DF_profiles[time][code][radius] / surface_area)
- star_surface_density[time].append(temp)
- lines = plt.plot(
- star_radius_DF_xs[time][code],
- star_surface_density[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.semilogy()
- plt.xlim(0, 14)
- plt.ylim(1e-1, 2e3)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- "$\mathrm{Newly\ Formed\ Stellar\ Surface\ Density\ (M_{\odot}/pc^2)}$",
- fontsize="large",
- )
- plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
-
- plt.subplot2grid((4, 1), (3, 0), rowspan=1)
- ave_profiles = np.mean(np.array(star_surface_density[time]), axis=0)
- for code in range(len(codes)):
- lines = plt.plot(
- star_radius_DF_xs[time][code],
- (star_surface_density[time][code] - ave_profiles) / ave_profiles,
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(-1.0, 3.0)
- plt.grid(True)
- plt.locator_params(axis="y", nbins=4)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Residual,}\/(\mathrm{\Sigma} - \mathrm{\bar \Sigma})/\mathrm{\bar \Sigma}$",
- fontsize="large",
- )
- if add_mean_fractional_dispersion == 1:
- mean_fractional_dispersion = (
- np.std(np.array(star_surface_density[time]), axis=0) / ave_profiles
- )[
- (mean_dispersion_radius_range[0] < star_radius_DF_xs[time][0])
- & (star_radius_DF_xs[time][0] < mean_dispersion_radius_range[1])
- ].mean()
- mean_fractional_dispersion_in_dex = np.log10(
- 1.0 + mean_fractional_dispersion
- )
- plt.annotate(
- "Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) "
- % (
- mean_dispersion_radius_range[0],
- mean_dispersion_radius_range[1],
- mean_fractional_dispersion,
- mean_fractional_dispersion_in_dex,
- ),
- xy=(0.02, 0.84),
- size=10,
- xycoords="axes fraction",
- )
- plt.savefig(
- "star_surface_density_radial_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
-
- if draw_star_radius_DF == 2 and time != 0:
- # plt.subplot(111, aspect=1)
- fig = plt.figure(figsize=(8, 8))
- gridspec.GridSpec(4, 1)
- plt.subplot2grid((4, 1), (0, 0), rowspan=3)
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- temp = []
- dr = 0.5 * (
- sfr_radius_DF_xs[time][code][1] - sfr_radius_DF_xs[time][code][0]
- )
- for radius in range(len(sfr_radius_DF_profiles[time][code])):
- surface_area = np.pi * (
- (sfr_radius_DF_xs[time][code][radius] + dr) ** 2
- - (sfr_radius_DF_xs[time][code][radius] - dr) ** 2
- )
- temp.append(
- sfr_radius_DF_profiles[time][code][radius] / surface_area
- )
- sfr_surface_density[time].append(temp)
- lines = plt.plot(
- sfr_radius_DF_xs[time][code],
- sfr_surface_density[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.semilogy()
- plt.xlim(0, 14)
- plt.ylim(1e-4, 1e1)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- "$\mathrm{Star\ Formation\ Rate\ Surface\ Density\ (M_{\odot}/yr/kpc^2)}$",
- fontsize="large",
- )
- plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
-
- plt.subplot2grid((4, 1), (3, 0), rowspan=1)
- ave_profiles = np.mean(np.array(sfr_surface_density[time]), axis=0)
- for code in range(len(codes)):
- lines = plt.plot(
- sfr_radius_DF_xs[time][code],
- (sfr_surface_density[time][code] - ave_profiles) / ave_profiles,
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 14)
- plt.ylim(-1.0, 3.0)
- plt.grid(True)
- plt.locator_params(axis="y", nbins=4)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Residual,}\/(\mathrm{\Sigma} - \mathrm{\bar \Sigma})/\mathrm{\bar \Sigma}$",
- fontsize="large",
- )
- # if add_mean_fractional_dispersion == 1:
- # mean_fractional_dispersion = (np.std(np.array(sfr_surface_density[time]), axis=0)/ave_profiles)[(mean_dispersion_radius_range[0] < sfr_radius_DF_xs[time][0]) & (sfr_radius_DF_xs[time][0] < mean_dispersion_radius_range[1])].mean()
- # mean_fractional_dispersion_in_dex = np.log10(1.+mean_fractional_dispersion)
- # plt.annotate("Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) " % (mean_dispersion_radius_range[0], mean_dispersion_radius_range[1], mean_fractional_dispersion, mean_fractional_dispersion_in_dex), xy=(0.02, 0.84), size=10, xycoords='axes fraction')
- plt.savefig(
- "sfr_surface_density_radial_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
-
- # Draw K-S plot; below assumes that surface_density (or radius_DF_profiles) and sfr_surface_density (or sfr_radius_DF_profiles) have the identical size (n_bins in ProfilePlot)
- # plt.subplot(111)
- fig = plt.figure(figsize=(8, 7))
- KS_fit_t1 = []
- KS_fit_t2 = []
- Bigiel_surface_density = []
- Bigiel_sfr_surface_density = []
- for code in range(len(codes)):
- # Remove bins where SFR surface density is zero
- KS_x = np.array(surface_density[time][code])
- KS_x[np.where(np.array(sfr_surface_density[time][code]) < 1e-10)] = 0
- KS_x = np.log10(KS_x)
- KS_y = np.log10(np.array(sfr_surface_density[time][code]))
- KS_x = KS_x[~np.isinf(KS_x)]
- KS_y = KS_y[~np.isinf(KS_y)]
- t1, t2 = np.polyfit(KS_x, KS_y, 1)
- KS_fit_t1.append(t1)
- KS_fit_t2.append(t2)
- plt.scatter(
- KS_x,
- KS_y,
- color=color_names[code],
- edgecolor=color_names[code],
- s=30,
- linewidth=0.7,
- marker=marker_names[code],
- alpha=0.8,
- )
- plt.xlim(0, 3)
- plt.ylim(-4, 1)
- plt.grid(True)
- plt.xlabel("$\mathrm{log[Gas\ Surface\ Density\ (M_{\odot}/pc^2)]}$")
- plt.ylabel(
- "$\mathrm{log[Star\ Formation\ Rate\ Surface\ Density\ (M_{\odot}/yr/kpc^2)]}$"
- )
- plt.legend(codes, loc=2, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
- t = np.arange(-2, 5, 0.01)
- # for line in import_text("./Bigieletal2008_Fig8_Contour.txt", " "):
- # Bigiel_surface_density.append(float(line[0]) + np.log10(1.36)) # for the factor 1.36, see Section 2.3.1 of Bigiel et al. 2008
- # Bigiel_sfr_surface_density.append(float(line[1]))
- # plt.fill(Bigiel_surface_density, Bigiel_sfr_surface_density, fill=True, color='b', alpha = 0.1, hatch='\\') # contour by Bigiel et al. 2008 (Fig 8), or Feldmann et al. 2012 (Fig 1)
- # plt.plot(Bigiel_surface_density, Bigiel_sfr_surface_density, 'b-', linewidth = 0.7, alpha = 0.4)
- plt.plot(
- t, 1.37 * t - 3.78, "k--", linewidth=2, alpha=0.7
- ) # observational fit by Kennicutt et al. 2007
- # plt.savefig("K-S_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
- for code in range(len(codes)):
- plt.plot(
- t,
- np.polyval([KS_fit_t1[code], KS_fit_t2[code]], t),
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- alpha=0.7,
- ) # linear fits
- plt.savefig(
- "K-S_with_fits_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
- if draw_height_DF == 1:
- plt.clf()
- # plt.subplot(111, aspect=0.5)
- fig = plt.figure(figsize=(8, 6))
- for code in range(len(codes)):
- lines = plt.plot(
- height_DF_xs[time][code],
- np.add.accumulate(height_DF_profiles[time][code]),
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.semilogy()
- plt.xlim(0, 1.4)
- plt.ylim(1e9, 2e10)
- plt.grid(True)
- plt.xlabel("$\mathrm{Vertical\ Height\ (kpc)}$", fontsize="large")
- plt.ylabel("$\mathrm{Mass\ (M_{\odot})}$", fontsize="large")
- plt.legend(codes, loc=4, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
- # plt.savefig("height_DF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
-
- plt.clf()
- # plt.subplot(111, aspect=0.8)
- fig = plt.figure(figsize=(8, 8))
- gridspec.GridSpec(4, 1)
- plt.subplot2grid((4, 1), (0, 0), rowspan=3)
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- temp = []
- dh = height_DF_xs[time][code][1] - height_DF_xs[time][code][0]
- for height in range(len(height_DF_profiles[time][code])):
- surface_area = (
- 2 * dh * 1e3 * figure_width * 1e3
- ) # surface_area = 2*d(height)*figure_width in pc^2
- temp.append(height_DF_profiles[time][code][height] / surface_area)
- height_surface_density[time].append(temp)
- lines = plt.plot(
- height_DF_xs[time][code],
- height_surface_density[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.semilogy()
- plt.xlim(0, 1.4)
- plt.ylim(1e-1, 2e3)
- plt.grid(True)
- plt.xlabel("$\mathrm{Vertical\ Height\ (kpc)}$", fontsize="large")
- plt.ylabel("$\mathrm{Surface\ Density\ (M_{\odot}/pc^2)}$", fontsize="large")
- plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
-
- plt.subplot2grid((4, 1), (3, 0), rowspan=1)
- ave_profiles = np.mean(np.array(height_surface_density[time]), axis=0)
- for code in range(len(codes)):
- lines = plt.plot(
- height_DF_xs[time][code],
- (height_surface_density[time][code] - ave_profiles) / ave_profiles,
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, 1.4)
- plt.ylim(-1.0, 3.0)
- plt.grid(True)
- plt.locator_params(axis="y", nbins=4)
- plt.xlabel("$\mathrm{Vertical\ Height\ (kpc)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Residual,}\/(\mathrm{\Sigma} - \mathrm{\bar \Sigma})/\mathrm{\bar \Sigma}$",
- fontsize="large",
- )
- if add_mean_fractional_dispersion == 1:
- mean_fractional_dispersion = (
- np.std(np.array(height_surface_density[time]), axis=0) / ave_profiles
- )[
- (mean_dispersion_height_range[0] < height_DF_xs[time][0])
- & (height_DF_xs[time][0] < mean_dispersion_height_range[1])
- ].mean()
- mean_fractional_dispersion_in_dex = np.log10(
- 1.0 + mean_fractional_dispersion
- )
- plt.annotate(
- "Mean fractional dispersion for %.1f < z < %.1f kpc = %.3f (%.3f dex) "
- % (
- mean_dispersion_height_range[0],
- mean_dispersion_height_range[1],
- mean_fractional_dispersion,
- mean_fractional_dispersion_in_dex,
- ),
- xy=(0.02, 0.84),
- size=10,
- xycoords="axes fraction",
- )
- plt.savefig(
- "gas_surface_density_vertical_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
- if draw_SFR >= 1 and time != 0:
- plt.clf()
- # plt.subplot(111)#, aspect=1e-7)
- fig = plt.figure(figsize=(8, 6))
- for code in range(len(codes)):
- lines = plt.plot(
- sfr_ts[time][code],
- sfr_cum_masses[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, times[time])
- plt.ylim(0, 2.5e9)
- plt.grid(True)
- plt.xlabel("$\mathrm{Time\ (Myr)}$", fontsize="large")
- plt.ylabel("$\mathrm{Stellar\ Mass\ (M_{\odot})}$", fontsize="large")
- plt.legend(codes, loc=2, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
- plt.savefig(
- "Stellar_mass_evolution_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
- # plt.subplot(111, aspect=55)
- fig = plt.figure(figsize=(8, 8))
- gridspec.GridSpec(4, 1)
- plt.subplot2grid((4, 1), (0, 0), rowspan=3)
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- for code in range(len(codes)):
- lines = plt.plot(
- sfr_ts[time][code],
- sfr_sfrs[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, times[time])
- plt.ylim(0, 8)
- plt.grid(True)
- plt.xlabel("$\mathrm{Time\ (Myr)}$", fontsize="large")
- plt.ylabel("$\mathrm{Star\ Formation\ Rate\ (M_{\odot}/yr)}$", fontsize="large")
- if draw_SFR == 2:
- pf_sfr_ref = load_dataset(
- 2,
- 1,
- 0,
- ["GIZMO"],
- [
- [
- "dummy",
- "/lustre/ki/pfs/mornkr/080112_CHaRGe/pfs-hyades/AGORA-DISK-repository/Grackle+SF+ThermalFbck/GIZMO/AGORA_disk_second_ps1.5/snapshot_temp_100_last",
- ]
- ],
- )
- pf_sfr_ref.unit_registry.add(
- "pccm", pf_sfr_ref.unit_registry.lut["pc"][0], length, "\\rm{pc}/(1+z)"
- )
- pf_sfr_ref.hubble_constant = 0.71
- pf_sfr_ref.omega_lambda = 0.73
- pf_sfr_ref.omega_matter = 0.27
- pf_sfr_ref.omega_curvature = 0.0
- sp_sfr_ref = pf_sfr_ref.sphere(center, (0.5 * figure_width, "kpc"))
- draw_SFR_mass_ref = sp_sfr_ref[
- (PartType_Star_to_use, "particle_mass")
- ].in_units("Msun")
- draw_SFR_ct_ref = sp_sfr_ref[
- (PartType_Star_to_use, FormationTimeType_to_use)
- ].in_units("Myr")
- sfr_ref = StarFormationRate(
- pf_sfr_ref,
- star_mass=draw_SFR_mass_ref,
- star_creation_time=draw_SFR_ct_ref,
- volume=sp_sfr_ref.volume(),
- bins=25,
- )
- lines = plt.plot(
- sfr_ref.time.in_units("Myr"),
- sfr_ref.Msol_yr,
- color="k",
- linestyle="--",
- marker="*",
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.legend(codes + ["GIZMO-ps2"], loc=2, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
- plt.subplot2grid((4, 1), (3, 0), rowspan=1)
- ave_profiles = np.mean(np.array(sfr_sfrs[time]), axis=0)
- for code in range(len(codes)):
- lines = plt.plot(
- sfr_ts[time][code],
- (sfr_sfrs[time][code].d - ave_profiles) / ave_profiles,
- color=color_names[code],
- linestyle="none",
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.xlim(0, times[time])
- plt.ylim(-1.0, 1.0)
- plt.grid(True)
- plt.xlabel("$\mathrm{Time\ (Myr)}$", fontsize="large")
- plt.ylabel(
- r"$\mathrm{Residual,}\/(\mathrm{\dot{\rho}_\star} - \mathrm{\bar{\dot{\rho}_\star}})/\mathrm{\bar{\dot{\rho}_\star}}$",
- fontsize="large",
- )
- if add_mean_fractional_dispersion == 1:
- mean_fractional_dispersion = (
- np.std(np.array(sfr_sfrs[time]), axis=0) / ave_profiles
- )[
- (mean_dispersion_time_range[0] < sfr_ts[time][0])
- & (sfr_ts[time][0] < mean_dispersion_time_range[1])
- ].mean()
- mean_fractional_dispersion_in_dex = np.log10(
- 1.0 + mean_fractional_dispersion
- )
- plt.annotate(
- "Mean fractional dispersion for %d < t < %d Myr = %.3f (%.3f dex) "
- % (
- mean_dispersion_time_range[0],
- mean_dispersion_time_range[1],
- mean_fractional_dispersion,
- mean_fractional_dispersion_in_dex,
- ),
- xy=(0.35, 0.08),
- size=10,
- xycoords="axes fraction",
- )
- plt.savefig(
- "SFR_%dMyr" % times[time], bbox_inches="tight", pad_inches=0.03, dpi=300
- )
- plt.clf()
- if draw_cut_through == 1:
- plt.clf()
- # plt.subplot(111, aspect=0.5)
- fig = plt.figure(figsize=(8, 6))
- for code in range(len(codes)):
- lines = plt.plot(
- cut_through_zs[time][code],
- cut_through_zvalues[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.semilogy()
- plt.xlim(-1.4, 1.4)
- plt.ylim(1e-26, 1e-22)
- plt.grid(True)
- plt.xlabel("$\mathrm{Vertical\ Height\ (kpc)}$", fontsize="large")
- plt.ylabel("$\mathrm{Density\ (g/cm^3)}$", fontsize="large")
- plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
- z = np.arange(-1.4, 1.4, 0.05)
- plt.plot(
- z,
- rho_agora_disk(0, np.abs(z)),
- linestyle="--",
- linewidth=2,
- color="k",
- alpha=0.7,
- )
- plt.savefig(
- "cut_through_z_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()
- # plt.subplot(111, aspect=0.5)
- fig = plt.figure(figsize=(8, 6))
- for code in range(len(codes)):
- lines = plt.plot(
- cut_through_xs[time][code],
- cut_through_xvalues[time][code],
- color=color_names[code],
- linestyle=linestyle_names[np.mod(code, len(linestyle_names))],
- marker=marker_names[code],
- markeredgecolor="none",
- linewidth=1.2,
- alpha=0.8,
- )
- plt.semilogy()
- plt.xlim(-14, 14)
- plt.ylim(1e-26, 1e-22)
- plt.grid(True)
- plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize="large")
- plt.ylabel("$\mathrm{Density\ (g/cm^3)}$", fontsize="large")
- plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
- leg = plt.gca().get_legend()
- ltext = leg.get_texts()
- plt.setp(ltext, fontsize="small")
- x = np.arange(-14, 14, 0.05)
- plt.plot(
- x,
- rho_agora_disk(np.abs(x), 0),
- linestyle="--",
- linewidth=2,
- color="k",
- alpha=0.7,
- )
- plt.savefig(
- "cut_through_x_%dMyr" % times[time],
- bbox_inches="tight",
- pad_inches=0.03,
- dpi=300,
- )
- plt.clf()