From af3f0226137e0d8813aee9554d8740cf9daf7d3a Mon Sep 17 00:00:00 2001
From: Anna Wirbel <anna@avaframe.org>
Date: Tue, 30 Jun 2026 13:00:50 +0200
Subject: [PATCH] add debug plot (com1 out3)

---
 avaframe/com1DFA/com1DFA.py        |  4 ++
 avaframe/out3Plot/outDebugPlots.py | 60 +++++++++++++++++++++++++++++-
 2 files changed, 63 insertions(+), 1 deletion(-)

diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py
index 41ce01f81..e539aa29a 100644
--- a/avaframe/com1DFA/com1DFA.py
+++ b/avaframe/com1DFA/com1DFA.py
@@ -1629,6 +1629,9 @@ def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", rel
         and not cfg.getboolean("initialiseParticlesFromFile")
         and len(relThField) == 0
     ):
+        if debugPlot:
+            xyParticlesAll = {"x": particles["x"], "y": particles["y"]}
+
         particles = geoTrans.checkParticlesInRelease(
             particles, releaseLine, cfg.getfloat("thresholdPointInPoly")
         )
@@ -1671,6 +1674,7 @@ def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", rel
 
     if debugPlot:
         debPlot.plotPartIni(particles, dem)
+        debPlot.plotParticlesRelease(particles, relRaster, releaseLine, dem, cfg, xyParticlesAll)
 
     # space for deposited particles TODO: extra dict for them?
     particles["stoppedParticles"] = {
diff --git a/avaframe/out3Plot/outDebugPlots.py b/avaframe/out3Plot/outDebugPlots.py
index 504f0744e..00e103508 100644
--- a/avaframe/out3Plot/outDebugPlots.py
+++ b/avaframe/out3Plot/outDebugPlots.py
@@ -362,7 +362,7 @@ def plotProfile(s, z, idsBetaPoint):
 def plotVolumeRelease(releaseLine, relThField, releaseLineField):
     """ create a plot of the release line raster, the relThField for release thickness,
         releaseLineField - combination of relThField and release line raster mask """
-        
+
     fig = plt.figure()
     ax1 = fig.add_subplot(131)
     ax2 = fig.add_subplot(132)
@@ -374,3 +374,61 @@ def plotVolumeRelease(releaseLine, relThField, releaseLineField):
     fig.colorbar(im1, ax=ax2)
     fig.colorbar(im2, ax=ax3)
     plt.show()
+
+
+def plotParticlesRelease(particles, relRaster, releaseLine, dem, cfg, xyParticlesAll):
+    """plot the release raster and release polygon outline and the particles that were placed and the ones that
+    are finally used for the computation after removing the ones that are outside the release polygon
+    - used in particleInitialisation to visualize initialisation process in com1DFA/initializeParticles
+
+    Parameters
+    ------------
+    particles: dict
+        final particles dict after initialization
+    relRaster: numpy ndarray
+        release raster
+    releaseLine: dict
+        release polygon info dict
+    dem: dict
+        dem info dict
+    cfg: configparser object
+        simulation configuration settings, requires rho
+    xyParticlesAll: dict
+        particles dict during initialization before removing those outside the release polygon
+
+    """
+
+    # compute volume of particles
+    volParticles = particles["mTot"] / cfg.getfloat("rho")
+
+    extentCellCenters, extentCellCorners = pU.createExtentMinMax(
+        relRaster, dem["originalHeader"], originLLCenter=True
+    )
+
+    # figure
+    fig, ax = plt.subplots(nrows=1, ncols=1)
+    # Minor ticks
+    ax.set_xticks(
+        np.arange(extentCellCorners[0], extentCellCorners[1], dem["originalHeader"]["cellsize"]), minor=True
+    )
+    ax.set_yticks(
+        np.arange(extentCellCorners[2], extentCellCorners[3], dem["originalHeader"]["cellsize"]), minor=True
+    )
+    # Gridlines based on minor ticks
+    ax.grid(which="minor", color="w", linestyle="-", linewidth=2)
+    ax.imshow(relRaster, extent=extentCellCorners, origin="lower")
+    # plot all particles before removing the ones outside of release polygon
+    ax.plot(
+        xyParticlesAll["x"] + dem["originalHeader"]["xllcenter"],
+        xyParticlesAll["y"] + dem["originalHeader"]["yllcenter"],
+        "+g",
+    )
+    # only particles that have not been removed
+    ax.plot(
+        particles["x"] + dem["originalHeader"]["xllcenter"],
+        particles["y"] + dem["originalHeader"]["yllcenter"],
+        "*r",
+    )
+    ax.plot(releaseLine["x"], releaseLine["y"], "-b")
+    ax.set_title("mass/rho: %.2fm3" % (volParticles))
+    plt.show()