Other scripts example

[Albertringwood]

Dear Dr. Dao Thank you so much for making the complete guide available, it will be very useful! As you are a great expert on Python and the program, I ask you if it is possible to insert sample scripts where you can view other quantities of the simulation products instead of the material id, especially entropy and temperature (or internal energy, as your profile picture shows). Is it possible to have scripts that can see these quantities and save the snapshots of the simulation in sequence? In this way it would be possible to see the giant impacts also according to other physical quantities. Thank you very much for your availability! Albert

[JingyaoDOU]

Hi Albert,

Perhaps the below code would help:

snapdir = "/usr/foo/"
for i in range(100):
    snapshotfile = snapdir+"snapshot_%04d.hdf5"%i
    data    = sw.load(snapshotfile)
    box_mid = 0.5 * data.metadata.boxsize[0] / R_earth
    pos = data.gas.coordinates - box_mid
    data.gas.densities.convert_to_mks()
    rho_mks = np.array(data.gas.densities)
    data.gas.internal_energies.convert_to_mks()
    u_mks   = np.array(data.gas.internal_energies) # internal energy in mks
    matid   = np.array(data.gas.material_ids)

    entropy = woma.eos.eos.A1_s_u_rho(u_mks, rho_mks, matid)
    T       = woma.eos.eos.A1_T_u_rho(u_mks, rho_mks, matid)

    fig, ax = plt.subplots(1,1,figsize=(6,6))
    ax.scatter(
    pos[:, 0],
    pos[:, 1],
    s=3.5 * size,
    c=entropy,
    edgecolors="none",
    marker=".",
    alpha=1,
    cmap='viridis',

)

    saveloc='/user/plotsave/'
    save = saveloc + "snapPLOT_%0d4.png"%i 
    fig.savefig(save,bbox_inches='tight',facecolor='w')

    plt.show()
    plt.cla()
    fig.clf() 
    plt.close()

[Albertringwood]

Dear Dr. Dao I managed to get the script working by editing it, in the figure is shown what I managed to get (temperature).

The problem is that, when I choose the limit of the axes, the produced figure becomes totally white. I don’t know why and I couldn’t solve the problem. Please, could you take a look at the script? Also, how can I add a label and a temperature range to display? Product snapshots always change color during collision because I think they always try to adapt. Thank you very much Albert script.txt

[JingyaoDOU]

Hi Albert,

The reason why it's white or empty might be because the center of mass of the system is draft and you are perhaps loading the area where there are actually no particles there. Instead of loading by the mask function, you could load all the particles all at once and relocate the center of mass of the system. I'm sorry my previous codes are misleading and use the mask function to load particles. I hope the below code will help. Please be aware that I'm setting internal unit of length to Earth Radius and mass to Earth mass when running SWIFTsim. The codes below add a colorbar and set the Tmin and Tmax to fix the Temperature range that the colormap covers. So you should first try to find what will be the max temperature of your simulation by tracking the T at the most violent snapshot during the first 1 or 2 hours of the impact. Then set the "Tmax" to this temperature.

loc='/foo/snapOUT_cooling_planet_npt51594_0d99893Earth.hdf5'

data    = sw.load(loc)
box_mid = 0.5 * data.metadata.boxsize[0]
pos = data.gas.coordinates - box_mid
data.gas.densities.convert_to_mks()
rho_mks = np.array(data.gas.densities)
data.gas.internal_energies.convert_to_mks()
u_mks   = np.array(data.gas.internal_energies) # internal energy in mks
matid   = np.array(data.gas.material_ids)
m       = np.array(data.gas.masses)

entropy = woma.eos.eos.A1_s_u_rho(u_mks, rho_mks, matid)
T       = woma.eos.eos.A1_T_u_rho(u_mks, rho_mks, matid)

pos_centerM = np.sum(pos * m[:, np.newaxis], axis=0) / np.sum(m)
pos -= pos_centerM

Tmin = 1000
Tmax = 6000

fig, ax = plt.subplots(1,1,figsize=(6,6))
cimg = ax.scatter(
pos[:, 0],
pos[:, 1],
s=3.5,
c=T,
edgecolors="none",
marker=".",
alpha=1,
cmap='viridis',
label='Temperature',
vmin = Tmin,
vmax = Tmax,

)
MAX = Tmax

cbaxes = inset_axes(ax,width="4%", height="60%",loc='right',bbox_to_anchor=(0.1,0,1,1),bbox_transform=ax.transAxes) 
cbar = fig.colorbar(cimg, ax=ax, cax=cbaxes)
cbar.ax.get_yaxis().labelpad = 12.5
cbar.ax.tick_params(labelsize=14,color='k')
cbar.outline.set_edgecolor('k')
cbar.ax.yaxis.set_ticks_position('right')
ax.text(1.15,0.7,r'T/k',fontsize=14)

ax.set_title("TEST")
ax.legend(scatterpoints=50)

plt.show()
plt.cla()
fig.clf() 
plt.close()

[Albertringwood]

Dear Dr Dao I have been able to plot the script as I desired:

Now I'm challenging problems about plotting my sim according to a plane (Zcut). The classic script doesn't work. What can I do? I attach there my script.

TE2.txt

The error that I get is:

` File "/opt/homebrew/lib/python3.10/site-packages/matplotlib/axes/_axes.py", line 4367, in _parse_scatter_color_args colors = mcolors.to_rgba_array(c) File "/opt/homebrew/lib/python3.10/site-packages/matplotlib/colors.py", line 487, in to_rgba_array rgba = np.array([to_rgba(cc) for cc in c]) File "/opt/homebrew/lib/python3.10/site-packages/matplotlib/colors.py", line 487, in rgba = np.array([to_rgba(cc) for cc in c]) File "/opt/homebrew/lib/python3.10/site-packages/matplotlib/colors.py", line 299, in to_rgba rgba = _to_rgba_no_colorcycle(c, alpha) File "/opt/homebrew/lib/python3.10/site-packages/matplotlib/colors.py", line 381, in _to_rgba_no_colorcycle raise ValueError(f"Invalid RGBA argument: {orig_c!r}") ValueError: Invalid RGBA argument: 1118.5159912109375

The above exception was the direct cause of the following exception:

Traceback (most recent call last): File "/Volumes/MacSwift/macswift/examples/IMPACTS/temperaturetest/TE2.py", line 126, in cimg = ax.scatter( File "/opt/homebrew/lib/python3.10/site-packages/matplotlib/init.py", line 1423, in inner return func(ax, *map(sanitize_sequence, args), **kwargs) File "/opt/homebrew/lib/python3.10/site-packages/matplotlib/axes/_axes.py", line 4530, in scatter self._parse_scatter_color_args( File "/opt/homebrew/lib/python3.10/site-packages/matplotlib/axes/_axes.py", line 4373, in _parse_scatter_color_args raise invalid_shape_exception(c.size, xsize) from err ValueError: 'c' argument has 1577062 elements, which is inconsistent with 'x' and 'y' with size 788558. `

I have no idea how to solve it

Thank you again Albert.

[Albertringwood]

In the end I managed to plot a section. The problem is that the values are totally wrong.

I don’t know why, but the values are like reversed. However, going forward with the simulation, the temeprature increases normally

The correct values should be these (thanks @francyrad for editing the original script)

Density values are also wrong, it's like they are reversed.

I don't have any idea why this happens. It may be due to the z rapresentation? I inspired from the new script in this guide and I wrote:

# Z<0 
` sel = np.where(pos[:, 2] < 100)[0]
matid = matid[sel]
id = id[sel]
m = m[sel]
pos = pos[sel]

# Sort in z order so higher particles are plotted on top
sort = np.argsort(pos[:, 2])
matid = matid[sort]
id = id[sort]
m = m[sort]
pos = pos[sort]`

If I write a number less than 1, it won't load all the particles. If I put a number a bigger than 1 but less than 10, it will only plot some snapshot, then I'll receive the error.

Whan i plotted without Z, values seemed corrects:

This is the script, I have no idea what to do honestly to plot the correct values...

help.txt

[JingyaoDOU]

Hi Alert,

Sorry for the late reply, I'm current not at school.

I think the reason is the when you set sel = np.where(pos[:, 2] < 100)[0] in the script, you will actually select basically all the particles in this way since most particles having z value less than 100 Earth Radius. I slightly changed your script, and marker the changed part with doller sign. I also add T = T[sel] and T = T[sort] since only partcial particles will be selected in this way. This should now plot the right value.

#!/usr/bin/env python
import woma
import swiftsimio as sw
import numpy as np
import unyt
import h5py
import matplotlib.pyplot as plt
import matplotlib
import matplotlib as mpl
from woma.misc import utils
from copy import copy
from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable
from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes
#import pylustrator
#pylustrator.start()

woma.load_eos_tables() # load the eos table for the calculation of entropy at later stage

snapshotfile = "example.hdf5"
#snapshotfile = "example_%04d.hdf5"%i
data    = sw.load(snapshotfile)
box_mid = 0.5 * data.metadata.boxsize[0].to(unyt.Rearth)
pos = data.gas.coordinates.to(unyt.Rearth) - box_mid
data.gas.densities.convert_to_mks()
rho_mks = np.array(data.gas.densities)
data.gas.internal_energies.convert_to_mks()
u_mks   = np.array(data.gas.internal_energies) # internal energy in mks
matid   = np.array(data.gas.material_ids)
data.gas.masses.convert_to_mks()
m       = np.array(data.gas.masses)
m = data.gas.masses
ax_lim = 5
x_min = box_mid - ax_lim * unyt.Rearth
x_max = box_mid + ax_lim * unyt.Rearth
load_region = [[x_min, x_max], [x_min, x_max], [x_min, box_mid]]
data.gas.pressures.convert_to_cgs()
p = np.array(data.gas.pressures)
entropy = woma.eos.eos.A1_s_u_rho(u_mks, rho_mks, matid)
T       = woma.eos.eos.A1_T_u_rho(u_mks, rho_mks, matid)

pos_centerM = np.sum(pos * m[:, np.newaxis], axis=0) / np.sum(m)
pos -= pos_centerM

#box_mid = 0.5 * data.metadata.boxsize[0]
#pos = data.gas.coordinates - box_mid
id = data.gas.particle_ids
#mat_id = data.gas.material_ids.value

#idoff = 200000000
#bodyoff = 100000000

# $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
sel = np.where(pos[:, 2] < .01)[0]  # half planet below with z<0                         #$
sel = np.where(np.logical_and(pos[:, 2] < .1,pos[:, 2] > -.1))[0]  #Thin slice cut       #$
# $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

matid = matid[sel]
id = id[sel]
m = m[sel]
T = T[sel]
pos = pos[sel]

# Sort in z order so higher particles are plotted on top
sort = np.argsort(pos[:, 2])
matid = matid[sort]
id = id[sort]
m = m[sort]
T = T[sort]
pos = pos[sort]

font_size = 20
params = {
"axes.labelsize": font_size,
"font.size": font_size,
"xtick.labelsize": font_size,
"ytick.labelsize": font_size,
"font.family": "serif",
}
matplotlib.rcParams.update(params)

# Material IDs ( = type_id * type_factor + unit_id )
type_factor = 100
type_idg = 0
type_Til = 1
type_HHe = 2
type_SESAME = 3
type_ANEOS = 4
id_body = 1000000
outtime   =100  #insert time in seconds of each snapshot output
txtsize   =15
font_size = 20

# Name and ID
Di_mat_id = {
"ANEOS_iron": type_ANEOS * type_factor + 1,
"ANEOS_iron_2": type_ANEOS * type_factor + 1 + id_body,
"ANEOS_forsterite": type_ANEOS * type_factor,
"ANEOS_forsterite_2": type_ANEOS * type_factor + id_body,
}

Tmax = 4000
Tmin = 1000

fig, ax = plt.subplots(1,1,figsize=(15,15))
cimg = ax.scatter(
pos[:, 0],
pos[:, 1],
s=10,
c=T,
edgecolors="none",
marker=".",
alpha=1,
cmap='hot',
label='Temperature [K]',
vmax = Tmax,
vmin = Tmin,

)
MAX = Tmax

cbaxes = inset_axes(ax,width="4%", height="100%",loc='right',bbox_to_anchor=(0.1,0,1,1),bbox_transform=ax.transAxes)
cbar = fig.colorbar(cimg, ax=ax, cax=cbaxes)
cbar.ax.get_yaxis().labelpad = 12.5
cbar.ax.tick_params(labelsize=20, colors='w') #Size of label's numbers
cbar.outline.set_edgecolor('w')
cbar.ax.yaxis.set_ticks_position('right')

ax.set_xlim(-ax_lim, ax_lim)
ax.set_xlim(ax_lim, -ax_lim)
ax.set_yticks(ax.get_xticks())
ax.set_ylim(ax_lim, -ax_lim)
ax.set_xlabel(r"x Position ($R_\oplus$)")
ax.set_ylabel(r"y Position ($R_\oplus$)")
ax.set_facecolor("black")
fig.patch.set_facecolor('black')
ax.legend(scatterpoints=50)
plt.subplots_adjust(right=0.85)

ax.yaxis.label.set_color('w')
ax.xaxis.label.set_color('w')
ax.tick_params(axis='x', colors='w',labelsize=20)
ax.tick_params(axis='y', colors='w',labelsize=20)

plt.show()
plt.cla()
fig.clf()
plt.close()

[Albertringwood]

Dear Dr. Dao I'm very happy to announce the the script works! :D Thank you very much for your time! Later I'll test other quantities to avoid other issues and I'll will close this one. Thank you again Albert

[Albertringwood]

The script works perfectly also for other quantities. I'm in debt with you! The only thing that I don't understand is why my legend moves while the simulation is going, for example:

Is there a way to make it stable? Albert

[JingyaoDOU]

Hi Albert,

Sorry for the super late reply, I did not notice your last message. I assume you already fix the issue but if not, the problem might be fixed by fixing the location of the label with

ax.legend(scatterpoints=50,loc = 'upper right')

Hope this should help.

[Albertringwood]

It didn't work :(

[JingyaoDOU]

How about trying

fig.legend(loc = 'upper right')

or

ax.legend(loc = 'upper right')

Please share the code if it's still not working.

[Albertringwood]

Sorry for my late response. I have been busy. The point is that python got a new upgrade and i cannot install both the classic woma and swiftsimio., so i'm waiting that swiftsim developers fix all :(. waiting for the solution of the problem, I noticed that Jacob has updated the original version of woma. do you think to integrate it with woma0? Let me know and thank you for your work

[JingyaoDOU]

Hi Albert,

Maybe you could try to use an older python version instead. I'm still using python 3.8 for most of my conda environments. I have not tried the lastest python, although I heard it should be really fast

Jacob has updated mainly the header of the eos tables and has no other major updates. I will find some time to try to merge woma0 with the last woma.

[Albertringwood]

I solved using "python3.10 script.py". Regarding the legend problem

ax.legend(loc = 'upper right')

doesn't work

instead

fig.legend(loc = 'upper right')

works, but put the legend outside the simulation box. It's not a problem however. The problem account when i want to put it in the "lower left", so it doesn't look bad alongside a title. The problem is that it creates me an upper and a lower right legend. the other problem is that when I try to plot the label logarithmically, it doesn't create barlines for me. This makes it difficult to appreciate the exact values. I searched all over the internet for a way to put them, but nothing. it's as if in some graphs it sets them automatically, in others, like this one, it doesn't. Would you know how to do it? Thanks a lot again for your help. I am attaching the images and the code.

THE CODE:

    fig, ax = plt.subplots(1,1,figsize=(15,15))
    cimg = ax.scatter(
    pos[:, 0],
    pos[:, 1],
    s=10,
    c=u_mks,
    edgecolors="none",
    marker=".",
    alpha=1,
    cmap='plasma',
    norm=mpl.colors.LogNorm(vmin=1000000, vmax=120000000),
    label='Specific Internal Energy [J/Kg]',
    #vmax = umax,
    #vmin = umin,

    )
    #MAX = umax

    cbaxes = inset_axes(ax,width="4%", height="100%",loc='right',bbox_to_anchor=(0.1,0,1,1),bbox_transform=ax.transAxes)
    cbar = fig.colorbar(cimg, ax=ax, cax=cbaxes)
    cbar.ax.get_yaxis().labelpad = 12.5
    cbar.ax.tick_params(labelsize=15, colors='w') #Size of label's numbers
    cbar.outline.set_edgecolor('w')
    cbar.ax.yaxis.set_ticks_position('right')
    cbar.ax.set_yscale("symlog")
    ticklabels = cbar.ax.get_ymajorticklabels()
    ticks = list(cbar.get_ticks())
    fig.legend(loc = 'lower right')
    fig.legend(scatterpoints=50)

    #ax.legend(scatterpoints=50,loc = 'upper right')
    #cbar.ax.yaxis.set_major_locator(mpl.ticker.LogLocator())
    #cbar.ax.yaxis.set_major_formatter(mpl.ticker.LogFormatter())

    ax.set_xlim(-ax_lim, ax_lim)
    ax.set_xlim(ax_lim, -ax_lim)
    ax.set_yticks(ax.get_xticks())
    ax.set_ylim(ax_lim, -ax_lim)
    ax.set_xlabel(r"x Position ($R_\oplus$)")
    ax.set_ylabel(r"y Position ($R_\oplus$)")
    ax.set_facecolor("black")
    fig.patch.set_facecolor('black')
    #ax.legend(scatterpoints=50)
    plt.subplots_adjust(right=0.85)

    ax.yaxis.label.set_color('w')
    ax.xaxis.label.set_color('w')
    ax.tick_params(axis='x', colors='w',labelsize=20)
    ax.tick_params(axis='y', colors='w',labelsize=20)

    saveloc='specific_internal_energy'
    save = "specific_internal_energy/log_%04d.png"%i
    fig.savefig(save, dpi=300)

WHAT I GET

WHAT I'D LKE TO GET (BARLINES)

Thank you so much for your work

[Francyrad]

@Albertringwood ChatGPT helped a lot to solve the issues. Here it's your edited full working script:

#!/usr/bin/env python
import woma
import swiftsimio as sw
import numpy as np
import unyt
import h5py
import matplotlib.pyplot as plt
import matplotlib
import matplotlib as mpl
from matplotlib.colors import LogNorm
from matplotlib.ticker import LogLocator, LogFormatter, FuncFormatter, FixedLocator # Import LogLocator and LogFormatter
from woma.misc import utils
from copy import copy
from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable
from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes

#import pylustrator
#pylustrator.start()

woma.load_eos_tables() # load the eos table for the calculation of entropy at later stage

internal_energy = "internal_energy"
if not os.path.exists(internal_energy):
    os.makedirs(internal_energy)

def log_formatter(x, pos):
    exponent = np.log10(x)
    if exponent.is_integer():
        return r'$10^{{{:.0f}}}$'.format(exponent)
    else:
        return ''

for i in range(1330, 1369):
    snapshotfile = "canup__%04d.hdf5"%i
    data    = sw.load(snapshotfile)
    box_mid = 0.5 * data.metadata.boxsize[0].to(unyt.Rearth)
    pos = data.gas.coordinates.to(unyt.Rearth) - box_mid
    data.gas.densities.convert_to_mks()
    rho_mks = np.array(data.gas.densities)
    data.gas.internal_energies.convert_to_mks()
    u_mks   = np.array(data.gas.internal_energies) # internal energy in mks
    matid   = np.array(data.gas.material_ids)
    data.gas.masses.convert_to_mks()
    m       = np.array(data.gas.masses)
    m = data.gas.masses
    ax_lim = 5
    x_min = box_mid - ax_lim * unyt.Rearth
    x_max = box_mid + ax_lim * unyt.Rearth
    load_region = [[x_min, x_max], [x_min, x_max], [x_min, box_mid]]
    data.gas.pressures.convert_to_cgs()
    p = np.array(data.gas.pressures)
    entropy = woma.eos.eos.A1_s_u_rho(u_mks, rho_mks, matid)
    T       = woma.eos.eos.A1_T_u_rho(u_mks, rho_mks, matid)

    pos_centerM = np.sum(pos * m[:, np.newaxis], axis=0) / np.sum(m)
    pos -= pos_centerM

    id = data.gas.particle_ids

    sel = np.where(pos[:, 2] < .01)[0]  # half planet below with z<0                         #$
    sel = np.where(np.logical_and(pos[:, 2] < .1,pos[:, 2] > -.1))[0]  #Thin slice cut       #$

    matid = matid[sel]
    id = id[sel]
    m = m[sel]
    u_mks = u_mks[sel]
    pos = pos[sel]

    # Sort in z order so higher particles are plotted on top
    sort = np.argsort(pos[:, 2])
    matid = matid[sort]
    id = id[sort]
    m = m[sort]
    u_mks = u_mks[sort]
    pos = pos[sort]

    font_size = 20
    params = {
    "axes.labelsize": font_size,
    "font.size": font_size,
    "xtick.labelsize": font_size,
    "ytick.labelsize": font_size,
    "font.family": "serif",
    }
    matplotlib.rcParams.update(params)

    # Material IDs ( = type_id * type_factor + unit_id )
    type_factor = 100
    type_idg = 0
    type_Til = 1
    type_HHe = 2
    type_SESAME = 3
    type_ANEOS = 4
    id_body = 1000000
    outtime   =100  #insert time in seconds of each snapshot output
    txtsize   =15
    font_size = 20

    # Name and ID
    Di_mat_id = {
    "ANEOS_iron": type_ANEOS * type_factor + 1,
    "ANEOS_iron_2": type_ANEOS * type_factor + 1 + id_body,
    "ANEOS_forsterite": type_ANEOS * type_factor,
    "ANEOS_forsterite_2": type_ANEOS * type_factor + id_body,
    }

    #umax = 120000000
    #umin = 1000000

    fig, ax = plt.subplots(1,1,figsize=(15,15))
    cimg = ax.scatter(
    pos[:, 0],
    pos[:, 1],
    s=10,
    c=u_mks,
    edgecolors="none",
    marker=".",
    alpha=1,
    cmap='plasma',
    norm=mpl.colors.LogNorm(vmin=1000000, vmax=120000000),
    label='Specific Internal Energy [J/Kg]',
    #vmax = umax,
    #vmin = umin,

    )
    #MAX = umax

    cbaxes = inset_axes(ax,width="4%", height="100%",loc='right',bbox_to_anchor=(0.1,0,1,1),bbox_transform=ax.transAxes)
    cbar = fig.colorbar(cimg, ax=ax, cax=cbaxes)
    cbar.ax.get_yaxis().labelpad = 12.5
    cbar.ax.tick_params(labelsize=15, colors='w')
    cbar.outline.set_edgecolor('w')

    cbar = fig.colorbar(cimg, ax=ax, cax=cbaxes)
    cbar.ax.get_yaxis().labelpad = 12.5
    cbar.ax.tick_params(labelsize=15, colors='w')
    cbar.outline.set_edgecolor('w')

    cbar.ax.yaxis.set_major_locator(LogLocator(subs='all'))  # Set tick positions using LogLocator with all logarithmic values
    cbar.ax.yaxis.set_major_formatter(FuncFormatter(log_formatter))  # Format ticks using FuncFormatter with the custom log_formatter function

    fig.legend(loc='lower right', bbox_to_anchor=(1, 0), scatterpoints=50)
    ax.set_title('title simulation. (C) Francesco Radica - Francyrad', color='w', fontsize=22)

    ax.text(0.9, 0.9, r'            {0:6.2f} h'.format(i * outtime / 3600),
        transform=ax.transAxes, ha='right', va='top', fontsize=txtsize, color='w')

    ax.set_xlim(-ax_lim, ax_lim)
    ax.set_yticks(ax.get_xticks())
    ax.set_ylim(-ax_lim, ax_lim)
    ax.set_xlabel(r"x Position ($R_\oplus$)")
    ax.set_ylabel(r"y Position ($R_\oplus$)")
    ax.set_facecolor("black")
    fig.patch.set_facecolor('black')
    plt.subplots_adjust(right=0.85)

    ax.yaxis.label.set_color('w')
    ax.xaxis.label.set_color('w')
    ax.tick_params(axis='x', colors='w',labelsize=20)
    ax.tick_params(axis='y', colors='w',labelsize=20)

    saveloc='specific_internal_energy'
    save = "specific_internal_energy/log_%04d.png"%i
    fig.savefig(save, dpi=300)

    #plt.show()
    #plt.cla()
    #fig.clf()
    #plt.close()

Thanks for sharing it, it will help a lot! Francesco

[Albertringwood]

It worked perfectly, thank you! I'll try to use ChatGPT for python things if it is so useful!

Thank you to both for all the help. @Francyrad please make a page with all the scripts when you will have time.

Albert