From 0bc8ea3d10c0cfa9a6749d74eccd7f9478d08932 Mon Sep 17 00:00:00 2001
From: antoinepado <sanxi.padovani@gmail.com>
Date: Sun, 25 Oct 2026 14:28:48 +0100
Subject: [PATCH 1/8] Fix bending magnet flux and sampling

---
 mcxtrace-comps/sources/Bending_magnet.comp | 367 +++++++++++++--------
 1 file changed, 228 insertions(+), 139 deletions(-)

diff --git a/mcxtrace-comps/sources/Bending_magnet.comp b/mcxtrace-comps/sources/Bending_magnet.comp
index 8be92132bc..69cc2d7165 100644
--- a/mcxtrace-comps/sources/Bending_magnet.comp
+++ b/mcxtrace-comps/sources/Bending_magnet.comp
@@ -13,6 +13,7 @@
 * Date: May, 2013.
 * Version: 1.0
 * Origin: DTU Physics
+* Modified by: Antoine Padovani, 25th Oct 2026 : fix flux and sampling
 *
 * Model of a bending magnet source
 * 
@@ -52,38 +53,62 @@ SETTING PARAMETERS (E0=0, dE=0, lambda0=0,dlambda=0, phase=0, randomphase=1, Ee=
 
 SHARE
 %{
-  //  %include "read_table-lib"
-
-  #ifndef MCCODE_BESSELKNU
-  #define MCCODE_BESSELKNU 1
-
-  #pragma acc routine seq
-  double
-  besselKnu (double nu, double x) {
-    const double h = 0.5;
-    double KK = 0, dK;
-    int r = 0;
-    const int maxiter = 1000;
-    KK = exp (-x) / 2.0;
-    dK = 1;
-    while (dK > DBL_EPSILON && r < maxiter) {
+//  %include "read_table-lib"
+
+#ifndef MCCODE_BESSELKNU
+#define MCCODE_BESSELKNU 1
+
+#pragma acc routine seq
+double besselKnu(double nu, double x){
+    const double h=0.5;
+    double KK=0,dK;
+    int r=0;
+    const int maxiter=1000;
+    KK=exp(-x)/2.0;
+    dK=1;
+    while (dK>DBL_EPSILON && r<maxiter){
       r++;
-      dK = exp (-x * cosh (r * h)) * cosh (nu * r * h);
-      KK += dK;
+      dK=exp(-x*cosh(r*h))*cosh(nu*r*h);
+      KK+=dK;
+    }
+#ifndef OPENACC
+    if (r>=maxiter) {
+      fprintf(stderr,"Warning: Maximum number of iterations exceeded in besselKnu(%g,%g).\n",nu,x);
     }
-    #ifndef OPENACC
-    if (r >= maxiter) {
-      fprintf (stderr, "Warning: Maximum number of iterations exceeded in besselKnu(%g,%g).\n", nu, x);
+#endif
+    KK*=h;
+    return KK;
+  }
+#endif /*MCCODE_BESSELKNU*/
+
+#ifndef M_SQRT1_2
+#define M_SQRT1_2 0.70710678118654752440
+#endif
+
+#ifndef MCCODE_INTEGBESSELKNU
+#define MCCODE_INTEGBESSELKNU 1
+
+  double IntegbesselKnu(double infbound){
+    //fprintf(stderr,"infbound (%g):.\n",infbound);
+    const double h=0.01;
+    double KK=0,dK;
+    int r=0;
+    const int maxiter=1000;
+    KK=0;
+    dK=1;
+    while (dK>DBL_EPSILON && r<maxiter){
+      r++;
+      dK=besselKnu(1.666666666666666666666666667,infbound+r*h);
+      KK+=dK;
     }
-    #endif
-    KK *= h;
+    //if (r>=maxiter) {
+    //  fprintf(stderr,"Warning (%s): Maximum number of iterations exceeded in IntegbesselKnu(%g).\n",NAME_CURRENT_COMP,infbound);
+    //}
+    KK*=h;
     return KK;
   }
-  #endif /*MCCODE_BESSELKNU*/
+#endif /*MCCODE_INTEGBESSELKNU*/
 
-  #ifndef M_SQRT1_2
-  #define M_SQRT1_2 0.70710678118654752440
-  #endif
 %}
 
 DECLARE
@@ -95,6 +120,16 @@ DECLARE
   double s1x;
   double s1y; /*beam's size at dist (convolution of sigex/sigey and igamma)*/
   double p0;
+
+  /*fine structure constant from CODATA*/
+  //const double alpha=7.2973525698e-3;
+  //const double MELE=9.10938291e-31;
+  
+  double ctemult;
+  double muc;
+  double cnstc;
+  double qe;    
+  
 %}
 
 
@@ -103,165 +138,219 @@ INITIALIZE
 
   // fprintf(stderr,"Warning (%s): Bending_magnet is an experimental component - testing is ongoing\n",NAME_CURRENT_COMP);
 
-  if (B <= 0 || Ee <= 0 || Ie <= 0) {
-    fprintf (stderr, "Error (%s): B, Ee, and Ie must all be >= 0. Found (%g %g %g). Aborting.\n", NAME_CURRENT_COMP, B, Ee, Ie);
-    exit (1);
+  if(B<=0 || Ee<=0 || Ie<=0 ){
+    fprintf(stderr, "Error (%s): B, Ee, and Ie must all be >= 0. Found (%g %g %g). Aborting.\n",NAME_CURRENT_COMP,B,Ee,Ie);
+    exit(1);
   }
 
-  if (sigex < 0 || sigey < 0) {
-    fprintf (stderr, "Error (%s): sigex and sigey must be > 0. Negative beam size isn't meaningful. Aborting.\n", NAME_CURRENT_COMP);
-    exit (1);
+  if (sigex <0 || sigey<0){
+    fprintf(stderr, "Error (%s): sigex and sigey must be > 0. Negative beam size isn't meaningful. Aborting.\n",NAME_CURRENT_COMP);
+    exit(1);
   }
-  if (dist <= 0) {
-    fprintf (stderr, "Error (%s): Target undefined.\n", NAME_CURRENT_COMP);
-    exit (1);
+  if (dist<=0){
+    fprintf(stderr,"Error (%s): Target undefined.\n",NAME_CURRENT_COMP);
+    exit(1);
   }
 
   /*compute gamma*/
-  gamma = (Ee * 1e9) / (MELECTRON / CELE * M_C * M_C); /*the extra CELE is to convert to eV*/
-  gamma2 = gamma * gamma;
-  igamma = 1.0 / gamma;
+  gamma=(Ee*1e9)/(MELECTRON/CELE*M_C*M_C);/*the extra CELE is to convert to eV*/
+  gamma2=gamma*gamma;
+  igamma=1.0/gamma;
+
+  ctemult=(3.0*ALPHA/(4.0*M_PI*M_PI))*(gamma2)*(Ie/CELE);
 
-  // printf("Bending_magnet (%s): gamma=%g, divergence is 1/gamma=%g rad.\n",NAME_CURRENT_COMP,gamma,igamma);
+  //printf("Bending_magnet (%s): gamma=%g, divergence is 1/gamma=%g rad.\n",NAME_CURRENT_COMP,gamma,igamma);
   /*compute characteristic energy in keV*/
-  double Ec = 0.665 * Ee * Ee * B;
-  // double Ec=1.5*gamma2*HBAR*CELE*B/MELECTRON *1e-3; /*check units on this one. The 1e-3 factor is because energy is assumed to be in keV*/
+  double Ec=0.665*Ee*Ee*B;
+  //double Ec=1.5*gamma2*HBAR*CELE*B/MELECTRON *1e-3; /*check units on this one. The 1e-3 factor is because energy is assumed to be in keV*/
   /*We normally do computations in k so use that for transfer*/
-  kc = E2K * Ec;
+  kc=E2K*Ec;
+
+  muc = kc*M_C/(2.0*M_PI);
+
+  //s1x=sqrt(sigex*sigex + igamma*igamma*dist*dist);
+  //s1y=sqrt(sigey*sigey + igamma*igamma*dist*dist);
+  if (focus_xw){
+  s1x=focus_xw/2.0;
+  }
+  if (focus_yh){
+  s1y=focus_yh/2.0;  
+  }
+  p0=1.0/mcget_ncount();
 
-  s1x = sqrt (sigex * sigex + igamma * igamma * dist * dist);
-  s1y = sqrt (sigey * sigey + igamma * igamma * dist * dist);
-  p0 = 1.0 / mcget_ncount ();
+  cnstc = erf(1.0/(sqrt(2.0)));
+  qe = 1.0/(2.0*dE);   
 %}
 
 
 TRACE
 %{
 
-  double xx, yy, x1, y1, z1;
-  double k, e, l;
-  double F1 = 1.0;
-  double dx, dy, dz;
+  double xx,yy,x1,y1,z1;
+  double k,e,l;
+  double F1=1.0;
+  double dx,dy,dz;
+  double dl;
 
+  double angle, psi, theta, xcap, oneplusxsquared, chi, K23, K13;
+  double qtheta;
+  double qpsi;
+    
   // initial source area
-  xx = randnorm ();
-  yy = randnorm ();
-  x = xx * sigex;
-  y = yy * sigey;
-  z = 0;
+  xx=randnorm();
+  yy=randnorm();
+  x=xx*sigex;
+  y=yy*sigey;
+  z=0;
 
   // Gaussian distribution at origin
-  p = p0; /*initial weight is p0*/
-  if (E0) {
-    if (!dE) {
-      e = E0;
-    } else {
-      e = randpm1 () * dE + E0;
+  p=p0;/*initial weight is p0*/
+  if (E0){
+    if(!dE){
+      e=E0;
+    }else {
+      //e=randpm1()*dE + E0;
+      e=randpm1()*dE + E0;
     }
-    k = E2K * e;
-  } else if (lambda0) {
-    if (!dlambda) {
-      l = lambda0;
-    } else {
-      l = randpm1 () * dlambda + lambda0;
+    k=E2K*e;
+  }else if (lambda0){
+    if (!dlambda){
+      l=lambda0;
+    }else{
+      l=randpm1()*dlambda + lambda0;
     }
-    k = (2 * M_PI / l);
+    k=(2*M_PI/l);
   }
+  
+  double oneoverk=1.0/k;
+  double mu = k*M_C/(2*M_PI);
+  double mumuc = mu/muc;
+      
   // targeted area calculation
-  if (focus_xw) {
-    if (!gauss_t) {
+  if (focus_xw){
+    if (!gauss_t){
       /*sample uniformly but adjust weight*/
-      x1 = randpm1 () * focus_xw / 2.0;
-      p *= exp (-(x1 * x1) / (2.0 * s1x * s1x));
-    } else {
+      x1=randpm1()*focus_xw/2.0;
+      theta = atan((x1-x)/dist);
+      p*=focus_xw; //sampling probability density is 1.0/focus_xw
+      p*= dist/(dist*dist + (x1-x)*(x1-x)); //take into account the fact that we are sampling a length, not an angle (and our flux is angular). So need to apply that correction also.
+    }else {
       do {
-        x1 = randnorm () * s1x;
-      } while (focus_xw != 0 && fabs (x1) > focus_xw / 2.0);
-      /*adjust for restricted sampling window*/
-      p *= erf (focus_xw * 0.5 * M_SQRT1_2 / s1x);
+        x1=randnorm()*s1x;
+      }while (focus_xw!=0 && fabs(x1)>focus_xw/2.0);
+      theta = atan((x1-x)/dist);
+      /*adjust for restricted sampling window*/            
+      p*=erf(focus_xw*0.5*M_SQRT1_2/s1x)*s1x*sqrt(2.0*M_PI)*exp((x1*x1)/(2*s1x*s1x)); // 1/sampling probability density
+      p*= dist/(dist*dist + (x1-x)*(x1-x));
     }
-  } else {
-    x1 = randnorm () * igamma;
+  }else{    
+    angle=2.5e-3; //angle is the angle psi and theta will be generated over. Range of (-angle,angle).
+    do {
+      theta = randnorm()*angle;
+    } while(theta>angle || theta<-angle);  
+
+    qtheta = 1.0/(sqrt(2*M_PI)*angle*cnstc)*exp(-(theta*theta)/(2*angle*angle)); //sampling probability density
+    p*=1.0/qtheta;
+    
+    x1=x+dist*tan(theta);
   }
-  if (focus_yh) {
-    if (!gauss_t) {
+  if (focus_yh){
+    if (!gauss_t){
       /*sample uniformly but adjust weight*/
-      y1 = randpm1 () * focus_yh / 2.0;
-      p *= exp (-(y1 * y1) / (2.0 * s1y * s1y));
-    } else {
+      y1=randpm1()*focus_yh/2.0;
+      psi = atan((y1-y)/dist);
+      p*=focus_yh; //sampling probability density is 1.0/focus_yh
+      p*= dist/(dist*dist + (y1-y)*(y1-y));
+    }else {
       do {
-        y1 = randnorm () * s1y;
-      } while (fabs (y1) > focus_yh / 2.0);
-      /*adjust for restricted sampling window*/
-      p *= erf (focus_yh * 0.5 * M_SQRT1_2 / s1y);
+        y1=randnorm()*s1y;
+      }while (fabs(y1)>focus_yh/2.0);
+      psi = atan((y1-y)/dist);
+      /*adjust for restricted sampling window*/      
+      p*=erf(focus_yh*0.5*M_SQRT1_2/s1y)*s1y*sqrt(2.0*M_PI)*exp((y1*y1)/(2*s1y*s1y)); // 1/sampling probability density
+      p*= dist/(dist*dist + (y1-y)*(y1-y));
     }
-  } else {
-    y1 = randnorm () * igamma;
+  }else{      
+    angle=2.5e-3; //angle is the angle psi and theta will be generated over. Range of (-angle,angle).
+    do {
+      psi = randnorm()*angle;
+    } while(psi>angle || psi<-angle);      
+    
+    qpsi = 1.0/(sqrt(2*M_PI)*angle*cnstc)*exp(-(psi*psi)/(2*angle*angle)); //sampling probability density
+    p*=1.0/qpsi;
+    
+    y1=y+dist*tan(psi);
   }
-  z1 = dist;
-  dx = x1 - x;
-  dy = y1 - y;
-  dz = sqrt (dx * dx + dy * dy + dist * dist);
+  
+  xcap = gamma*psi;
+  oneplusxsquared = 1+xcap*xcap;
 
-  kx = (k * dx) / dz;
-  ky = (k * dy) / dz;
-  kz = (k * dist) / dz;
+  chi = mumuc*pow(oneplusxsquared,1.5)*(1.0/2.0);
+  
+  K23 = besselKnu(0.666666666666666666666666667,chi);
+  K13 = besselKnu(0.333333333333333333333333334,chi);
+  
+  z1=z+dist;
+  dx=x1-x;
+  dy=y1-y;  
+  dz=z1-z;
+    
+  dl = sqrt(dx*dx+dy*dy+dz*dz);
+    
+  kx=(k*dx)/dl;
+  ky=(k*dy)/dl;
+  kz=(k*dz)/dl;  
+  
+  p*=ctemult*mumuc*mumuc*(oneplusxsquared*oneplusxsquared)*(K23*K23+(xcap*xcap/(oneplusxsquared))*K13*K13);
 
-  /*spectral strength of radiation is given by Patterson*/
-  double k_kc = k / kc;
-  double K2_3 = besselKnu (0.666666666666666666666666667, k_kc * 0.5);
-  p *= 1.33e13 * Ee * Ee * Ie * k_kc * k_kc * K2_3 * K2_3;
-  // p*=ALPHA/(M_PI*M_PI)*gamma2*Ie/CELE* 1e-4 * 0.75 *k_kc*k_kc*K2_3*K2_3;
+  p*= 1.0/(qe*e);  // energy sampling density 1.0/qe and 1/e coming from the physical flux
 
   /*randomly pick phase*/
-  if (randomphase) {
-    phi = rand01 () * 2 * M_PI;
-  } else {
-    phi = phase;
+  if (randomphase){
+    phi=rand01()*2*M_PI;
+  }else{
+    phi=phase;
   }
 
   /*set polarization vector*/
-  Ex = 0;
-  Ey = 0;
-  Ez = 0;
+  Ex=0;Ey=0;Ez=0;
 %}
 
 MCDISPLAY
 %{
+  
 
-
-  double radius, D;
-  radius = 3.3 * Ee / B;
-  ;
-  D = 0.1;
-  double x0, x1, z0, z1;
-  const double phimin = -2.0 * DEG2RAD, phimax = 2.0 * DEG2RAD;
-  double phi = phimin, dphi;
-  dphi = (phimax - phimin) / 32;
-  while (phi < phimax) {
-    x0 = radius * (1 - cos (phi));
-    x1 = radius * (1 - cos (phi + dphi));
-    z0 = radius * sin (phi);
-    z1 = radius * sin (phi + dphi);
-    line (x0, 0.0, z0, x1, 0.0, z1);
-    phi += dphi;
+  double radius,D;
+  radius=3.3*Ee/B;;
+  D=0.1;
+  double x0,x1,z0,z1;
+  const double phimin=-2.0*DEG2RAD, phimax=2.0*DEG2RAD;
+  double phi=phimin,dphi;  
+  dphi=(phimax-phimin)/32;
+  while (phi<phimax){
+    x0=radius*(1-cos(phi));
+    x1=radius*(1-cos(phi+dphi));
+    z0=radius*sin(phi);
+    z1=radius*sin(phi+dphi);
+    line(x0,0.0,z0,x1,0.0,z1);
+    phi+=dphi;
   }
 
-  line (0.0, 0.0, 0.0, D* sin (igamma), 0.0, D);
-  line (0.0, 0.0, 0.0, -D* sin (igamma), 0.0, D);
-  line (0.0, 0.0, 0.0, 0.0, D* sin (igamma), D);
-  line (0.0, 0.0, 0.0, 0.0, -D* sin (igamma), D);
-
-  phi = -igamma;
-  dphi = 2.0 * igamma / 32;
-  while (phi < igamma) {
-    x0 = D * sin (phi);
-    x1 = D * sin (phi + dphi);
-    z0 = D * cos (phi);
-    z1 = D * cos (phi + dphi);
-    line (x0, 0.0, z0, x1, 0.0, z1);
-    line (0.0, x0, z0, 0.0, x1, z1);
-    phi += dphi;
+  line(0.0,0.0,0.0, D*sin(igamma), 0.0, D);
+  line(0.0,0.0,0.0,-D*sin(igamma), 0.0, D);
+  line(0.0,0.0,0.0, 0.0, D*sin(igamma), D);
+  line(0.0,0.0,0.0, 0.0,-D*sin(igamma), D);
+  
+  phi =-igamma;
+  dphi= 2.0*igamma/32;
+  while(phi<igamma){
+    x0=D*sin(phi);
+    x1=D*sin(phi+dphi);
+    z0=D*cos(phi);
+    z1=D*cos(phi+dphi);
+    line(x0,0.0,z0,x1,0.0,z1);
+    line(0.0,x0,z0,0.0,x1,z1);
+    phi+=dphi;
   }
 %}
 

From 6db1bbba77a554d0a813a1ea856961037be4e75c Mon Sep 17 00:00:00 2001
From: Emmanuel FARHI <emmanuel.farhi@synchrotron-soleil.fr>
Date: Mon, 26 Oct 2026 21:35:49 +0100
Subject: [PATCH 2/8] Bending_magnet: E0,dE option only. More tests.

I suggest a few changes:
- The 'lamdba' branch is broken: I removed it.
- I added a test over dE as it is needed for `qe` init.
- The Bessel integral function is not used. Removed.
---
 mcxtrace-comps/sources/Bending_magnet.comp | 74 +++++-----------------
 1 file changed, 17 insertions(+), 57 deletions(-)

diff --git a/mcxtrace-comps/sources/Bending_magnet.comp b/mcxtrace-comps/sources/Bending_magnet.comp
index 69cc2d7165..f5e36c4779 100644
--- a/mcxtrace-comps/sources/Bending_magnet.comp
+++ b/mcxtrace-comps/sources/Bending_magnet.comp
@@ -39,21 +39,18 @@
 * gauss_t:  [0/1] If 0 the target window will be sampled uniformly and the weight adjusted accordingly, otherwise we will use a gaussian sampling scheme
 * E0:       [keV] Center of emitted energy spectrum (overrides lambda0)
 * dE:       [keV] Half-width of emitted energy spectrum
-* lambda0:  [AA]  Center of emitted wavelength spectrum
-* dlambda:  [AA]  Half-width of emitted wavelength spectrum
 * 
 * %End
 **********************/
 
 DEFINE COMPONENT Bending_magnet
 
-SETTING PARAMETERS (E0=0, dE=0, lambda0=0,dlambda=0, phase=0, randomphase=1, Ee=2.4, Ie=0.4, B=1.6, sigey=0, sigex=0, focus_xw=0, focus_yh=0, dist=1, gauss_t=0)
+SETTING PARAMETERS (E0=0, dE=0, phase=0, randomphase=1, Ee=2.4, Ie=0.4, B=1.6, sigey=0, sigex=0, focus_xw=0, focus_yh=0, dist=1, gauss_t=0)
 
-/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */ 
+/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */
 
 SHARE
 %{
-//  %include "read_table-lib"
 
 #ifndef MCCODE_BESSELKNU
 #define MCCODE_BESSELKNU 1
@@ -85,30 +82,6 @@ double besselKnu(double nu, double x){
 #define M_SQRT1_2 0.70710678118654752440
 #endif
 
-#ifndef MCCODE_INTEGBESSELKNU
-#define MCCODE_INTEGBESSELKNU 1
-
-  double IntegbesselKnu(double infbound){
-    //fprintf(stderr,"infbound (%g):.\n",infbound);
-    const double h=0.01;
-    double KK=0,dK;
-    int r=0;
-    const int maxiter=1000;
-    KK=0;
-    dK=1;
-    while (dK>DBL_EPSILON && r<maxiter){
-      r++;
-      dK=besselKnu(1.666666666666666666666666667,infbound+r*h);
-      KK+=dK;
-    }
-    //if (r>=maxiter) {
-    //  fprintf(stderr,"Warning (%s): Maximum number of iterations exceeded in IntegbesselKnu(%g).\n",NAME_CURRENT_COMP,infbound);
-    //}
-    KK*=h;
-    return KK;
-  }
-#endif /*MCCODE_INTEGBESSELKNU*/
-
 %}
 
 DECLARE
@@ -136,13 +109,14 @@ DECLARE
 INITIALIZE
 %{
 
-  // fprintf(stderr,"Warning (%s): Bending_magnet is an experimental component - testing is ongoing\n",NAME_CURRENT_COMP);
-
   if(B<=0 || Ee<=0 || Ie<=0 ){
-    fprintf(stderr, "Error (%s): B, Ee, and Ie must all be >= 0. Found (%g %g %g). Aborting.\n",NAME_CURRENT_COMP,B,Ee,Ie);
+    fprintf(stderr, "Error (%s): B, Ee, and Ie must all be > 0. Found (%g %g %g). Aborting.\n",NAME_CURRENT_COMP,B,Ee,Ie);
     exit(1);
   }
-
+  if (E0 <= 0 || dE <= 0) {
+    fprintf (stderr, "Error (%s): dE and E0 need to be set. Aborting.\n", NAME_CURRENT_COMP);
+    exit (1);
+  }
   if (sigex <0 || sigey<0){
     fprintf(stderr, "Error (%s): sigex and sigey must be > 0. Negative beam size isn't meaningful. Aborting.\n",NAME_CURRENT_COMP);
     exit(1);
@@ -171,15 +145,15 @@ INITIALIZE
   //s1x=sqrt(sigex*sigex + igamma*igamma*dist*dist);
   //s1y=sqrt(sigey*sigey + igamma*igamma*dist*dist);
   if (focus_xw){
-  s1x=focus_xw/2.0;
+    s1x=focus_xw/2.0;
   }
   if (focus_yh){
-  s1y=focus_yh/2.0;  
+    s1y=focus_yh/2.0;  
   }
   p0=1.0/mcget_ncount();
 
   cnstc = erf(1.0/(sqrt(2.0)));
-  qe = 1.0/(2.0*dE);   
+  qe = 1.0/(2.0*dE);    // requires dE > 0, no lambda handling here
 %}
 
 
@@ -205,27 +179,14 @@ TRACE
 
   // Gaussian distribution at origin
   p=p0;/*initial weight is p0*/
-  if (E0){
-    if(!dE){
-      e=E0;
-    }else {
-      //e=randpm1()*dE + E0;
-      e=randpm1()*dE + E0;
-    }
-    k=E2K*e;
-  }else if (lambda0){
-    if (!dlambda){
-      l=lambda0;
-    }else{
-      l=randpm1()*dlambda + lambda0;
-    }
-    k=(2*M_PI/l);
-  }
+  // E0 and dE are defined
+  e=randpm1()*dE + E0;
+  k=E2K*e;
   
   double oneoverk=1.0/k;
   double mu = k*M_C/(2*M_PI);
   double mumuc = mu/muc;
-      
+
   // targeted area calculation
   if (focus_xw){
     if (!gauss_t){
@@ -313,15 +274,14 @@ TRACE
   }
 
   /*set polarization vector*/
-  Ex=0;Ey=0;Ez=0;
+  Ex=Ey=Ez=0;
 %}
 
 MCDISPLAY
 %{
-  
 
   double radius,D;
-  radius=3.3*Ee/B;;
+  radius=3.3*Ee/B;
   D=0.1;
   double x0,x1,z0,z1;
   const double phimin=-2.0*DEG2RAD, phimax=2.0*DEG2RAD;
@@ -340,7 +300,7 @@ MCDISPLAY
   line(0.0,0.0,0.0,-D*sin(igamma), 0.0, D);
   line(0.0,0.0,0.0, 0.0, D*sin(igamma), D);
   line(0.0,0.0,0.0, 0.0,-D*sin(igamma), D);
-  
+
   phi =-igamma;
   dphi= 2.0*igamma/32;
   while(phi<igamma){

From 71fdf3c71894dac0d57af1b3dd7f535da597bb7f Mon Sep 17 00:00:00 2001
From: Emmanuel FARHI <emmanuel.farhi@synchrotron-soleil.fr>
Date: Tue, 27 Oct 2026 22:08:52 +0100
Subject: [PATCH 3/8] Bending_magnet: restore 'lambda' mode

The lambda mode is restored.
Checks for lambda0-dlambda> 0 and E0-dE>0
---
 mcxtrace-comps/sources/Bending_magnet.comp | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/mcxtrace-comps/sources/Bending_magnet.comp b/mcxtrace-comps/sources/Bending_magnet.comp
index f5e36c4779..016b5c0a77 100644
--- a/mcxtrace-comps/sources/Bending_magnet.comp
+++ b/mcxtrace-comps/sources/Bending_magnet.comp
@@ -39,13 +39,16 @@
 * gauss_t:  [0/1] If 0 the target window will be sampled uniformly and the weight adjusted accordingly, otherwise we will use a gaussian sampling scheme
 * E0:       [keV] Center of emitted energy spectrum (overrides lambda0)
 * dE:       [keV] Half-width of emitted energy spectrum
+* lambda0:  [AA]  Center of emitted wavelength spectrum
+* dlambda:  [AA]  Half-width of emitted wavelength spectrum
 * 
 * %End
 **********************/
 
 DEFINE COMPONENT Bending_magnet
 
-SETTING PARAMETERS (E0=0, dE=0, phase=0, randomphase=1, Ee=2.4, Ie=0.4, B=1.6, sigey=0, sigex=0, focus_xw=0, focus_yh=0, dist=1, gauss_t=0)
+SETTING PARAMETERS (E0=0, dE=0, lambda0=0,dlambda=0, 
+  phase=0, randomphase=1, Ee=2.4, Ie=0.4, B=1.6, sigey=0, sigex=0, focus_xw=0, focus_yh=0, dist=1, gauss_t=0)
 
 /* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */
 
@@ -113,7 +116,14 @@ INITIALIZE
     fprintf(stderr, "Error (%s): B, Ee, and Ie must all be > 0. Found (%g %g %g). Aborting.\n",NAME_CURRENT_COMP,B,Ee,Ie);
     exit(1);
   }
-  if (E0 <= 0 || dE <= 0) {
+  if (lambda0 > 0 && dlambda > 0 && lambda0-dlambda > 0) {
+    // compute equivalent E0 and dE
+    double E1 = (2 * M_PI / (lambda0+dlambda)) * K2E;
+    double E2 = (2 * M_PI / (lambda0-dlambda)) * K2E;
+    E0 = fabs(E1 - E2);
+    dE = (E1 + E2)/2;
+  }
+  if (E0 <= 0 || dE <= 0 || E0-dE <= 0) {
     fprintf (stderr, "Error (%s): dE and E0 need to be set. Aborting.\n", NAME_CURRENT_COMP);
     exit (1);
   }

From 309a706fca123a0578af8f61c82f76268779057a Mon Sep 17 00:00:00 2001
From: Emmanuel FARHI <emmanuel.farhi@synchrotron-soleil.fr>
Date: Wed, 28 Oct 2026 07:28:21 +0100
Subject: [PATCH 4/8] Bending_magnet: fix lambda0 -> E0

---
 mcxtrace-comps/sources/Bending_magnet.comp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mcxtrace-comps/sources/Bending_magnet.comp b/mcxtrace-comps/sources/Bending_magnet.comp
index 016b5c0a77..1f414162a5 100644
--- a/mcxtrace-comps/sources/Bending_magnet.comp
+++ b/mcxtrace-comps/sources/Bending_magnet.comp
@@ -120,8 +120,8 @@ INITIALIZE
     // compute equivalent E0 and dE
     double E1 = (2 * M_PI / (lambda0+dlambda)) * K2E;
     double E2 = (2 * M_PI / (lambda0-dlambda)) * K2E;
-    E0 = fabs(E1 - E2);
-    dE = (E1 + E2)/2;
+    E0 =     (E1 + E2)/2;
+    dE = fabs(E1 - E2)/2;
   }
   if (E0 <= 0 || dE <= 0 || E0-dE <= 0) {
     fprintf (stderr, "Error (%s): dE and E0 need to be set. Aborting.\n", NAME_CURRENT_COMP);

From 1c7352fff0c00b51c0237fefddae6cb32a3b74cc Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Thu, 29 Oct 2026 18:25:15 +0100
Subject: [PATCH 5/8] McXtrace examples: update %Example references for new
 Bending_magnet flux

The fixed Bending_magnet now returns photons/s integrated over the sampled
solid angle and energy band, so reference intensities change by design.
Values from mxrun (mccode-dev, macOS arm64, CPU) with this PR's component,
1e6-1e7 rays. Test_BM SOURCE=1 (Wiggler) and the %Scan references are
unaffected and left as-is.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr  | 2 +-
 mcxtrace-comps/examples/BARC/Czerny_Turner/README.md            | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md           | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr  | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/README.md            | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/SOLEIL_MARS.instr    | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/README.md            | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/SOLEIL_ROCK.instr    | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/README.md           | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/SOLEIL_SAMBA.instr  | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/README.md          | 2 +-
 .../examples/SOLEIL/SOLEIL_SIRIUS/SOLEIL_SIRIUS.instr           | 2 +-
 .../examples/Tests_optics/Test_grating_reflect/README.md        | 2 +-
 .../Test_grating_reflect/Test_grating_reflect.instr             | 2 +-
 mcxtrace-comps/examples/Tests_sources/Test_BM/README.md         | 2 +-
 mcxtrace-comps/examples/Tests_sources/Test_BM/Test_BM.instr     | 2 +-
 16 files changed, 16 insertions(+), 16 deletions(-)

diff --git a/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr b/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr
index 5ceec336f9..d232538c20 100644
--- a/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr
+++ b/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr
@@ -19,7 +19,7 @@
 * It can be found here: https://inis.iaea.org/collection/NCLCollectionStore/_Public/19/019/19019134.pdf
 * Example: Do a scan (scan parameter =1) from x_screw=20 to 103 mm. The calculated monitor "Wavelength(Ang) as a function of x_screw(mm)" will be a linear function. The monochromator functions properly.
 *
-* %Example: x_screw=60 Detector: w_monitor_I=1.10881e-06
+* %Example: x_screw=60 Detector: w_monitor_I=0.0723
 * %Scan: mxrun Czerny_Turner.instr -N21 x_screw=20,103 scan=1 -n1e6 Detector: Wavelength_I={
 *   2007.08,2406.85,2820.95,3242.17,3659.02,4075.98,4485.76,4901.23,
 *   5324.84,5733.95,6156.24,6569.68,6979.66,7395.58,7805.61,8222.2,
diff --git a/mcxtrace-comps/examples/BARC/Czerny_Turner/README.md b/mcxtrace-comps/examples/BARC/Czerny_Turner/README.md
index fce75b8183..2a5b25cf87 100644
--- a/mcxtrace-comps/examples/BARC/Czerny_Turner/README.md
+++ b/mcxtrace-comps/examples/BARC/Czerny_Turner/README.md
@@ -22,7 +22,7 @@ Example: Do a scan (scan parameter =1) from x_screw=20 to 103 mm. The calculated
 
 ## Examples
 
-- **Test: x_screw=60 Detector: w_monitor_I=1.10881e-06**
+- **Test: x_screw=60 Detector: w_monitor_I=0.0723**
 - **TestScan: mxrun Czerny_Turner.instr -N21 x_screw=20,103 scan=1 -n1e6 Detector: Wavelength_I={21 values}**
 
 ## Input parameters
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md
index 5eed2352af..f08dca5bb5 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md
@@ -20,7 +20,7 @@ You may DISCO_scan the grating monochromator with e.g.:
 
 - **TestScan: mxrun SOLEIL_DISCO.instr -N21 DISCO_x_screw=20,103 DISCO_scan=1 Detector: Wavelength_I={21 values}**
 - 
-- **Test: DISCO_grazing_angle_one=22.5 -n 1e5 Detector: DISCO_second_HFP_I=2.08292e+11**
+- **Test: DISCO_grazing_angle_one=22.5 -n 1e5 Detector: DISCO_second_HFP_I=1.448e+16**
 
 ## Input parameters
 
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr
index 2174c71802..28a3884d9a 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr
@@ -18,7 +18,7 @@
 *   5306.77,5738.39,6151.91,6565.99,6979.47,7395.13,7786.17,8199.52,
 *   8631.53,9051.46,9483.21,9913.98,10299.5 }
 *
-* %Example: DISCO_grazing_angle_one=22.5 -n 1e5 Detector: DISCO_second_HFP_I=2.08292e+11
+* %Example: DISCO_grazing_angle_one=22.5 -n 1e5 Detector: DISCO_second_HFP_I=1.448e+16
 *
 * %Parameters
 * DISCO_grazing_angle_one: [deg] Mirrors M1/M2 rotation angle
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/README.md
index a7eb00befa..24bb970fd8 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/README.md
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/README.md
@@ -28,7 +28,7 @@ This model implements the XRD station for powders.
 
 ## Examples
 
-- **Test: MARS_E0=16.99 Detector: MARS_detector_diffraction_I=2.56e+09**
+- **Test: MARS_E0=16.99 Detector: MARS_detector_diffraction_I=1.334e+08**
 
 ## Input parameters
 
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/SOLEIL_MARS.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/SOLEIL_MARS.instr
index c422ef5893..ebbdefccb1 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/SOLEIL_MARS.instr
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/SOLEIL_MARS.instr
@@ -24,7 +24,7 @@
 *
 * This model implements the XRD station for powders.
 *
-* %Example: MARS_E0=16.99 Detector: MARS_detector_diffraction_I=2.56e+09
+* %Example: MARS_E0=16.99 Detector: MARS_detector_diffraction_I=1.334e+08
 *
 * %Parameters
 * MARS_E0: [keV]     Central energy to be emitted by the source
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/README.md
index 276519732e..c7bf570a60 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/README.md
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/README.md
@@ -48,7 +48,7 @@ Pretty energy repartition monitor result: mxrun SOLEIL_ROCK.instr ROCK_E0=17.000
 
 ## Examples
 
-- **Test: ROCK_E0=15.918 Detector: ROCK_fluo_monitor_I=2e+09**
+- **Test: ROCK_E0=15.918 Detector: ROCK_fluo_monitor_I=1.29e+11**
 
 ## Input parameters
 
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/SOLEIL_ROCK.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/SOLEIL_ROCK.instr
index 7faf684e6c..398adb2de3 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/SOLEIL_ROCK.instr
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/SOLEIL_ROCK.instr
@@ -44,7 +44,7 @@
 * Manganese and chrome ROCK_scan: mxrun SOLEIL_ROCK.instr ROCK_E0=5.700,6.800 ROCK_scan=1 ROCK_sample_file=MnCr ROCK_cc=2 -N100
 * Pretty energy repartition monitor result: mxrun SOLEIL_ROCK.instr ROCK_E0=17.000 ROCK_cc=2 -n1e8
 *
-* %Example: ROCK_E0=15.918 Detector: ROCK_fluo_monitor_I=2e+09
+* %Example: ROCK_E0=15.918 Detector: ROCK_fluo_monitor_I=1.29e+11
 *
 * %Parameters
 * ROCK_E0:               [keV]    Energy to hit, i.e. selected by the channel-cut monochromator
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/README.md
index 0e88d31f6f..ca94ec984f 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/README.md
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/README.md
@@ -34,7 +34,7 @@ Position | Element
 
 ## Examples
 
-- **Test: SAMBA_E0=15 Detector: SAMBA_det_fluo_I=1.23968e+10**
+- **Test: SAMBA_E0=15 Detector: SAMBA_det_fluo_I=5.92e+10**
 
 ## Input parameters
 
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/SOLEIL_SAMBA.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/SOLEIL_SAMBA.instr
index a179c44abb..fe6e9b84d5 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/SOLEIL_SAMBA.instr
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/SOLEIL_SAMBA.instr
@@ -34,7 +34,7 @@
 * 30       | Sample Stage (XAS configuration)
 * 35       | Sample Stage (SEXAFS configuration)
 *
-* %Example: SAMBA_E0=15 Detector: SAMBA_det_fluo_I=1.23968e+10
+* %Example: SAMBA_E0=15 Detector: SAMBA_det_fluo_I=5.92e+10
 *
 * %Parameters
 * SAMBA_E0:               [keV]  Central energy
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/README.md
index 008b20f8b3..faa6f108fe 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/README.md
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/README.md
@@ -30,7 +30,7 @@ Position | Element
 
 ## Examples
 
-- **Test: SIRIUS_E0=8 Detector: SIRIUS_def_loc_I=2.14522e+10**
+- **Test: SIRIUS_E0=8 Detector: SIRIUS_def_loc_I=6.51e+09**
 
 ## Input parameters
 
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/SOLEIL_SIRIUS.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/SOLEIL_SIRIUS.instr
index fcbb84bbd0..01d5e14c9d 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/SOLEIL_SIRIUS.instr
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/SOLEIL_SIRIUS.instr
@@ -28,7 +28,7 @@
 * 30       | Sample Stage (XAS configuration)
 * 35       | Sample Stage (SEXAFS configuration)
 *
-* %Example: SIRIUS_E0=8 Detector: SIRIUS_def_loc_I=2.14522e+10
+* %Example: SIRIUS_E0=8 Detector: SIRIUS_def_loc_I=6.51e+09
 *
 * %Parameters
 * SIRIUS_E0:               [keV]  Central energy (1.4 - 13 keV)
diff --git a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md
index 29095ae97c..330c0ef250 100644
--- a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md
+++ b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md
@@ -17,7 +17,7 @@ Simply a bending magnet illuminating a reflection grating.
 
 ## Examples
 
-- **Test: E0=1 Detector: psd_monitor_I=3.79378e+06**
+- **Test: E0=1 Detector: psd_monitor_I=3.47e+09**
 
 ## Input parameters
 
diff --git a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr
index 685f61520e..00332db8f8 100644
--- a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr
+++ b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr
@@ -15,7 +15,7 @@
 *
 * Simply a bending magnet illuminating a reflection grating.
 *
-* %Example: E0=1 Detector: psd_monitor_I=3.79378e+06
+* %Example: E0=1 Detector: psd_monitor_I=3.47e+09
 *
 * %Parameters
 * E0: [keV] Source's center of emitted energy spectrum
diff --git a/mcxtrace-comps/examples/Tests_sources/Test_BM/README.md b/mcxtrace-comps/examples/Tests_sources/Test_BM/README.md
index 8a662d9e62..c503862be0 100644
--- a/mcxtrace-comps/examples/Tests_sources/Test_BM/README.md
+++ b/mcxtrace-comps/examples/Tests_sources/Test_BM/README.md
@@ -17,7 +17,7 @@ This is a simple test-instrument for the bending magnet component.
 
 ## Examples
 
-- **Test: Test_BM.instr SOURCE=0 e0=12.5 de=0.5 Detector: emon_I=2.93138e+13**
+- **Test: Test_BM.instr SOURCE=0 e0=12.5 de=0.5 Detector: emon_I=9.146e+14**
 - **Test: Test_BM.instr SOURCE=1 e0=12.5 de=0.5 Detector: emon_I=1.28549e+15**
 
 ## Input parameters
diff --git a/mcxtrace-comps/examples/Tests_sources/Test_BM/Test_BM.instr b/mcxtrace-comps/examples/Tests_sources/Test_BM/Test_BM.instr
index 8133248629..9abadf6208 100644
--- a/mcxtrace-comps/examples/Tests_sources/Test_BM/Test_BM.instr
+++ b/mcxtrace-comps/examples/Tests_sources/Test_BM/Test_BM.instr
@@ -16,7 +16,7 @@
 * %Description
 * This is a simple test-instrument for the bending magnet component. 
 *
-* %Example: Test_BM.instr SOURCE=0 e0=12.5 de=0.5 Detector: emon_I=2.93138e+13
+* %Example: Test_BM.instr SOURCE=0 e0=12.5 de=0.5 Detector: emon_I=9.146e+14
 * %Example: Test_BM.instr SOURCE=1 e0=12.5 de=0.5 Detector: emon_I=1.28549e+15
 *
 * %Parameters

From 18107d40192126d9739a774d0be9e3a7d9e6ad29 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Thu, 29 Oct 2026 20:01:26 +0100
Subject: [PATCH 6/8] McXtrace examples: refine Bending_magnet %Example
 references

Use the mean of 5 seeds at 1e7 rays (mxrun, mccode-dev, macOS arm64 CPU)
instead of single runs. Largest corrections: SOLEIL_ROCK 1.29e11 -> 1.61e11,
SOLEIL_SAMBA 5.92e10 -> 5.69e10. ROCK, SAMBA and SIRIUS show seed-to-seed
scatter (7-18% at 1e7) well above the reported per-run error.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr  | 2 +-
 mcxtrace-comps/examples/BARC/Czerny_Turner/README.md            | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md           | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr  | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/README.md            | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/SOLEIL_MARS.instr    | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/README.md            | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/SOLEIL_ROCK.instr    | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/README.md           | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/SOLEIL_SAMBA.instr  | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/README.md          | 2 +-
 .../examples/SOLEIL/SOLEIL_SIRIUS/SOLEIL_SIRIUS.instr           | 2 +-
 .../examples/Tests_optics/Test_grating_reflect/README.md        | 2 +-
 .../Test_grating_reflect/Test_grating_reflect.instr             | 2 +-
 mcxtrace-comps/examples/Tests_sources/Test_BM/README.md         | 2 +-
 mcxtrace-comps/examples/Tests_sources/Test_BM/Test_BM.instr     | 2 +-
 16 files changed, 16 insertions(+), 16 deletions(-)

diff --git a/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr b/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr
index d232538c20..2a606bb777 100644
--- a/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr
+++ b/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr
@@ -19,7 +19,7 @@
 * It can be found here: https://inis.iaea.org/collection/NCLCollectionStore/_Public/19/019/19019134.pdf
 * Example: Do a scan (scan parameter =1) from x_screw=20 to 103 mm. The calculated monitor "Wavelength(Ang) as a function of x_screw(mm)" will be a linear function. The monochromator functions properly.
 *
-* %Example: x_screw=60 Detector: w_monitor_I=0.0723
+* %Example: x_screw=60 Detector: w_monitor_I=0.0721
 * %Scan: mxrun Czerny_Turner.instr -N21 x_screw=20,103 scan=1 -n1e6 Detector: Wavelength_I={
 *   2007.08,2406.85,2820.95,3242.17,3659.02,4075.98,4485.76,4901.23,
 *   5324.84,5733.95,6156.24,6569.68,6979.66,7395.58,7805.61,8222.2,
diff --git a/mcxtrace-comps/examples/BARC/Czerny_Turner/README.md b/mcxtrace-comps/examples/BARC/Czerny_Turner/README.md
index 2a5b25cf87..b8458b2f9c 100644
--- a/mcxtrace-comps/examples/BARC/Czerny_Turner/README.md
+++ b/mcxtrace-comps/examples/BARC/Czerny_Turner/README.md
@@ -22,7 +22,7 @@ Example: Do a scan (scan parameter =1) from x_screw=20 to 103 mm. The calculated
 
 ## Examples
 
-- **Test: x_screw=60 Detector: w_monitor_I=0.0723**
+- **Test: x_screw=60 Detector: w_monitor_I=0.0721**
 - **TestScan: mxrun Czerny_Turner.instr -N21 x_screw=20,103 scan=1 -n1e6 Detector: Wavelength_I={21 values}**
 
 ## Input parameters
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md
index f08dca5bb5..dcd76ba85b 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md
@@ -20,7 +20,7 @@ You may DISCO_scan the grating monochromator with e.g.:
 
 - **TestScan: mxrun SOLEIL_DISCO.instr -N21 DISCO_x_screw=20,103 DISCO_scan=1 Detector: Wavelength_I={21 values}**
 - 
-- **Test: DISCO_grazing_angle_one=22.5 -n 1e5 Detector: DISCO_second_HFP_I=1.448e+16**
+- **Test: DISCO_grazing_angle_one=22.5 -n 1e5 Detector: DISCO_second_HFP_I=1.445e+16**
 
 ## Input parameters
 
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr
index 28a3884d9a..b4ca4b54b4 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr
@@ -18,7 +18,7 @@
 *   5306.77,5738.39,6151.91,6565.99,6979.47,7395.13,7786.17,8199.52,
 *   8631.53,9051.46,9483.21,9913.98,10299.5 }
 *
-* %Example: DISCO_grazing_angle_one=22.5 -n 1e5 Detector: DISCO_second_HFP_I=1.448e+16
+* %Example: DISCO_grazing_angle_one=22.5 -n 1e5 Detector: DISCO_second_HFP_I=1.445e+16
 *
 * %Parameters
 * DISCO_grazing_angle_one: [deg] Mirrors M1/M2 rotation angle
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/README.md
index 24bb970fd8..aa2adae939 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/README.md
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/README.md
@@ -28,7 +28,7 @@ This model implements the XRD station for powders.
 
 ## Examples
 
-- **Test: MARS_E0=16.99 Detector: MARS_detector_diffraction_I=1.334e+08**
+- **Test: MARS_E0=16.99 Detector: MARS_detector_diffraction_I=1.341e+08**
 
 ## Input parameters
 
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/SOLEIL_MARS.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/SOLEIL_MARS.instr
index ebbdefccb1..8984370585 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/SOLEIL_MARS.instr
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/SOLEIL_MARS.instr
@@ -24,7 +24,7 @@
 *
 * This model implements the XRD station for powders.
 *
-* %Example: MARS_E0=16.99 Detector: MARS_detector_diffraction_I=1.334e+08
+* %Example: MARS_E0=16.99 Detector: MARS_detector_diffraction_I=1.341e+08
 *
 * %Parameters
 * MARS_E0: [keV]     Central energy to be emitted by the source
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/README.md
index c7bf570a60..04c4c8f2ea 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/README.md
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/README.md
@@ -48,7 +48,7 @@ Pretty energy repartition monitor result: mxrun SOLEIL_ROCK.instr ROCK_E0=17.000
 
 ## Examples
 
-- **Test: ROCK_E0=15.918 Detector: ROCK_fluo_monitor_I=1.29e+11**
+- **Test: ROCK_E0=15.918 Detector: ROCK_fluo_monitor_I=1.61e+11**
 
 ## Input parameters
 
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/SOLEIL_ROCK.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/SOLEIL_ROCK.instr
index 398adb2de3..7cd24dd0c8 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/SOLEIL_ROCK.instr
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/SOLEIL_ROCK.instr
@@ -44,7 +44,7 @@
 * Manganese and chrome ROCK_scan: mxrun SOLEIL_ROCK.instr ROCK_E0=5.700,6.800 ROCK_scan=1 ROCK_sample_file=MnCr ROCK_cc=2 -N100
 * Pretty energy repartition monitor result: mxrun SOLEIL_ROCK.instr ROCK_E0=17.000 ROCK_cc=2 -n1e8
 *
-* %Example: ROCK_E0=15.918 Detector: ROCK_fluo_monitor_I=1.29e+11
+* %Example: ROCK_E0=15.918 Detector: ROCK_fluo_monitor_I=1.61e+11
 *
 * %Parameters
 * ROCK_E0:               [keV]    Energy to hit, i.e. selected by the channel-cut monochromator
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/README.md
index ca94ec984f..cb3d87ca31 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/README.md
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/README.md
@@ -34,7 +34,7 @@ Position | Element
 
 ## Examples
 
-- **Test: SAMBA_E0=15 Detector: SAMBA_det_fluo_I=5.92e+10**
+- **Test: SAMBA_E0=15 Detector: SAMBA_det_fluo_I=5.69e+10**
 
 ## Input parameters
 
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/SOLEIL_SAMBA.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/SOLEIL_SAMBA.instr
index fe6e9b84d5..f3e96f61a1 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/SOLEIL_SAMBA.instr
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/SOLEIL_SAMBA.instr
@@ -34,7 +34,7 @@
 * 30       | Sample Stage (XAS configuration)
 * 35       | Sample Stage (SEXAFS configuration)
 *
-* %Example: SAMBA_E0=15 Detector: SAMBA_det_fluo_I=5.92e+10
+* %Example: SAMBA_E0=15 Detector: SAMBA_det_fluo_I=5.69e+10
 *
 * %Parameters
 * SAMBA_E0:               [keV]  Central energy
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/README.md
index faa6f108fe..4ebd4abbea 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/README.md
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/README.md
@@ -30,7 +30,7 @@ Position | Element
 
 ## Examples
 
-- **Test: SIRIUS_E0=8 Detector: SIRIUS_def_loc_I=6.51e+09**
+- **Test: SIRIUS_E0=8 Detector: SIRIUS_def_loc_I=6.69e+09**
 
 ## Input parameters
 
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/SOLEIL_SIRIUS.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/SOLEIL_SIRIUS.instr
index 01d5e14c9d..2ad3ecd670 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/SOLEIL_SIRIUS.instr
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/SOLEIL_SIRIUS.instr
@@ -28,7 +28,7 @@
 * 30       | Sample Stage (XAS configuration)
 * 35       | Sample Stage (SEXAFS configuration)
 *
-* %Example: SIRIUS_E0=8 Detector: SIRIUS_def_loc_I=6.51e+09
+* %Example: SIRIUS_E0=8 Detector: SIRIUS_def_loc_I=6.69e+09
 *
 * %Parameters
 * SIRIUS_E0:               [keV]  Central energy (1.4 - 13 keV)
diff --git a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md
index 330c0ef250..25074a06d0 100644
--- a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md
+++ b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md
@@ -17,7 +17,7 @@ Simply a bending magnet illuminating a reflection grating.
 
 ## Examples
 
-- **Test: E0=1 Detector: psd_monitor_I=3.47e+09**
+- **Test: E0=1 Detector: psd_monitor_I=3.53e+09**
 
 ## Input parameters
 
diff --git a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr
index 00332db8f8..10c59582b9 100644
--- a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr
+++ b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr
@@ -15,7 +15,7 @@
 *
 * Simply a bending magnet illuminating a reflection grating.
 *
-* %Example: E0=1 Detector: psd_monitor_I=3.47e+09
+* %Example: E0=1 Detector: psd_monitor_I=3.53e+09
 *
 * %Parameters
 * E0: [keV] Source's center of emitted energy spectrum
diff --git a/mcxtrace-comps/examples/Tests_sources/Test_BM/README.md b/mcxtrace-comps/examples/Tests_sources/Test_BM/README.md
index c503862be0..28987cad2b 100644
--- a/mcxtrace-comps/examples/Tests_sources/Test_BM/README.md
+++ b/mcxtrace-comps/examples/Tests_sources/Test_BM/README.md
@@ -17,7 +17,7 @@ This is a simple test-instrument for the bending magnet component.
 
 ## Examples
 
-- **Test: Test_BM.instr SOURCE=0 e0=12.5 de=0.5 Detector: emon_I=9.146e+14**
+- **Test: Test_BM.instr SOURCE=0 e0=12.5 de=0.5 Detector: emon_I=9.212e+14**
 - **Test: Test_BM.instr SOURCE=1 e0=12.5 de=0.5 Detector: emon_I=1.28549e+15**
 
 ## Input parameters
diff --git a/mcxtrace-comps/examples/Tests_sources/Test_BM/Test_BM.instr b/mcxtrace-comps/examples/Tests_sources/Test_BM/Test_BM.instr
index 9abadf6208..29bd12c2c4 100644
--- a/mcxtrace-comps/examples/Tests_sources/Test_BM/Test_BM.instr
+++ b/mcxtrace-comps/examples/Tests_sources/Test_BM/Test_BM.instr
@@ -16,7 +16,7 @@
 * %Description
 * This is a simple test-instrument for the bending magnet component. 
 *
-* %Example: Test_BM.instr SOURCE=0 e0=12.5 de=0.5 Detector: emon_I=9.146e+14
+* %Example: Test_BM.instr SOURCE=0 e0=12.5 de=0.5 Detector: emon_I=9.212e+14
 * %Example: Test_BM.instr SOURCE=1 e0=12.5 de=0.5 Detector: emon_I=1.28549e+15
 *
 * %Parameters

From 8645e0e11ab1438cfba18eb8c325f26a54da4a98 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Thu, 29 Oct 2026 21:04:15 +0100
Subject: [PATCH 7/8] McXtrace examples: raise ncount for noisy Bending_magnet
 tests

SOLEIL_DISCO_2: the %Scan ran at the mctest default of 1e5 rays per
point. With the new Bending_magnet sampling +-2.5 mrad (was +-1/gamma),
2 of 21 points got no rays on w_monitor, so Wavelength=0. Use -n1e6,
as Czerny_Turner already does. Reproduced and checked locally with
seed 1000 and 2 MPI ranks (as in CI): all 21 points within 1.5%.

Test_grating_reflect: 25% seed-to-seed scatter at 1e6 (CI seed gave
65% of the 5-seed mean). Use -n 1e7: CI seed gives 89%.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md           | 2 +-
 mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr  | 2 +-
 .../examples/Tests_optics/Test_grating_reflect/README.md        | 2 +-
 .../Test_grating_reflect/Test_grating_reflect.instr             | 2 +-
 4 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md
index dcd76ba85b..7943269940 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md
@@ -18,7 +18,7 @@ You may DISCO_scan the grating monochromator with e.g.:
 
 ## Examples
 
-- **TestScan: mxrun SOLEIL_DISCO.instr -N21 DISCO_x_screw=20,103 DISCO_scan=1 Detector: Wavelength_I={21 values}**
+- **TestScan: mxrun SOLEIL_DISCO.instr -N21 DISCO_x_screw=20,103 DISCO_scan=1 -n1e6 Detector: Wavelength_I={21 values}**
 - 
 - **Test: DISCO_grazing_angle_one=22.5 -n 1e5 Detector: DISCO_second_HFP_I=1.445e+16**
 
diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr
index b4ca4b54b4..43a787d493 100644
--- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr
+++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr
@@ -13,7 +13,7 @@
 * %Description
 * DISCO beamline description here: https://hal.archives-ouvertes.fr/hal-01479318/document
 * You may DISCO_scan the grating monochromator with e.g.:
-* %Scan: mxrun SOLEIL_DISCO.instr -N21 DISCO_x_screw=20,103 DISCO_scan=1 Detector: Wavelength_I={
+* %Scan: mxrun SOLEIL_DISCO.instr -N21 DISCO_x_screw=20,103 DISCO_scan=1 -n1e6 Detector: Wavelength_I={
 *   2019.39,2365.97,2837.41,3242.59,3655.92,4070.23,4487.58,4921.15,
 *   5306.77,5738.39,6151.91,6565.99,6979.47,7395.13,7786.17,8199.52,
 *   8631.53,9051.46,9483.21,9913.98,10299.5 }
diff --git a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md
index 25074a06d0..143cc28999 100644
--- a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md
+++ b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md
@@ -17,7 +17,7 @@ Simply a bending magnet illuminating a reflection grating.
 
 ## Examples
 
-- **Test: E0=1 Detector: psd_monitor_I=3.53e+09**
+- **Test: E0=1 -n 1e7 Detector: psd_monitor_I=3.53e+09**
 
 ## Input parameters
 
diff --git a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr
index 10c59582b9..44005c09e7 100644
--- a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr
+++ b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr
@@ -15,7 +15,7 @@
 *
 * Simply a bending magnet illuminating a reflection grating.
 *
-* %Example: E0=1 Detector: psd_monitor_I=3.53e+09
+* %Example: E0=1 -n 1e7 Detector: psd_monitor_I=3.53e+09
 *
 * %Parameters
 * E0: [keV] Source's center of emitted energy spectrum

From fa469fae63806b21bd204bf2d1a39cd51dd3ffb6 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Thu, 29 Oct 2026 21:21:39 +0100
Subject: [PATCH 8/8] Bending_magnet: doc/cleanup fixes from review

- E0 now overrides lambda0, as documented (lambda mode only if E0==0).
- %D: reference Kim (X-Ray Data Booklet eq. 1), state that p is photons/s
  into the sampled solid angle and [E0-dE,E0+dE], that dE>0 is required,
  and that without focus the angles are sampled within +-2.5 mrad.
- No-focus mode: sample theta/psi uniformly in +-2.5 mrad instead of a
  gaussian truncated at 1 sigma (nearly flat anyway). Same expectation,
  but RNG use changes, so seeded results shift within statistics.
- Remove unused variables and stale/commented code.

CI-style check (seed 1000, 2 MPI ranks): all %Example tests within 14%
of the references, both %Scan tests within 1.2% on all points.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 mcxtrace-comps/sources/Bending_magnet.comp | 65 ++++++++--------------
 1 file changed, 24 insertions(+), 41 deletions(-)

diff --git a/mcxtrace-comps/sources/Bending_magnet.comp b/mcxtrace-comps/sources/Bending_magnet.comp
index 1f414162a5..2b5e1d58c6 100644
--- a/mcxtrace-comps/sources/Bending_magnet.comp
+++ b/mcxtrace-comps/sources/Bending_magnet.comp
@@ -18,7 +18,15 @@
 * Model of a bending magnet source
 * 
 * %Description
-* A source model based on the derivation from B.D. Patterson, Am. J. Phys. 79, 1046 (2011); doi: 10.1119/1.3614033
+* A source model using the full angular distribution of synchrotron radiation from a bending magnet,
+* eq. (1) of K.-J. Kim, X-Ray Data Booklet, Section 2.1 (https://xdb.lbl.gov/Section2/Sec_2-1.html).
+* In the horizontal plane (psi=0) this reduces to eq. (28) of B.D. Patterson, Am. J. Phys. 79, 1046 (2011); doi: 10.1119/1.3614033
+*
+* The ray weight p is the flux in [photons/s] emitted into the sampled solid angle and the energy
+* band [E0-dE, E0+dE] (i.e. NOT per mrad^2 or per 0.1% bandwidth). The energy is sampled uniformly, so dE>0
+* is required. In lambda0/dlambda mode the corresponding energy interval is sampled uniformly in energy.
+* With focus_xw/focus_yh the target window is sampled; without them, the horizontal/vertical emission angle
+* is sampled uniformly within +-2.5 mrad, so the total flux is then the flux into a 5 mrad horizontal fan.
 *
 * Example: Bending_magnet(
 *   E0 = 14, dE = 7, Ee = 2.75,
@@ -94,18 +102,11 @@ DECLARE
   double igamma;
   double kc; /*characteristic wavenumber of radiation from bending magnet*/
   double s1x;
-  double s1y; /*beam's size at dist (convolution of sigex/sigey and igamma)*/
+  double s1y; /*sigma of the gaussian sampling of the target window (gauss_t)*/
   double p0;
-
-  /*fine structure constant from CODATA*/
-  //const double alpha=7.2973525698e-3;
-  //const double MELE=9.10938291e-31;
-  
   double ctemult;
   double muc;
-  double cnstc;
-  double qe;    
-  
+  double qe;
 %}
 
 
@@ -116,7 +117,7 @@ INITIALIZE
     fprintf(stderr, "Error (%s): B, Ee, and Ie must all be > 0. Found (%g %g %g). Aborting.\n",NAME_CURRENT_COMP,B,Ee,Ie);
     exit(1);
   }
-  if (lambda0 > 0 && dlambda > 0 && lambda0-dlambda > 0) {
+  if (!E0 && lambda0 > 0 && dlambda > 0 && lambda0-dlambda > 0) {
     // compute equivalent E0 and dE
     double E1 = (2 * M_PI / (lambda0+dlambda)) * K2E;
     double E2 = (2 * M_PI / (lambda0-dlambda)) * K2E;
@@ -124,7 +125,7 @@ INITIALIZE
     dE = fabs(E1 - E2)/2;
   }
   if (E0 <= 0 || dE <= 0 || E0-dE <= 0) {
-    fprintf (stderr, "Error (%s): dE and E0 need to be set. Aborting.\n", NAME_CURRENT_COMP);
+    fprintf (stderr, "Error (%s): E0 and dE>0 (or lambda0 and dlambda>0) need to be set, with E0-dE>0. Aborting.\n", NAME_CURRENT_COMP);
     exit (1);
   }
   if (sigex <0 || sigey<0){
@@ -152,8 +153,6 @@ INITIALIZE
 
   muc = kc*M_C/(2.0*M_PI);
 
-  //s1x=sqrt(sigex*sigex + igamma*igamma*dist*dist);
-  //s1y=sqrt(sigey*sigey + igamma*igamma*dist*dist);
   if (focus_xw){
     s1x=focus_xw/2.0;
   }
@@ -162,8 +161,7 @@ INITIALIZE
   }
   p0=1.0/mcget_ncount();
 
-  cnstc = erf(1.0/(sqrt(2.0)));
-  qe = 1.0/(2.0*dE);    // requires dE > 0, no lambda handling here
+  qe = 1.0/(2.0*dE);    // energy sampling density, dE > 0 checked above
 %}
 
 
@@ -171,15 +169,13 @@ TRACE
 %{
 
   double xx,yy,x1,y1,z1;
-  double k,e,l;
-  double F1=1.0;
+  double k,e;
   double dx,dy,dz;
   double dl;
 
-  double angle, psi, theta, xcap, oneplusxsquared, chi, K23, K13;
-  double qtheta;
-  double qpsi;
-    
+  const double angle=2.5e-3; /*without focus, psi and theta are sampled uniformly in (-angle,angle)*/
+  double psi, theta, xcap, oneplusxsquared, chi, K23, K13;
+
   // initial source area
   xx=randnorm();
   yy=randnorm();
@@ -192,8 +188,7 @@ TRACE
   // E0 and dE are defined
   e=randpm1()*dE + E0;
   k=E2K*e;
-  
-  double oneoverk=1.0/k;
+
   double mu = k*M_C/(2*M_PI);
   double mumuc = mu/muc;
 
@@ -208,21 +203,15 @@ TRACE
     }else {
       do {
         x1=randnorm()*s1x;
-      }while (focus_xw!=0 && fabs(x1)>focus_xw/2.0);
+      }while (fabs(x1)>focus_xw/2.0);
       theta = atan((x1-x)/dist);
       /*adjust for restricted sampling window*/            
       p*=erf(focus_xw*0.5*M_SQRT1_2/s1x)*s1x*sqrt(2.0*M_PI)*exp((x1*x1)/(2*s1x*s1x)); // 1/sampling probability density
       p*= dist/(dist*dist + (x1-x)*(x1-x));
     }
   }else{    
-    angle=2.5e-3; //angle is the angle psi and theta will be generated over. Range of (-angle,angle).
-    do {
-      theta = randnorm()*angle;
-    } while(theta>angle || theta<-angle);  
-
-    qtheta = 1.0/(sqrt(2*M_PI)*angle*cnstc)*exp(-(theta*theta)/(2*angle*angle)); //sampling probability density
-    p*=1.0/qtheta;
-    
+    theta = randpm1()*angle;
+    p*=2*angle; // 1/sampling probability density
     x1=x+dist*tan(theta);
   }
   if (focus_yh){
@@ -242,14 +231,8 @@ TRACE
       p*= dist/(dist*dist + (y1-y)*(y1-y));
     }
   }else{      
-    angle=2.5e-3; //angle is the angle psi and theta will be generated over. Range of (-angle,angle).
-    do {
-      psi = randnorm()*angle;
-    } while(psi>angle || psi<-angle);      
-    
-    qpsi = 1.0/(sqrt(2*M_PI)*angle*cnstc)*exp(-(psi*psi)/(2*angle*angle)); //sampling probability density
-    p*=1.0/qpsi;
-    
+    psi = randpm1()*angle;
+    p*=2*angle; // 1/sampling probability density
     y1=y+dist*tan(psi);
   }