HighResQuant.py

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from gadget import *
from gadget_subfind import *
import const as c
from sys import argv, maxsize
import numpy as np
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from matplotlib.legend_handler import HandlerTuple
import cmocean as cmo
from cgm import ionisation, avgWeighted, colDensDist, avg, sumBins, coolingTime, starBirth
from statistics import median
from matplotlib.ticker import AutoMinorLocator
import csv
from matplotlib.markers import MarkerStyle
import seaborn as sns
import pandas as pd
from scipy.ndimage import gaussian_filter1d

numpy.set_printoptions(threshold=False)

snap = np.array([int(sys.argv[1])])
halo_d = ['level3_cgm_1e11', 'level4_cgm']

lcool_med_all, lcool_upper_all, lcool_lower_all, lcool_radii_all	=	[], [], [], []
TInflowMed, TInflowUpper, TInflowLower, TInflowRadii = [], [], [], []
TOutflowMed, TOutflowUpper, TOutflowLower, TOutflowRadii = [], [], [], []
VelInflow, VelOutflow = [], []
virialMass, coolRate, temp, tempcool, density, metallicity, metallicitycool, velocity, avgrefinement, uprefinement, lorefinement = [], [], [], [], [], [], [], [], [], [], []
vortin_med, vortin_radii = [], []
vortout_med, vortout_radii = [], []
haloes, lCool, mass, volume  = [], [], [], []
TempInflow, TempOutflow = [], []
MassInflow, MassOutflow = [], []
coolmass, SFR = [], []
cFracHIRad, cFracCIVRad, cFracSiIIRad, cFracMgIIRad, cFracOVIRad, colDensCovRad  = [], [], [], [], [], []
stellarBirthRedshift, stellarInitMass	= [], []
stelMass, virialRadius2, redshift = [], [], []
H_med, H_lower, H_upper, H_radii			= [], [], [], []
frac_cold, birth, birthx = [], [], []
t_med, t_upper, t_lower, t_upper2, t_lower2, t_radii																											= [], [], [], [], [], []
d_med, d_upper, d_lower, d_upper2, d_lower2, d_radii																											= [], [], [], [], [], []
m_med, m_upper, m_lower, m_upper2, m_lower2, m_radii																											= [], [], [], [], [], []
b_med, b_upper, b_lower, b_upper2, b_lower2, b_radii																											= [], [], [], [], [], []
v_med, v_upper, v_lower, v_upper2, v_lower2, v_radii																											= [], [], [], [], [], []

coolingmass = []
#A few cheeky constants
m_p: float			= 1.67e-24
k_b: float			= 1.38e-16
kpcincm: float		= c.parsec*1e3
G: float			= 6.67e-8
H_0: float			= 2.17e-18

for k in range(len(snap)):
	lcool_med, lcool_upper, lcool_lower, lcool_radii	=	[], [], [], []
	TInflowMed, TInflowUpper, TInflowLower, TInflowRadii = [], [], [], []
	for i in range(len(halo_d)):
		if halo_d[i] == 'level3_cgm_1e11':
			halo = ['h10_standard', 'h10_500pc', 'h10_250pc']
		if halo_d[i] == 'level4_cgm':
			halo = ['h5_standard', 'h5_1kpc', 'h5_500pc']
		for j in range(len(halo)):
			print(snap[k])
			if halo[j] == 'h10_250pc' or halo[j] == 'h10_500pc' or halo[j] == 'h5_500pc':
				sf = load_subfind(snap[k], dir = '/home/tango/spxfv/surge/' + halo_d[i] + '/%s/output/' % halo[j])
				load    = gadget_readsnap(snap[k], snappath = '/home/tango/spxfv/surge/' + halo_d[i] + '/%s/output/' % halo[j], loadonlytype=[0,4], lazy_load=True, subfind = sf)
			elif halo[j] == 'h10_standard':
				sf = load_subfind(snap[k], dir = '/home/cosmos/spxfv/Auriga/' + halo_d[i] + '/%s/output/' % halo[j])
				load    = gadget_readsnap(snap[k], snappath = '/home/cosmos/spxfv/Auriga/' + halo_d[i] + '/%s/output/' % halo[j], loadonlytype=[0,4], lazy_load=True, subfind = sf)
			else:
				sf = load_subfind(snap[k], dir = '/home/cosmos_spxfv/Auriga/' + halo_d[i] + '/%s/output/' % halo[j])
				load    = gadget_readsnap(snap[k], snappath = '/home/cosmos_spxfv/Auriga/' + halo_d[i] + '/%s/output/' % halo[j], loadonlytype=[0,4], lazy_load=True, subfind = sf)
			print(load.redshift)
			load.calc_sf_indizes(sf)
			load.select_halo(sf, do_rotation = True)
			haloes.append(halo[j])
			redshift.append(load.redshift)

			haloid = 0
			virialRadius = sf.data['frc2'][0]*1e3
			cent = sf.data['fpos'][haloid,:]
#			load.pos -= np.array(cent)
#			load.vel -= np.array(sf.data['fvel'][haloid,:])
			virialMass.append(sf.data['fmc2'][0]*1e10), stelMass.append(sf.data['fmty'][0,4]*1e10), SFR.append(sf.data['ssfr'][sf.data['ffsh'][0]])
			load.data['vort'] = np.sqrt((load.vort**2).sum(axis=1))

			#Investigating the birth of stars
			birthRedshift, birthMass = starBirth(load)
			stellarBirthRedshift.append(birthRedshift), stellarInitMass.append(birthMass)

			virialRadius2.append(virialRadius)

			load.pos *= 1e3
			load.vol *= 1e9 # Mpc^3 to kpc^3
			load.mass *= 1e10 # 1e10 Msun to Msun
			avgrefinement.append(median((load.vol*1e9)**(1/3)))#, uprefinement.append(np.percentile((load.vol)**(1/3), 84)), lorefinement.append(np.percentile((load.vol)**(1/3), 16))
			load.data['dist'] = np.sqrt((load.pos**2).sum(axis=1))

			normRadius = load.data['dist']/virialRadius

#			load.data['pos'] -= cent
#			load.data['pos'] *= 1e3
#			load.data['vel'] -= sf.data['fvel'][haloid,:]
#			load.data['vel'] *= 1e3 #convert to m/s

			#element number   0     1      2      3      4      5      6      7      8      9      10     11     12      13		14
			elements       = ['H',  'He',  'C',   'N',   'O',   'Ne',  'Mg',  'Si',  'Fe',  'Y',   'Sr',  'Zr',  'Ba',   'Pb',	'S']
			elements_Z     = [1,    2,     6,     7,     8,     10,    12,    14,    26,    39,    38,    40,    56,     82, 	16]
			elements_mass  = [1.01, 4.00,  12.01, 14.01, 16.00, 20.18, 24.30, 28.08, 55.85, 88.91, 87.62, 91.22, 137.33, 207.2, 32.06]
			elements_solar = [12.0, 10.93, 8.43,  7.83,  8.69,  7.93,  7.60,  7.51,  7.50,  2.21,  2.87,  2.58,  2.18,   1.75,	]

			Zsolar = 0.0127
			omegabaryon0 = 0.048

			print(load.data['vort'])

			rhocrit = 3. * (load.omega0 * (1+load.redshift)**3. + load.omegalambda) * (load.hubbleparam * 100*1e5/(c.parsec*1e6))**2. / ( 8. * pi * c.G)
			rhomean = 3. * (load.omega0 * (1+load.redshift)**3.) * (load.hubbleparam * 100*1e5/(c.parsec*1e6))**2. / ( 8. * pi * c.G)

			meanweight = sum(load.gmet[:,0:9][:len(load.ne)], axis = 1) / ( sum(load.gmet[:,0:9][:len(load.ne)]/elements_mass[0:9], axis = 1) + load.ne*load.gmet[:,0][:len(load.ne)] )
			Tfac = 1. / meanweight * (1.0 / (5./3.-1.)) * c.KB / c.amu * 1e10 * c.msol / 1.989e53 #converts to internal to temperature

			#Calculations of physical properties

			Tfac				= 1. / meanweight * (1.0 / (5./3.-1.)) * c.KB / c.amu * 1e10 * c.msol / 1.989e53

			load.bfld			= np.multiply(load.bfld, np.multiply(c.bfac,1e6)) #Converts to microG
			load.data['bfld']	= abs(np.sqrt((load.bfld**2).sum(axis=1))) #Magnitude of bfield vectors, microGaus
			load.data['dense']	= load.rho/(np.multiply(c.parsec, 1e6)) ** 3 * np.multiply(c.msol, 1e10)									#Density of cells in g/cm^3
			load.data['T']		= np.divide(load.u, Tfac)														#Temperature in K
			load.data['gz']		= abs(np.divide(load.gz, 0.0127))													#Metallicity normalised to solar metallicity
			load.data['thermP']	= np.divide(np.multiply(load.data['dense'], load.data['T']), np.multiply(meanweight, m_p))									#Thermal pressure in Pascals/k_b
			load.data['bfldP']	= ((load.data['bfld']/1e6)**2)/(8*np.pi*k_b)										#Magnetic pressure also in Pascals/k_b
			load.data['totP']	= np.add(load.data['thermP'], load.data['bfldP'])										#Total pressure
			load.data['n_H']	= np.divide(np.multiply(load.data['dense'], load.gmet[:,0][:len(load.ne)]), np.multiply(elements_mass[0], c.amu))				#Number density of hydrogen
			load.data['pRat']	= np.divide(load.data['thermP'], load.data['bfldP'])							#Pressure ratio
			load.data['gcol']	= load.gcol
			load.data['pos_m']	= np.multiply(load.pos, 3.086e16)
			load.data['sr']		= (load.vol*1e9)**(1/3)

			mag = np.sqrt(((load.pos[:,0])**2) + ((load.pos[:,1])**2) + ((load.pos[:,2])**2)) #Magnitude of the position of each cell from the centre of the halo.
			load.data['velocity'] = (load.data['vel'][:,0]*load.pos[:,0]+load.pos[:,1]*load.pos[:,1]+load.vel[:,2]*load.pos[:,2])/mag
			load.data['velocity'][:len(load.ne)] *= load.data['dense']

			#"where" statements
			wvortout, = np.where((load.sfr[:len(load.ne)] == 0) & (normRadius[:len(load.ne)] >= 0) & (normRadius[:len(load.ne)] <= 1) & (load.data['vort'][:len(load.ne)]>0) & (load.halo[:len(load.ne)] == 0) & (load.subhalo[:len(load.ne)] <= 0))
			w, = np.where((load.sfr[:len(load.ne)] == 0) & (normRadius[:len(load.ne)] >= 0.25) & (normRadius[:len(load.ne)] <= 1) & (load.halo[:len(load.ne)] == 0) & (load.subhalo[:len(load.ne)] <= 0))
			wAll, = np.where((load.sfr[:len(load.ne)] == 0) & (normRadius[:len(load.ne)] >= 0.0) & (normRadius[:len(load.ne)] <= 1) & (load.halo[:len(load.ne)] == 0) & (load.subhalo[:len(load.ne)] <= 0))
			fraccoldgas, = np.where(((load.data['T']>=10**4)) & (load.data['T']<=10**4.5) & (load.sfr[:len(load.ne)] == 0) & (normRadius[:len(load.ne)] >= 0.25) & (normRadius[:len(load.ne)] <= 1) & (load.halo[:len(load.ne)] == 0) & (load.subhalo[:len(load.ne)] <= 0))

			heatinggas, = np.where((load.gcol[:len(load.ne)]>0) & (load.sfr[:len(load.ne)] == 0) & (normRadius[:len(load.ne)] >= 0.25) & (normRadius[:len(load.ne)] <= 1) & (load.halo[:len(load.ne)] == 0) & (load.subhalo[:len(load.ne)] <= 0))
			coolinggas, = np.where((load.gcol[:len(load.ne)]<0) & (load.sfr[:len(load.ne)] == 0) & (normRadius[:len(load.ne)] >= 0.25) & (normRadius[:len(load.ne)] <= 1) & (load.halo[:len(load.ne)] == 0) & (load.subhalo[:len(load.ne)] <= 0))
			coolinggasall, = np.where((load.gcol[:len(load.ne)]<0) & (load.sfr[:len(load.ne)] == 0) & (normRadius[:len(load.ne)] >= 0) & (normRadius[:len(load.ne)] <= 1) & (load.halo[:len(load.ne)] == 0) & (load.subhalo[:len(load.ne)] <= 0))

			intCGM, = np.where((load.sfr[:len(load.ne)] == 0) & (normRadius[:len(load.ne)] >= 0.25) & (normRadius[:len(load.ne)] <= 1) & (load.halo[:len(load.ne)] == 0) & (load.subhalo[:len(load.ne)] <= 0))
			Inflow = np.where((load.sfr[:len(load.ne)] == 0) & (normRadius[:len(load.ne)] >= 0.25) & (normRadius[:len(load.ne)] <= 1) & (load.halo[:len(load.ne)] == 0) & (load.subhalo[:len(load.ne)] <= 0) & (load.data['velocity'][:len(load.ne)] < 0))
			Outflow = np.where((load.sfr[:len(load.ne)] == 0) & (normRadius[:len(load.ne)] >= 0.25) & (normRadius[:len(load.ne)] <= 1) & (load.halo[:len(load.ne)] == 0) & (load.subhalo[:len(load.ne)] <= 0) & (load.data['velocity'][:len(load.ne)] > 0))

			#Calculations for plotting
#			te_med, te_upper, te_lower, te_radii	= avgWeighted(normRadius[w], load.data['T'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, load.mass[w])
#			me_med, me_upper, me_lower, me_radii	= avgWeighted(normRadius[w], load.data['gz'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, load.mass[w])
#			de_med, de_upper, de_lower, de_radii	= avgWeighted(normRadius[w], load.data['n_H'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, load.vol[w])
			be_med, be_upper, be_lower, be_radii	= avgWeighted(normRadius[w], load.data['bfld'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, load.mass[w])
			vorteout_med, vorteout_upper, vorteout_lower, vorteout_radii	= avgWeighted(normRadius[wvortout], load.data['vort'][wvortout], min(normRadius[wvortout]), max(normRadius[wvortout]), min(normRadius[wvortout]), 100, load.mass[wvortout])
			birth_sum, birth_x						= sumBins(birthRedshift, birthMass, min(birthRedshift), max(birthRedshift), min(birthRedshift), 200)
			ve_med, ve_upper, ve_lower, ve_radii	= avg(normRadius[wAll], load.data['sr'][wAll], min(normRadius[wAll]), 2, 50) #, min(normRadius[w]))
			birth.append(birth_sum), birthx.append(birth_x)

#			t_med.append(te_med), t_upper.append(te_upper), t_lower.append(te_lower), t_radii.append(te_radii)
#			d_med.append(de_med), d_upper.append(de_upper), d_lower.append(de_lower), d_radii.append(de_radii)
#			m_med.append(me_med), m_upper.append(me_upper), m_lower.append(me_lower), m_radii.append(me_radii)
			b_med.append(be_med), b_upper.append(be_upper), b_lower.append(be_lower), b_radii.append(be_radii)
			v_med.append(ve_med), v_upper.append(ve_upper), v_lower.append(ve_lower), v_radii.append(ve_radii)
			TempInflow.append(load.data['T'][Inflow]), TempOutflow.append(load.data['T'][Outflow])
			MassInflow.append(load.mass[Inflow]), MassOutflow.append(load.mass[Outflow])
			VelInflow.append(load.data['velocity'][Inflow]), VelOutflow.append(load.data['velocity'][Outflow])
			vortout_med.append(vorteout_med), vortout_radii.append(vorteout_radii)

			load.data['Tcool'] = load.data['T'][fraccoldgas]
			coolmass.append(load.mass[fraccoldgas])
			frac_cold.append(len(load.mass[fraccoldgas])/len(load.mass[w]))

			coolingmass.append(load.mass[coolinggas])
			temp.append(load.data['T'][w])
			tempcool.append(load.data['T'][coolinggas])
			density.append(load.data['n_H'][w])
			metallicity.append(load.data['gz'][w])
			metallicitycool.append(load.data['gz'][coolinggas])
			velocity.append(load.data['velocity'][w])
			mass.append(load.mass[w])
			volume.append(load.vol[w])

			lcool, tcool = coolingTime(load, normRadius, 0.25, 1)
			lcool2, tcool2 = coolingTime(load, normRadius, 0.0, 1)
			lCool.append(lcool)

			TinflowMed, TinflowUpper, TinflowLower, TinflowRadii = avgWeighted(normRadius[Inflow], load.data['T'][Inflow], min(normRadius[Inflow]), max(normRadius[Inflow]), min(normRadius[Inflow]), 100, load.mass[Inflow])
			ToutflowMed, ToutflowUpper, ToutflowLower, ToutflowRadii = avgWeighted(normRadius[Outflow], load.data['T'][Outflow], min(normRadius[Outflow]), max(normRadius[Outflow]), min(normRadius[Outflow]), 100, load.mass[Outflow])
			lcool_med, lcool_upper, lcool_lower, lcool_radii = avgWeighted(normRadius[coolinggasall], lcool2, min(normRadius[coolinggasall]), max(normRadius[coolinggasall]), min(normRadius[coolinggasall]), 100, load.mass[coolinggasall])
			lcool_med_all.append(lcool_med), lcool_upper_all.append(lcool_upper), lcool_lower_all.append(lcool_lower), lcool_radii_all.append(lcool_radii)
			TInflowMed.append(TinflowMed), TInflowUpper.append(TinflowUpper), TInflowLower.append(TinflowLower), TInflowRadii.append(TinflowRadii)
			TOutflowMed.append(ToutflowMed), TOutflowUpper.append(ToutflowUpper), TOutflowLower.append(ToutflowLower), TOutflowRadii.append(ToutflowRadii)
			print(len(mass[i]), len(temp[i]))
			print(len(TempInflow[i]), len(MassInflow[i]))

linestyles = ['-', '--', '-.']
colours = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', 'tab:brown']

countcool0, bin_edgescool0, _ = plt.hist(np.log10(lCool[0]), range = [0,8],  weights = np.log10(coolingmass[0]), density = True, bins = 50, color = colours[0], linestyle = linestyles[0], label = 'SMR')
countcool1, bin_edgescool1, _ = plt.hist(np.log10(lCool[1]), range = [0,8],  weights = np.log10(coolingmass[1]), density = True, bins = 50, color = colours[1], linestyle = linestyles[1], label = '~500pc')
countcool2, bin_edgescool2, _ = plt.hist(np.log10(lCool[2]), range = [0,8],  weights = np.log10(coolingmass[2]), density = True, bins = 50, color = colours[2], linestyle = linestyles[2], label = '~250pc')
countcool3, bin_edgescool3, _ = plt.hist(np.log10(lCool[3]), range = [0,8],  weights = np.log10(coolingmass[3]), density = True, bins = 50, color = colours[3], linestyle = linestyles[0])
countcool4, bin_edgescool4, _ = plt.hist(np.log10(lCool[4]), range = [0,8],  weights = np.log10(coolingmass[4]), density = True, bins = 50, color = colours[4], linestyle = linestyles[1])
countcool5, bin_edgescool5, _ = plt.hist(np.log10(lCool[5]), range = [0,8],  weights = np.log10(coolingmass[5]), density = True, bins = 50, color = colours[5], linestyle = linestyles[2])

counttemp0, bin_edgestemp0, _ = plt.hist(np.log10(tempcool[0]), range = [4,6.5], weights = np.log10(coolingmass[0]), density = True, bins = 50, color = colours[0], linestyle = linestyles[0], label = 'SMR')
counttemp1, bin_edgestemp1, _ = plt.hist(np.log10(tempcool[1]), range = [4,6.5], weights = np.log10(coolingmass[1]), density = True, bins = 50, color = colours[1], linestyle = linestyles[1], label = '~500pc')
counttemp2, bin_edgestemp2, _ = plt.hist(np.log10(tempcool[2]), range = [4,6.5], weights = np.log10(coolingmass[2]), density = True, bins = 50, color = colours[2], linestyle = linestyles[2], label = '~250pc')
counttemp3, bin_edgestemp3, _ = plt.hist(np.log10(tempcool[3]), range = [4,6.5], weights = np.log10(coolingmass[3]), density = True, bins = 50, color = colours[3], linestyle = linestyles[0])
counttemp4, bin_edgestemp4, _ = plt.hist(np.log10(tempcool[4]), range = [4,6.5], weights = np.log10(coolingmass[4]), density = True, bins = 50, color = colours[4], linestyle = linestyles[1])
counttemp5, bin_edgestemp5, _ = plt.hist(np.log10(tempcool[5]), range = [4,6.5], weights = np.log10(coolingmass[5]), density = True, bins = 50, color = colours[5], linestyle = linestyles[2])

countmet0, bin_edgesmet0, _ = plt.hist(np.log10(metallicitycool[0]), range = [-5,1.0], weights = np.log10(coolingmass[0]), density = True, bins = 50, color = colours[0], linestyle = linestyles[0], label = 'SMR')
countmet1, bin_edgesmet1, _ = plt.hist(np.log10(metallicitycool[1]), range = [-5,1.0], weights = np.log10(coolingmass[1]), density = True, bins = 50, color = colours[1], linestyle = linestyles[1], label = '~500pc')
countmet2, bin_edgesmet2, _ = plt.hist(np.log10(metallicitycool[2]), range = [-5,1.0], weights = np.log10(coolingmass[2]), density = True, bins = 50, color = colours[2], linestyle = linestyles[2], label = '~250pc')
countmet3, bin_edgesmet3, _ = plt.hist(np.log10(metallicitycool[3]), range = [-5,1.0], weights = np.log10(coolingmass[3]), density = True, bins = 50, color = colours[3], linestyle = linestyles[0])
countmet4, bin_edgesmet4, _ = plt.hist(np.log10(metallicitycool[4]), range = [-5,1.0], weights = np.log10(coolingmass[4]), density = True, bins = 50, color = colours[4], linestyle = linestyles[1])
countmet5, bin_edgesmet5, _ = plt.hist(np.log10(metallicitycool[5]), range = [-5,1.0], weights = np.log10(coolingmass[5]), density = True, bins = 50, color = colours[5], linestyle = linestyles[2])

bin_centrescool0 = (bin_edgescool0[:-1] + bin_edgescool0[1:]) / 2
bin_centrescool1 = (bin_edgescool1[:-1] + bin_edgescool1[1:]) / 2
bin_centrescool2 = (bin_edgescool2[:-1] + bin_edgescool2[1:]) / 2
bin_centrescool3 = (bin_edgescool3[:-1] + bin_edgescool3[1:]) / 2
bin_centrescool4 = (bin_edgescool4[:-1] + bin_edgescool4[1:]) / 2
bin_centrescool5 = (bin_edgescool5[:-1] + bin_edgescool5[1:]) / 2

bin_centrestemp0 = (bin_edgestemp0[:-1] + bin_edgestemp0[1:]) / 2
bin_centrestemp1 = (bin_edgestemp1[:-1] + bin_edgestemp1[1:]) / 2
bin_centrestemp2 = (bin_edgestemp2[:-1] + bin_edgestemp2[1:]) / 2
bin_centrestemp3 = (bin_edgestemp3[:-1] + bin_edgestemp3[1:]) / 2
bin_centrestemp4 = (bin_edgestemp4[:-1] + bin_edgestemp4[1:]) / 2
bin_centrestemp5 = (bin_edgestemp5[:-1] + bin_edgestemp5[1:]) / 2

bin_centresmet0 = (bin_edgesmet0[:-1] + bin_edgesmet0[1:]) / 2
bin_centresmet1 = (bin_edgesmet1[:-1] + bin_edgesmet1[1:]) / 2
bin_centresmet2 = (bin_edgesmet2[:-1] + bin_edgesmet2[1:]) / 2
bin_centresmet3 = (bin_edgesmet3[:-1] + bin_edgesmet3[1:]) / 2
bin_centresmet4 = (bin_edgesmet4[:-1] + bin_edgesmet4[1:]) / 2
bin_centresmet5 = (bin_edgesmet5[:-1] + bin_edgesmet5[1:]) / 2

smooth_countcool0 = gaussian_filter1d(countcool0, sigma=1.0)
smooth_countcool1 = gaussian_filter1d(countcool1, sigma=1.0)
smooth_countcool2 = gaussian_filter1d(countcool2, sigma=1.0)
smooth_countcool3 = gaussian_filter1d(countcool3, sigma=1.0)
smooth_countcool4 = gaussian_filter1d(countcool4, sigma=1.0)
smooth_countcool5 = gaussian_filter1d(countcool5, sigma=1.0)

smooth_counttemp0 = gaussian_filter1d(counttemp0, sigma=1.0)
smooth_counttemp1 = gaussian_filter1d(counttemp1, sigma=1.0)
smooth_counttemp2 = gaussian_filter1d(counttemp2, sigma=1.0)
smooth_counttemp3 = gaussian_filter1d(counttemp3, sigma=1.0)
smooth_counttemp4 = gaussian_filter1d(counttemp4, sigma=1.0)
smooth_counttemp5 = gaussian_filter1d(counttemp5, sigma=1.0)

smooth_countmet0 = gaussian_filter1d(countmet0, sigma=1.0)
smooth_countmet1 = gaussian_filter1d(countmet1, sigma=1.0)
smooth_countmet2 = gaussian_filter1d(countmet2, sigma=1.0)
smooth_countmet3 = gaussian_filter1d(countmet3, sigma=1.0)
smooth_countmet4 = gaussian_filter1d(countmet4, sigma=1.0)
smooth_countmet5 = gaussian_filter1d(countmet5, sigma=1.0)

fig, axs = plt.subplots(nrows=2, ncols=3, figsize = (12,8), sharex = 'col', sharey = 'row')

dwarfSMR = axs[0,0].plot(bin_centrescool0, smooth_countcool0, color = colours[0], linestyle = linestyles[0], label = 'SMR')
dwarf500 = axs[0,0].plot(bin_centrescool1, smooth_countcool1, color = colours[1], linestyle = linestyles[1], label = '~500pc')
dwarf250 = axs[0,0].plot(bin_centrescool2, smooth_countcool2, color = colours[2], linestyle = linestyles[2], label = '~250pc')
MWSMR = axs[1,0].plot(bin_centrescool3, smooth_countcool3, color = colours[3], linestyle = linestyles[0], label = 'SMR')
MW1k = axs[1,0].plot(bin_centrescool4, smooth_countcool4, color = colours[4], linestyle = linestyles[1], label = '~1kpc')
MW500 = axs[1,0].plot(bin_centrescool5, smooth_countcool5, color = colours[5], linestyle = linestyles[2], label = '~500pc')

axs[0,1].plot(bin_centrestemp0, smooth_counttemp0, color = colours[0], linestyle = linestyles[0])
axs[0,1].plot(bin_centrestemp1, smooth_counttemp1, color = colours[1], linestyle = linestyles[1])
axs[0,1].plot(bin_centrestemp2, smooth_counttemp2, color = colours[2], linestyle = linestyles[2])
axs[1,1].plot(bin_centrestemp3, smooth_counttemp3, color = colours[3], linestyle = linestyles[0])
axs[1,1].plot(bin_centrestemp4, smooth_counttemp4, color = colours[4], linestyle = linestyles[1])
axs[1,1].plot(bin_centrestemp5, smooth_counttemp5, color = colours[5], linestyle = linestyles[2])

axs[0,2].plot(bin_centresmet0, smooth_countmet0, color = colours[0], linestyle = linestyles[0])
axs[0,2].plot(bin_centresmet1, smooth_countmet1, color = colours[1], linestyle = linestyles[1])
axs[0,2].plot(bin_centresmet2, smooth_countmet2, color = colours[2], linestyle = linestyles[2])
axs[1,2].plot(bin_centresmet3, smooth_countmet3, color = colours[3], linestyle = linestyles[0])
axs[1,2].plot(bin_centresmet4, smooth_countmet4, color = colours[4], linestyle = linestyles[1])
axs[1,2].plot(bin_centresmet5, smooth_countmet5, color = colours[5], linestyle = linestyles[2])

axs[0,0].axvline(np.log10(avgrefinement[1]), c = colours[1], linestyle = linestyles[1])
axs[0,0].axvline(np.log10(avgrefinement[2]), c = colours[2], linestyle = linestyles[2])
axs[1,0].axvline(np.log10(avgrefinement[4]), c = colours[4], linestyle = linestyles[1])
axs[1,0].axvline(np.log10(avgrefinement[5]), c = colours[5], linestyle = linestyles[2])

axs[-1,0].set_xticks([0,2,4,6,8])
axs[-1,1].set_xticks([4,4.5,5,5.5,6,6.5])
axs[-1,2].set_xticks([-5,-4,-3,-2,-1,0,1])

axs[1,0].set_xlabel(r'Log$_{10}l_{\mathrm{cool}}$ [pc]')
axs[1,1].set_xlabel(r'Log$_{10}$T [K]')
axs[1,2].set_xlabel(r'Log$_{10}$Z/Z$_{\odot}$')
axs[0,0].set_ylabel(r'PDF ($10^{11}$M$_{\odot}$)')
axs[1,0].set_ylabel(r'PDF ($10^{12}$M$_{\odot}$)')

axs[0,0].legend()
axs[1,0].legend()

axs[-1,0].set_xlim(0,8)
axs[-1,1].set_xlim(4,6.5)
axs[-1,2].set_xlim(-5,1)

for ax in axs.flat:
	ax.tick_params(axis = 'both', which = 'both', bottom = True, left = True, right = True, top = True, direction = 'in')
	ax.xaxis.set_minor_locator(AutoMinorLocator(4))
	ax.yaxis.set_minor_locator(AutoMinorLocator(4))
	ax.grid()
	ax.label_outer()

plt.subplots_adjust(hspace = 0, wspace = 0.2)
plt.savefig(f'./Figures/Paper3Plots/cooling_length_soimething_Radii_{snap[k]}.pdf', dpi = 300, transparent = True, bbox_inches = 'tight')
plt.close("fig")

plt.figure(2)
fig, axs = plt.subplots(nrows=1, ncols=2, figsize = (8,4), sharey = 'row')

axs[0].plot(v_radii[0], np.log10(v_med[0]), c='tab:blue', label = 'SMR')
axs[0].fill_between(v_radii[0], np.log10(v_lower[0]), np.log10(v_upper[0]), color = 'tab:blue', alpha = 0.3)
axs[0].plot(v_radii[1], np.log10(v_med[1]), c='tab:orange', label = '~500pc', linestyle = '--')
axs[0].fill_between(v_radii[1], np.log10(v_lower[1]), np.log10(v_upper[1]), color = 'tab:red', alpha = 0.3)
axs[0].plot(v_radii[2], np.log10(v_med[2]), c='tab:green', label = '~250pc', linestyle = '-.')
axs[0].fill_between(v_radii[2], np.log10(v_lower[2]), np.log10(v_upper[2]), color = 'tab:green', alpha = 0.3)
axs[1].plot(v_radii[3], np.log10(v_med[3]), c='tab:red', label = 'SMR')
axs[1].fill_between(v_radii[3], np.log10(v_lower[3]), np.log10(v_upper[3]), color = 'tab:red', alpha = 0.3)
axs[1].plot(v_radii[4], np.log10(v_med[4]), c='tab:purple', label = '~1kpc', linestyle = '--')
axs[1].fill_between(v_radii[4], np.log10(v_lower[4]), np.log10(v_upper[4]), color = 'tab:purple', alpha = 0.3)
axs[1].plot(v_radii[5], np.log10(v_med[5]), c='tab:brown', label = '~500pc', linestyle = '-.')
axs[1].fill_between(v_radii[5], np.log10(v_lower[5]), np.log10(v_upper[5]), color = 'tab:brown', alpha = 0.3)

axs[0].set_ylabel(r'Log$_{10}$(Cell Size) [pc]')

axs[0].set_title(r'$10^{11}$M$_{\odot}$')
axs[1].set_title(r'$10^{12}$M$_{\odot}$')

axs[0].legend()
axs[1].legend()

axs[0].grid()
axs[1].grid()

for ax in axs.flat:
	ax.label_outer()
	ax.tick_params(axis = 'both', which = 'both', bottom = True, left = True, right = True, top = True, direction = 'in')
	ax.xaxis.set_minor_locator(AutoMinorLocator(4))
	ax.yaxis.set_minor_locator(AutoMinorLocator(4))
	ax.set_xlabel(r'R/R$_{200c}$')
	ax.set_xlim(0,1)

axs[0].set_xticks([0.0,0.2,0.4,0.6,0.8])
axs[1].set_xticks([0.2,0.4,0.6,0.8,1.0])

plt.subplots_adjust(hspace = 0, wspace = 0)
plt.savefig(f'./Figures/Paper3Plots/spatial_refinement_radial_dependence_{snap[k]}.pdf', dpi = 300, transparent = True, bbox_inches = 'tight')
plt.close("all")

refinement = ['SMR', r'Medium\\Refinement', r'High\\Refinement']

countin0, bin_edgesin0, _ = plt.hist(np.log10(TempInflow[0]), range = [3.9,6.1], weights = np.log10(MassInflow[0]), density = True, bins = 50, color = colours[0], linestyle = linestyles[0], label = 'SMR')
countin1, bin_edgesin1, _ = plt.hist(np.log10(TempInflow[1]), range = [3.9,6.1], weights = np.log10(MassInflow[1]), density = True, bins = 50, color = colours[1], linestyle = linestyles[1], label = '~500pc')
countin2, bin_edgesin2, _ = plt.hist(np.log10(TempInflow[2]), range = [3.9,6.1], weights = np.log10(MassInflow[2]), density = True, bins = 50, color = colours[2], linestyle = linestyles[2], label = '~250pc')
countout0, bin_edgesout0, _ = plt.hist(np.log10(TempOutflow[0]), range = [3.9,6.1], weights = np.log10(MassOutflow[0]), density = True, bins = 50, color = colours[0], linestyle = linestyles[0], label = 'SMR')
countout1, bin_edgesout1, _ = plt.hist(np.log10(TempOutflow[1]), range = [3.9,6.1], weights = np.log10(MassOutflow[1]), density = True, bins = 50, color = colours[1], linestyle = linestyles[1], label = '~1kpc')
countout2, bin_edgesout2, _ = plt.hist(np.log10(TempOutflow[2]), range = [3.9,6.1], weights = np.log10(MassOutflow[2]), density = True, bins = 50, color = colours[2], linestyle = linestyles[2], label = '~500pc')

countin3, bin_edgesin3, _ = plt.hist(np.log10(TempInflow[3]), range = [3.9,6.1], weights = np.log10(MassInflow[3]), density = True, bins = 50, color = colours[3], linestyle = linestyles[0])
countin4, bin_edgesin4, _ = plt.hist(np.log10(TempInflow[4]), range = [3.9,6.1], weights = np.log10(MassInflow[4]), density = True, bins = 50, color = colours[4], linestyle = linestyles[1])
countin5, bin_edgesin5, _ = plt.hist(np.log10(TempInflow[5]), range = [3.9,6.1], weights = np.log10(MassInflow[5]), density = True, bins = 50, color = colours[5], linestyle = linestyles[2])
countout3, bin_edgesout3, _ = plt.hist(np.log10(TempOutflow[3]), range = [3.9,6.1], weights = np.log10(MassOutflow[3]), density = True, bins = 50, color = colours[3], linestyle = linestyles[0])
countout4, bin_edgesout4, _ = plt.hist(np.log10(TempOutflow[4]), range = [3.9,6.1], weights = np.log10(MassOutflow[4]), density = True, bins = 50, color = colours[4], linestyle = linestyles[1])
countout5, bin_edgesout5, _ = plt.hist(np.log10(TempOutflow[5]), range = [3.9,6.1], weights = np.log10(MassOutflow[5]), density = True, bins = 50, color = colours[5], linestyle = linestyles[2])

bin_centresin0 = (bin_edgesin0[:-1] + bin_edgesin0[1:]) / 2
bin_centresin1 = (bin_edgesin1[:-1] + bin_edgesin1[1:]) / 2
bin_centresin2 = (bin_edgesin2[:-1] + bin_edgesin2[1:]) / 2
bin_centresin3 = (bin_edgesin3[:-1] + bin_edgesin3[1:]) / 2
bin_centresin4 = (bin_edgesin4[:-1] + bin_edgesin4[1:]) / 2
bin_centresin5 = (bin_edgesin5[:-1] + bin_edgesin5[1:]) / 2
bin_centresout0 = (bin_edgesout0[:-1] + bin_edgesout0[1:]) / 2
bin_centresout1 = (bin_edgesout1[:-1] + bin_edgesout1[1:]) / 2
bin_centresout2 = (bin_edgesout2[:-1] + bin_edgesout2[1:]) / 2
bin_centresout3 = (bin_edgesout3[:-1] + bin_edgesout3[1:]) / 2
bin_centresout4 = (bin_edgesout4[:-1] + bin_edgesout4[1:]) / 2
bin_centresout5 = (bin_edgesout5[:-1] + bin_edgesout5[1:]) / 2

smooth_countin0 = gaussian_filter1d(countin0, sigma=1.0)
smooth_countin1 = gaussian_filter1d(countin1, sigma=1.0)
smooth_countin2 = gaussian_filter1d(countin2, sigma=1.0)
smooth_countout0 = gaussian_filter1d(countout0, sigma=1.0)
smooth_countout1 = gaussian_filter1d(countout1, sigma=1.0)
smooth_countout2 = gaussian_filter1d(countout2, sigma=1.0)
smooth_countin3 = gaussian_filter1d(countin3, sigma=1.0)
smooth_countin4 = gaussian_filter1d(countin4, sigma=1.0)
smooth_countin5 = gaussian_filter1d(countin5, sigma=1.0)
smooth_countout3 = gaussian_filter1d(countout3, sigma=1.0)
smooth_countout4 = gaussian_filter1d(countout4, sigma=1.0)
smooth_countout5 = gaussian_filter1d(countout5, sigma=1.0)

plt.figure(3)
fig, axs = plt.subplots(nrows=2, ncols=2, figsize = (8,8), sharey = 'row', sharex = 'col')
'''
axs[0,0].hist(np.log10(TempInflow[0]), range = [3.9,6.1], weights = np.log10(MassInflow[0]), density = True, bins = 50, histtype = 'step', color = colours[0], linestyle = linestyles[0], label = 'SMR')
axs[0,0].hist(np.log10(TempInflow[1]), range = [3.9,6.1], weights = np.log10(MassInflow[1]), density = True, bins = 50, histtype = 'step', color = colours[1], linestyle = linestyles[1], label = '~500pc')
axs[0,0].hist(np.log10(TempInflow[2]), range = [3.9,6.1], weights = np.log10(MassInflow[2]), density = True, bins = 50, histtype = 'step', color = colours[2], linestyle = linestyles[2], label = '~250pc')
axs[0,1].hist(np.log10(TempOutflow[0]), range = [3.9,6.1], weights = np.log10(MassOutflow[0]), density = True, bins = 50, histtype = 'step', color = colours[0], linestyle = linestyles[0], label = 'SMR')
axs[0,1].hist(np.log10(TempOutflow[1]), range = [3.9,6.1], weights = np.log10(MassOutflow[1]), density = True, bins = 50, histtype = 'step', color = colours[1], linestyle = linestyles[1], label = '~1kpc')
axs[0,1].hist(np.log10(TempOutflow[2]), range = [3.9,6.1], weights = np.log10(MassOutflow[2]), density = True, bins = 50, histtype = 'step', color = colours[2], linestyle = linestyles[2], label = '~500pc')

axs[1,0].hist(np.log10(TempInflow[3]), range = [3.9,6.1], weights = np.log10(MassInflow[3]), density = True, bins = 50, histtype = 'step', color = colours[3], linestyle = linestyles[0])
axs[1,0].hist(np.log10(TempInflow[4]), range = [3.9,6.1], weights = np.log10(MassInflow[4]), density = True, bins = 50, histtype = 'step', color = colours[4], linestyle = linestyles[1])
axs[1,0].hist(np.log10(TempInflow[5]), range = [3.9,6.1], weights = np.log10(MassInflow[5]), density = True, bins = 50, histtype = 'step', color = colours[5], linestyle = linestyles[2])
axs[1,1].hist(np.log10(TempOutflow[3]), range = [3.9,6.1], weights = np.log10(MassOutflow[3]), density = True, bins = 50, histtype = 'step', color = colours[3], linestyle = linestyles[0])
axs[1,1].hist(np.log10(TempOutflow[4]), range = [3.9,6.1], weights = np.log10(MassOutflow[4]), density = True, bins = 50, histtype = 'step', color = colours[4], linestyle = linestyles[1])
axs[1,1].hist(np.log10(TempOutflow[5]), range = [3.9,6.1], weights = np.log10(MassOutflow[5]), density = True, bins = 50, histtype = 'step', color = colours[5], linestyle = linestyles[2])
'''
axs[0,0].plot(bin_centresin0, smooth_countin0, color = colours[0], linestyle = linestyles[0], label = 'SMR')
axs[0,0].plot(bin_centresin1, smooth_countin1, color = colours[1], linestyle = linestyles[1], label = '~500pc')
axs[0,0].plot(bin_centresin2, smooth_countin2, color = colours[2], linestyle = linestyles[2], label = '~250pc')
axs[0,1].plot(bin_centresout0, smooth_countout0, color = colours[0], linestyle = linestyles[0])
axs[0,1].plot(bin_centresout1, smooth_countout1, color = colours[1], linestyle = linestyles[1])
axs[0,1].plot(bin_centresout2, smooth_countout2, color = colours[2], linestyle = linestyles[2])

axs[1,0].plot(bin_centresin3, smooth_countin3, color = colours[3], linestyle = linestyles[0], label = 'SMR')
axs[1,0].plot(bin_centresin4, smooth_countin4, color = colours[4], linestyle = linestyles[1], label = '~1kpc')
axs[1,0].plot(bin_centresin5, smooth_countin5, color = colours[5], linestyle = linestyles[2], label = '~500pc')
axs[1,1].plot(bin_centresout3, smooth_countout3, color = colours[3], linestyle = linestyles[0])
axs[1,1].plot(bin_centresout4, smooth_countout4, color = colours[4], linestyle = linestyles[1])
axs[1,1].plot(bin_centresout5, smooth_countout5, color = colours[5], linestyle = linestyles[2])

axs[0,0].set_ylabel(r'PDF ($10^{11}$M$_{\odot}$)')
axs[1,0].set_ylabel(r'PDF ($10^{12}$M$_{\odot}$)')
axs[1,0].set_xlabel(r'Temperature [K]')
axs[1,1].set_xlabel(r'Temperature [K]')

axs[0,0].set_title(r'Inflow')
axs[0,1].set_title(r'Outflow')

axs[0,0].legend()
axs[1,0].legend()

axs[0,0].set_ylim(0,2.5)
axs[1,0].set_ylim(0,2.5)

axs[1,0].set_xlim(4.0,6.5)
axs[1,1].set_xlim(4.0,6.5)

axs[0,0].grid()
axs[1,0].grid()
axs[0,1].grid()
axs[1,1].grid()

for ax in axs.flat:
	ax.label_outer()
	ax.tick_params(axis = 'both', which = 'both', bottom = True, left = True, right = True, top = True, direction = 'in')
	ax.xaxis.set_minor_locator(AutoMinorLocator(4))
	ax.yaxis.set_minor_locator(AutoMinorLocator(4))

axs[0,0].set_yticks([0.5,1,1.5,2,2.5])
axs[1,0].set_yticks([0,0.5,1,1.5,2])

#plt.tick_params(axis = 'both', which = 'both', bottom = True, left = True, right = True, top = True, direction = 'in')

plt.subplots_adjust(hspace = 0, wspace = 0)
plt.savefig(f'./Figures/Paper3Plots/inflow_outflow_temp_{snap[k]}.pdf', dpi = 300, transparent = True, bbox_inches = 'tight')
plt.close("all")

print(lcool_med)
plt.figure(4)
fig, axs = plt.subplots(nrows=2, ncols=1, figsize = (5,10), sharey = 'row')

axs[0].plot(lcool_radii_all[0], np.log10(lcool_med_all[0]), c = 'tab:blue', label = 'SMR')
axs[0].plot(lcool_radii_all[1], np.log10(lcool_med_all[1]), c = 'tab:orange', label = r'${\sim}500$pc', linestyle = '--')
axs[0].plot(lcool_radii_all[2], np.log10(lcool_med_all[2]), c = 'tab:green', label = r'${\sim}250$pc', linestyle = '-.')
axs[1].plot(lcool_radii_all[3], np.log10(lcool_med_all[3]), c = 'tab:red', label = 'SMR')
axs[1].plot(lcool_radii_all[4], np.log10(lcool_med_all[4]), c = 'tab:purple', label = r'${\sim}1$kpc', linestyle = '--')
axs[1].plot(lcool_radii_all[5], np.log10(lcool_med_all[5]), c = 'tab:brown', label = r'${\sim}500$pc', linestyle = '-.')

axs[0].fill_between([0,0.25], 0, 6.8, alpha = 0.2, color = 'grey')
axs[1].fill_between([0,0.25], 0, 6, alpha = 0.2, color = 'grey')

axs[0].set_yticks([1,2,3,4,5,6])
axs[1].set_yticks([0,1,2,3,4,5])

axs[0].set_ylabel(r'Cooling length [pc], ($10^{11}$M$_{\odot}$)')
axs[1].set_ylabel(r'Cooling length [pc], ($10^{12}$M$_{\odot}$)')
axs[1].set_xlabel(r'R/R$_{\mathrm{200c}}$')

axs[0].legend()
axs[1].legend()

axs[0].set_ylim(0,6.8)
axs[1].set_ylim(0,6)

for ax in axs.flat:
	ax.label_outer()
	ax.tick_params(axis = 'both', which = 'both', bottom = True, left = True, right = True, top = True, direction = 'in')
	ax.xaxis.set_minor_locator(AutoMinorLocator(4))
	ax.yaxis.set_minor_locator(AutoMinorLocator(4))
	ax.grid()

plt.subplots_adjust(hspace = 0, wspace = 0)
plt.savefig(f'./Figures/Paper3Plots/cooling_length_radial_profile_{snap[k]}.pdf', dpi = 300, transparent = True, bbox_inches = 'tight')
plt.close("all")

plt.figure(5)
fig, axs = plt.subplots(nrows=1, ncols=2, figsize = (10,5), sharey = 'row')

axs[0].plot(b_radii[0], np.log10(b_med[0]), c = 'tab:blue', label = 'SMR')
axs[0].plot(b_radii[1], np.log10(b_med[1]), c = 'tab:orange', label = r'${\sim}500$pc')
axs[0].plot(b_radii[2], np.log10(b_med[2]), c = 'tab:green', label = r'${\sim}250$pc')
axs[1].plot(b_radii[3], np.log10(b_med[3]), c = 'tab:red', label = 'SMR')
axs[1].plot(b_radii[4], np.log10(b_med[4]), c = 'tab:purple', label = r'${\sim}1$kpc')
axs[1].plot(b_radii[5], np.log10(b_med[5]), c = 'tab:brown', label = r'${\sim}500$pc')

axs[0].set_ylabel(r'Magnetic Field Strength, B [$\mu$G]')
axs[1].set_ylabel(r'Magnetic Field Strength, B [$\mu$G]')
axs[0].set_xlabel(r'R/R$_{\mathrm{200c}}$')
axs[1].set_xlabel(r'R/R$_{\mathrm{200c}}$')

axs[0].legend()
axs[1].legend()

for ax in axs.flat:
	ax.label_outer()
	ax.tick_params(axis = 'both', which = 'both', bottom = True, left = True, right = True, top = True, direction = 'in')
	ax.xaxis.set_minor_locator(AutoMinorLocator(4))
	ax.yaxis.set_minor_locator(AutoMinorLocator(4))
	ax.grid()

plt.subplots_adjust(hspace = 0, wspace = 0)
plt.savefig(f'./Figures/Paper3Plots/Magnetic_Field_radial_profile_{snap[k]}.pdf', dpi = 300, transparent = True, bbox_inches = 'tight')
plt.close("all")

plt.figure(6)
fig, axs = plt.subplots(nrows=1, ncols=2, figsize = (10,5), sharey = 'row')

axs[0].plot(vortout_radii[0], np.log10(vortout_med[0]), c = 'tab:blue', label = 'SMR')
axs[0].plot(vortout_radii[1], np.log10(vortout_med[1]), c = 'tab:orange', label = r'${\sim}500$pc')
axs[0].plot(vortout_radii[2], np.log10(vortout_med[2]), c = 'tab:green', label = r'${\sim}250$pc')
axs[1].plot(vortout_radii[3], np.log10(vortout_med[3]), c = 'tab:red', label = 'SMR')
axs[1].plot(vortout_radii[4], np.log10(vortout_med[4]), c = 'tab:purple', label = r'${\sim}1$kpc')
axs[1].plot(vortout_radii[5], np.log10(vortout_med[5]), c = 'tab:brown', label = r'${\sim}500$pc')

axs[0].set_ylabel(r'Vorticity, $\omega$')
axs[1].set_ylabel(r'Vorticity, $\omega$')
axs[0].set_xlabel(r'R/R$_{\mathrm{200c}}$')
axs[1].set_xlabel(r'R/R$_{\mathrm{200c}}$')

axs[0].legend()
axs[1].legend()

for ax in axs.flat:
	ax.label_outer()
	ax.tick_params(axis = 'both', which = 'both', bottom = True, left = True, right = True, top = True, direction = 'in')
	ax.xaxis.set_minor_locator(AutoMinorLocator(4))
	ax.yaxis.set_minor_locator(AutoMinorLocator(4))
	ax.grid()

plt.subplots_adjust(hspace = 0, wspace = 0)
plt.savefig(f'./Figures/Paper3Plots/Vorticity_radial_profile_{snap[k]}.pdf', dpi = 300, transparent = True, bbox_inches = 'tight')
plt.close("all")

plt.figure(7)
fig, axs = plt.subplots(nrows=2, ncols=2, figsize = (8,8), sharey = True)

axs[0,0].plot(TInflowRadii[0], np.log10(TInflowMed[0]),color = colours[0], linestyle = linestyles[0], label = 'SMR')
axs[0,0].plot(TInflowRadii[1], np.log10(TInflowMed[1]),color = colours[1], linestyle = linestyles[1], label = '~500pc')
axs[0,0].plot(TInflowRadii[2], np.log10(TInflowMed[2]),color = colours[2], linestyle = linestyles[2], label = '~250pc')
axs[0,1].plot(TOutflowRadii[0], np.log10(TOutflowMed[0]), color = colours[0], linestyle = linestyles[0], label = 'SMR')
axs[0,1].plot(TOutflowRadii[1], np.log10(TOutflowMed[1]), color = colours[1], linestyle = linestyles[1], label = '~1kpc')
axs[0,1].plot(TOutflowRadii[2], np.log10(TOutflowMed[2]), color = colours[2], linestyle = linestyles[2], label = '~500pc')

axs[1,0].plot(TInflowRadii[3], np.log10(TInflowMed[3]),color = colours[3], linestyle = linestyles[0])
axs[1,0].plot(TInflowRadii[4], np.log10(TInflowMed[4]),color = colours[4], linestyle = linestyles[1])
axs[1,0].plot(TInflowRadii[5], np.log10(TInflowMed[5]),color = colours[5], linestyle = linestyles[2])
axs[1,1].plot(TOutflowRadii[3], np.log10(TOutflowMed[3]), color = colours[3], linestyle = linestyles[0])
axs[1,1].plot(TOutflowRadii[4], np.log10(TOutflowMed[4]), color = colours[4], linestyle = linestyles[1])
axs[1,1].plot(TOutflowRadii[5], np.log10(TOutflowMed[5]), color = colours[5], linestyle = linestyles[2])

axs[0,0].fill_between(TInflowRadii[0], np.log10(TInflowLower[0]), np.log10(TInflowUpper[0]),color = colours[0], linestyle = linestyles[0], alpha = 0.2)
axs[0,0].fill_between(TInflowRadii[1], np.log10(TInflowLower[1]), np.log10(TInflowUpper[1]),color = colours[1], linestyle = linestyles[1], alpha = 0.2)
axs[0,0].fill_between(TInflowRadii[2], np.log10(TInflowLower[2]), np.log10(TInflowUpper[2]),color = colours[2], linestyle = linestyles[2], alpha = 0.2)
axs[0,1].fill_between(TOutflowRadii[0], np.log10(TOutflowLower[0]), np.log10(TOutflowUpper[0]), color = colours[0], linestyle = linestyles[0], alpha = 0.2)
axs[0,1].fill_between(TOutflowRadii[1], np.log10(TOutflowLower[1]), np.log10(TOutflowUpper[1]), color = colours[1], linestyle = linestyles[1], alpha = 0.2)
axs[0,1].fill_between(TOutflowRadii[2], np.log10(TOutflowLower[2]), np.log10(TOutflowUpper[2]), color = colours[2], linestyle = linestyles[2], alpha = 0.2)

axs[1,0].fill_between(TInflowRadii[3], np.log10(TInflowLower[3]), np.log10(TInflowUpper[3]),color = colours[3], linestyle = linestyles[0], alpha = 0.2)
axs[1,0].fill_between(TInflowRadii[4], np.log10(TInflowLower[4]), np.log10(TInflowUpper[4]),color = colours[4], linestyle = linestyles[1], alpha = 0.2)
axs[1,0].fill_between(TInflowRadii[5], np.log10(TInflowLower[5]), np.log10(TInflowUpper[5]),color = colours[5], linestyle = linestyles[2], alpha = 0.2)
axs[1,1].fill_between(TOutflowRadii[3], np.log10(TOutflowLower[3]), np.log10(TOutflowUpper[3]), color = colours[3], linestyle = linestyles[0], alpha = 0.2)
axs[1,1].fill_between(TOutflowRadii[4], np.log10(TOutflowLower[4]), np.log10(TOutflowUpper[4]), color = colours[4], linestyle = linestyles[1], alpha = 0.2)
axs[1,1].fill_between(TOutflowRadii[5], np.log10(TOutflowLower[5]), np.log10(TOutflowUpper[5]), color = colours[5], linestyle = linestyles[2], alpha = 0.2)

axs[0,0].set_xlabel(r'R/R$_{200c}$')
axs[1,0].set_xlabel(r'R/R$_{200c}$')
axs[1,0].set_ylabel(r'Temperature [K]')
axs[1,1].set_ylabel(r'Temperature [K]')

axs[0,0].set_title(r'Inflow')
axs[0,1].set_title(r'Outflow')

axs[0,0].legend()
axs[1,0].legend()

axs[0,0].set_ylim(4,6)
axs[1,0].set_ylim(4,6)

axs[0,0].grid()
axs[1,0].grid()
axs[0,1].grid()
axs[1,1].grid()

for ax in axs.flat:
	ax.label_outer()
	ax.tick_params(axis = 'both', which = 'both', bottom = True, left = True, right = True, top = True, direction = 'in')
	ax.xaxis.set_minor_locator(AutoMinorLocator(4))
	ax.yaxis.set_minor_locator(AutoMinorLocator(4))

axs[0,0].set_yticks([4.5, 5, 5.5, 6])
axs[1,0].set_yticks([4, 4.5, 5, 5.5])

#plt.tick_params(axis = 'both', which = 'both', bottom = True, left = True, right = True, top = True, direction = 'in')

plt.subplots_adjust(hspace = 0, wspace = 0)
plt.savefig(f'./Figures/Paper3Plots/inflow_outflow_temp_radialprofile_{snap[k]}.pdf', dpi = 300, transparent = True, bbox_inches = 'tight')
plt.close("all")

plt.figure(8)
fig, axs = plt.subplots(nrows=2, ncols=3, figsize = (9,6))

axs[0,0].hist2d(np.log10(abs(VelInflow[0])), np.log10(TempInflow[0]), weights = MassInflow[0], density = True, bins = (50,50))
axs[0,1].hist2d(np.log10(abs(VelInflow[1])), np.log10(TempInflow[1]), weights = MassInflow[1], density = True, bins = (50,50))
axs[0,2].hist2d(np.log10(abs(VelInflow[2])), np.log10(TempInflow[2]), weights = MassInflow[2], density = True, bins = (50,50))

axs[1,0].hist2d(np.log10(VeloOutflow[3]), np.log10(TempOutflow[3]), weights = MassOutflow[3], density = True, bins = (50,50))
axs[1,1].hist2d(np.log10(VeloOutflow[4]), np.log10(TempOutflow[4]), weights = MassOutflow[4], density = True, bins = (50,50))
axs[1,2].hist2d(np.log10(VeloOutflow[5]), np.log10(TempOutflow[5]), weights = MassOutflow[5], density = True, bins = (50,50))

axs[0,0].set_xticks([])
axs[0,1].set_xticks([])

axs[0,0].set_ylabel(r'Temperature (Inflow) [K]')
axs[1,0].set_ylabel(r'Temperature (Outflow) [K]')
axs[1,0].set_xlabel(r'Cooling Length [pc]')
axs[1,1].set_xlabel(r'Cooling Length [pc]')
axs[1,2].set_xlabel(r'Cooling Length [pc]')

for ax in axs.flat:
	ax.label_outer()
	ax.tick_params(axis = 'both', which = 'both', bottom = True, left = True, right = True, top = True, direction = 'in')
	ax.xaxis.set_minor_locator(AutoMinorLocator(4))
	ax.yaxis.set_minor_locator(AutoMinorLocator(4))

plt.subplots_adjust(hspace = 0, wspace = 0)
plt.savefig(f'./Figures/Paper3Plots/vel_temp_2dhist_{snap[k]}.pdf', dpi = 300, transparent = True, bbox_inches = 'tight')
plt.close("all")