bugfix convert snapshots to density

sbmy_control/analysis/power_spectrum.py

@@ -1,7 +1,7 @@
 import numpy as np
 import os
 
-kmin = 1e-1
+nfirst_kmodes = 10
 kmax = 2e0
 Nk = 50
 AliasingCorr=False
@@ -24,12 +24,17 @@ def crop_field(field, Ncrop):
         field.L2 = field.N2*d2
 
 
-def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk, G=None):
+def get_power_spectrum(field, nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, G=None):
     from pysbmy.power import PowerSpectrum
     from pysbmy.fft import FourierGrid
     from pysbmy.correlations import get_autocorrelation
 
     if G is None:
+        default_k_modes = PowerSpectrum(field.L0,field.L1,field.L2,int(field.N0),int(field.N1),int(field.N2),).FourierGrid.k_modes
+        if kmax is None:
+            kmax = default_k_modes[-1]
+        if kmax > default_k_modes[-1]:
+            raise ValueError(f"kmax ({kmax}) is larger than the maximum k mode available in the field ({default_k_modes[-1]}).")
         G = FourierGrid(
             field.L0,
             field.L1,
@@ -37,8 +42,8 @@ def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk, G=None):
             int(field.N0),
             int(field.N1),
             int(field.N2),
-            k_modes=np.concat([PowerSpectrum(field.L0,field.L1,field.L2,int(field.N0),int(field.N1),int(field.N2),).FourierGrid.k_modes[:10],np.logspace(
-                np.log10(kmin),
+            k_modes=np.concat([default_k_modes[:nfirst_kmodes],np.logspace(
+                np.log10(default_k_modes[nfirst_kmodes]),
                 np.log10(kmax),
                 Nk,
             )]),
@@ -51,13 +56,18 @@ def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk, G=None):
     return G, k, Pk
 
 
-def get_cross_correlations(field_A, field_B, kmin=kmin, kmax=kmax, Nk=Nk, G=None):
+def get_cross_correlations(field_A, field_B, nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, G=None):
     from pysbmy.power import PowerSpectrum
     from pysbmy.fft import FourierGrid
     from pysbmy.correlations import get_crosscorrelation
 
 
     if G is None:
+        default_k_modes = PowerSpectrum(field_A.L0,field_A.L1,field_A.L2,int(field_A.N0),int(field_A.N1),int(field_A.N2),).FourierGrid.k_modes
+        if kmax is None:
+            kmax = default_k_modes[-1]
+        if kmax > default_k_modes[-1]:
+            raise ValueError(f"kmax ({kmax}) is larger than the maximum k mode available in the field ({default_k_modes[-1]}).")
         G = FourierGrid(
             field_A.L0,
             field_A.L1,
@@ -65,8 +75,8 @@ def get_cross_correlations(field_A, field_B, kmin=kmin, kmax=kmax, Nk=Nk, G=None
             int(field_A.N0),
             int(field_A.N1),
             int(field_A.N2),
-            k_modes=np.concat([PowerSpectrum(field_A.L0,field_A.L1,field_A.L2,int(field_A.N0),int(field_A.N1),int(field_A.N2),).FourierGrid.k_modes[:10],np.logspace(
-                np.log10(kmin),
+            k_modes=np.concat([default_k_modes[:nfirst_kmodes],np.logspace(
+                np.log10(default_k_modes[nfirst_kmodes]),
                 np.log10(kmax),
                 Nk,
             )]),
@@ -104,7 +114,7 @@ def plot_power_spectra(filenames,
                        yticks = np.linspace(0.9,1.1,11),
                        bound1=0.01,
                        bound2=0.02,
-                       kmin=kmin,
+                       nfirst_kmodes=nfirst_kmodes,
                        kmax=kmax,
                        Nk=Nk,
                        figsize=(8,4),
@@ -139,7 +149,7 @@ def plot_power_spectra(filenames,
                                                   color=colors[i],
                                                   linestyle=linestyles[i],
                                                   marker=markers[i],),
-                                   power_args=dict(kmin=kmin,
+                                   power_args=dict(nfirst_kmodes=nfirst_kmodes,
                                                    kmax=kmax,
                                                    Nk=Nk),
                                     )
@@ -178,7 +188,7 @@ def plot_cross_correlations(filenames_A,
                             yticks = np.linspace(0.99,1.001,12),
                             bound1=0.001,
                             bound2=0.002,
-                            kmin=kmin,
+                            nfirst_kmodes=nfirst_kmodes,
                             kmax=kmax,
                             Nk=Nk,
                             figsize=(8,4),
@@ -217,7 +227,7 @@ def plot_cross_correlations(filenames_A,
                                                        color=colors[i],
                                                        linestyle=linestyles[i],
                                                        marker=markers[i],),
-                                       power_args=dict(kmin=kmin,
+                                       power_args=dict(nfirst_kmodes=nfirst_kmodes,
                                                         kmax=kmax,
                                                         Nk=Nk),
                                         )   
@@ -283,8 +293,15 @@ def console_main():
     parser.add_argument('-yrp', '--ylim_power', type=float, nargs=2, default=[0.9,1.1], help='Y-axis limits.')
     parser.add_argument('-yrc', '--ylim_corr', type=float, nargs=2, default=[0.99,1.001], help='Y-axis limits.')
     parser.add_argument('--crop', type=int, default=None, help='Remove the outter N pixels of the fields.')
+    parser.add_argument('--nfirst_kmodes', type=float, default=10, help='First k modes from all available k modes to be used in the power spectrum before the geomspace part.')
+    parser.add_argument('--kmax', type=float, default=None, help='Maximum k value for the power spectrum geomspace part. If None, it will be set to the maximum k mode available in the field.')
+    parser.add_argument('--Nk', type=int, default=50, help='Number of k values for the power spectrum geomspace part.')
 
     args = parser.parse_args()
+    
+    nfirst_kmodes = args.nfirst_kmodes
+    kmax = args.kmax
+    Nk = args.Nk
 
     if not args.power_spectrum and not args.cross_correlation:
         print('You must choose between power_spectrum and cross_correlation.')
@@ -298,7 +315,7 @@ def console_main():
         from pysbmy.field import read_field
         F_ref = read_field(args.directory+args.reference)
         crop_field(F_ref, args.crop)
-        G, _, Pk_ref = get_power_spectrum(F_ref, kmin=kmin, kmax=kmax, Nk=Nk)
+        G, _, Pk_ref = get_power_spectrum(F_ref, nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk)
     else:
         Pk_ref = None
         G = None
@@ -325,7 +342,7 @@ def console_main():
                            yticks = yticks_power,
                            bound1=0.01,
                            bound2=0.02,
-                           kmin=kmin,
+                           nfirst_kmodes=nfirst_kmodes,
                            kmax=kmax,
                            Nk=Nk,
                            ax=axes[0],
@@ -344,7 +361,7 @@ def console_main():
                                 yticks = yticks_corr,
                                 bound1=0.001,
                                 bound2=0.002,
-                                kmin=kmin,
+                                nfirst_kmodes=nfirst_kmodes,
                                 kmax=kmax,
                                 Nk=Nk,
                                 ax=axes[1],
@@ -371,7 +388,7 @@ def console_main():
                                     yticks = yticks_power,
                                      bound1=0.01,
                                      bound2=0.02,
-                                     kmin=kmin,
+                                     nfirst_kmodes=nfirst_kmodes,
                                      kmax=kmax,
                                      Nk=Nk,
                                     Ncrop=args.crop,
@@ -392,7 +409,7 @@ def console_main():
                                         yticks = yticks_corr,
                                           bound1=0.001,
                                           bound2=0.002,
-                                          kmin=kmin,
+                                          nfirst_kmodes=nfirst_kmodes,
                                           kmax=kmax,
                                           Nk=Nk,
                                         Ncrop=args.crop,

sbmy_control/low_level.py

@@ -64,10 +64,14 @@ def stdout_redirector(stream):
     try:
         # Create a temporary file and redirect stdout to it
         tfile = tempfile.TemporaryFile(mode="w+b")
+
         _redirect_stdout(tfile.fileno())
+
         # Yield to caller, then redirect stdout back to the saved fd
         yield
+
         _redirect_stdout(saved_stdout_fd)
+
         # Copy contents of temporary file to the given stream
         tfile.flush()
         tfile.seek(0, io.SEEK_SET)
@@ -215,9 +219,9 @@ def get_progress_from_logfile(filename):
     with open(filename, "r") as f:
         lines = f.readlines()
     for line in lines:
-        if " operation " in line:
+        if "Step " in line:
             try:
-                splitted_line = line.split(" operation ")[1]
+                splitted_line = line.split("Step ")[1]
                 splitted_line = splitted_line.split(":")[0]
                 current_operation = int(splitted_line.split("/")[0])
                 total_operations = int(splitted_line.split("/")[1])
@@ -237,7 +241,7 @@ def progress_bar_from_logfile(filename:str, desc:str="", verbose:int=1, **kwargs
     from tqdm import tqdm
     from time import sleep
     k=0
-    limit=600
+    limit=6000
     update_interval=0.2
     wait_until_file_exists(filename, verbose=verbose, limit=limit)
     current_operation, total_operations = get_progress_from_logfile(filename)

sbmy_control/main.py

@@ -199,25 +199,26 @@ def main(parsed_args):
 
 
 def check_consistency(card_dict, mode):
-    ## Check consistency of EvolutionMode and ModulePMCOLA
-    if mode == "PM" or mode == "allPM":
-        if card_dict["EvolutionMode"] != 1:
-            raise ValueError(f"EvolutionMode is not 1: EvolutionMode={card_dict["EvolutionMode"]}")
-        if card_dict["ModulePMCOLA"] != 1:
-            raise ValueError(f"ModulePMCOLA is not 1: ModulePMCOLA={card_dict["ModulePMCOLA"]}")
-    elif mode == "tCOLA" or mode == "alltCOLA":
-        if card_dict["EvolutionMode"] != 2:
-            raise ValueError(f"EvolutionMode is not 2: EvolutionMode={card_dict["EvolutionMode"]}")
-        if card_dict["ModulePMCOLA"] != 1:
-            raise ValueError(f"ModulePMCOLA is not 1: ModulePMCOLA={card_dict["ModulePMCOLA"]}")
-    elif mode == "LPT":
-        if card_dict["ModulePMCOLA"] !=0:
-            raise ValueError(f"ModulePMCOLA is not 0: ModulePMCOLA={card_dict["ModulePMCOLA"]}")
-    elif mode == "pre_sCOLA" or mode == "post_sCOLA" or mode == "sCOLA" or mode=="allsCOLA":
-        if card_dict["EvolutionMode"] != 3:
-            raise ValueError(f"EvolutionMode is not 3: EvolutionMode={card_dict["EvolutionMode"]}")
-        if card_dict["ModulePMCOLA"] != 1:
-            raise ValueError(f"ModulePMCOLA is not 1: ModulePMCOLA={card_dict["ModulePMCOLA"]}")
+    # ## Check consistency of EvolutionMode and ModulePMCOLA
+    # if mode == "PM" or mode == "allPM":
+    #     if card_dict["EvolutionMode"] != 1:
+    #         raise ValueError(f"EvolutionMode is not 1: EvolutionMode={card_dict["EvolutionMode"]}")
+    #     if card_dict["ModulePMCOLA"] != 1:
+    #         raise ValueError(f"ModulePMCOLA is not 1: ModulePMCOLA={card_dict["ModulePMCOLA"]}")
+    # elif mode == "tCOLA" or mode == "alltCOLA":
+    #     if card_dict["EvolutionMode"] != 2:
+    #         raise ValueError(f"EvolutionMode is not 2: EvolutionMode={card_dict["EvolutionMode"]}")
+    #     if card_dict["ModulePMCOLA"] != 1:
+    #         raise ValueError(f"ModulePMCOLA is not 1: ModulePMCOLA={card_dict["ModulePMCOLA"]}")
+    # elif mode == "LPT":
+    #     if card_dict["ModulePMCOLA"] !=0:
+    #         raise ValueError(f"ModulePMCOLA is not 0: ModulePMCOLA={card_dict["ModulePMCOLA"]}")
+    # elif mode == "pre_sCOLA" or mode == "post_sCOLA" or mode == "sCOLA" or mode=="allsCOLA":
+    #     if card_dict["EvolutionMode"] != 3:
+    #         raise ValueError(f"EvolutionMode is not 3: EvolutionMode={card_dict["EvolutionMode"]}")
+    #     if card_dict["ModulePMCOLA"] != 1:
+    #         raise ValueError(f"ModulePMCOLA is not 1: ModulePMCOLA={card_dict["ModulePMCOLA"]}")
+    pass
 
 
 def console_main():

sbmy_control/monofonic.py

@@ -32,7 +32,7 @@ def main_monofonic(parsed_args):
         print_message("Running monofonic in slurm mode.", 1, "monofonic", verbose=parsed_args.verbose)
         slurm_dict=parse_arguments_slurm(parsed_args)
         main_dict=parse_arguments_main(parsed_args)
-        slurm_script = slurm_dict["scripts"]+"monofonic.sh"
+        slurm_script = slurm_dict["scripts"]+f"monofonic_{main_dict['simname']}.sh"
 
         if not isfile(slurm_script) or parsed_args.force:
             print_message(f"SLURM script {slurm_script} does not exist (or forced). Creating it.", 2, "monofonic", verbose=parsed_args.verbose)

sbmy_control/parameters_card.py

@@ -47,6 +47,15 @@ def register_arguments_card(parser:ArgumentParser):
     parser.add_argument("--OutputFCsDensity", type=str, default=None, help="Output FCs density file.")
     parser.add_argument("--OutputFCsSnapshot", type=str, default=None, help="Output FCs snapshot file.")
     parser.add_argument("--OutputRngStateLPT", type=str, default=None, help="Output RNG state file.")
+    ## Tests with phiBCs and density
+    parser.add_argument("--WriteGravPot", type=bool, default=False, help="Write gravitational potential.")
+    parser.add_argument("--OutputGravitationalPotentialBase", type=str, default=None, help="Output gravitational potential base.")
+    parser.add_argument("--MeshGravPot", type=int, default=None, help="Mesh for gravitational potential.")
+    parser.add_argument("--WriteDensity", type=bool, default=False, help="Write density.")
+    parser.add_argument("--OutputDensityBase", type=str, default=None, help="Output density base.")
+    parser.add_argument("--MeshDensity", type=int, default=None, help="Mesh for density.")
+    parser.add_argument("--LoadPhiBCs", type=bool, default=False, help="Load phiBCs.")
+    parser.add_argument("--InputPhiBCsBase", type=str, default=None, help="Input phiBCs file base.")
     parser.add_argument("--WriteReferenceFrame", type=bool, default=False, help="Write reference frame (COCA).")
     parser.add_argument("--OutputMomentaBase", type=str, default=None, help="Output momenta base (COCA).")
     parser.add_argument("--ReadReferenceFrame", type=bool, default=False, help="Read reference frame (COCA).")
@@ -126,6 +135,16 @@ def parse_arguments_card(parsed_args):
             OutputFCsDensity=parsed_args.OutputFCsDensity,
             OutputFCsSnapshot=parsed_args.OutputFCsSnapshot,
             OutputRngStateLPT=parsed_args.OutputRngStateLPT,
+            ## Tests with phiBCs and density
+            WriteGravPot=parsed_args.WriteGravPot,
+            OutputGravitationalPotentialBase=parsed_args.OutputGravitationalPotentialBase,
+            MeshGravPot=parsed_args.MeshGravPot,
+            WriteDensity=parsed_args.WriteDensity,
+            OutputDensityBase=parsed_args.OutputDensityBase,
+            MeshDensity=parsed_args.MeshDensity,
+            LoadPhiBCs=parsed_args.LoadPhiBCs,
+            InputPhiBCsBase=parsed_args.InputPhiBCsBase,
+            ## Cosmological parameters
             h=cosmo_dict["h"],
             Omega_m=cosmo_dict["Omega_m"],
             Omega_b=cosmo_dict["Omega_b"],
@@ -211,6 +230,17 @@ def parse_arguments_card(parsed_args):
         card_dict["OutputFCsSnapshot"] = main_dict["resultdir"]+"final_particles_"+main_dict["simname"]+".gadget3"
     if card_dict["OutputRngStateLPT"] is None:
         card_dict["OutputRngStateLPT"] = main_dict["workdir"]+"rng_state.h5"
+    ## Tests with phiBCs and density
+    if card_dict["OutputGravitationalPotentialBase"] is None:
+        card_dict["OutputGravitationalPotentialBase"] = main_dict["workdir"]+"gravpot_"+main_dict["simname"]
+    if card_dict["MeshGravPot"] is None:
+        card_dict["MeshGravPot"] = card_dict["N_PM_mesh"]
+    if card_dict["OutputDensityBase"] is None:
+        card_dict["OutputDensityBase"] = main_dict["workdir"]+"density_"+main_dict["simname"]
+    if card_dict["MeshDensity"] is None:
+        card_dict["MeshDensity"] = card_dict["N_PM_mesh"]
+    if card_dict["InputPhiBCsBase"] is None:
+        card_dict["InputPhiBCsBase"] = main_dict["workdir"]+"gravpot_tCOLA"
     if card_dict["OutputMomentaBase"] is None:
         card_dict["OutputMomentaBase"] = main_dict["workdir"]+"momenta_"+main_dict["simname"]+"_"
     if card_dict["InputMomentaBase"] is None:
@@ -333,6 +363,15 @@ def create_parameter_card_dict(
                     OutputFCsDensity:str = 'fcs_density.h5',
                     OutputFCsSnapshot:str = 'fcs_particles.gadget3',
 
+                    ## Tests with phiBCs and density
+                    WriteGravPot:bool = True,
+                    OutputGravitationalPotentialBase:str = 'gravitational_potential.h5',
+                    MeshGravPot:int = 128,
+                    WriteDensity:bool = False,
+                    OutputDensityBase:str = 'density.h5',
+                    MeshDensity:int = 128,
+                    LoadPhiBCs:bool = False,
+                    InputPhiBCsBase:str = 'gravitational_potential.h5',
                     ## COCA parameters
                     WriteReferenceFrame:bool = False,
                     OutputMomentaBase:str = 'momenta_',
@@ -403,6 +442,46 @@ def create_parameter_card_dict(
             OutputLPTPotential1=OutputLPTPotential1,
             OutputLPTPotential2=OutputLPTPotential2,
             OutputTilesBase=OutputTilesBase,
+
+            # Tests with phiBCs and density
+            WriteGravPot=int(WriteGravPot),
+            OutputGravitationalPotentialBase=OutputGravitationalPotentialBase,
+            MeshGravPot=MeshGravPot,
+            WriteDensity=int(WriteDensity),
+            OutputDensityBase=OutputDensityBase,
+            MeshDensity=MeshDensity,
+            LoadPhiBCs=int(LoadPhiBCs),
+            InputPhiBCsBase=InputPhiBCsBase,
+            
+            
+            # P3M default parameters
+            n_Tiles                      = 60,
+            # Number of tiles in each dimension for the P3M code
+            nPairsForceDiagnostic               =0,
+            # Number of particle-hole pairs to be used for the force diagnostic
+            nBinsForceDiagnostic               =0,
+            maxTrialsForceDiagnostic             =0,
+            OutputForceDiagnostic ='',
+            PrintOutputTimestepsLog             = 0,
+            OutputTimestepsLog               = "",
+            fac_dyn_custom                 = 0.03,
+            da_max_early_custom               = 0.0013,
+            fac_H_custom                 =  0.05,
+            fac_bend                    = 0.3125,
+            sub_bend1                    = 0.012,
+            sub_bend2                    = 0.014,
+            fac_p3m_fit                   = 0.15,
+            da_max_late_custom               = 0.02,
+            # Polynomial coefficients for the P3M time step limiter. Only used if use_p3m_fit=1
+            use_p3m_fit                   = 1,
+            p3m_fit_coeff0                 = -3.19544528,
+            p3m_fit_coeff1                 = -0.01948223,
+            p3m_fit_coeff2                 = 0.5180241,
+            p3m_fit_coeff3                 = -1.08372045,
+            p3m_fit_coeff4                 = -0.71808255,
+            p3m_fit_coeff5                 = -0.14370312,
+            RunForceDiagnostic               = 0,
+
             h=h,
             Omega_m=Omega_m,
             Omega_b=Omega_b,

sbmy_control/parameters_monofonic.py

@@ -95,8 +95,6 @@ def get_config_from_dict(monofonic_dict):
             "Omega_L": monofonic_dict["Omega_q"],
             "H0": monofonic_dict["h"]*100.,
             "n_s": monofonic_dict["n_s"],
-            "sigma_8": monofonic_dict["sigma8"],
-            "A_s": monofonic_dict["A_s"],
             "Tcmb": monofonic_dict["Tcmb"],
             "k_p": monofonic_dict["k_p"],
             "N_ur": monofonic_dict["N_ur"],
@@ -110,7 +108,14 @@ def get_config_from_dict(monofonic_dict):
             "transfer": "CLASS",
             "ztarget": 2.5,
             "ZeroRadiation": False
-        } 
+        }
+
+        if "A_s" in monofonic_dict.keys() and monofonic_dict["A_s"] is not None:
+            config["cosmology"]["A_s"] = monofonic_dict["A_s"]
+        elif "sigma8" in monofonic_dict.keys() and monofonic_dict["sigma8"] is not None:
+            config["cosmology"]["sigma8"] = monofonic_dict["sigma8"]
+        else:
+            raise KeyError("Either A_s or sigma8 must be provided in the cosmology parameters.")
 
     config["random"] = {
         "generator": "NGENIC",

sbmy_control/scripts/convert_snapshot_to_density.py

@@ -1,9 +1,7 @@
-from pysbmy.density import get_density_pm_snapshot
-from pysbmy.snapshot import read_snapshot
 import argparse
 
 
-def convert_snapshot_to_density(snapshot_path, output_path, N=None, corner=(0.0, 0.0, 0.0)):
+def convert_snapshot_to_density(snapshot_path, output_path, N=None, corner=(0.0, 0.0, 0.0), time=1.0):
     """
     Convert a snapshot to a density field.
 
@@ -18,6 +16,9 @@ def convert_snapshot_to_density(snapshot_path, output_path, N=None, corner=(0.0,
     corner : tuple of float
         Corner of the box (x, y, z).
     """
+    from pysbmy.density import get_density_pm_snapshot
+    from pysbmy.snapshot import read_snapshot
+
     # Read the snapshot
     print("Reading snapshot...")
     snap = read_snapshot(snapshot_path)
@@ -28,12 +29,86 @@ def convert_snapshot_to_density(snapshot_path, output_path, N=None, corner=(0.0,
     # Calculate density
     print("Calculating density...")
     F = get_density_pm_snapshot(snap, N, N, N, corner[0], corner[1], corner[2])
+    F.time = time  # Set the time for the field
 
     # Write density to file
     print("Writing density...")
     F.write(output_path)
     print("Density written to", output_path)
     print("Done.")
+    
+
+def convert_snapshots_to_density(snapshot_base_path, output_path, N=None, corner=(0.0, 0.0, 0.0), time=1.0):
+    """
+    Convert multiple snapshots to density fields.
+
+    Parameters
+    ----------
+    snapshot_base_path : str
+        Base path for the snapshot files.
+    output_path : str
+        Path to the output density file.
+    N : int
+        Size of the density field grid (N x N x N).
+    corner : tuple of float
+        Corner of the box (x, y, z).
+    """
+    from pysbmy.density import get_density_pm, density_to_delta
+    from pysbmy.snapshot import read_snapshot
+    from pysbmy.field import Field
+    import numpy as np
+    
+    # Get all snapshot files with path "snapshot_base_path*"
+    import glob
+    
+    snapshot_files = glob.glob(snapshot_base_path + "*")
+    if not snapshot_files:
+        raise FileNotFoundError(f"No snapshot files found at {snapshot_base_path}")
+    
+    A = None
+    L0 = L1 = L2 = None
+    d0 = d1 = d2 = None
+    
+    for snapshot_path in snapshot_files:
+        # Read the snapshot
+        print(f"Reading snapshot {snapshot_path}...")
+        snap = read_snapshot(snapshot_path)
+        
+        if N is None:
+            N = snap.Np0
+        
+        if L0 is None:
+            # Get the grid parameters from the snapshot
+            L0 = snap.L0
+            L1 = snap.L1
+            L2 = snap.L2
+            d0 = L0 / N
+            d1 = L1 / N
+            d2 = L2 / N
+        else:
+            # Ensure the grid parameters match the first snapshot
+            if not (snap.L0 == L0 and snap.L1 == L1 and snap.L2 == L2):
+                raise ValueError(f"All snapshots must have the same grid parameters. Got {snap.L0}, {snap.L1}, {snap.L2} but expected {L0}, {L1}, {L2}.")
+        
+        if A is None:
+            # Initialize the density field
+            A = np.zeros((N, N, N), dtype=np.float32)
+
+        # Calculate density
+        print("Calculating density for this snapshot...")
+        A = get_density_pm(snap.pos, A, d0=d0, d1=d1, d2=d2)
+    
+    # Convert density to delta field
+    print("Converting density to delta field...")
+    A = density_to_delta(A)
+    
+    # Create a Field object
+    F = Field(L0=L0,L1=L1,L2=L2,corner0=corner[0],corner1=corner[1],corner2=corner[2],rank=1,N0=N,N1=N,N2=N,time=time,data=A)
+    # Write density to file
+    print("Writing density...")
+    F.write(output_path)
+    print("Density written to", output_path)
+    print("Done.")
 
 
 def console_main():
@@ -70,13 +145,23 @@ def console_main():
 
     args = parser.parse_args()
 
-    convert_snapshot_to_density(
-        snapshot_path=args.snapshot,
-        output_path=args.output,
-        N=args.N,
-        corner=args.corner,
-    )
+    if args.snapshot.endswith("h5") or args.snapshot.endswith("hdf5") or args.snapshot.endswith("gadget3"):
+        # If the snapshot is a single file, convert it to density
+        convert_snapshot_to_density(
+            snapshot_path=args.snapshot,
+            output_path=args.output,
+            N=args.N,
+            corner=args.corner,
+        )
+    else:
+        # If the snapshot is a base path, convert all snapshots to density
+        convert_snapshots_to_density(
+            snapshot_base_path=args.snapshot,
+            output_path=args.output,
+            N=args.N,
+            corner=args.corner,
+        )
 
 
 if __name__ == "__main__":
-    console_main()
+    console_main()

sbmy_control/simbelmyne.py

@@ -16,7 +16,7 @@ def main_simbelmyne(parsed_args):
 
         if isfile(log_file): # Remove the preexisting log file to allow for the progress_bar to be run normally
             from os import remove
-            oremove(log_file)
+            remove(log_file)
 
         command_args = ["simbelmyne", paramfile, log_file]