"""
.. module:: plot_mc
:synopsis: Plotting routines for molecular cloud |prodimo| models (0D chemistry).
.. moduleauthor:: Ch. Rab
"""
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import prodimopy.utils as putils
from prodimopy.read import Chemistry
from prodimopy.plot import reacToStr, scale_figs, spnToLatex
# That crashes in the bash on Mac OS X.. maybe because of installation
[docs]
class PlotMcModel(object):
"""
Plot routines for a single molecular cloud |prodimo| model (0D chemistry).
.. todo::
- Redesign this. it is not very useful to have this dokwargs legend etc. routines copied all the time.
"""
def __init__(
self,
pdf: mpl.backends.backend_pdf.PdfPages | None = None,
colors=None,
styles=None,
markers=None,
fs_legend=6, # TODO: init it with the system default legend fontsize
ncol_legend=0,
):
"""
Parameters
----------
pdf :
a PdfPages object used to save the plot or None to have output on screen.
colors : list
a list of matplotlib colors used for different models. (default: None)
styles : list
a list of matplotlib styles used for different models. (default: None)
markers : list
a list of matplotlib markers used for different models. (default: None)
Attributes
----------
"""
self.pdf = pdf
if colors is None:
self.colors = [
"b",
"g",
"r",
"c",
"m",
"y",
"k",
"sienna",
"lime",
"pink",
"DarkOrange",
"Olive",
]
else:
self.colors = colors
if styles is None:
self.styles = ["-"] * len(self.colors)
else:
self.styles = styles
if markers is None:
self.markers = [None] * len(self.colors)
else:
self.markers = markers
self.fs_legend = fs_legend
self.ncol_legend = ncol_legend
def _dokwargs(self, ax, **kwargs):
"""
Handles the passed kwargs elements (assumes that defaults are already set)
"""
if "ylim" in kwargs:
ax.set_ylim(kwargs["ylim"])
if "xlim" in kwargs:
ax.set_xlim(kwargs["xlim"])
if "xlog" in kwargs:
if kwargs["xlog"]:
ax.semilogx()
if "ylog" in kwargs:
if kwargs["ylog"]:
ax.semilogy()
if "grid" in kwargs: # add a grid
ax.grid(kwargs["grid"])
def _legend(self, ax):
"""
plots the legend, deals with multiple columns
"""
handles, labels = ax.get_legend_handles_labels()
ncol = 1
if self.ncol_legend > 1 and len(labels) > self.ncol_legend:
ncol = int(len(labels) / self.ncol_legend)
ax.legend(
handles,
labels,
loc="best",
fancybox=False,
framealpha=0.75,
ncol=ncol,
fontsize=self.fs_legend,
)
def _closefig(self, fig):
"""
Save and close the plot (Figure).
If `self.pdf` is None than nothing is done and the figure is returned
#set the transparent attribute (used rcParam savefig.transparent)
"""
# trans=mpl.rcParams['savefig.transparent']
if self.pdf is not None:
self.pdf.savefig(figure=fig, transparent=False)
plt.close(fig)
return None
else:
return fig
def _initfig(
self, ax: mpl.axes.Axes = None, **kwargs
) -> tuple[mpl.figure.Figure, mpl.axes.Axes]:
"""
Inits Figure and Axes object.
If an Axes object is passed, it is returned together with the Figure object.
This is for a single plot (i.e. only one panel).
Returns
-------
The figure and axes object.
"""
if ax is not None:
fig = ax.get_figure()
else:
fig, ax = plt.subplots(1, 1, figsize=self._sfigs(**kwargs))
return fig, ax
def _sfigs(self, **kwargs):
"""
Scale the figure size from matplotlibrc by the factors given in the
array sfigs (in kwargs) the first element is for the width the second for
the height
"""
if "sfigs" in kwargs:
fac = kwargs["sfigs"]
return scale_figs(fac)
else:
return scale_figs([1.0, 1.0])
[docs]
def plot_species(self, model, spnames, xunit="years", ax=None, **kwargs):
"""
Plot the species abundances as function of age (time).
One can also pass further arguments such as `xlim (Array)`, `ylim (Array)`, `xlog`, `ylog`, `grid`.
Parameters
----------
model : :class:`prodimopy.read_mc.Data_mc`
a :class:`prodimopy.read_mc.Data_mc` object used to save the plot.
spnames : list or str
The names of the species that should be plotted. Can also be a single
string.
ax : :class:`matplotlib.axes.Axes`
The axes object to plot the data on. If None a new figure is created.
xunit : str
The unit for the x-axis (age). Default is "years". Options are `"years"` and `"days"`
.. todo::
the x axis (ages) is always plotted in log-scale so the zero age (
initial abundance) is not shown. Needs more flexibility.
"""
if isinstance(spnames, str):
spnames = [spnames]
if ax is not None:
fig = ax.get_figure()
else:
fig, ax = plt.subplots(1, 1)
i = 0
xmax = 1.0e-99
xmin = 1.0e99
for spname in spnames:
if spname in model.species:
if xunit == "days":
x = model.ages * 365.25
else:
x = model.ages
ax.plot(
x,
model.abundances[:, model.species.index(spname)],
self.styles[i],
marker=self.markers[i],
color=self.colors[i],
label=spname,
)
xmax = max([max(x), xmax])
# exclude zero because of log plot
# nozeros=
xmin = min([min(x[x > 0.0]), xmin])
i = i + 1
else:
print("Species " + spname + " not found.")
# no data to plot just return
if i == 0:
return
ax.set_xlim(xmin, xmax)
if kwargs.get("xlog", True):
ax.semilogx()
if kwargs.get("ylog", True):
ax.semilogy()
self._dokwargs(ax, **kwargs)
ax.set_xlabel(xunit)
ax.set_ylabel(r" species abundance")
self._legend(ax)
if "title" in kwargs and kwargs["title"] is not None:
ax.set_title(kwargs["title"])
return self._closefig(fig)
[docs]
def plot_reac(
self,
rtype: str,
chemanas: list[Chemistry],
chemana_steadystate: Chemistry | None = None,
ax=None,
**kwargs,
):
"""
Plots the rates for the most important reactions at the point ``ix``, ``iz`` for the
given chemanalysis list, which usually come from a time-dependent |prodimo| disk model (i.e time steps).
That ages (x-axis) are taken from the model connected the the chemanalysis objects.
Parameters
----------
rtype :
``f`` for formation reactions, ``d`` for destruction reactions.
chemanas :
A list of chemanalysis objects (e.g. for each
age in a time-dependent model). All have to be for the same species.
chemana_steadystate :
The chemistry analysis of an equivalent steady-state model (i.e. same grid, etc.) for the
same species. Optional
"""
ix = 0
iz = 0
# Fixme: Would make sense to have this as a parameter
nmaxreac = 3
if rtype == "d":
tit = "dest. reactions for "
else:
tit = "form. reactions for "
totrates = list()
reacidx = list()
# FIXME: attaching a model to chemana in mc mode is not so usefull
ages = chemanas[0].model.ages
for chemana in chemanas:
if rtype == "d":
totrates.append(chemana.totdrate[ix, iz])
reacs = chemana.gridd[ix, iz, 0]
else:
totrates.append(chemana.totfrate[ix, iz])
reacs = chemana.gridf[ix, iz, 0]
# figure out what reactions indices I want to plot
for i in range(nmaxreac):
if i < len(reacs):
reacidx.append(reacs[i])
reacidx = np.unique(np.array(reacidx))
reactions = list()
for ridx in reacidx:
rates = list()
for chemana in chemanas:
if rtype == "d":
rate = chemana.gridd[ix, iz, 1][chemana.gridd[ix, iz, 0] == ridx]
else:
rate = chemana.gridf[ix, iz, 1][chemana.gridf[ix, iz, 0] == ridx]
if len(rate) == 0:
rates.append(0.0)
else:
rates.append(rate[0])
reactions.append(rates)
fig, ax = self._initfig(ax, **kwargs)
ax.plot(ages, totrates, label="totrate", color="black", linewidth="3")
for i in range(len(reacidx)):
# This assume that all chemanas have the same chemnet in the background, otherwise it would not make much sense
ax.plot(
ages,
reactions[i],
label=reacToStr(chemanas[0].model.chemnet.reactions[reacidx[i] - 1]),
)
# for the steady state model
if chemana_steadystate is not None:
reacidxsss = list()
ratesss = list()
for i in range(nmaxreac):
if rtype == "d":
if i < len(chemana_steadystate.gridd[ix, iz, 0]):
reacidxsss.append(chemana_steadystate.gridd[ix, iz, 0][i])
ratesss.append(chemana_steadystate.gridd[ix, iz, 1][i])
else:
if i < len(chemana_steadystate.gridf[ix, iz, 0]):
reacidxsss.append(chemana_steadystate.gridf[ix, iz, 0][i])
ratesss.append(chemana_steadystate.gridf[ix, iz, 1][i])
colors = np.arange(0.1, 0.9, 0.7 / nmaxreac)
for i, (reacidx, rate) in enumerate(zip(reacidxsss, ratesss)):
ax.scatter(
ages[-1] * 1.05,
rate,
label=reacToStr(chemana_steadystate.model.chemnet.reactions[reacidx - 1]),
marker="<",
color=str(colors[i]),
s=20 / np.log(i + 2),
)
ax.legend(
bbox_to_anchor=(1.01, 1.0), loc="upper left", prop={"family": "monospace", "size": 5}
)
if kwargs.get("xlog", True):
ax.semilogx()
if kwargs.get("ylog", True):
ax.semilogy()
ax.set_xlabel("age [yr]")
ax.set_ylabel(r"rate $[cm^{-3}\,s^{-1}]$")
ax.set_title(tit + " ix={:5d}, iz={:5d}".format(ix, iz))
self._dokwargs(ax, **kwargs)
return self._closefig(fig)
[docs]
class PlotMcModels(object):
"""
Plot routines for a molecular cloud |prodimo| models (0D chemistry).
"""
def __init__(
self,
pdf: mpl.backends.backend_pdf.PdfPages | None = None,
colors=None,
styles=None,
markers=None,
fs_legend=6, # TODO: init it with the system default legend fontsize
ncol_legend=0,
):
"""
Parameters
----------
name :
a `PdfPages` object used to save the plot.
colors : list
a list of matplotlib colors used for different models. (default: None)
styles : list
a list of matplotlib styles used for different models. (default: None)
markers : list
a list of matplotlib markers used for different models. (default: None)
Attributes
----------
"""
self.pdf = pdf
if colors is None:
self.colors = [
"b",
"g",
"r",
"c",
"m",
"y",
"k",
"sienna",
"lime",
"pink",
"DarkOrange",
"Olive",
]
else:
self.colors = colors
if styles is None:
self.styles = ["-"] * len(self.colors)
else:
self.styles = styles
if markers is None:
self.markers = [None] * len(self.colors)
else:
self.markers = markers
self.fs_legend = fs_legend
self.ncol_legend = ncol_legend
def _dokwargs(self, ax, **kwargs):
"""
Handles the passed kwargs elements (assumes that defaults are already set)
"""
if "ylim" in kwargs:
ax.set_ylim(kwargs["ylim"])
if "xlim" in kwargs:
ax.set_xlim(kwargs["xlim"])
if "xlog" in kwargs:
if kwargs["xlog"]:
ax.semilogx()
if "ylog" in kwargs:
if kwargs["ylog"]:
ax.semilogy()
if "grid" in kwargs: # add a grid
ax.grid(kwargs["grid"])
def _legend(self, ax):
"""
plots the legend, deals with multiple columns
"""
handles, labels = ax.get_legend_handles_labels()
ncol = 1
if self.ncol_legend > 1 and len(labels) > self.ncol_legend:
ncol = int(len(labels) / self.ncol_legend)
ax.legend(
handles,
labels,
loc="best",
fancybox=False,
framealpha=0.75,
ncol=ncol,
fontsize=self.fs_legend,
)
def _closefig(self, fig):
"""
Save and close the plot (Figure).
If self.pdf is None than nothing is done and the figure is returned
#set the transparent attribut (used rcParam savefig.transparent)
"""
# trans=mpl.rcParams['savefig.transparent']
if self.pdf is not None:
self.pdf.savefig(figure=fig, transparent=False)
plt.close(fig)
return None
else:
return fig
[docs]
def plot_species(self, models, spname, ice=False, ax=None, **kwargs):
"""
Plot the given species (spname) for all the given models.
One can also pass further arguments such as `xlim (Array)`, `ylim (Array)`, `xlog`, `ylog`, `grid`.
Parameters
----------
models : array_like(:class:`prodimopy.read_mc.Data_mc`,ndim=1)
list of molecular cloud models
spname : str
The name fo the species to plot
"""
if ax is not None:
fig = ax.get_figure()
else:
fig, ax = plt.subplots(1, 1)
i = 0
xmax = 1.0e-99
xmin = 1.0e99
ymax = 1.0e-99
ymin = 1.0e99
for model in models:
if spname in model.species:
abun = model.abundances[:, model.species.index(spname)]
(line,) = ax.plot(
model.ages,
abun,
self.styles[i],
marker=self.markers[i],
color=self.colors[i],
label=model.name,
)
ymax = max([max(abun), ymax])
ymin = min([min(abun[model.ages > 0.0]), ymin])
if ice and (spname + "#" in model.species):
abun = model.abundances[:, model.species.index(spname + "#")]
(line,) = ax.plot(
model.ages, abun, "--", marker=self.markers[i], color=line.get_color()
)
ymax = max([max(abun), ymax])
ymin = min([min(abun[model.ages > 0.0]), ymin])
xmax = max([max(model.ages), xmax])
# exclude zero because of log plot
# nozeros=
xmin = min([min(model.ages[model.ages > 0.0]), xmin])
i = i + 1
# no data to plot just return
if i == 0:
print("Species " + spname + " not found.")
return
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin * 0.9, ymax * 1.1)
ax.semilogx()
ax.semilogy()
self._dokwargs(ax, **kwargs)
ax.set_xlabel(r"years")
ylabel = r" $\mathrm{\epsilon(" + spnToLatex(spname) + ")}$"
if ice and (spname + "#" in model.species):
ylabel += r", $\mathrm{\epsilon(" + spnToLatex(spname + "#") + ") dashed}$"
ax.set_ylabel(ylabel)
self._legend(ax)
if "title" in kwargs and kwargs["title"] is not None:
ax.set_title(kwargs["title"])
return self._closefig(fig)
[docs]
def plot_species_diff(self, models, spname, ax=None, **kwargs):
"""
Plot the difference of the species abundance relative to the last model.
The relative diffferencie is definded as abs(modelval/modelrefval-1.0).
One can also pass further arguments such as `xlim (Array)`, `ylim (Array)`, `xlog`, `ylog`, `grid`.
Parameters
----------
models : array_like(:class:`prodimopy.read_mc.Data_mc`,ndim=1)
list of molecular cloud models
spname : str
The name fo the species to plot
"""
if ax is not None:
fig = ax.get_figure()
else:
fig, ax = plt.subplots(1, 1)
i = 0
xmax = 1.0e-99
xmin = 1.0e99
ymax = 1.0e-99
ymin = 1.0e99
for model in models[:-1]:
if spname in model.species:
# FIXME: add small number to avoid problems with zero
diff = np.absolute(
(model.abundances[1:, model.species.index(spname)] + 1.0e-200)
/ (models[-1].abundances[1:, models[-1].species.index(spname)] + 1.0e-200)
- 1.0
)
# print(diff)
ax.plot(
model.ages[1:],
diff,
self.styles[i],
marker=self.markers[i],
color=self.colors[i],
label=model.name,
)
xmax = max([max(model.ages), xmax])
xmin = min([min(model.ages[1:]), xmin])
ymax = max([max(diff), ymax])
ymin = min([min(diff), ymin])
i = i + 1
# no data to plot just return
if i == 0:
print("Species " + spname + " not found.")
return
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin * 0.9, ymax * 1.1)
ax.semilogx()
ax.semilogy()
self._dokwargs(ax, **kwargs)
ax.set_xlabel(r"years")
ax.set_ylabel(r" rel. diff $\mathrm{\epsilon(" + spnToLatex(spname) + ")}$")
self._legend(ax)
# if "title" in kwargs and kwargs["title"] != None:
# ax.set_title(kwargs["title"])
return self._closefig(fig)
[docs]
def plot_abunratio(self, models, spname1, spname2, **kwargs):
"""
Plots the ratio of species spname1 over species spname2 for all the passed models.
One can also pass further arguments such as `xlim (Array)`, `ylim (Array)`, `xlog`, `ylog`, `grid`.
Parameters
----------
models : array_like(:class:`prodimopy.read_mc.Data_mc`,ndim=1)
list of molecular cloud models
spname1 : str
The name of the first species used for the ratio.
spname2 : str
The name of the second species used for the ratio (sp1/sp2)
"""
fig, ax = plt.subplots(1, 1)
i = 0
xmax = 1.0e-99
xmin = 1.0e99
for model in models:
if spname1 in model.species and spname2 in model.species:
ratio = (
model.abundances[:, model.species.index(spname1)]
/ model.abundances[:, model.species.index(spname2)]
)
ax.plot(
model.ages,
ratio,
self.styles[i],
marker=self.markers[i],
color=self.colors[i],
label=model.name,
)
xmax = max([max(model.ages), xmax])
# exclude zero because of log plot
# nozeros=
xmin = min([min(model.ages[model.ages > 0.0]), xmin])
i = i + 1
if i == 0:
print("Species " + spname1 + "and/or " + spname2 + " not found.")
return
ax.set_xlim(xmin, xmax)
ax.semilogx()
ax.semilogy()
self._dokwargs(ax, **kwargs)
ax.set_xlabel(r"years")
ax.set_ylabel(
r" $\mathrm{\epsilon("
+ spnToLatex(spname1)
+ ")/"
+ r"\epsilon("
+ spnToLatex(spname2)
+ ")}$"
)
self._legend(ax)
# if "title" in kwargs and kwargs["title"] != None:
# ax.set_title(kwargs["title"])
return self._closefig(fig)
[docs]
def load_style(style="prodimopy"):
"""
Simple wrapper that calls :func:`~prodimopy.utils.load_mplstyle` with the given style.
Parameters
----------
style : str
The name of the style to load. Default: "prodimopy"
"""
putils.load_mplstyle(style)