Commit a6d2859c authored by Maxime Rey's avatar Maxime Rey
Browse files

New function to plot theoretical SFR_ff.

parent 65dabd7c
# Just a file to store some basic functions.
# Just a file to store some basic or one time use functions.
# h5 related functions: h5_save, h5_read and h5_clean.
"""
Basic & one-time use functions.
==============================
- h5_save, h5_read, h5_clean, h5_list, h5rename
- compute_aout
- minimalmap
- SFR_ff
"""
###################################################################################################################################
######################################################## h5 base functions ########################################################
###################################################################################################################################
def h5_save(data,key,filename):
"""
Save data with name key in filename if it does not already exists.
......@@ -67,7 +75,7 @@ def h5_list(genpath,folder):
def h5rename(genpath, folders, timesteps, filename):
import h5py
for folder in folders:
RamsesDir = genpath+folders
RamsesDir = genpath+folder
filepath = f'{RamsesDir}/ratadat/{filename}.h5'
thefile = h5py.File(filepath,'r+')
......@@ -83,6 +91,9 @@ def h5rename(genpath, folders, timesteps, filename):
del thefile[oldname]
###################################################################################################################################
##################################################### one-time use functions #####################################################
###################################################################################################################################
def compute_aout(dt=10, z_init=15, nbins=10, zmin=None, checktimes=False):
""" Generate sequence of aexp for evenly-timed (dt in Myrs) ramses outputs in the ramses namelist format. JB & modifs from MR."""
from astropy import cosmology
......@@ -103,6 +114,7 @@ def compute_aout(dt=10, z_init=15, nbins=10, zmin=None, checktimes=False):
print((cosmo.age(1/a-1)).to(u.Myr))
return aout
def minimalmap():
""" Plot whatever quickly with this."""
from minirats.utils.py import cellutils as cu
......@@ -134,3 +146,72 @@ def minimalmap():
nx,ny = cu.py_cell_utils.get_map_nxny(lmax,0,1,0,1)
map2plot,_ = cu.py_cell_utils.make_map_new(lmax,True,0,1,0,1,0,1,cell_l,np.ones_like(cells[:,-1]),cell_pos[:,0],cell_pos[:,1],cell_pos[:,2],cell_l,nx,ny)
im = plt.imshow(map2plot.T,interpolation='nearest',origin='lower',norm=None, extent=(0,1,0,1))
def SFR_ff(what, thresh=1e-3, nblines=25):
from matplotlib import colors, cm, pyplot as plt
from ratatouille import plotutils as put
from scipy import special
import numpy as np
"""
Plots SFR_ff as a function of alpha_vir and the Mach number for Kim, Kretschmer and the ratio of both.
You might want to use "%matplotlib widget" if in a notebook to interact with the figure. Use "%matplotlib inline" to go back to normal.
"""
nbpoints = 100
avir_rg = np.linspace(0.01,100,nbpoints)
avir_rg = np.logspace(-2,2,nbpoints)
# Mach_rg = np.linspace(1,100,100)
Mach_rg = np.logspace(0,2,100)
avir, Mach = np.meshgrid(avir_rg, Mach_rg)
sig2 = np.log(1+0.4**2*Mach**2)
scritKi = np.log(0.067/0.33**2*avir*Mach**2)
scritKr = np.log(avir*(1+2*Mach**4/(1+Mach**2)))
SFRKi = 0.5/2/0.57*np.exp(3/8*sig2) * (1 + special.erf((sig2-scritKi)/np.sqrt(2*sig2)))
SFRKr = 1/2*np.exp(3/8*sig2) * (1 + special.erf((sig2-scritKr)/np.sqrt(2*sig2)))
fig = plt.figure(figsize=(10,10))
ax = plt.axes(projection ='3d')
if what=='Kimm':
SFR_ff = SFRKi.copy()
if thresh:
SFR_ff[SFRKi<thresh] = np.nan
colors = cm.inferno(norm(SFR_ff))
norm = colors.LogNorm(np.nanmin(SFR_ff), np.nanmax(SFR_ff))
ax.set_zlabel(r'$SFR_\mathrm{ff}$')
ax.set_xlim([-2,2])
ax.set_ylim([0,2])
ax.set_zlim([-3,2])
elif what=='Kretschmer':
SFR_ff = SFRKr.copy()
if thresh:
SFR_ff[SFRKi<thresh] = np.nan
norm = colors.LogNorm(np.nanmin(SFR_ff), np.nanmax(SFR_ff))
colors = cm.inferno(norm(SFR_ff))
ax.set_zlabel(r'$SFR_\mathrm{ff}$')
ax.set_xlim([-2,2])
ax.set_ylim([0,2])
ax.set_zlim([-3,2])
elif what=='ratio':
SFR_ff = SFRKr/SFRKi
if thresh:
SFR_ff[(SFRKi>thresh)&(SFRKr>thresh)] = np.nan
norm = colors.Normalize(np.nanmin(SFR_ff), np.nanmax(SFR_ff))
colors = cm.inferno(norm(SFR_ff))
ax.set_zlabel(r'$SFR_\mathrm{ff, Kr} / SFR_\mathrm{ff, Ki}$')
if nblines:
rcount, ccount = nblines, nblines
else:
rcount, ccount, _ = colors.shape
map = ax.plot_surface(np.log10(avir), np.log10(Mach), np.log10(SFR_ff), rcount=rcount, ccount=ccount, facecolors=colors, shade=False)
map.set_facecolor((0,0,0,0))
ax.xaxis.set_major_formatter(plt.FuncFormatter(put.axis_format_10pow))
ax.yaxis.set_major_formatter(plt.FuncFormatter(put.axis_format_10pow))
ax.zaxis.set_major_formatter(plt.FuncFormatter(put.axis_format_10pow))
ax.set_ylabel(r'$\mathcal{M}$')
ax.set_xlabel(r'$\alpha_\mathrm{vir}$')
plt.show()
\ No newline at end of file
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