Commit 2491d87e authored by Maxime Rey's avatar Maxime Rey
Browse files

New function: 'rad_col2', new version to compute the radial column density with the rascas method.

parent 1ff43e82
......@@ -643,7 +643,7 @@ def rad_prof(genpath, folders, labels, timesteps, var_str='H', overwhat='V', r_m
print('\t Computing... \t timestep', timestep)
# Get radius of cells
_, cells, cell_pos, _ = ras.extract_cells(RamsesDir, timestep, factor=factor, saveinfile=saveinfile)
cell_dx, rho, _, _, vx, vy, vz, xHII, xHeII, xHeIII, _, _, mass, _ = cells
cell_dx, rho, _, _, vx, vy, vz, xHII, _, _, _, _, mass, _ = cells
if var_str=='H':
var = rho*g2Msun
......@@ -1046,6 +1046,78 @@ def rad_col(genpath, folders, timesteps, labels, r_max=150, lmax=18, var_str='HI
plt.tight_layout()
def rad_col2(genpath, folders, timesteps, ions=None, nbins=1000, xlims=None, ylims=None, logx=True, savefig=False):
""" Plots radial column densities using results from rascas. """
from ratatouille import readNsave as ras
from ratatouille import plotutils as put
import matplotlib.pyplot as plt
import numpy as np
import warnings
warnings.filterwarnings("ignore")
# Read ion to plot from first file and assume it's the same for the others if none given.
with open(f'{genpath}{folders[0]}/rascas/{timesteps[0]:05d}/ColDens.dat') as f:
ions_file = f.readline().split('of')[1].strip().split(' ')
# If no ions specified, read all the ions the file.
prop_cycle = plt.rcParams['axes.prop_cycle']
colors = prop_cycle.by_key()['color']
for iion, ion in enumerate(ions if ions else ions_file):
# If not going through all ions, associate the correct index to the ion.
if ions:
iion = ions_file.index(ion)
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(15,9))
# Load column densities and impact parameters into arrays
for ifold, folder in enumerate(folders):
Nion = []
b_kpc = []
for timestep in timesteps:
print(f'{folder}: timestep {timestep}', end='\r')
my_data = np.genfromtxt(f'{genpath+folder}/rascas/{timestep:05d}/ColDens.dat', comments="#", delimiter='')
info = ras.extract_info(timestep, genpath+folder, saveinfile=True)
cu2kpc = info['unit_l'] * info['boxlen'] / 3.086e21
b_kpc.append(my_data[:,0]*cu2kpc)
Nion.append(my_data[:,iion+1])
Nion = np.array(Nion )
b_kpc = np.array(b_kpc)
# Median
min_b, max_b = np.min(b_kpc), np.max(b_kpc)
if logx:
bins = np.logspace(np.log10(min_b), np.log10(max_b), num=nbins)
ax.set_xscale('log')
else:
bins = np.linspace(min_b, max_b, num=nbins)
x_axis = (bins[:-1]+bins[1:])/2
idx = np.digitize(b_kpc,bins, right=False) # Returns an index for each bin. Those below min = 0 and those higher are len(bin)
med_dens = [np.median(Nion[idx==k]) for k in range(1,nbins)] # Remove those with index 0 and k goes to len(bin)-1 -> keep only values between bins.
ax.plot(x_axis, med_dens, label=' '.join(folder.split('_')[1:])[:-1], color=colors[ifold])
# Scatter
perc_min, perc_max = 15.9, 84.1
ax.fill_between(x_axis,[np.nanpercentile(Nion[idx==k],perc_min) for k in range(1,nbins)], [np.nanpercentile(Nion[idx==k],perc_max) for k in range(1,nbins)], alpha=0.25, color=colors[ifold])
put.plot_dexter(ion, path2files='/scratch/Cral/mrey/outputs/data/Ions/', ax=ax)
put.plot_data(ion, path2files='/scratch/Cral/mrey/outputs/data/Ions2/', ax=ax, logx=logx)
ax.legend(frameon=False)
plt.rcParams['font.size'] = 13
plt.yscale('log')
if not xlims:
plt.xlim([min_b, max_b])
else:
plt.xlim(xlims)
plt.ylim(ylims)
plt.xlabel('Impact parameter [kpc]')
plt.ylabel(rf'$N_{{{str(ion)}}}\rm\ [cm^{{-2}}]$')
plt.tight_layout()
plt.show()
if savefig:
fig.savefig(f'N{ion}_VS_R.svg')
def check_folder(genpath, folders, timesteps, r_max=150, lmax=20, split_flows=True, var_str='HI'):
"""
Check the lines saved in an hdf5 file. Done to check the particularly slow rad_col routine.
......@@ -3219,10 +3291,10 @@ def plot_dexter(ion, path2files='data/Ions/', ax=None):
ax.errorbar(impact_param[nonlim], coldens[nonlim], yerr=[yerrm[nonlim], yerrp[nonlim]], fmt='o', capsize=5, color=col, label=fileuh[len(ion)+1:-4])
for indix, TruFal in enumerate(limsup):
if TruFal:
plt.errorbar(impact_param[indix], coldens[indix], yerr=0, uplims=True, color=col)
plt.errorbar(impact_param[indix], coldens[indix], yerr=0, uplims=True, color='grey', capsize=4)
for indix, TruFal in enumerate(liminf):
if TruFal:
plt.errorbar(impact_param[indix], coldens[indix], yerr=0, lolims=True, color=col)
plt.errorbar(impact_param[indix], coldens[indix], yerr=0, lolims=True, color='grey', capsize=4)
def plot_data(ion, path2files='data/Ions2/', ax=None, logx=True):
......@@ -3262,10 +3334,10 @@ def plot_data(ion, path2files='data/Ions2/', ax=None, logx=True):
ax.errorbar(impact_param[nonlim], 10**coldens[nonlim], yerr=yerr, xerr=xerr[nonlim], fmt='o', capsize=5, color=col, label=str_dat)
for indix, TruFal in enumerate(limsup):
if TruFal:
plt.errorbar(impact_param[indix], 10**coldens[indix], yerr=0, uplims=True, color=col)
plt.errorbar(impact_param[indix], 10**coldens[indix], yerr=0, uplims=True, color='grey', capsize=4)
for indix, TruFal in enumerate(liminf):
if TruFal:
plt.errorbar(impact_param[indix], 10**coldens[indix], yerr=0, lolims=True, color=col)
plt.errorbar(impact_param[indix], 10**coldens[indix], yerr=0, lolims=True, color='grey', capsize=4)
def runtime(RamsesDir, logstr='logs'):
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment