Commit 224b7b9a authored by Maxime Rey's avatar Maxime Rey
Browse files

BUGFIX: 'rad_col' was only plotting outflow cells ! Due to combined mistake...

BUGFIX: 'rad_col' was only plotting outflow cells ! Due to combined mistake between jibs and a wrong string in the routine... all previous rad_col.h5 files should be deleted.
parent b715b5f6
......@@ -1282,7 +1282,7 @@ def rad_Zfrac(genpath, folders, labels, timesteps, r_max=150, factor=2, nbins=10
if savefig:
fig.savefig(f'rad_Zprof.svg')
def rad_col(genpath, folders, timesteps, labels, r_max=150, lmax=20, var_str='HI', split_flows=False, nbins = 1000, factor = 2,
def rad_col(genpath, folders, timesteps, labels, r_max=150, lmax=20, var_str='HI', split_flows=False, nbins=1000, factor=2, rmsat=False,
plot=True, xlims=None, ylims=None, saveinfile=True, no_err=True, logx=True):
"""
Slow routine.
......@@ -1322,8 +1322,7 @@ def rad_col(genpath, folders, timesteps, labels, r_max=150, lmax=20, var_str='HI
print('Folder', folder)
# Initialisation
all_varbin, array_x_tstep, array_xyz = [], [], []
all_varout, all_varin = [], []
all_varbin, all_varout, all_varin = [], [], []
# Create file to store data.
filepath = f'{RamsesDir}/ratadat/'
......@@ -1346,6 +1345,8 @@ def rad_col(genpath, folders, timesteps, labels, r_max=150, lmax=20, var_str='HI
if split_flows:
name_out = f'out_{var_str}_{timestep:05d}_{r_max}_{lmax}'
name_in = f'in_{var_str}_{timestep:05d}_{r_max}_{lmax}'
if rmsat:
name_out, name_in = name_out+'_nosat', name_in+'_nosat'
try:
for proj in ['x', 'y', 'z']:
all_varout.append(file[f'{name_out}_{proj}'][:])
......@@ -1356,11 +1357,11 @@ def rad_col(genpath, folders, timesteps, labels, r_max=150, lmax=20, var_str='HI
pass
else:
base_name = f'out_{var_str}_{timestep:05d}_{r_max}_{lmax}'
if rmsat:
base_name = base_name+'_nosat'
try:
for proj in ['x', 'y', 'z']:
all_varbin.append(file[f'{base_name}_{proj}'][:])
array_xyz. append(file[f'{base_name}_{proj}'][:])
array_x_tstep. append(file[f'{base_name}_x'][:])
continue
except KeyError:
pass
......@@ -1369,7 +1370,7 @@ def rad_col(genpath, folders, timesteps, labels, r_max=150, lmax=20, var_str='HI
info = ras.extract_info(timestep, RamsesDir, saveinfile=saveinfile)
cu2cm = info['unit_l'] * info['boxlen']
cu2kpc = cu2cm / kpc2cm
center, r200 = ras.get_r200_cen(RamsesDir, timestep)
center, _ = ras.get_r200_cen(RamsesDir, timestep)
_, rad = ras.get_rad(RamsesDir, timestep, factor)
rad_kpc = rad*cu2kpc
if (rad_kpc<r_max):
......@@ -1377,14 +1378,18 @@ def rad_col(genpath, folders, timesteps, labels, r_max=150, lmax=20, var_str='HI
# Load data and select satellites
_, cells, cell_pos, cell_l = ras.extract_cells(RamsesDir,timestep,factor=factor,saveinfile=saveinfile,verbose=False)
sat = is_sat(RamsesDir,timestep, factor=factor)
nsat = np.logical_not(sat)
if rmsat:
sat = is_sat(RamsesDir,timestep, factor=factor)
nsat = np.logical_not(sat)
# Select cells within r_max
center_cm = [cen*cu2cm for cen in center]
vec_r = np.array([cell_pos[:,0]-center_cm[0], cell_pos[:,1]-center_cm[1], cell_pos[:,2]-center_cm[2]]).transpose()
norm_r = np.linalg.norm(vec_r,2, axis=1)
# Retsrict to selection
cells, cell_pos, cell_l = [truc[nsat & (norm_r<r_max_cm)] for truc in cells], cell_pos[nsat & (norm_r<r_max_cm)], cell_l[nsat & (norm_r<r_max_cm)]
if rmsat:
cells, cell_pos, cell_l = [truc[nsat & (norm_r<r_max_cm)] for truc in cells], cell_pos[nsat & (norm_r<r_max_cm)], cell_l[nsat & (norm_r<r_max_cm)]
else:
cells, cell_pos, cell_l = [truc[(norm_r<r_max_cm)] for truc in cells], cell_pos[(norm_r<r_max_cm)], cell_l[(norm_r<r_max_cm)]
_, _, nH, _, vx, vy, vz, xHII, _, _, _, _, _, _ = cells
if var_str=='H':
......@@ -1395,10 +1400,16 @@ def rad_col(genpath, folders, timesteps, labels, r_max=150, lmax=20, var_str='HI
try:
var = ras.extract_ions(RamsesDir, timestep, var_str, factor)
var.reshape((len(var)))
var = var[nsat & (norm_r<r_max_cm)]
if rmsat:
var = var[nsat & (norm_r<r_max_cm)]
else:
var = var[(norm_r<r_max_cm)]
except:
raise ValueError(f'The string {var_str} does not seem to have been computed beforehand.')
vec_r, norm_r = vec_r[nsat & (norm_r<r_max_cm)], norm_r[nsat & (norm_r<r_max_cm)]
if rmsat:
vec_r, norm_r = vec_r[nsat & (norm_r<r_max_cm)], norm_r[nsat & (norm_r<r_max_cm)]
else:
vec_r, norm_r = vec_r[(norm_r<r_max_cm)], norm_r[(norm_r<r_max_cm)]
# Make arrays to split infow and outflows
if split_flows:
......@@ -1456,23 +1467,19 @@ def rad_col(genpath, folders, timesteps, labels, r_max=150, lmax=20, var_str='HI
file.create_dataset(strinflow,data=varnrmlz_in)
else:
# All
strallflow = f'all_{timestep:05d}_{proj}'
strallflow = f'{base_name}_{proj}'
try:
varnrmlz_all = file[strallflow][:]
except KeyError:
map2d,_ = cu.py_cell_utils.make_map_new(lmax,False,xmin,xmax,ymin,ymax,zmin,zmax,var,np.ones_like(var),xloc,yloc,zloc,cell_l,nx,ny)
map2d = map2d*(zmax-zmin)*cu2cm # Now it's a column density.
var_bin, _ = np.histogram(norm, range=[0,r_max], weights=map2d, bins=nbins) # Bin radially
nrmlz, _ = np.histogram(norm, weights=None, bins=nbins) # nrmlz
nrmlz, _ = np.histogram(norm, weights=None, bins=nbins) # nrmlz
varnrmlz_all = var_bin/nrmlz
all_varbin.append(varnrmlz_all)
temparr_xyz.append(varnrmlz_all)
file.create_dataset(strallflow,data=varnrmlz_all)
if proj=='x': # timestep std on a single proj.
array_x_tstep.append(varnrmlz_all)
# End proj loop
if not split_flows:
array_xyz.append(np.std(temparr_xyz, axis=0)) # standard deviation of the projs on ONE tstep
# End timestep loop
# Plot everything
......@@ -1483,25 +1490,11 @@ def rad_col(genpath, folders, timesteps, labels, r_max=150, lmax=20, var_str='HI
axes[0].plot(x_axis, med_varbin, label=labels[index], color=colors[index])
axes[0].fill_between(x_axis,np.percentile(all_varbin,perc_min, axis=0), np.percentile(all_varbin,perc_max, axis=0), color=colors[index], alpha=0.25)
# Impact of projection and time on deviation: which is the main responsible for variation ?
if not no_err:
err_x_tstep = np.std(array_x_tstep, axis=0) # median over timesteps for a single projection axis
err_xyz = np.std(array_xyz, axis=0) # median of the standard deviation due to projections
axes[1].plot(x_axis,err_xyz/err_x_tstep, colors[index]) # std proj/std tstep
axes[1].set_ylabel(fr'$\left( \sigma_{{x}}\ \left/\ \overline{{\sigma_{{x,y,z}}}}\ \right.\right)_{{{len_tstep}\ snaps}}$')
else:
med_varout = np.median(all_varout, axis=0)
med_varin = np.median(all_varin, axis=0)
axes[0].plot(x_axis, med_varout, label=labels[index], color=colors[index])
axes[0].plot(x_axis, med_varin, ls='--', color=colors[index])
if not no_err:
permin_out, permax_out = np.percentile(all_varout,perc_min, axis=0), np.percentile(all_varout,perc_max, axis=0)
permin_in, permax_in = np.percentile(all_varin,perc_min, axis=0), np.percentile(all_varin,perc_max, axis=0)
# Show ratio of percentiles in the second axis
axes[1].plot(x_axis,permax_out/permin_out, colors[index])
axes[2].plot(x_axis,permax_in/permin_in, colors[index], ls='--')
# axes[1].set_ylabel(fr'$\left( P_{{{perc_max}}}\ \left/\ P_{{{perc_min}}}\ \right.\right)_{{out}}$')
axes[1].set_ylabel(fr'$\left( \frac{{ P_{{{perc_max}}} }}{{ P_{{{perc_min}}} }} \right)_{{out}}$')
axes[2].set_ylabel(fr'$\left( \frac{{ P_{{{perc_max}}} }}{{ P_{{{perc_min}}} }} \right)_{{in}}$')
if plot:
plot_dexter(var_str, path2files='data/Ions/', ax=axes[0])
......
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