HighResColDensQuant.py

import matplotlib
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, cum_avg, cum_sum
from statistics import median
from matplotlib.ticker import AutoMinorLocator
import csv
import time
import pandas as pd

plt.switch_backend('agg')
start = time.time()
snap = np.arange(77,128,1)
snap = np.array([int(sys.argv[1])])
halo_d = ['level3_cgm_1e11', 'level4_cgm']
halo_d = ['level3_cgm_1e11']
n_Hx, n_H_edgex, n_Hix, n_CIVx, n_SiIIx, n_Ovix, n_MgIIx, n_Hiex, n_CIVex, n_SiIIex, n_Oviex, n_MgIIex, Tx, metx, T_edgex, met_edgex, velx, volx = [], [], [], [], [], [], [], [], [], [], [], [], [], [], [], [], [], []
n_Hy, n_H_edgey, n_Hiy, n_CIVy, n_SiIIy, n_Oviy, n_MgIIy, n_Hiey, n_CIVey, n_SiIIey, n_Oviey, n_MgIIey, Ty, mety, T_edgey, met_edgey, vely, voly = [], [], [], [], [], [], [], [], [], [], [], [], [], [], [], [], [], []
densgrid, dens_edgegrid, n_Hgrid, n_H_edgegrid, n_Higrid, n_CIVgrid, n_SiIIgrid, n_Ovigrid, n_MgIIgrid, n_Hiegrid, n_CIVegrid, n_SiIIegrid, n_Oviegrid, n_MgIIegrid, Tgrid, metgrid, T_edgegrid, met_edgegrid, velgrid, volgrid = [], [], [], [], [], [], [], [], [], [], [], [], [], [], [], [], [], [], [], []
x1, x2 = [], []
virialMass, coolRate, temp = [], [], []

temp_HI, tupper_HI, tlower_HI, dense_HI, dupper_HI, dlower_HI, met_HI, mupper_HI, mlower_HI, vel_HI, vupper_HI, vlower_HI = [], [], [], [], [], [], [], [], [], [], [], []

temp_CIV, temp_OVI, temp_SiII, temp_MgII = [], [], [], []
tupper_CIV, tupper_OVI, tupper_SiII, tupper_MgII = [], [], [], []
tlower_CIV, tlower_OVI, tlower_SiII, tlower_MgII = [], [], [], []

dense_CIV, dense_OVI, dense_SiII, dense_MgII = [], [], [], []
dupper_CIV, dupper_OVI, dupper_SiII, dupper_MgII = [], [], [], []
dlower_CIV, dlower_OVI, dlower_SiII, dlower_MgII = [], [], [], []

met_CIV, met_OVI, met_SiII, met_MgII = [], [], [], []
mupper_CIV, mupper_OVI, mupper_SiII, mupper_MgII = [], [], [], []
mlower_CIV, mlower_OVI, mlower_SiII, mlower_MgII = [], [], [], []

vel_CIV, vel_OVI, vel_SiII, vel_MgII = [], [], [], []
vupper_CIV, vupper_OVI, vupper_SiII, vupper_MgII = [], [], [], []
vlower_CIV, vlower_OVI, vlower_SiII, vlower_MgII = [], [], [], []

HI_mass_med, HI_mass_upp, HI_mass_low, HI_mass_rad = [], [], [], []
n_HI, n_HI_radii, n_HI_low, n_HI_upp = [], [], [], []
n_CIV, n_CIV_radii, n_CIV_low, n_CIV_upp  = [], [], [], []
n_CI, n_CI_radii, n_CI_low, n_CI_upp  = [], [], [], []
n_CII, n_CII_radii, n_CII_low, n_CII_upp  = [], [], [], []
n_CIII, n_CIII_radii, n_CIII_low, n_CIII_upp  = [], [], [], []
n_CV, n_CV_radii, n_CV_low, n_CV_upp  = [], [], [], []
n_CVI, n_CVI_radii, n_CVI_low, n_CVI_upp  = [], [], [], []
n_MgI, n_MgI_radii = [], []
n_HI_sr, n_HI_sr_radii = [], []
n_OVI_sr, n_OVI_sr_radii = [], []
n_CIV_sr, n_CIV_sr_radii = [], []
n_OVI, n_OVI_radii, n_OVI_low, n_OVI_upp  = [], [], [], []
n_SiII, n_SiII_radii, n_SiII_low, n_SiII_upp  = [], [], [], []
haloes, lCool, coolmass = [], [], []

cFracHI13, cFracHI14, cFracHI15, cFracHI16, cFracHI17  = [], [], [], [], []
cFracCII13, cFracCI14, cFracCI15, cFracCI16, cFracCI17  = [], [], [], [], []
cFracSII13, cFracSII14, cFracSII15, cFracSII16, cFracSII17  = [], [], [], [], []
cFracSiII13, cFracSiI14, cFracSiI15, cFracSiI16, cFracSiI17  = [], [], [], [], []
cFracOIII13, cFracOIII14, cFracOIII15, cFracOIII16, cFracOIII17  = [], [], [], [], []

cFracHIRad, cFracCIVRad, cFracCIRad, cFracCIIRad, cFracCIIIRad, cFracCVRad, cFracCVIRad, cFracSiIIRad, cFracMgIIRad, cFracOVIRad  = [], [], [], [], [], [], [], [], [], []
colDensCovRad = np.arange(13,21,0.01)
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	= [], [], [], [], [], []
coolingmass = []
HI, CIV, OVI, SiII, SR = [], [], [], [], []

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

fHI = open('./highResSims/HIz.txt', 'w')
fCIV = open('./highResSims/CIVz.txt', 'w')
fOVI = open('./highResSims/OVIz.txt', 'w')
writerHI = csv.writer(fHI)
writerCIV = csv.writer(fCIV)
writerOVI = csv.writer(fOVI)
header = ['h10_standard', 'h10_500pc', 'h10_250pc', 'h5_standard', 'h5_1kpc', 'h5_500pc']
writerHI.writerow(header)
writerCIV.writerow(header)
writerOVI.writerow(header)

def col_dens_dist(ion):
	global ionw	#Please ignore

	#Initialise 2-D arrays for the radius and column density (column density is added to a new array which this function spits out)
	s = (len(load.data[str(ion) + 'proj']['x']), len(load.data[str(ion) + 'proj']['y']))

	load.data[str(ion) + 'rad'], load.data[str(ion) + 'ionw']	= np.zeros(s), np.zeros(s)

	#Iterate over x and y coordinates generated by get_Aslice. This creates a magnitude for each pixel
	for n in range(len(load.data[str(ion) + 'proj']['x'])):
		for m in range(len(load.data[str(ion) + 'proj']['y'])):
			load.data[str(ion) + 'rad'][n, m]		= np.sqrt(np.add(load.data[str(ion) + 'proj']['x'][n]**2, load.data[str(ion) + 'proj']['y'][m]**2))

	#Iterate over the length of x and y coordinates again, this time to generate a new array of column density for each pixel in order for the
	#column density values to match the position of the new array
	for n in range(len(load.data[str(ion) + 'proj']['x'])-1):
		for m in range(len(load.data[str(ion) + 'proj']['y'])-1):
			load.data[str(ion) + 'ionw'][n, m]		= load.data[str(ion) + '0'][n,m]

	load.data[str(ion) + 'rad'], load.data[str(ion) + 'ionw']	= load.data[str(ion) + 'rad'].flatten(), load.data[str(ion) + 'ionw'].flatten()

	load.data[str(ion) + 'rad'].tolist()
	load.data[str(ion) + 'rad'] = np.divide(load.data[str(ion) + 'rad'], virialRadius)	#Normalise new radius array to virial radius of halo

	#Remove any nan values from column density array
	nan_array = np.isnan(load.data[str(ion) + 'ionw'])
	not_nan_array = ~nan_array

	load.data[str(ion) + 'ionw']	= load.data[str(ion) + 'ionw'][not_nan_array]
	load.data[str(ion) + 'ionw'].tolist()

	ionw, = np.where((load.data[str(ion) + 'rad'] <= 1))	#Create conditional where whatever the name of the radius variable is only including the virial radius

	return np.log10(load.data[str(ion) + 'ionw'][ionw], where = load.data[str(ion) + 'ionw'][ionw] != 0), load.data[str(ion) + 'rad'][ionw], load.data[str(ion) + 'ionw'][ionw]

for k in range(len(snap)):
	preHIz, preCIVz, preOVIz = [], [], []
	for i in range(len(halo_d)):
		if halo_d[i] == 'level3_cgm_1e11':
			halo = ['h10_standard', 'h10_500pc', 'h10_250pc']
			pixres = 0.2
		if halo_d[i] == 'level4_cgm':
			halo = ['h5_standard', 'h5_1kpc', 'h5_500pc']
			pixres = 0.4
		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)

			load.calc_sf_indizes(sf)
			load.select_halo(sf, do_rotation = True)

			haloes.append(halo[j])
			redshift.append(load.redshift)
			print(load.redshift)
			haloid = 0
			virialRadius = sf.data['frc2'][0]*1e3
			cent = sf.data['fpos'][haloid,:]
			virialMass.append(sf.data['fmc2'][0]*1e10)
			stelMass.append(sf.data['fmty'][0,4]*1e10)

			#Investigating the birth of stars
			stellarAge	= load.age[load.age > 0]	#Time when stars formed. load.age > 0 added to exclude wind phase gas cells. This is in terms of the scale factor
			whereStars,	= np.where(load.age > 0)
			birthRedshift = (1/stellarAge) - 1
			birthMass = np.multiply(load.data['gima'][whereStars], 1e10)
			stellarBirthRedshift.append(birthRedshift)
			stellarInitMass.append(birthMass)

			virialRadius2.append(virialRadius)
			df = h5py.File('fg2009_ss_hr.h5')

			load.pos *= 1e3
			load.vol *= 1e9 # Mpc^3 to kpc^3
			load.mass *= 1e10 # 1e10 Msun to Msun

			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

			load.data['volume'] = load.vol

			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)

			theta = np.linspace(0, 2*pi, 100)

			r = virialRadius
			x1.append(r*np.sin(theta))
			x2.append(r*np.cos(theta))

			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

			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))
			#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']

			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.3) & (normRadius[:len(load.ne)] <= 1) & (load.halo[:len(load.ne)] == 0) & (load.subhalo[:len(load.ne)] <= 0))
			load.data['Tcool'] = load.data['T'][fraccoldgas]
			coolmass.append(load.mass[fraccoldgas])
			coolinggas, = np.where((load.gcol<0) & (load.sfr[:len(load.ne)] == 0) & (normRadius[:len(load.ne)] >= 0.3) & (normRadius[:len(load.ne)] <= 1) & (load.halo[:len(load.ne)] == 0) & (load.subhalo[:len(load.ne)] <= 0))
			frac_cold.append(len(load.mass[fraccoldgas])/len(load.mass[w]))

			meanweightion = 1.2
			tcool = load.u[coolinggas] * 1e10 * load.data['dense'][coolinggas] / (abs(load.gcol[coolinggas]) * load.data['n_H'][coolinggas]**2)
			csound = sqrt(5./3. * k_b * load.data['T'][coolinggas] / ((load.nh[coolinggas] * 0.6 + (1-load.nh[coolinggas]) * meanweightion) * m_p))

			coolingmass.append(load.mass[coolinggas])
			
			lcool = tcool * csound / kpcincm
			lCool.append(lcool)
			'''
			boxsize		= 2*virialRadius
			boxloscol	= 2*virialRadius
			if halo_d[i] == 'level3_cgm_1e11':
				pixres = 0.5
			elif halo_d[i] == 'level4_cgm':
				pixres = 1
			pixres = 3
			pixreslos	= 0.5 #0.2 #0.1
			Mpcincm		= c.parsec*1e6
			convert		= pixreslos*Mpcincm/1e3
			imgcent		= [0,0,0]
			edge		= [2,0] #[0,1]
			face		= [2,1]

			load.data['n_HI']								= ((load.data['dense'][:len(load.ne)] / (elements_mass[0]*c.amu)) * load.gmet[:,0][:len(load.ne)]) * load.data['nh'][:len(load.ne)]
			start = time.time()
			load.data['n_SiII'], SiII_mass, SiII_mass_frac	= ionisation('Si', 7, 1, load.redshift, df, elements_mass, load = load)
			load.data['n_CIV'], CIV_mass, CIV_mass_frac		= ionisation('C', 2, 3, load.redshift, df, elements_mass, load = load)
			load.data['n_OVI'], OVI_mass, OVI_mass_frac		= ionisation('O', 4, 5, load.redshift, df, elements_mass, load = load)
			end = time.time()
			print(f'Elapsed time {end-start}')
			HI_mass = load.data['mass'][:len(load.ne)]*load.gmet[:,0][:len(load.ne)]*load.data['nh'][:len(load.ne)]

#			HImassmed, HImassrad = avg(normRadius[w][:len(load.ne)], HI_mass[w][:len(load.ne)], min(normRadius[w][:len(load.ne)]), max(normRadius[w][:len(load.ne)]))

#			HI_mass_med.append(HImassmed), HI_mass_rad.append(HImassrad)
			
			load.data['HIproj']		= load.get_Aslice("n_HI", box = [boxsize,boxsize], center = imgcent, nx = int(boxsize/pixres), ny = int(boxsize/pixres), nz = int(boxloscol/pixreslos), boxz = boxloscol, axes = edge, proj = True, numthreads = 16)
			load.data['CIVproj']	= load.get_Aslice("n_CIV", box = [boxsize,boxsize], center = imgcent, nx = int(boxsize/pixres), ny = int(boxsize/pixres), nz = int(boxloscol/pixreslos), boxz = boxloscol, axes = edge, proj = True, numthreads = 16)
			load.data['OVIproj']	= load.get_Aslice("n_OVI", box = [boxsize,boxsize], center = imgcent, nx = int(boxsize/pixres), ny = int(boxsize/pixres), nz = int(boxloscol/pixreslos), boxz = boxloscol, axes = edge, proj = True, numthreads = 16)
			load.data['SiIIproj']	= load.get_Aslice("n_SiII", box = [boxsize,boxsize], center = imgcent, nx = int(boxsize/pixres), ny = int(boxsize/pixres), nz = int(boxloscol/pixreslos), boxz = boxloscol, axes = edge, proj = True, numthreads = 16)
			load.data['srproj']		= load.get_Aslice("sr", box = [boxsize,boxsize], center = imgcent, nx = int(boxsize/pixres), ny = int(boxsize/pixres), nz = int(boxloscol/pixreslos), boxz = boxloscol, axes = edge, proj = True, numthreads = 16)
			
			load.data['HI0']	= load.data['HIproj']['grid']*convert
			load.data['HI']		= load.data['HIproj']['grid'].flatten()*convert

			load.data['CIV0']	= load.data['CIVproj']['grid']*convert
			load.data['CIV']	= load.data['CIVproj']['grid'].flatten()*convert

			load.data['OVI0']	= load.data['OVIproj']['grid']*convert
			load.data['OVI']	= load.data['OVIproj']['grid'].flatten()*convert

			load.data['SiII0']	= load.data['SiIIproj']['grid']*convert
			load.data['SiII']	= load.data['SiIIproj']['grid'].flatten()*convert

			load.data['sr']		= load.data['srproj']['grid'].flatten()

			load.data['HIionwE'], load.data['HIrad'], _		= col_dens_dist('HI')
			load.data['CIVionwE'], load.data['CIVrad'], _	= col_dens_dist('CIV')
			load.data['OVIionwE'], load.data['OVIrad'], _	= col_dens_dist('OVI')
			load.data['SiIIionwE'], load.data['SiIIrad'], _	= col_dens_dist('SiII')

			print(load.data['srproj']['grid'][0])

			HI.append(load.data['HI0'].flatten()), CIV.append(load.data['CIV0'].flatten()), OVI.append(load.data['OVI0'].flatten()), SiII.append(load.data['SiII0'].flatten()), SR.append(load.data['sr'])
			
			HI_med, HI_upper, HI_lower, HI_radii			= avg(load.data['sr'], load.data['HI0'].flatten(), min(load.data['sr']), max(load.data['sr']), 100)
			CIV_med, CIV_upper, CIV_lower, CIV_radii		= avg(load.data['sr'], load.data['CIV0'].flatten(), min(load.data['sr']), max(load.data['sr']), 100)
			OVI_med, OVI_upper, OVI_lower, OVI_radii		= avg(load.data['sr'], load.data['OVI0'].flatten(), min(load.data['sr']), max(load.data['sr']), 100)
			SiII_med, SiII_upper, SiII_lower, SiII_radii	= avg(load.data['sr'], load.data['SiII0'].flatten(), min(load.data['sr']), max(load.data['sr']), 100)

			n_HI_low.append(HI_lower), n_HI_upp.append(HI_upper), n_HI.append(HI_med), n_HI_radii.append(HI_radii)
			n_CIV_low.append(CIV_lower), n_CIV_upp.append(CIV_upper), n_CIV.append(CIV_med), n_CIV_radii.append(CIV_radii)
			n_OVI_low.append(OVI_lower), n_OVI_upp.append(OVI_upper), n_OVI.append(OVI_med), n_OVI_radii.append(OVI_radii)
			n_SiII_low.append(SiII_lower), n_SiII_upp.append(SiII_upper), n_SiII.append(SiII_med), n_SiII_radii.append(SiII_radii)
			
			cFracHI, cFracCIV, cFracCI, cFracCII, cFracCIII, cFracCV, cFracCVI, cFracSiII, cFracOVI, cFracMgII, colDensCov  = [], [], [], [], [], [], [], [], [], [], []

			for covfrac in colDensCovRad:
				cFracHI.append(np.divide(len(load.data['HIionwE'][(load.data['HIionwE'] >= covfrac)]), len(load.data['HIionwE'])))
				cFracCIV.append(np.divide(len(load.data['CIVionwE'][(load.data['CIVionwE'] >= covfrac)]), len(load.data['CIVionwE'])))
				cFracOVI.append(np.divide(len(load.data['OVIionwE'][(load.data['OVIionwE'] >= covfrac)]), len(load.data['OVIionwE'])))
				cFracSiII.append(np.divide(len(load.data['SiIIionwE'][(load.data['SiIIionwE'] >= covfrac)]), len(load.data['SiIIionwE'])))

			cFracHIRad.append(cFracHI)
			cFracCIVRad.append(cFracCIV)
			cFracOVIRad.append(cFracOVI)
			cFracSiIIRad.append(cFracSiII)
			'''
			temp_HI, tupper_HI, tlower_HI, tradii_HI = avgWeighted(normRadius[w], load.data['T'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, HI_mass[w])
			dense_HI, dupper_HI, dlower_HI, dradii_HI = avgWeighted(normRadius[w], load.data['n_H'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, HI_mass[w])
			met_HI, mupper_HI, mlower_HI, mradii_HI = avgWeighted(normRadius[w], load.data['gz'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, HI_mass[w])
			vel_HI, vupper_HI, vlower_HI, vradii_HI = avgWeighted(normRadius[w], load.data['velocity'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, HI_mass[w])

			temp_SiII, tupper_SiII, tlower_SiII, tradii_SiII = avgWeighted(normRadius[w], load.data['T'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, SiII_mass[w])
			dense_SiII, dupper_SiII, dlower_SiII, dradii_SiII = avgWeighted(normRadius[w], load.data['n_H'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, SiII_mass[w])
			met_SiII, mupper_SiII, mlower_SiII, mradii_SiII = avgWeighted(normRadius[w], load.data['gz'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, SiII_mass[w])
			vel_SiII, vupper_SiII, vlower_SiII, vradii_SiII = avgWeighted(normRadius[w], load.data['velocity'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, SiII_mass[w])

			temp_CIV, tupper_CIV, tlower_CIV, tradii_CIV = avgWeighted(normRadius[w], load.data['T'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, CIV_mass[w])
			dense_CIV, dupper_CIV, dlower_CIV, dradii_CIV = avgWeighted(normRadius[w], load.data['n_H'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, CIV_mass[w])
			met_CIV, mupper_CIV, mlower_CIV, mradii_CIV = avgWeighted(normRadius[w], load.data['gz'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, CIV_mass[w])
			vel_CIV, vupper_CIV, vlower_CIV, vradii_CIV = avgWeighted(normRadius[w], load.data['velocity'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, CIV_mass[w])

			temp_OVI, tupper_OVI, tlower_OVI, tradii_OVI = avgWeighted(normRadius[w], load.data['T'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, OVI_mass[w])
			dense_OVI, dupper_OVI, dlower_OVI, dradii_OVI = avgWeighted(normRadius[w], load.data['n_H'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, OVI_mass[w])
			met_OVI, mupper_OVI, mlower_OVI, mradii_OVI = avgWeighted(normRadius[w], load.data['gz'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, OVI_mass[w])
			vel_OVI, vupper_OVI, vlower_OVI, vradii_OVI = avgWeighted(normRadius[w], load.data['velocity'][w], min(normRadius[w]), max(normRadius[w]), min(normRadius[w]), 100, OVI_mass[w])
			'''
			preHIz.append(median(load.data['HIionwE'][(load.data['HIrad']>0.3) & (load.data['HIrad']<1)])), preCIVz.append(median(load.data['CIVionwE'][(load.data['CIVrad']>0.3) & (load.data['CIVrad']<1)])), preOVIz.append(median(load.data['OVIionwE'][(load.data['OVIrad']>0.3) & (load.data['OVIrad']<1)])), 

#	writerHI.writerow(preHIz), writerCIV.writerow(preCIVz), writerOVI.writerow(preOVIz)
#	preHIz.clear(), preCIVz.clear(), preOVIz.clear()

cmaprho = plt.cm.get_cmap('viridis')
cmapT = plt.cm.get_cmap('Spectral_r')
cmapmet = cmo.cm.thermal
	
colours = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', 'tab:brown']
labels = ['SMR', '500pc', '250pc']
labels2 = ['SMR', '1kpc', '500pc']

linestyles = ['-', '--', '-.', '-', '--', '-.']
print(n_HI[0][1:], n_HI_radii[0][1:])

fig, axs = plt.subplots(2, 2, figsize = (8,8), sharex = 'col', sharey = 'row')
#axs[0,0].hist2d(np.log10(SR[0]), np.log10(HI[0]), bins = (60,60), norm = matplotlib.colors.LogNorm())
#axs[0,1].hist2d(np.log10(SR[0]), np.log10(CIV[0]), bins = (60,60), norm = matplotlib.colors.LogNorm())
#axs[1,0].hist2d(np.log10(SR[0]), np.log10(OVI[0]), bins = (60,60), norm = matplotlib.colors.LogNorm())
#axs[1,1].hist2d(np.log10(SR[0]), np.log10(SiII[0]), bins = (60,60), norm = matplotlib.colors.LogNorm())
for i in range(3):
	axs[0,0].plot(np.log10(n_HI_radii[i][1:]), np.log10(n_HI[i][1:]))
	axs[0,1].plot(np.log10(n_CIV_radii[i][1:]), np.log10(n_CIV[i][1:]))
	axs[1,0].plot(np.log10(n_OVI_radii[i][1:]), np.log10(n_OVI[i][1:]))
	axs[1,1].plot(np.log10(n_SiII_radii[i][1:]), np.log10(n_SiII[i][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.0, wspace = 0.0)
plt.savefig(f'./Figures/Paper3Plots/BIGGY_test.pdf', dpi = 300, transparent = True, bbox_inches = 'tight')
plt.close('fig')
'''
plt.figure(2)
fig, axs = plt.subplots(4, len(haloes), figsize = (2*len(haloes),9), sharex = 'col', sharey = 'row')
for i in range(len(haloes)):

	axs[0,i].plot(te_radii[i], te_med[i], c = str(color_cycle[i]), linestyle = linestyles[i])
	axs[1,i].plot(de_radii[i], de_med[i], c = str(color_cycle[i]), linestyle = linestyles[i])
	axs[2,i].plot(me_radii[i], me_med[i], c = str(color_cycle[i]), linestyle = linestyles[i])
	axs[3,i].plot(ve_radii_in[i], ve_med_in[i], c = str(color_cycle[i]), label = 'Inflows', linestyle = linestyles[i])
	axs[3,i].plot(ve_radii_out[i], ve_med_out[i], ls=':', c = str(color_cycle[i]), label = 'Outflows')

	axs[0,i].set_ylim(4,6.5)
	axs[1,i].set_ylim(-6,-1)
	axs[2,i].set_ylim(-3,1)
	axs[3,i].set_ylim(0.9,2.1)
	axs[0,i].set_xlim(0,1)
	axs[1,i].set_xlim(0,1)
	axs[2,i].set_xlim(0,1)
	axs[3,i].set_xlim(0,1)
	axs[3,i].set_xlabel(r'R/R$_{\mathrm{200c}}$')
#	axs[0,i].set_title(mass[i])
	axs[0,i].grid()
	axs[1,i].grid()
	axs[2,i].grid()
	axs[3,i].grid()

for ax in axs.flat:
	ax.label_outer()
	ax.set_xticks([0.25,0.5,0.75])
	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_title(r'$10^{11}$M$_{\odot}$,SMR')
axs[0,1].set_title(r'$10^{11}$M$_{\odot}$,${\sim}500$pc')
axs[0,2].set_title(r'$10^{11}$M$_{\odot}$,${\sim}250$pc')
axs[0,3].set_title(r'$10^{12}$M$_{\odot}$,SMR')
axs[0,4].set_title(r'$10^{12}$M$_{\odot}$,${\sim}1$kpc')
axs[0,5].set_title(r'$10^{12}$M$_{\odot}$,${\sim}500$pc')

axs[3,-1].legend()
axs[0,0].legend()
axs[0,0].set_ylabel(r'Log$_{10}$T [K]')
axs[1,0].set_ylabel(r'Log$_{10}$n$_{\mathrm{H}}$ [$\mathrm{cm}^{-3}$]')
axs[2,0].set_ylabel(r'Log$_{10}$Z/Z$_{\odot}$')
axs[3,0].set_ylabel(r'Log$_{10}$V$_{\mathrm{rad}}$ [km s$^{-1}$]')

axs[0,0].set_yticks([4.5,5,5.5,6.0,6.5])
axs[1,0].set_yticks([-5,-4,-3,-2])
axs[2,0].set_yticks([-2,-1,0])
axs[3,0].set_yticks([1,1.5,2])

axs[3,0].set_xticks([0.00,0.25,0.5,0.75])
axs[3,-1].set_xticks([0.25,0.5,0.75,1])

plt.subplots_adjust(hspace = 0.0, wspace = 0.0)
plt.savefig(f'./Figures/Paper3Plots/collated_physical_properties_medians_high_res.pdf', dpi = 300, transparent = True, bbox_inches = 'tight')
plt.close('all')
'''
plt.figure(13)
fig, axs = plt.subplots(nrows=2, ncols=4, figsize = (10,6), sharey = 'row')
for i in range(3):
	axs[0,0].plot(colDensCovRad, cFracHIRad[i], c = colours[i], label = labels[i])
	axs[1,0].plot(colDensCovRad, cFracHIRad[i+3], c = colours[i+3], label = labels2[i])
	axs[0,1].plot(colDensCovRad, cFracSiIIRad[i], c = colours[i], label = labels[i])
	axs[1,1].plot(colDensCovRad, cFracSiIIRad[i+3], c = colours[i+3], label = labels2[i])
	axs[0,2].plot(colDensCovRad, cFracCIVRad[i], c = colours[i], label = labels[i])
	axs[1,2].plot(colDensCovRad, cFracCIVRad[i+3], c = colours[i+3], label = labels2[i])
	axs[0,3].plot(colDensCovRad, cFracOVIRad[i], c = colours[i], label = labels[i])
	axs[1,3].plot(colDensCovRad, cFracOVIRad[i+3], c = colours[i+3], label = labels2[i])

	axs[0,i].set_ylabel(r'Covering Fraction, $f_{\mathrm{N}}$')
	axs[1,i].set_ylabel(r'Covering Fraction, $f_{\mathrm{N}}$')

axs[0,-1].legend(fontsize = '10', loc = 'upper right')
axs[1,-1].legend(fontsize = '10', loc = 'upper right')

axs[-1,0].set_xlabel(r'Log$_{10}$N$_{\mathrm{HI}}$')
axs[-1,1].set_xlabel(r'Log$_{10}$N$_{\mathrm{SiII}}$')
axs[-1,2].set_xlabel(r'Log$_{10}$N$_{\mathrm{CIV}}$')
axs[-1,3].set_xlabel(r'Log$_{10}$N$_{\mathrm{OVI}}$')

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.label_outer()
	ax.grid()
	ax.set_xlim([13,21])
	ax.set_xticks([15,17,19])

axs[-1,0].set_xticks([13,15,17,19])
axs[-1,-1].set_xticks([15,17,19,21])

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

end = time.time
print(f'Elapsed time = {end-start}/60')