diff --git a/examples/GEAR/AgoraDisk/cleanupSwift.py b/examples/GEAR/AgoraDisk/cleanupSwift.py
index e3c364081fc82986968727695e1142690781bb67..31e6ed8ab0ce6650bad861ac39462cab8e4a0f4b 100644
--- a/examples/GEAR/AgoraDisk/cleanupSwift.py
+++ b/examples/GEAR/AgoraDisk/cleanupSwift.py
@@ -13,9 +13,10 @@ copyfile(filename, out)
 
 f = File(out)
 
-name = "PartType0/ElementAbundance"
-if name in f:
-    del f[name]
+for i in range(6):
+    name = "PartType{}/ElementAbundance".format(i)
+    if name in f:
+        del f[name]
 
 for i in range(NPartType):
     name = "PartType%i" % i
@@ -25,4 +26,10 @@ for i in range(NPartType):
     grp = f[name + "/SmoothingLength"]
     grp[:] *= 1.823
 
+cosmo = f["Cosmology"].attrs
+head = f["Header"].attrs
+head["OmegaLambda"] = cosmo["Omega_lambda"]
+head["Omega0"] = cosmo["Omega_b"]
+head["HubbleParam"] = cosmo["H0 [internal units]"]
+
 f.close()
diff --git a/examples/GEAR/AgoraDisk/getSolution.sh b/examples/GEAR/AgoraDisk/getSolution.sh
index 41fbab800098bfaa381ca2e83cab0210c8dcf924..37b7f8031b90a5612d8fd3e88b1ecffbc02451a6 100755
--- a/examples/GEAR/AgoraDisk/getSolution.sh
+++ b/examples/GEAR/AgoraDisk/getSolution.sh
@@ -1,36 +1,12 @@
 #!/bin/bash
 
-OPTIND=1
-
-with_cooling=0
-
-function show_help {
-    echo "Valid options are:"
-    echo "\t -h \t Show this help"
-    echo "\t -C \t Download solution with cooling"
-}
-
-while getopts "h?C" opt; do
-    case "$opt" in
-	h|\?)
-	    show_help
-	    exit
-	    ;;
-	C)
-	    with_cooling=1
-	    ;;
-    esac
-done
-
 # cleanup work space
 rm snapshot_0000.hdf5 snapshot_0500.hdf5
 
-if [ $with_cooling -eq 1 ]; then
-    wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/with_cooling/snapshot_0000.hdf5
-    wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/with_cooling/snapshot_0500.hdf5
-else
-    wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/without_cooling/snapshot_0000.hdf5
-    wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/without_cooling/snapshot_0500.hdf5
-fi
+wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/with_cooling/snapshot_0000.hdf5
+wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/with_cooling/snapshot_0500.hdf5
+
+# wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/without_cooling/snapshot_0000.hdf5
+# wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/without_cooling/snapshot_0500.hdf5
 
 
diff --git a/examples/GEAR/AgoraDisk/plotSolution.py b/examples/GEAR/AgoraDisk/plotSolution.py
index fd5ed8a5e9fb1a1475202633595699d682edf281..8db2c98dbca87104bb119845edcbb934af08195e 100644
--- a/examples/GEAR/AgoraDisk/plotSolution.py
+++ b/examples/GEAR/AgoraDisk/plotSolution.py
@@ -1,12 +1,11 @@
-#!/usr/bin/env python2
+#!/usr/bin/env python3
 #######################################################################
 #
 #  UNIFIED ANALYSIS SCRIPT FOR DISK SIMULATION FOR THE AGORA PROJECT
 #
 #  FOR SCRIPT HISTORY SEE VERSION CONTROL CHANGELOG
 #
-#  Note: This script requires yt-3.2 or yt-dev. Older versions may
-#        yield incorrect results.
+#  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/)
@@ -37,35 +36,32 @@ 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
 
-if "-C" in sys.argv:
-    with_cooling = 1
-else:
-    with_cooling = 0
-    
 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                    = 0#0           # 0/1         = OFF/ON
-draw_star_clump_stats            = 0#0           # 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                     = 0###0         # 0/1/2       = OFF/ON/ON with local K-S plot using patches
-draw_PDF                         = 0###0         # 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                 = 0##0          # 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            = 0##0          # 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              = 0##0          # 0/1/2/3     = OFF/ON/ON with 1D profile/ON with analytic ftn subtracted
+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                   = 0##0          # 0/1         = OFF/ON
-draw_density_DF                  = 0#0           # 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                   = 0#0           # 0/1         = OFF/ON
-draw_star_radius_DF              = 0#0           # 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                   = 0#0           # 0/1         = OFF/ON
-draw_SFR                         = 0#0           # 0/1/2       = OFF/ON/ON with extra line with GIZMO-ps2
+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   = 0#0           # 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?"
@@ -113,10 +109,10 @@ marker_names             = ['s', 'o', 'p', 'v', '^', '<', '>', 'h', '*']
 #             [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"],
+filenames = [[["./agora_disk_IC.hdf5", "./agora_disk_500Myr.hdf5"], # Sim-noSFF
               ["./snapshot_0000.hdf5", "./snapshot_0500.hdf5"]],
-             [["./agora_disk_IC.hdf5", "./agora_disk_500Myr.hdf5"],
-              ["./snapshot_0000.hdf5", "./snapshot_0500.hdf5"]]]
+             [["./agora_disk_IC.hdf5", "./agora_disk_500Myr.hdf5"], # Sim-SFF (with star formation and feedback)
+              ["./snapshot_0000.hdf5", "./snapshot_0500.hdf5"]]] # I did not check the order, they can be switched
 
 # codes = ["SWIFT"]
 # filenames = [[["./agora_disk_0000.hdf5", "./agora_disk_0050.hdf5"]],
@@ -243,6 +239,69 @@ def load_dataset(dataset_num, time, code, codes, filenames_entry):
                           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             = []
@@ -513,14 +572,23 @@ for time in range(len(times)):
                         MassType_to_use = "Masses"
                 elif codes[code] == "SWIFT":
                         PartType_Gas_to_use = "PartType0"
-                        PartType_Star_to_use = "PartType2"
-                        PartType_StarBeforeFiltered_to_use = "PartType2"
+                        PartType_Star_to_use = "PartType4"
+                        PartType_StarBeforeFiltered_to_use = "PartType4"
+                        if time != 0:
+                            def _FormationTime(field, data):
+                                return pf.arr(data["PartType4", "BirthTime"].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
@@ -670,25 +738,16 @@ for time in range(len(times)):
                         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, "particle_position_cylindrical_z")])
+                                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)
-                        def _Density_2(field, data):
-                                return data[(PartType_Gas_to_use, "Density")].in_units('g/cm**3')
-                        pf.add_field((PartType_Gas_to_use, "Density_2"), 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.add_field((PartType_Gas_to_use, "Temperature_2"), 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.add_field((PartType_Gas_to_use, "Mass_2"), function=_Mass_2, take_log=True, particle_type=False, display_name="Mass", units="Msun")
                         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=False, 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, "particle_position_cylindrical_radius")], data[(PartType_Gas_to_use, "particle_position_cylindrical_z_abs")])
+                                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")
 
                 ####################################
@@ -701,7 +760,7 @@ for time in range(len(times)):
                 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], 'code_length') # naive usage such as YTArray([cen3[1], cen3[2], cen3[3]]) doesn't work somehow for ART-II data
+                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:
@@ -1109,13 +1168,18 @@ for time in range(len(times)):
                                 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_2"), (PartType_Gas_to_use, "Temperature_2"), (PartType_Gas_to_use, "Mass_2"), weight_field=None, fontsize=12, x_bins=300, y_bins=300)
-                                p3.set_zlim((PartType_Gas_to_use, "Mass_2"), 1e3, 1e8)
-                                p3.set_cmap((PartType_Gas_to_use, "Mass_2"), my_cmap2)
-                                plot3 = p3.plots[(PartType_Gas_to_use, "Mass_2")]
+                                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
@@ -1153,11 +1217,11 @@ for time in range(len(times)):
                                 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, "particle_position_cylindrical_radius")].take_log = False
+                                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, "particle_position_cylindrical_radius"), (PartType_Gas_to_use, "particle_velocity_cylindrical_theta"), \
+                                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("particle_position_cylindrical_radius", 'kpc')
+                                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)
@@ -1189,11 +1253,11 @@ for time in range(len(times)):
                                         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, "particle_position_cylindrical_radius"), (PartType_Gas_to_use, "particle_velocity_cylindrical_theta"), \
+                                        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("particle_position_cylindrical_radius", False)
-                                        p5.set_unit("particle_position_cylindrical_radius", 'kpc')
+                                        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)
@@ -1244,8 +1308,8 @@ for time in range(len(times)):
                                                 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, "particle_position_cylindrical_radius")].in_units("kpc") >= (radius - dr)) & \
-                                                                               (data[(PartType_Gas_to_use, "particle_position_cylindrical_radius")].in_units("kpc") < (radius + dr)))
+                                                        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")
@@ -1253,8 +1317,8 @@ for time in range(len(times)):
                                                 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, "particle_position_cylindrical_radius")].in_units("kpc") >= (radius - dr)) & \
-                                                                               (data[(PartType_Gas_to_use, "particle_position_cylindrical_radius")].in_units("kpc") < (radius + dr)))
+                                                        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")
@@ -1265,10 +1329,10 @@ for time in range(len(times)):
                                         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, "particle_position_cylindrical_radius"), (PartType_Gas_to_use, "particle_velocity_minus_local_rotational_velocity_squared"), \
+                                        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("particle_position_cylindrical_radius", False)
-                                        p55.set_unit("particle_position_cylindrical_radius", 'kpc')
+                                        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)
@@ -1290,10 +1354,10 @@ for time in range(len(times)):
                                                 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, "particle_position_cylindrical_radius"), (PartType_Gas_to_use, "particle_velocity_z_squared"), \
+                                        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("particle_position_cylindrical_radius", False)
-                                        p551.set_unit("particle_position_cylindrical_radius", 'kpc')
+                                        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)
@@ -1306,12 +1370,12 @@ for time in range(len(times)):
                 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, "particle_position_cylindrical_radius")].take_log = False
+                        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, "particle_position_cylindrical_radius"), (PartType_Star_to_use, "particle_velocity_cylindrical_theta"), \
+                        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("particle_position_cylindrical_radius", 'kpc')
+                        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
@@ -1333,11 +1397,11 @@ for time in range(len(times)):
 
                         # Add 1D profile line if requested
                         if draw_star_pos_vel_PDF >= 2 and time != 0:
-                                p51 = ProfilePlot(sp, (PartType_Star_to_use, "particle_position_cylindrical_radius"), (PartType_Star_to_use, "particle_velocity_cylindrical_theta"), \
+                                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, "particle_position_cylindrical_radius"), False)
-                                p51.set_unit("particle_position_cylindrical_radius", 'kpc')
+                                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)
@@ -1351,8 +1415,8 @@ for time in range(len(times)):
                                         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, "particle_position_cylindrical_radius")].in_units("kpc") >= (radius - dr)) & \
-                                                                       (data[(PartType_StarBeforeFiltered_to_use, "particle_position_cylindrical_radius")].in_units("kpc") < (radius + dr)))
+                                                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")
@@ -1360,8 +1424,8 @@ for time in range(len(times)):
                                         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, "particle_position_cylindrical_radius")].in_units("kpc") >= (radius - dr)) & \
-                                                                       (data[(PartType_StarBeforeFiltered_to_use, "particle_position_cylindrical_radius")].in_units("kpc") < (radius + dr)))
+                                                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")
@@ -1373,10 +1437,10 @@ for time in range(len(times)):
                                              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, "particle_position_cylindrical_radius"), (PartType_Star_to_use, "particle_velocity_minus_local_rotational_velocity_squared"), \
+                                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, "particle_position_cylindrical_radius"), False)
-                                p515.set_unit("particle_position_cylindrical_radius", 'kpc')
+                                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)
@@ -1389,10 +1453,10 @@ for time in range(len(times)):
                                              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, "particle_position_cylindrical_radius"), (PartType_Star_to_use, "particle_velocity_z_squared"), \
+                                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, "particle_position_cylindrical_radius"), False)
-                                p5151.set_unit("particle_position_cylindrical_radius", 'kpc')
+                                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)
@@ -1426,21 +1490,21 @@ for time in range(len(times)):
                         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, "particle_position_cylindrical_radius"), (PartType_Gas_to_use, "particle_position_cylindrical_z_abs"), (PartType_Gas_to_use, "Density_2"), weight_field=(PartType_Gas_to_use, "Mass_2"), fontsize=12, x_bins=200, y_bins=200)
-                                        p55.set_zlim((PartType_Gas_to_use, "Density_2"), 1e-26, 1e-21)
-                                        p55.set_cmap((PartType_Gas_to_use, "Density_2"), 'algae')
-                                        p55.set_log("particle_position_cylindrical_radius", False)
+                                        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("particle_position_cylindrical_radius", 'kpc')
+                                        p55.set_unit("radius", 'kpc')
                                         p55.set_unit("particle_position_cylindrical_z_abs", 'kpc')
-                                        plot55 = p55.plots[(PartType_Gas_to_use, "Density_2")]
+                                        plot55 = p55.plots[(PartType_Gas_to_use, "density")]
                                 elif draw_rad_height_PDF == 3:
-                                        p55 = PhasePlot(sp, (PartType_Gas_to_use, "particle_position_cylindrical_radius"), (PartType_Gas_to_use, "particle_position_cylindrical_z_abs"), (PartType_Gas_to_use, "Density_2_minus_analytic"), weight_field=(PartType_Gas_to_use, "Mass_2"), fontsize=12, x_bins=200, y_bins=200)
+                                        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("particle_position_cylindrical_radius", False)
+                                        p55.set_log("radius", False)
                                         p55.set_log("particle_position_cylindrical_z_abs", False)
-                                        p55.set_unit("particle_position_cylindrical_radius", 'kpc')
+                                        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")]
 
@@ -1469,11 +1533,11 @@ for time in range(len(times)):
                                         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, "particle_position_cylindrical_radius"), (PartType_Gas_to_use, "particle_position_cylindrical_z_abs"), \
+                                        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("particle_position_cylindrical_radius", False)
+                                        p56.set_log("radius", False)
                                         p56.set_log("particle_position_cylindrical_z_abs", False)
-                                        p56.set_unit("particle_position_cylindrical_radius", 'kpc')
+                                        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)
@@ -1501,7 +1565,7 @@ for time in range(len(times)):
                         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 = PhasePlot(sp, (PartType_Gas_to_use, "Density_2"), (PartType_Gas_to_use, "Temperature_2"), (PartType_Gas_to_use, "Metallicity_2"), weight_field=(PartType_Gas_to_use, "Mass_2"), fontsize=12, x_bins=300, y_bins=300)
+                                p3 = PhasePlot(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")]
@@ -1530,12 +1594,11 @@ for time in range(len(times)):
                                 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_2"),  (PartType_Gas_to_use, "Mass_2"), weight_field=None, n_bins=50, x_log=True, accumulation=False)
-#				p6 = ProfilePlot(sp, (PartType_Gas_to_use, "Density_2"),  (PartType_Gas_to_use, "Mass_2"), weight_field=None, n_bins=50, x_log=True, accumulation=True)
-                                p6.set_log("Mass_2", True)
-                                p6.set_xlim(1e-29, 1e-21)
+                                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]["Mass_2"].in_units('Msun').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:
@@ -1559,15 +1622,15 @@ for time in range(len(times)):
                                         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_2"), function=_Density_2, take_log=True, particle_type=False, display_name="Density", 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, "Mass_2"), function=_Mass_2, take_log=True, particle_type=False, display_name="Mass", units="Msun")
-                                                p61 = ProfilePlot(sp_1st, (PartType_Gas_to_use, "Density_2"),  (PartType_Gas_to_use, "Mass_2"), weight_field=None, n_bins=50, x_log=True, accumulation=False)
-                                                p61.set_log("Mass_2", True)
+                                                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]["Mass_2"].in_units('Msun').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
@@ -1585,13 +1648,14 @@ for time in range(len(times)):
                                 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, "particle_position_cylindrical_radius"),  (PartType_Gas_to_use, "Mass_2"), weight_field=None, n_bins=50, x_log=False, accumulation=False)
-                                p7.set_log("Mass_2", True)
-                                p7.set_log("particle_position_cylindrical_radius", False)
-                                p7.set_unit("particle_position_cylindrical_radius", 'kpc')
+
+                                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]["Mass_2"].in_units('Msun').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:
@@ -1599,9 +1663,9 @@ for time in range(len(times)):
                         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, "particle_position_cylindrical_radius")].take_log = False
-                        p71 = ProfilePlot(sp, (PartType_Star_to_use, "particle_position_cylindrical_radius"),  (PartType_Star_to_use, "particle_mass"), weight_field=None, n_bins=50, x_log=False, accumulation=False)
-                        p71.set_unit("particle_position_cylindrical_radius", 'kpc')
+                        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)
@@ -1619,10 +1683,10 @@ for time in range(len(times)):
 
                                 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, "particle_position_cylindrical_radius")].take_log = False
-                                p72 = ProfilePlot(sp, (PartType_Star_to_use, "particle_position_cylindrical_radius"), (PartType_Star_to_use, "particle_mass_young_stars_star_radius_DF"), \
+                                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("particle_position_cylindrical_radius", 'kpc')
+                                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
@@ -1640,13 +1704,13 @@ for time in range(len(times)):
                                 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, "Mass_2"), weight_field=None, n_bins=10, x_log=False, accumulation=False)
-                                p8.set_log("Mass_2", True)
+                                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]["Mass_2"].in_units('Msun').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:
@@ -1698,9 +1762,9 @@ for time in range(len(times)):
         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 >= 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']])
@@ -1755,7 +1819,7 @@ for time in range(len(times)):
                         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.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))
@@ -1790,7 +1854,7 @@ for time in range(len(times)):
                         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.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))
@@ -1817,7 +1881,7 @@ for time in range(len(times)):
                         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.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))
@@ -1870,9 +1934,12 @@ for time in range(len(times)):
                                 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]))
-                                KS_x = KS_x[~np.isinf(KS_x)]
-                                KS_y = KS_y[~np.isinf(KS_y)]
-                                if codes[code] == "CHANGA":
+                                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
@@ -1893,10 +1960,10 @@ for time in range(len(times)):
                         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)
+                        # 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
@@ -1909,13 +1976,13 @@ for time in range(len(times)):
                         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_2"), function=_Density_2, take_log=True, particle_type=False, display_name="Density", 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_2"), function=_Temperature_2, take_log=True, particle_type=False, display_name="Temperature", 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, "Mass_2"), function=_Mass_2, take_log=True, particle_type=False, display_name="Mass", 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"))
@@ -1930,16 +1997,16 @@ for time in range(len(times)):
                         # 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_2"), (PartType_Gas_to_use, "Temperature_2"), weight_field=(PartType_Gas_to_use, "Mass_2"), n_bins=30)
+                        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_2", 10, 1e7)
+                        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_2"].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)
+                #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)
@@ -2050,10 +2117,10 @@ for time in range(len(times)):
                                 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.+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.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)
+                # 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)
@@ -2163,10 +2230,10 @@ for time in range(len(times)):
                                 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.+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.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)
+                # 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)
@@ -2292,7 +2359,7 @@ for time in range(len(times)):
                 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.savefig("radius_DF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
 
                 plt.clf()
 #		plt.subplot(111, aspect=1)
@@ -2350,7 +2417,7 @@ for time in range(len(times)):
                 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.savefig("star_radius_DF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
 
                 plt.clf()
 #		plt.subplot(111, aspect=1)
@@ -2465,13 +2532,13 @@ for time in range(len(times)):
                         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)
+                        # 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)
+                        # 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)
@@ -2492,7 +2559,7 @@ for time in range(len(times)):
                 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.savefig("height_DF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
 
                 plt.clf()
 #		plt.subplot(111, aspect=0.8)
diff --git a/examples/GEAR/AgoraDisk/run.sh b/examples/GEAR/AgoraDisk/run.sh
index 5b85be7df875cee69b513e36a327b1469e35b60e..8284098a1a71093a4de269ddd0189c8f047c90f9 100755
--- a/examples/GEAR/AgoraDisk/run.sh
+++ b/examples/GEAR/AgoraDisk/run.sh
@@ -5,31 +5,21 @@
 # currently only the low resolution is available
 sim=low
 
-# enable cooling or not
-cooling=0
-
 # make run.sh fail if a subcommand fails
 set -e
 
-# define flags
-flag=
-if [ $cooling -eq 1 ]
-then
-    flag=-C
-fi
-
 # Generate the initial conditions if they are not present.
 if [ ! -e $sim.hdf5 ]
 then
     echo "Fetching initial glass file for the Sedov blast example..."
-    ./getIC.sh $sim.hdf5
+    ./getIC.sh $sim
 fi
 
 # Get the Grackle cooling table
 if [ ! -e CloudyData_UVB=HM2012.h5 ]
 then
     echo "Fetching the Cloudy tables required by Grackle..."
-    ../getCoolingTable.sh
+    ../../Cooling/getGrackleCoolingTable.sh 
 fi
 
 # copy the initial conditions
@@ -38,7 +28,7 @@ cp $sim.hdf5 agora_disk.hdf5
 python3 changeType.py agora_disk.hdf5
 
 # Run SWIFT
-#../../swift $flag --hydro --self-gravity --threads=4 agora_disk.yml 2>&1 | tee output.log
+../../swift --cooling --hydro --self-gravity --threads=4 agora_disk.yml 2>&1 | tee output.log
 
 
 echo "Changing smoothing length to be Gadget compatible"
@@ -46,7 +36,8 @@ python3 cleanupSwift.py agora_disk_0000.hdf5 agora_disk_IC.hdf5
 python3 cleanupSwift.py agora_disk_0050.hdf5 agora_disk_500Myr.hdf5
 
 echo "Fetching GEAR solution..."
-./getSolution.sh $flag
+./getSolution.sh
+
 
 echo "Plotting..."
-python3 plotSolution.py $flag
+python3 plotSolution.py