Skip to content
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
93 changes: 60 additions & 33 deletions kippenhahn/kippenhahn.py
Original file line number Diff line number Diff line change
Expand Up @@ -135,10 +135,10 @@ def InterpolateOneProfile(profile, NY, Yaxis, Ymin, Ymax, Variable):
Variable (str) -- the variable to be interpolated
Possible values: eps_nuc, velocity, entropy, total_energy, j_rot,
eps_recombination, ionization_energy, energy, potential_plus_kinetic,
extra_heat, v_div_vesc, v_div_csound, pressure, temperature, density,
extra_heat, v_div_v_escape, v_div_csound, pressure, temperature, density,
tau, opacity, gamma1, gamma3, dq, L_div_Ledd, Lrad_div_Ledd. t_thermal, t_dynamical,
t_dynamical_down, t_thermal_div_t_dynamical, omega_div_omega_crit, omega
super_ad, vconv, vconv_div_vesc, conv_vel_div_csound, total_energy_plus_vconv2
super_ad, vconv, vconv_div_v_escape, conv_vel_div_csound, total_energy_plus_vconv2
t_thermal_div_t_expansion, E_kinetic_div_E_thermal


Expand Down Expand Up @@ -205,12 +205,12 @@ def InterpolateOneProfile(profile, NY, Yaxis, Ymin, Ymax, Variable):
except Exception:
raise ValueError("Column 'total_energy' and/or 'energy' is missing from the profile files")

if (not "v_div_vesc" in profile.dtype.names and (Variable == 'v_div_vesc' )):
if (not "v_div_v_escape" in profile.dtype.names and (Variable == 'v_div_v_escape' )):
G = const.G.to('cm3/(g s2)').value
Msun = const.M_sun.to('g').value
Rsun = const.R_sun.to('cm').value
try:
profile = numpy.lib.recfunctions.append_fields(profile,'v_div_vesc',
profile = numpy.lib.recfunctions.append_fields(profile,'v_div_v_escape',
data = profile['velocity']/np.sqrt(2.*G*(profile['mass']*Msun)/(profile['radius']*Rsun)), asrecarray=True)
except Exception:
raise ValueError("Column 'radius' and/or 'velocity' is missing from the profile files")
Expand Down Expand Up @@ -325,12 +325,12 @@ def InterpolateOneProfile(profile, NY, Yaxis, Ymin, Ymax, Variable):
except Exception:
raise ValueError("Column 'conv_vel' is missing from the profile files")

if (not "vconv_div_vesc" in profile.dtype.names and (Variable == 'vconv_div_vesc' )):
if (not "vconv_div_v_escape" in profile.dtype.names and (Variable == 'vconv_div_v_escape' )):
G = const.G.to('cm3/(g s2)').value
Msun = const.M_sun.to('g').value
Rsun = const.R_sun.to('cm').value
try:
profile = numpy.lib.recfunctions.append_fields(profile,'vconv_div_vesc',
profile = numpy.lib.recfunctions.append_fields(profile,'vconv_div_v_escape',
data = profile['conv_vel']/np.sqrt(2.*G*(profile['mass']*Msun)/(profile['radius']*Rsun)), asrecarray=True)
except Exception:
raise ValueError("Column 'radius' and/or 'conv_vel' is missing from the profile files")
Expand Down Expand Up @@ -431,6 +431,11 @@ def __init__(self, **kwargs):
masses_TML (list[float]) -- mass locations (in Msun) to trace in radius (default [])
xvals_TML (list[float]) -- xvals, between 0 and 1, where labels of lines that trace constant mass are places.
Must be same size as masses_TML. If left empty, labels are placed automatically (default [])
fig -- Figure object, if not given by the user the object will be created inside the routine (default None)
ax -- Axes object, if not given by the user the object will be created inside the routine (default None)
cmap_label -- Label for the third quantity ploted as the color map.
This must be given in case the user needs to plot a quantity that is not pre-defined in the ploting code. (default None)
log_cmap -- Define if the quantity ploted using the colormap will be shown in log10 scale or not (default True)
"""

self._param = {'data_path':"./", 'NX':1024, 'NY':1024, 'Yaxis':'mass', 'Xaxis':'star_age',
Expand All @@ -439,7 +444,8 @@ def __init__(self, **kwargs):
'onscreen':False, 'parallel':True, 'abundances':False, 'log_abundances':True, 'czones':False,
'signed_log_cmap':True, 'orbit':False, 'tau10':False, 'tau100':False, 'Nprofiles_to_plot':10,
'profiles_to_plot':[], 'masses_TML':[], 'xvals_TML':[], 'Xaxis_min':None, 'Xaxis_max':None, 'Xaxis_label':None,
'Yaxis_label':None, 'Yaxis_min':None, 'Yaxis_max':None, 'cmap_min':None, 'cmap_max':None, 'cmap_label':None}
'Yaxis_label':None, 'Yaxis_min':None, 'Yaxis_max':None, 'cmap_min':None, 'cmap_max':None, 'cmap_label':None,
'ax': None, 'fig' : None, 'log_cmap' : True}

for key in kwargs:
if (key in self._param):
Expand Down Expand Up @@ -596,6 +602,19 @@ def cmap_max(self):
@property
def cmap_label(self):
return self._param['cmap_label']

@property
def ax(self):
return self._param['ax']

@property
def fig(self):
return self._param['fig']

@property
def log_cmap(self):
return self._param['log_cmap']


def help(self):
#TODO: add a list of all parameters, the default values and the possible option, add a list of functions, and an example
Expand All @@ -614,7 +633,7 @@ def CheckParameters(self):
Yaxis -- mass, radius, q, log_mass, log_radius, log_q
Variable -- eps_nuc, velocity, entropy, total_energy, j_rot,
eps_recombination, ionization_energy, energy, potential_plus_kinetic,
extra_heat, v_div_vesc, v_div_csound, pressure, temperature, density,
extra_heat, v_div_v_escape, v_div_csound, pressure, temperature, density,
tau, opacity, gamma1, dq,L_div_Ledd, Lrad_div_Ledd, t_thermal, t_dynamical,
t_dynamical_down, t_thermal_div_t_dynamical, omega_div_omega_crit, omega

Expand All @@ -628,14 +647,17 @@ def CheckParameters(self):
if not (self._param['Xaxis'] in ['model_number', 'star_age', 'inv_star_age', 'log_model_number', 'log_star_age',
'log_inv_star_age']):
raise ValueError(self._param['Xaxis']+"not a valid option for parameter Xaxis")
if not (self._param['Variable'] in ['eps_nuc', 'velocity', 'entropy', 'total_energy', 'j_rot', 'eps_recombination'
, 'ionization_energy', 'energy', 'potential_plus_kinetic', 'extra_heat', 'v_div_vesc',
'v_div_csound', 'pressure', 'temperature', 'density', 'tau', 'opacity', 'gamma1', 'gamma3', 'dq',
'L_div_Ledd', 'Lrad_div_Ledd', 't_thermal', 't_dynamical', 't_dynamical_down', 't_thermal_div_t_dynamical',
'omega_div_omega_crit', 'omega','super_ad', 'vconv', 'vconv_div_vesc', 'conv_vel_div_csound',
't_thermal_div_t_expansion','E_kinetic_div_E_thermal','total_energy_plus_vconv2']):
raise ValueError(self._param['Variable']+"not a valid option for parameter Variable")


"""
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do you start a second docstring here? Btw. never use docstrings to multi line comments, this is very dirty coding!

if (self._param["cmap_label"] is None):
if (not (self._param['Variable'] in ['eps_nuc', 'velocity', 'entropy', 'total_energy', 'j_rot', 'eps_recombination'
, 'ionization_energy', 'energy', 'potential_plus_kinetic', 'extra_heat', 'v_div_v_escape',
'v_div_csound', 'pressure', 'temperature', 'density', 'tau', 'opacity', 'gamma1', 'gamma3', 'dq',
'L_div_Ledd', 'Lrad_div_Ledd', 't_thermal', 't_dynamical', 't_dynamical_down', 't_thermal_div_t_dynamical',
'omega_div_omega_crit', 'omega','super_ad', 'vconv', 'vconv_div_v_escape', 'conv_vel_div_csound',
't_thermal_div_t_expansion','E_kinetic_div_E_thermal','total_energy_plus_vconv2'])):
raise ValueError(self._param['Variable']+" not a valid option for parameter Variable, you must deffine cmap_label.")
"""

return

Expand Down Expand Up @@ -912,7 +934,7 @@ def plot(self, ax=None):
cmap_label = "log(specific potential+kinetic energy [erg/gr])"
elif self._param['Variable'] == "extra_heat":
cmap_label = "log(specific injected energy rate [erg/s/gr])"
elif self._param['Variable'] == "v_div_vesc":
elif self._param['Variable'] == "v_div_v_escape":
cmap_label = "log($v/v_{esc}$)"
elif self._param['Variable'] == "v_div_csound":
cmap_label = "log($v/v_{sound}$)"
Expand Down Expand Up @@ -952,7 +974,7 @@ def plot(self, ax=None):
cmap_label = "log($\\nabla - \\nabla_{ad}$)"
elif self._param['Variable'] == "vconv":
cmap_label = "log($v_{conv}$ [cm/s])"
elif self._param['Variable'] == "vconv_div_vesc":
elif self._param['Variable'] == "vconv_div_v_escape":
cmap_label = "log($v_{conv}/v_{esc}$)"
elif self._param['Variable'] == "conv_vel_div_csound":
cmap_label = "log($v_{conv}/v_{sound}$)"
Expand All @@ -968,16 +990,17 @@ def plot(self, ax=None):



if self._param['signed_log_cmap'] and (self._cmap_label is not None):
if self._param['signed_log_cmap'] and (cmap_label is not None):
cmap_label = "sign x log(max(1,abs(" + cmap_label[4:]+"))"


if ax is None:
if self._param['ax'] is None:
fig1 = plt.figure()
ax1 = fig1.add_subplot(111)
fig1.subplots_adjust(top=0.80, left=0.12, right=0.9, bottom=0.12)
else:
ax1 = ax
ax1 = self._param['ax']
fig1 = self._param['fig']
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

After this change, shouldn't we remove the argument ax from the plot function? If it is still used somewhere else you should replace it there with the new axis object stored in the kippenhahn object.



ax1.set_xlabel(Xlabel,fontsize=self._param['font_large'])
Expand All @@ -991,7 +1014,10 @@ def plot(self, ax=None):
if self._param['signed_log_cmap']:
data_to_plot = np.sign(np.transpose(self._data)) * np.log10(np.maximum(1.,np.abs(np.transpose(self._data))))
else:
data_to_plot = np.log10(np.transpose(self._data))
if self._param['log_cmap']:
data_to_plot = np.log10(np.abs(np.transpose(self._data)))
else:
data_to_plot = np.transpose(self._data)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just a matter of style, should you better use the structure:

if self._param['signed_log_cmap']:
    ...
elif self._param['log_cmap']:
    ...
else:
    ...

instead of

if self._param['signed_log_cmap']:
    ...
else:
    if self._param['log_cmap']:
        ...
    else:
        ...


#Define the minimum and maximum of the color scale
# When using signed_log_cmap, ignore cmap_dynamic_range
Expand Down Expand Up @@ -1345,14 +1371,15 @@ def plot(self, ax=None):


#Re-enforcing the calculated limits for X and Y axis. Without this there may be a white band on the right of the plot.
ax1.set_xlim([self._Xaxis_min,self._Xaxis_max])
ax1.set_ylim([self._Yaxis_min,self._Yaxis_max])
if self._param['abundances']:
ax2.set_xlim([self._Xaxis_min,self._Xaxis_max])
if self._param['log_abundances']:
ax2.set_ylim([1e-5,1.])
else:
ax2.set_ylim([0.,1.])
if ax is not None:
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you sure you don't want to check for self._param['ax'] here?
Even a deeper thought: shouldn't we check, whether self._Xaxis_min, self._Xaxis_max, ... are given? Usually, I'd let the user overwrite old limits, when new ones are given.

ax1.set_xlim([self._Xaxis_min,self._Xaxis_max])
ax1.set_ylim([self._Yaxis_min,self._Yaxis_max])
if self._param['abundances']:
ax2.set_xlim([self._Xaxis_min,self._Xaxis_max])
if self._param['log_abundances']:
ax2.set_ylim([1e-5,1.])
else:
ax2.set_ylim([0.,1.])


# fig1.tight_layout()
Expand Down Expand Up @@ -1409,16 +1436,16 @@ def Demo(self):
#Options for Xaxis: 'model_number', 'star_age', 'inv_star_age', 'log_model_number', 'log_star_age', 'log_inv_star_age'
#Options for Yaxis: 'mass', 'radius', 'q', 'log_mass', 'log_radius', 'log_q'
#Options for Variable: "eps_nuc", "velocity", "entropy", "total_energy", "j_rot", "eps_recombination", "ionization_energy",
# "energy", "potential_plus_kinetic", "extra_heat", "v_div_vesc", "v_div_csound"
# "energy", "potential_plus_kinetic", "extra_heat", "v_div_v_escape", "v_div_csound"
# "pressure", "temperature", "density", "tau", "opacity", "gamma1", "dq"
# "super_ad", "vconv", "vconv_div_vesc", "conv_vel_div_csound", "total_energy_plus_vconv2"
# "super_ad", "vconv", "vconv_div_v_escape", "conv_vel_div_csound", "total_energy_plus_vconv2"
# "omega_div_omega_crit", "omega"



data_path = "/data/disk1/fragkos/repos/CE_mesa/working/LOGS/"
a = kippenhahn(data_path=data_path, parallel=False, abundances=False, log_abundances = True, Yaxis='radius', Xaxis="star_age",
czones=True, Variable='v_div_vesc', orbit=True, masses_TML = [4., 6., 8., 10.,12.,14., 20., 25., 28.])
czones=True, Variable='v_div_v_escape', orbit=True, masses_TML = [4., 6., 8., 10.,12.,14., 20., 25., 28.])
a.SetParameters(onscreen=True, cmap = 'jet', cmap_dynamic_range=5, signed_log_cmap=False, tau10=True, tau100=True)


Expand Down