From f2184d2310a085600681a47e8a4ab32c8b3aaf82 Mon Sep 17 00:00:00 2001
From: salman <zaher.salman@psi.ch>
Date: Tue, 14 Apr 2026 17:28:23 +0200
Subject: [PATCH] Clarify TrimSP helper routines

---
 fortran/trimspNL.F | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/fortran/trimspNL.F b/fortran/trimspNL.F
index 699dde7..6ec41bc 100644
--- a/fortran/trimspNL.F
+++ b/fortran/trimspNL.F
@@ -3933,6 +3933,10 @@ C
       REAL*8 X1SY,X2SY,X3SY,X4SY,X5SY,X6SY,X1S,X2S,X3S,X4S,X5S,X6S,Y
       LOGICAL EQUAL
 
+C     Convert accumulated raw sums into per-particle moments by dividing
+C     through the population size Y. The higher-level MOMENTS/MOMENTN
+C     routines use these normalized values to build means, spreads and
+C     uncertainties for the reported observables.
       IF(EQUAL(Y,0.D0))GOTO 10
       X1SY=X1S/Y
       X2SY=X2S/Y
@@ -3950,6 +3954,9 @@ C
      &     ,SPHI(N)
       REAL*8 SRAT,CX2,CY2,CZ2,UNIT
 
+C     Rotate the current direction cosines by the sampled scattering
+C     angles (psi,phi). The output direction is explicitly renormalized
+C     because repeated collisions amplify small floating-point drift.
       DO IV=1,N
          SRAT=SPSI(IV)/SINE(IV)
          CX2=CPSI(IV)*COSX(IV)+SPSI(IV)*SINE(IV)*CPHI(IV)
@@ -3972,6 +3979,8 @@ C     MAKE SURE COSZ.NE.0.
       INTEGER n,I
 C     
 C     FETCH A NEW VELOCITY FROM A MAXWELLIAN FLUX AT A SURFACE
+C     Vectorized version used when multiple projectiles are initialized
+C     together with thermal surface velocities.
 C     
       REAL*8 FG(2*N),FFG(N),E(N),COSX(N),COSY(N),COSZ(N),SINE(N)
 C     DIMENSIOM E(N),COSX(N),COSY(N),COSZ(N),SINE(N)
@@ -3999,6 +4008,7 @@ C     DIMENSIOM E(N),COSX(N),COSY(N),COSZ(N),SINE(N)
       SUBROUTINE VELOC(E,COSX,COSY,COSZ,SINE)
 C
 C     FETCH A NEW VELOCITY FROM A MAXWELLIAN FLUX AT A SURFACE
+C     Scalar version that reuses Gaussian deviates generated by FGAUSS.
 C     
       IMPLICIT NONE
       INTEGER INIV1,INIV3
@@ -4041,6 +4051,8 @@ C     ZA=0.
 C     ZS=1.
 C     
 C     IT IS THE BOX-MUELLER METHOD
+C     FG returns symmetric normal deviates, while FFG stores the positive
+C     half-space distribution used for Maxwellian flux sampling.
 C     
       IMPLICIT NONE
       INTEGER IND,IND2,IANZ,JJ
@@ -4216,6 +4228,9 @@ C     why is it needed??
       INTEGER I,N,J,INC
       REAL*8 ARRAY(N),TARGT
 
+C     Return the first index whose value is strictly greater than TARGT.
+C     This is used to map a continuous coordinate/depth onto the current
+C     layer or histogram bin in monotonic arrays.
       J=1
       IF(INC.LT.0) J=N*(-INC)
       DO I=1,N
@@ -4233,6 +4248,8 @@ C     why is it needed??
       REAL*8 p1,p2,p3,pi
       real*4 random(2)
       DATA pi/3.14159265358979D0/
+C     Sample a Gaussian-distributed incident energy around E0 using a
+C     Box-Mueller transform. Epar is the per-projectile starting energy.
       call ranlux(random, 2)
       p1  = Esig*DSQRT(-2.D0*DLOG(1.D0-DBLE(random(1)))) 
       p2  = 2.D0*pi*DBLE(random(2))
@@ -4250,6 +4267,8 @@ C     why is it needed??
       REAL*8 p1,p2,p3,pi
       real*4 random(2)
       DATA pi/3.14159265358979D0/
+C     Sample a Gaussian-distributed incidence angle around ALPHA and
+C     return both the angle in radians (ALFA) and its sine/cosine.
       call ranlux(random, 2)	
       p1  = ALPHASIG*DSQRT(-2.D0*DLOG(1.D0-DBLE(random(1)))) 
       p2  = 2.D0*pi*DBLE(random(2))
@@ -4269,6 +4288,8 @@ C     why is it needed??
       REAL*8 F1,F2
       REAL*8 TINY
       PARAMETER  (TINY = 1.0D-10)
+C     Centralized floating-point comparison used throughout the transport
+C     code to avoid exact equality checks on derived real-valued quantities.
       IF (DABS(F1-F2).LE.TINY) THEN
          EQUAL = .TRUE.
       ELSE