From 4db2b3a6f1f70f6b98eaa9d5b7203995b4d74dfe Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Mon, 6 Jul 2026 11:07:49 +0200
Subject: [PATCH 01/17] Drop prms union in reflectivity-lib.h

---
 mcxtrace-comps/share/reflectivity-lib.h | 16 +++++++---------
 1 file changed, 7 insertions(+), 9 deletions(-)

diff --git a/mcxtrace-comps/share/reflectivity-lib.h b/mcxtrace-comps/share/reflectivity-lib.h
index 353659bda..0b0c6aa54 100644
--- a/mcxtrace-comps/share/reflectivity-lib.h
+++ b/mcxtrace-comps/share/reflectivity-lib.h
@@ -48,15 +48,13 @@ typedef struct t_reflec_eth_prmtc{
 } t_reflec_theta_e_prmtc;
 typedef struct reflec_T {
   enum reflec_Type type;
-  union {
-    struct t_reflec_bare rb;
-    struct t_reflec_coating rc;
-    struct t_reflec_q_prmtc rqpm;
-    struct t_reflec_parratt rp;
-    struct t_reflec_constant rconst;
-    struct t_reflec_kinematic rk;
-    struct t_reflec_eth_prmtc rethpm;
-  } prms;
+  struct t_reflec_bare rb;
+  struct t_reflec_coating rc;
+  struct t_reflec_q_prmtc rqpm;
+  struct t_reflec_parratt rp;
+  struct t_reflec_constant rconst;
+  struct t_reflec_kinematic rk;
+  struct t_reflec_eth_prmtc rethpm;
 } t_Reflec;
 /* reflectivity-lib.c */
 int reflec_Init(t_Reflec *R, enum reflec_Type type, char *file, void *pars);

From 29e6b0833987225b8368e8075a94c0992a3d267f Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Mon, 6 Jul 2026 11:09:14 +0200
Subject: [PATCH 02/17] Removal of prms union in reflectivity-lib.c, change to
 R=1.0 default, PARRATT allowed but not on GPU

---
 mcxtrace-comps/share/reflectivity-lib.c | 145 ++++++++++++------------
 1 file changed, 75 insertions(+), 70 deletions(-)

diff --git a/mcxtrace-comps/share/reflectivity-lib.c b/mcxtrace-comps/share/reflectivity-lib.c
index 9930c4a26..c6ece75c2 100644
--- a/mcxtrace-comps/share/reflectivity-lib.c
+++ b/mcxtrace-comps/share/reflectivity-lib.c
@@ -37,7 +37,7 @@ int reflec_Init(t_Reflec *R, enum reflec_Type typ, char *file, void *pars){
           if (pars){
 	    reflec_Init_const(R,((double *)pars)[0]);
           }else{
-              reflec_Init_const(R,0);
+              reflec_Init_const(R,1.0);
 	  }
 	  break;
       }
@@ -48,10 +48,10 @@ int reflec_Init(t_Reflec *R, enum reflec_Type typ, char *file, void *pars){
           if(pars){
               reflec_Init_parratt(R, N,((double **) pars)[1], ((double **) pars)[2], ((double **) pars)[3]);        
           }else{
-              fprintf(stderr,"Warning: %s: No parameters specified to Parratt reflectivity algortihm. Setting R=0.\n",
+              fprintf(stderr,"Warning: %s: No parameters specified to Parratt reflectivity algortihm. Setting R=1.0.\n",
                       REFLIBNAME);
               R->type=CONSTANT;
-              R->prms.rconst.R=0;
+              R->rconst.R=1.0;
           }
           break;
       }
@@ -60,41 +60,44 @@ int reflec_Init(t_Reflec *R, enum reflec_Type typ, char *file, void *pars){
           if(pars){
 	    reflec_Init_kinematic(R, (int) ((double *)pars)[0],((double *)pars)[1],((double *)pars)[2],((double *)pars)[3]);
           }else{
-            reflec_Init_kinematic(R, (int) 0, 0.0, 0.0, 0.0);
+              fprintf(stderr,"Warning: %s: No parameters specified to KINEMATIC reflectivity algortihm. Setting R=1.0.\n",
+                      REFLIBNAME);
+              R->type=CONSTANT;
+              R->rconst.R=1.0;
           }
           break;
       }
     case BARE:
       {
-          stracpy(R->prms.rb.matrl,file,255);
-          reflec_Init_File(R,R->prms.rb.matrl);
+          stracpy(R->rb.matrl,file,255);
+          reflec_Init_File(R,R->rb.matrl);
           break;
       }
     case COATING:
       {
-          stracpy(R->prms.rc.matrl,file,255);
-          reflec_Init_File(R,R->prms.rc.matrl);
+          stracpy(R->rc.matrl,file,255);
+          reflec_Init_File(R,R->rc.matrl);
           break;
       }
     case Q_PARAMETRIC:
       {
-          stracpy(R->prms.rqpm.fname,file,255);
-          reflec_Init_File(R,R->prms.rqpm.fname);
+          stracpy(R->rqpm.fname,file,255);
+          reflec_Init_File(R,R->rqpm.fname);
           break;
       }
     case ETH_PARAMETRIC:
       {
-          stracpy(R->prms.rethpm.fname,file,255);
-          reflec_Init_File(R,R->prms.rethpm.fname);
+          stracpy(R->rethpm.fname,file,255);
+          reflec_Init_File(R,R->rethpm.fname);
           break;
       }
     default:
       {
-        fprintf(stderr,"Error: %s: Undetermined reflectivity parameterization type. Setting R=0.\n",
+        fprintf(stderr,"Error: %s: Undetermined reflectivity parameterization type. ABORT.\n",
                 REFLIBNAME);
 	free(R);
 	R=NULL;
-	return -1;
+	exit(-1);
       }
   }
   return 0;
@@ -102,24 +105,24 @@ int reflec_Init(t_Reflec *R, enum reflec_Type typ, char *file, void *pars){
 
 int reflec_Init_parratt(t_Reflec *R, int N, double *d, double *delta, double *beta){
     R->type=PARRATT;
-    R->prms.rp.N=N;
-    R->prms.rp.d=d;
-    R->prms.rp.delta=delta;
-    R->prms.rp.beta=beta;
+    R->rp.N=N;
+    R->rp.d=d;
+    R->rp.delta=delta;
+    R->rp.beta=beta;
     return 0;
 }
 
 int reflec_Init_kinematic(t_Reflec *R, int N, double Gamma, double Lambda, double rhoAB){
     R->type=KINEMATIC;
-    R->prms.rk.N=N;
-    R->prms.rk.Gamma=Gamma;
-    R->prms.rk.Lambda=Lambda;
-    R->prms.rk.rho_AB=rhoAB;
+    R->rk.N=N;
+    R->rk.Gamma=Gamma;
+    R->rk.Lambda=Lambda;
+    R->rk.rho_AB=rhoAB;
     return 0;
 }
  int reflec_Init_const(t_Reflec *R, double R0){
    R->type=CONSTANT;
-   R->prms.rconst.R=R0;
+   R->rconst.R=R0;
    return 0;
  }
 
@@ -140,7 +143,7 @@ int reflec_Init_File(t_Reflec *R, char *filename){
       fprintf(stderr,"Warning: %s: No reflectivity file given. Surface is opaque.\n",
               REFLIBNAME);
       R->type=CONSTANT;
-      R->prms.rconst.R=0;
+      R->rconst.R=0;
       return 0;
     }
 
@@ -152,7 +155,7 @@ int reflec_Init_File(t_Reflec *R, char *filename){
     switch(R->type){
       case CONSTANT:
         {
-          R->prms.rconst.R=Table_Index(*table, 0, 0);
+          R->rconst.R=Table_Index(*table, 0, 0);
           break;
         }
 
@@ -164,8 +167,8 @@ int reflec_Init_File(t_Reflec *R, char *filename){
                     REFLIBNAME, filename);
             exit(-1);
           }
-          stracpy(R->prms.rb.matrl,header_parsed[0],255);
-          R->prms.rb.d=strtod(header_parsed[1], NULL);
+          stracpy(R->rb.matrl,header_parsed[0],255);
+          R->rb.d=strtod(header_parsed[1], NULL);
           break;
         }
 
@@ -177,22 +180,22 @@ int reflec_Init_File(t_Reflec *R, char *filename){
                     REFLIBNAME,filename);
             exit(-1);
           } else {
-            stracpy(R->prms.rc.matrl,header_parsed[0],255);
-            R->prms.rc.T = table;
-            R->prms.rc.d = strtod(header_parsed[4], NULL);
-            R->prms.rc.rho=strtod(header_parsed[3],NULL);
-            R->prms.rc.Z=strtod(header_parsed[1],NULL);
-            R->prms.rc.At=strtod(header_parsed[2],NULL);
+            stracpy(R->rc.matrl,header_parsed[0],255);
+            R->rc.T = table;
+            R->rc.d = strtod(header_parsed[4], NULL);
+            R->rc.rho=strtod(header_parsed[3],NULL);
+            R->rc.Z=strtod(header_parsed[1],NULL);
+            R->rc.At=strtod(header_parsed[2],NULL);
           }
           break;
         }
 
       case Q_PARAMETRIC:
         {
-          stracpy(R->prms.rqpm.fname,filename,255);
-          R->prms.rqpm.T=table;
-          R->prms.rqpm.qmin=Table_Index(*table,0,0);
-          R->prms.rqpm.qmax=Table_Index(*table,table->rows,0);
+          stracpy(R->rqpm.fname,filename,255);
+          R->rqpm.T=table;
+          R->rqpm.qmin=Table_Index(*table,0,0);
+          R->rqpm.qmax=Table_Index(*table,table->rows,0);
           break;
         }
 
@@ -205,13 +208,13 @@ int reflec_Init_File(t_Reflec *R, char *filename){
             exit(-1);
           }
           unsigned long N=strtol(header_parsed[0], NULL, 10);
-          R->prms.rp.N = N;
-          R->prms.rp.d = calloc(N,sizeof(double));
-          *(R->prms.rp.d) = strtod(header_parsed[1], NULL);
-          R->prms.rp.delta = calloc(N,sizeof(double));
-          *(R->prms.rp.delta) = strtod(header_parsed[2], NULL);
-          R->prms.rp.beta = calloc(N,sizeof(double));
-          *(R->prms.rp.beta) = strtod(header_parsed[3], NULL);
+          R->rp.N = N;
+          R->rp.d = calloc(N,sizeof(double));
+          *(R->rp.d) = strtod(header_parsed[1], NULL);
+          R->rp.delta = calloc(N,sizeof(double));
+          *(R->rp.delta) = strtod(header_parsed[2], NULL);
+          R->rp.beta = calloc(N,sizeof(double));
+          *(R->rp.beta) = strtod(header_parsed[3], NULL);
           break;
         }
 
@@ -224,17 +227,17 @@ int reflec_Init_File(t_Reflec *R, char *filename){
             exit(-1);
           }
 
-          R->prms.rk.N = strtol(header_parsed[0], NULL, 10);
-          R->prms.rk.Gamma = strtod(header_parsed[2], NULL);
-          R->prms.rk.Lambda = strtod(header_parsed[3], NULL);
-          R->prms.rk.rho_AB = strtod(header_parsed[4], NULL);
+          R->rk.N = strtol(header_parsed[0], NULL, 10);
+          R->rk.Gamma = strtod(header_parsed[2], NULL);
+          R->rk.Lambda = strtod(header_parsed[3], NULL);
+          R->rk.rho_AB = strtod(header_parsed[4], NULL);
           break;
         }
 
       case ETH_PARAMETRIC:
         {
-          stracpy(R->prms.rethpm.fname,filename,255);
-          R->prms.rethpm.T=table;
+          stracpy(R->rethpm.fname,filename,255);
+          R->rethpm.T=table;
 
           /*parse header for E_min E_max etc.*/
           char **header_parsed = Table_ParseHeader(table->header,"e_min=","e_max=","theta_min=","theta_max=",NULL);
@@ -243,12 +246,12 @@ int reflec_Init_File(t_Reflec *R, char *filename){
                     REFLIBNAME,filename); //1619
             exit(-1);
           }
-          R->prms.rethpm.emin=strtod(header_parsed[0],NULL);
-          R->prms.rethpm.emax=strtod(header_parsed[1],NULL);
-          R->prms.rethpm.thetamin=strtod(header_parsed[2],NULL);
-          R->prms.rethpm.thetamax=strtod(header_parsed[3],NULL);
-          int rows = R->prms.rethpm.T->rows;
-          int cols = R->prms.rethpm.T->columns;
+          R->rethpm.emin=strtod(header_parsed[0],NULL);
+          R->rethpm.emax=strtod(header_parsed[1],NULL);
+          R->rethpm.thetamin=strtod(header_parsed[2],NULL);
+          R->rethpm.thetamax=strtod(header_parsed[3],NULL);
+          int rows = R->rethpm.T->rows;
+          int cols = R->rethpm.T->columns;
           if(rows == 0){ //implies cols == 0 as well
             fprintf(stderr,"Error: %s: File %s contains no table.\n",
                     REFLIBNAME, filename);
@@ -257,16 +260,16 @@ int reflec_Init_File(t_Reflec *R, char *filename){
             if(rows == 1){
               printf("Warning: %s: File %s contains only a single row. Setting e_step = 0.\n",
                      REFLIBNAME, filename);
-              R->prms.rethpm.estep=0;
+              R->rethpm.estep=0;
             } else {
-              R->prms.rethpm.estep=(R->prms.rethpm.emax - R->prms.rethpm.emin)/(rows-1);
+              R->rethpm.estep=(R->rethpm.emax - R->rethpm.emin)/(rows-1);
             }
             if(cols == 1){
               printf("Warning: %s: File %s contains only a single column. Setting theta_step = 0.\n",
                      REFLIBNAME, filename);
-              R->prms.rethpm.thetastep=0;
+              R->rethpm.thetastep=0;
             } else {
-              R->prms.rethpm.thetastep=(R->prms.rethpm.thetamax - R->prms.rethpm.thetamin)/(cols-1);
+              R->rethpm.thetastep=(R->rethpm.thetamax - R->rethpm.thetamin)/(cols-1);
             }
           }
           break;
@@ -333,7 +336,7 @@ enum reflec_Type get_table_reflec_type(t_Table *t){
 
 /*This section contains the functions that compute the actual reflectivity*/
 cdouble reflec_coating(t_Reflec r_handle, double q, double g, double k){
-    struct t_reflec_coating *ptr=&(r_handle.prms.rc);
+    struct t_reflec_coating *ptr=&(r_handle.rc);
     /*adjust p according to reflectivity*/
     double Qc,f1,f2,na,e;
     /*length of wavevector transfer may be compute from s and n_ above*/
@@ -374,7 +377,7 @@ cdouble reflec_bare(t_Reflec r_handle, double q, double g){
 
 cdouble reflec_kinematic(t_Reflec r_handle, double q, double g){
     cdouble r1,rN,num_factor,den_factor;
-    struct t_reflec_kinematic *ptr=&(r_handle.prms.rk);
+    struct t_reflec_kinematic *ptr=&(r_handle.rk);
     double zeta=ptr->Lambda*q/(2*M_PI);
     double beta=0;
     double r1_scale;
@@ -396,7 +399,7 @@ cdouble reflec_kinematic(t_Reflec r_handle, double q, double g){
 
 cdouble reflec_q_prmtc(t_Reflec r_handle, double q, double g){
     double r;
-    struct t_reflec_q_prmtc *ptr=&(r_handle.prms.rqpm);
+    struct t_reflec_q_prmtc *ptr=&(r_handle.rqpm);
     if (ptr->T->columns>2){
         double c=(ptr->T->columns-2)*g+1;
         r=Table_Value2d(*(ptr->T),ptr->T->rows*q/(ptr->qmax-ptr->qmin),c);
@@ -408,7 +411,7 @@ cdouble reflec_q_prmtc(t_Reflec r_handle, double q, double g){
 
 cdouble reflec_eth_prmtc(t_Reflec r_handle, double g, double e, double th){
     double r,ec,thc;
-    struct t_reflec_eth_prmtc *ptr=&(r_handle.prms.rethpm);
+    struct t_reflec_eth_prmtc *ptr=&(r_handle.rethpm);
     ec=(ptr->T->rows-1) * (e-ptr->emin)/(ptr->emax - ptr->emin);
     thc=(ptr->T->columns-1)*(th-ptr->thetamin)/(ptr->thetamax - ptr->thetamin);
     r=Table_Value2d(*(ptr->T),ec,thc);
@@ -421,7 +424,7 @@ cdouble reflec_parratt(t_Reflec r_handle, double q, double g, double k){
     double k2=k*k;
     cdouble qinf;
     cdouble qpd,rd,p;
-    t_reflec_parratt *pp=&(r_handle.prms.rp);
+    t_reflec_parratt *pp=&(r_handle.rp);
 
     qp=cplx(q,0.0);
     /*q*q - 8*k2* *(pp->delta) + I*8*k2* *(pp->beta)*/
@@ -477,7 +480,7 @@ cdouble refleccq( t_Reflec r_handle, double q, double g, double k, double theta)
     switch(r_handle.type){
       case CONSTANT:
         {
-          r=cplx(r_handle.prms.rconst.R,0.0);
+          r=cplx(r_handle.rconst.R,0.0);
           break;
         }
       case BARE:
@@ -528,7 +531,7 @@ double reflecq( t_Reflec r_handle, double q, double g, double k, double theta){
     switch(r_handle.type){
       case CONSTANT:
         {
-          r=fabs(r_handle.prms.rconst.R);
+          r=1.0; //fabs(r_handle.rconst.R);
           break;
         }
     case BARE:
@@ -541,18 +544,20 @@ double reflecq( t_Reflec r_handle, double q, double g, double k, double theta){
           cdouble rp;
           rp=reflec_coating(r_handle,q,g,k);
           r= sqrt(creal(cmul(rp, conj(rp))));
-          break;
+	  break;
         }
     case Q_PARAMETRIC:
         {
           r=cabs(reflec_q_prmtc(r_handle,q,g));
           break;
         }
-        /*case PARRATT:
+    #ifndef OPENACC
+    case PARRATT:
         {
           r=cabs(reflec_parratt(r_handle,q,g,k));
           break;
-	  }*/
+	}
+    #endif
     case ETH_PARAMETRIC:
         {
             r=cabs(reflec_eth_prmtc(r_handle,g,k*K2E,theta));

From 35733adde2efb19385b468768f86053dee1b1310 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Mon, 6 Jul 2026 11:10:48 +0200
Subject: [PATCH 03/17] Drop unused EE, switch to var-naming like
 Mirror_curved, correct INIT

---
 .../optics/Multilayer_elliptic.comp           | 38 +++++++++----------
 1 file changed, 19 insertions(+), 19 deletions(-)

diff --git a/mcxtrace-comps/optics/Multilayer_elliptic.comp b/mcxtrace-comps/optics/Multilayer_elliptic.comp
index 062497430..0ce4fc559 100644
--- a/mcxtrace-comps/optics/Multilayer_elliptic.comp
+++ b/mcxtrace-comps/optics/Multilayer_elliptic.comp
@@ -67,7 +67,7 @@ SETTING PARAMETERS (string coating="Ref_W_B4C.txt",
 
 SHARE 
 %{
-  #include <complex.h>
+  %include "mccode-complex-lib"
   %include "read_table-lib"
   %include "reflectivity-lib"
   /*something that would be relevant for ALL elliptical mirrors*/
@@ -132,7 +132,7 @@ INITIALIZE
       status = reflec_Init (&re, COATING_UNDEFINED, coating, NULL);
     } else {
       /*assume a constant reflectivity*/
-      status = reflec_Init_const (&re, R0);
+      status = reflec_Init (&re, CONSTANT, NULL, &(R0));
     }
   } else {
     kinematical = 1;
@@ -144,7 +144,7 @@ INITIALIZE
 
 TRACE
 %{
-  double K, vink;
+  double k, s;
   double x_el, y_el, z_el;    // beginning coordinates transformed into the ellipse system
   double kx_el, ky_el, kz_el; // kvector transformed into the ellipse system, hence
 
@@ -154,7 +154,7 @@ TRACE
   double kxn, kyn, kzn; // reflected ray's kvector
 
   /* get the photon's coordinates and kvector in the ellipse frame */
-  K = sqrt (kx * kx + ky * ky + kz * kz);
+  k = sqrt (kx * kx + ky * ky + kz * kz);
 
   CoordTransME (&x_el, &y_el, &z_el, x, y, z, Z0, Y0, xi);
   CoordTransME (&kx_el, &ky_el, &kz_el, kx, ky, kz, 0, 0, xi);
@@ -193,10 +193,10 @@ TRACE
           nz = 1.0;
         }
         NORM (nx, ny, nz);
-        vink = scalar_prod (nx, ny, nz, kx_el, ky_el, kz_el);
-        kxn = kx_el - 2.0 * vink * nx;
-        kyn = ky_el - 2.0 * vink * ny;
-        kzn = kz_el - 2.0 * vink * nz;
+        s = scalar_prod (nx, ny, nz, kx_el, ky_el, kz_el);
+        kxn = kx_el - 2.0 * s * nx;
+        kyn = ky_el - 2.0 * s * ny;
+        kzn = kz_el - 2.0 * s * nz;
         NORM (kxn, kyn, kzn);
 
         double kxo, kyo, kzo;
@@ -204,13 +204,13 @@ TRACE
         kyo = ky, kzo = kz;
         CoordTransEM (&kx, &ky, &kz, kxn, kyn, kzn, 0, 0, xi);
 
-        kx = K * kx;
-        ky = K * ky;
-        kz = K * kz;
+        kx = k * kx;
+        ky = k * ky;
+        kz = k * kz;
 
-        double QQ, EE, Ref;
+        double QQ;
+        cdouble Ref;
         QQ = sqrt ((kx - kxo) * (kx - kxo) + (ky - kyo) * (ky - kyo) + (kz - kzo) * (kz - kzo));
-        EE = K * K2E;
 
         if (kinematical) {
           /*
@@ -219,17 +219,17 @@ TRACE
            * r1(zeta) = 2 i r_0 \rho_{AB} \left(\frac{\Lambda^2 \Gamma}{\zeta}\right) \frac{\sin\left(\pi\Gamma\zeta\right)}{\pi\Gamma\zeta);
            */
           Ref = reflecq (re, QQ, 0, 0, 0);
-          if (Ref > 1) {
-            /*Reflectivity can't be >1*/
-            Ref = 1.0;
-          }
         } else {
           /*interpolate in table*/
           Ref = reflecq (re, QQ, 0, 0, 0);
         }
-
+        double pmul = sqrt (creal (cmul(Ref, conj (Ref))));
+        /*Reflectivity can't be >1*/
+        if (pmul>1) {
+          pmul = 1.0;
+        }
         /* apply reflectivity */
-        p *= Ref;
+        p *= pmul;
         SCATTER;
 
       } else {

From ebd6f2cf515480c4e5652564949c6b2e05c2e167 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Mon, 6 Jul 2026 11:18:59 +0200
Subject: [PATCH 04/17] Test Multilayer_elliptic in 3 cases and pothole just in
 1...

---
 .../examples/Tests_optics/Test_Mirrors/Test_Mirrors.instr | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/mcxtrace-comps/examples/Tests_optics/Test_Mirrors/Test_Mirrors.instr b/mcxtrace-comps/examples/Tests_optics/Test_Mirrors/Test_Mirrors.instr
index 64242f131..453359176 100644
--- a/mcxtrace-comps/examples/Tests_optics/Test_Mirrors/Test_Mirrors.instr
+++ b/mcxtrace-comps/examples/Tests_optics/Test_Mirrors/Test_Mirrors.instr
@@ -64,8 +64,8 @@ INITIALIZE %{
     case 2: printf("Using Mirror_parabolic\n"); break; // OK when rotated 180 deg around z
     case 3: printf("Using Mirror_curved\n"); break; // OK
     case 4: printf("Using Mirror_toroid\n"); break;         // OK
-    case 5: printf("Using Mirror_toroid_pothole\n"); break; // OK
-    case 6: printf("Using Multilayer_elliptic\n"); break; // not OK
+    case 5: printf("Using Multilayer_elliptic\n"); break; // not OK
+    case 6: printf("Using Mirror_toroid_pothole\n"); break; // OK
     default: exit(fprintf(stderr, "Unknown Mirror index\n"));
   }
   /* If "None" is given it in fact means "" */
@@ -128,13 +128,13 @@ COMPONENT mt6_multilayer = Multilayer_elliptic( // default in XZ
     coating=coating, theta = gamma,
     s1 = L, s2 = L, length = 0.2, width = 0.05, R0 = 1,
     Emin=E0-1, Emax=E0+1, Estep=0.05)
-WHEN (index==6) 
+WHEN (index==5)
 AT(0,0,0) RELATIVE mirror_stage
 ROTATED (0, 0, 0) RELATIVE mirror_stage
 
 COMPONENT mt5_toroid_pot = Mirror_toroid_pothole( // default in XZ
     radius=0.1, radius_o=radius, xwidth=5e-2, zdepth=2e-1,R0=1, coating="") /* Does not support standard coatings it seems */
-WHEN (index==5) 
+WHEN (index==6)
 AT(0,0,0) RELATIVE mirror_stage
 ROTATED (0, 0, 180) RELATIVE mirror_stage
 

From e80838474ac6caa0e09f4ddf13236162a0214bea Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Mon, 6 Jul 2026 11:30:25 +0200
Subject: [PATCH 05/17] Fix comp instance naming in agreement with index var

---
 .../examples/Tests_optics/Test_Mirrors/Test_Mirrors.instr     | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mcxtrace-comps/examples/Tests_optics/Test_Mirrors/Test_Mirrors.instr b/mcxtrace-comps/examples/Tests_optics/Test_Mirrors/Test_Mirrors.instr
index 453359176..cd1da3618 100644
--- a/mcxtrace-comps/examples/Tests_optics/Test_Mirrors/Test_Mirrors.instr
+++ b/mcxtrace-comps/examples/Tests_optics/Test_Mirrors/Test_Mirrors.instr
@@ -124,7 +124,7 @@ WHEN (index==4)
 AT(0,0,0) RELATIVE mirror_stage
 ROTATED (0, 0, 180) RELATIVE mirror_stage
 
-COMPONENT mt6_multilayer = Multilayer_elliptic( // default in XZ
+COMPONENT mt5_multilayer = Multilayer_elliptic( // default in XZ
     coating=coating, theta = gamma,
     s1 = L, s2 = L, length = 0.2, width = 0.05, R0 = 1,
     Emin=E0-1, Emax=E0+1, Estep=0.05)
@@ -132,7 +132,7 @@ WHEN (index==5)
 AT(0,0,0) RELATIVE mirror_stage
 ROTATED (0, 0, 0) RELATIVE mirror_stage
 
-COMPONENT mt5_toroid_pot = Mirror_toroid_pothole( // default in XZ
+COMPONENT mt6_toroid_pot = Mirror_toroid_pothole( // default in XZ
     radius=0.1, radius_o=radius, xwidth=5e-2, zdepth=2e-1,R0=1, coating="") /* Does not support standard coatings it seems */
 WHEN (index==6)
 AT(0,0,0) RELATIVE mirror_stage

From 97030f65e08fd9bd710539a544f8cd5db9a274d6 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Mon, 6 Jul 2026 11:39:43 +0200
Subject: [PATCH 06/17] Switch back to a simple double here

---
 mcxtrace-comps/optics/Multilayer_elliptic.comp | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/mcxtrace-comps/optics/Multilayer_elliptic.comp b/mcxtrace-comps/optics/Multilayer_elliptic.comp
index 0ce4fc559..c80d9944c 100644
--- a/mcxtrace-comps/optics/Multilayer_elliptic.comp
+++ b/mcxtrace-comps/optics/Multilayer_elliptic.comp
@@ -209,7 +209,7 @@ TRACE
         kz = k * kz;
 
         double QQ;
-        cdouble Ref;
+	double Ref=1.0;
         QQ = sqrt ((kx - kxo) * (kx - kxo) + (ky - kyo) * (ky - kyo) + (kz - kzo) * (kz - kzo));
 
         if (kinematical) {
@@ -223,13 +223,12 @@ TRACE
           /*interpolate in table*/
           Ref = reflecq (re, QQ, 0, 0, 0);
         }
-        double pmul = sqrt (creal (cmul(Ref, conj (Ref))));
-        /*Reflectivity can't be >1*/
-        if (pmul>1) {
-          pmul = 1.0;
+	/*Reflectivity can't be >1*/
+        if (Ref>1) {
+          Ref = 1.0;
         }
         /* apply reflectivity */
-        p *= pmul;
+        p *= Ref;
         SCATTER;
 
       } else {

From d7b49ec0ebf855dcd05c2283927fc84e98987a8e Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Mon, 6 Jul 2026 11:45:51 +0200
Subject: [PATCH 07/17] Use scalar product here for similarity with
 Mirror_curved

---
 mcxtrace-comps/optics/Multilayer_elliptic.comp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mcxtrace-comps/optics/Multilayer_elliptic.comp b/mcxtrace-comps/optics/Multilayer_elliptic.comp
index c80d9944c..19d3ed9d3 100644
--- a/mcxtrace-comps/optics/Multilayer_elliptic.comp
+++ b/mcxtrace-comps/optics/Multilayer_elliptic.comp
@@ -154,7 +154,7 @@ TRACE
   double kxn, kyn, kzn; // reflected ray's kvector
 
   /* get the photon's coordinates and kvector in the ellipse frame */
-  k = sqrt (kx * kx + ky * ky + kz * kz);
+  k = sqrt (scalar_prod (kx, ky, kz, kx, ky, kz));
 
   CoordTransME (&x_el, &y_el, &z_el, x, y, z, Z0, Y0, xi);
   CoordTransME (&kx_el, &ky_el, &kz_el, kx, ky, kz, 0, 0, xi);

From e26e3e4c4513b6174efc4875216611cf5eade369 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Mon, 6 Jul 2026 12:52:31 +0200
Subject: [PATCH 08/17] Hack around constant reflectivity-case in
 Multilayer_elliptic

---
 .../optics/Multilayer_elliptic.comp           | 39 ++++++++++---------
 1 file changed, 21 insertions(+), 18 deletions(-)

diff --git a/mcxtrace-comps/optics/Multilayer_elliptic.comp b/mcxtrace-comps/optics/Multilayer_elliptic.comp
index 19d3ed9d3..29ef643d3 100644
--- a/mcxtrace-comps/optics/Multilayer_elliptic.comp
+++ b/mcxtrace-comps/optics/Multilayer_elliptic.comp
@@ -101,6 +101,7 @@ DECLARE
   double cost0;
   int kinematical;
   t_Reflec re;
+  int isconst;
 %}
 
 INITIALIZE
@@ -108,7 +109,7 @@ INITIALIZE
   /* calculation of the elliptical parameters according to the input mirror parameters:
   ellipse major axis a/2, minor axis b/2, M-magnification factor, Z0&Y0 - position of the mirror centre in the elliptical coordinate system.*/
   double Theta = DEG2RAD * theta;
-
+  isconst=0;
   if (xwidth)
     width = xwidth;
   if (zdepth)
@@ -132,6 +133,7 @@ INITIALIZE
       status = reflec_Init (&re, COATING_UNDEFINED, coating, NULL);
     } else {
       /*assume a constant reflectivity*/
+      isconst=1;
       status = reflec_Init (&re, CONSTANT, NULL, &(R0));
     }
   } else {
@@ -209,24 +211,25 @@ TRACE
         kz = k * kz;
 
         double QQ;
-	double Ref=1.0;
+	double Ref=R0;
         QQ = sqrt ((kx - kxo) * (kx - kxo) + (ky - kyo) * (ky - kyo) + (kz - kzo) * (kz - kzo));
-
-        if (kinematical) {
-          /*
-           * \Lambda: thickness of bilayer - following notation in Als-Nielsen/McMorrow
-           * \Gamma:  \Gamma*\Lambda thickness of high electron density material.
-           * r1(zeta) = 2 i r_0 \rho_{AB} \left(\frac{\Lambda^2 \Gamma}{\zeta}\right) \frac{\sin\left(\pi\Gamma\zeta\right)}{\pi\Gamma\zeta);
-           */
-          Ref = reflecq (re, QQ, 0, 0, 0);
-        } else {
-          /*interpolate in table*/
-          Ref = reflecq (re, QQ, 0, 0, 0);
-        }
-	/*Reflectivity can't be >1*/
-        if (Ref>1) {
-          Ref = 1.0;
-        }
+	if (!isconst) {
+          if (kinematical) {
+            /*
+             * \Lambda: thickness of bilayer - following notation in Als-Nielsen/McMorrow
+             * \Gamma:  \Gamma*\Lambda thickness of high electron density material.
+             * r1(zeta) = 2 i r_0 \rho_{AB} \left(\frac{\Lambda^2 \Gamma}{\zeta}\right) \frac{\sin\left(\pi\Gamma\zeta\right)}{\pi\Gamma\zeta);
+             */
+            Ref = reflecq (re, QQ, 0, 0, 0);
+          } else {
+            /*interpolate in table*/
+            Ref = reflecq (re, QQ, 0, 0, 0);
+          }
+	  /*Reflectivity can't be >1*/
+          if (Ref>1) {
+            Ref = 1.0;
+          }
+	}
         /* apply reflectivity */
         p *= Ref;
         SCATTER;

From 7e5db7a2f773b874a20c345c8e5b50c530d29072 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Mon, 6 Jul 2026 12:54:39 +0200
Subject: [PATCH 09/17] Let these mirrors identify string NULL as 'no file'

---
 mcxtrace-comps/optics/Mirror_curved.comp   | 2 +-
 mcxtrace-comps/optics/Mirror_elliptic.comp | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/mcxtrace-comps/optics/Mirror_curved.comp b/mcxtrace-comps/optics/Mirror_curved.comp
index 0f8821c09..f5e0d6c27 100644
--- a/mcxtrace-comps/optics/Mirror_curved.comp
+++ b/mcxtrace-comps/optics/Mirror_curved.comp
@@ -56,7 +56,7 @@ INITIALIZE
 %{
   int status;
 
-  if (coating && strlen (coating)) {
+  if (coating && strlen (coating) && strcmp (coating, "NULL")) {
     status = reflec_Init (&re, COATING_UNDEFINED, coating, NULL);
   } else {
     /*assume a constant reflectivity*/
diff --git a/mcxtrace-comps/optics/Mirror_elliptic.comp b/mcxtrace-comps/optics/Mirror_elliptic.comp
index f68174c22..ee56051d7 100644
--- a/mcxtrace-comps/optics/Mirror_elliptic.comp
+++ b/mcxtrace-comps/optics/Mirror_elliptic.comp
@@ -74,7 +74,7 @@ INITIALIZE
   if (radius)
     x_a = y_b = z_c = radius;
 
-  if (coating && strlen (coating)) {
+  if (coating && strlen (coating) && strcmp (coating, "NULL")) {
     status = reflec_Init (&re, COATING_UNDEFINED, coating, NULL);
   } else {
     /*assume a constant reflectivity*/

From c35083688c58ed2ae9d67e2dc003075ca1cc51b0 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Mon, 6 Jul 2026 12:55:29 +0200
Subject: [PATCH 10/17] Interpret "NULL" as no file

---
 mcxtrace-comps/contrib/Mirror_toroid_pothole.comp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mcxtrace-comps/contrib/Mirror_toroid_pothole.comp b/mcxtrace-comps/contrib/Mirror_toroid_pothole.comp
index c41a7f116..b3837692e 100644
--- a/mcxtrace-comps/contrib/Mirror_toroid_pothole.comp
+++ b/mcxtrace-comps/contrib/Mirror_toroid_pothole.comp
@@ -62,7 +62,7 @@ DECLARE
 
 INITIALIZE
 %{
-  if (coating && strlen (coating)) {
+  if (coating && strlen (coating) && strcmp (coating, "NULL")) {
     char** header_parsed;
     t_Table* tp = &reflec_table;
     /* read 1st block data from file into tp */

From 904c1d7e394ca851179400355def44da9e67f31e Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Mon, 6 Jul 2026 12:56:40 +0200
Subject: [PATCH 11/17] Use new default mirror 'NULL'

---
 .../examples/Tests_optics/Test_Mirrors/Test_Mirrors.instr     | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/mcxtrace-comps/examples/Tests_optics/Test_Mirrors/Test_Mirrors.instr b/mcxtrace-comps/examples/Tests_optics/Test_Mirrors/Test_Mirrors.instr
index cd1da3618..8f8ab2bd3 100644
--- a/mcxtrace-comps/examples/Tests_optics/Test_Mirrors/Test_Mirrors.instr
+++ b/mcxtrace-comps/examples/Tests_optics/Test_Mirrors/Test_Mirrors.instr
@@ -55,7 +55,7 @@
 *******************************************************************************/
 
 /* Change name of instrument and input parameters with default values */
-DEFINE INSTRUMENT Test_Mirrors(gamma=5, int index=1, L=2, radius=1000, E0=12.5, string coating="None")
+DEFINE INSTRUMENT Test_Mirrors(gamma=5, int index=1, L=2, radius=1000, E0=12.5, string coating="NULL")
 
 INITIALIZE %{
   switch (index) {
@@ -68,8 +68,6 @@ INITIALIZE %{
     case 6: printf("Using Mirror_toroid_pothole\n"); break; // OK
     default: exit(fprintf(stderr, "Unknown Mirror index\n"));
   }
-  /* If "None" is given it in fact means "" */
-  if (!strcmp(coating,"None")) sprintf(coating,"");
 %}
 
 TRACE

From cabe2a52b13970fa7a6081476a76aaf2b8781bf8 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Mon, 6 Jul 2026 13:08:06 +0200
Subject: [PATCH 12/17] Drop using explicit filenames in this instr, use comp
 naming

---
 .../examples/MAXIV/MAXIV_FemtoMAX/MAXIV_FemtoMAX.instr        | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mcxtrace-comps/examples/MAXIV/MAXIV_FemtoMAX/MAXIV_FemtoMAX.instr b/mcxtrace-comps/examples/MAXIV/MAXIV_FemtoMAX/MAXIV_FemtoMAX.instr
index d17345eb8..2dd4efbc2 100644
--- a/mcxtrace-comps/examples/MAXIV/MAXIV_FemtoMAX/MAXIV_FemtoMAX.instr
+++ b/mcxtrace-comps/examples/MAXIV/MAXIV_FemtoMAX/MAXIV_FemtoMAX.instr
@@ -133,7 +133,7 @@ COMPONENT slit1 = Slit(
 AT(0,0,0.9355) RELATIVE a1
 
 COMPONENT popin_yag1 = PSD_monitor(
-    xwidth=1e-3, yheight=1e-3,filename="popin_yag1")
+    xwidth=1e-3, yheight=1e-3)
 AT(0,0, 0.3355) RELATIVE slit1
 
 
@@ -165,7 +165,7 @@ COMPONENT slit2 = Slit(
 	xwidth=dxs2, yheight=dys2)
 AT(0,0,1.2) RELATIVE focus_mirror_exit1
 
-COMPONENT popin_yag2 = COPY(popin_yag1)(yheight=1e-2,filename="popin_yag2")
+COMPONENT popin_yag2 = COPY(popin_yag1)(yheight=1e-2)
 AT(0,0,0.4) RELATIVE slit2
 
 COMPONENT mono = Bragg_crystal(

From dbd80670ae4215de09713754cd3ad90066850e52 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Mon, 6 Jul 2026 13:53:43 +0200
Subject: [PATCH 13/17] Switch these includes to %include "mccode-complex-lib"

---
 mcxtrace-comps/optics/Grating_reflect.comp | 2 +-
 mcxtrace-comps/optics/Mirror_elliptic.comp | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/mcxtrace-comps/optics/Grating_reflect.comp b/mcxtrace-comps/optics/Grating_reflect.comp
index 26e041d46..5aa7fc5ea 100644
--- a/mcxtrace-comps/optics/Grating_reflect.comp
+++ b/mcxtrace-comps/optics/Grating_reflect.comp
@@ -56,7 +56,7 @@ SETTING PARAMETERS (d_phi=1, R0=1, rho_l=800,
 
 SHARE
 %{  
-  #include <complex.h>
+  %include "mccode-complex-lib"
   %include "read_table-lib"
 %}
 
diff --git a/mcxtrace-comps/optics/Mirror_elliptic.comp b/mcxtrace-comps/optics/Mirror_elliptic.comp
index ee56051d7..27a2bf74e 100644
--- a/mcxtrace-comps/optics/Mirror_elliptic.comp
+++ b/mcxtrace-comps/optics/Mirror_elliptic.comp
@@ -54,7 +54,7 @@ DEPENDENCY "-std=c99"
 SHARE
 %{
   %include "perfect_crystals-lib"
-  #include <complex.h>
+  %include "mccode-complex-lib"
   %include "read_table-lib"
   %include "reflectivity-lib"
 %}

From e36c792091a4cd87b5ec7348e1255cdbf094e5b2 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Mon, 6 Jul 2026 13:56:55 +0200
Subject: [PATCH 14/17] Apply terms present in other mirrors

---
 mcxtrace-comps/optics/Multilayer_elliptic.comp | 13 ++-----------
 1 file changed, 2 insertions(+), 11 deletions(-)

diff --git a/mcxtrace-comps/optics/Multilayer_elliptic.comp b/mcxtrace-comps/optics/Multilayer_elliptic.comp
index 29ef643d3..fe798060e 100644
--- a/mcxtrace-comps/optics/Multilayer_elliptic.comp
+++ b/mcxtrace-comps/optics/Multilayer_elliptic.comp
@@ -214,17 +214,8 @@ TRACE
 	double Ref=R0;
         QQ = sqrt ((kx - kxo) * (kx - kxo) + (ky - kyo) * (ky - kyo) + (kz - kzo) * (kz - kzo));
 	if (!isconst) {
-          if (kinematical) {
-            /*
-             * \Lambda: thickness of bilayer - following notation in Als-Nielsen/McMorrow
-             * \Gamma:  \Gamma*\Lambda thickness of high electron density material.
-             * r1(zeta) = 2 i r_0 \rho_{AB} \left(\frac{\Lambda^2 \Gamma}{\zeta}\right) \frac{\sin\left(\pi\Gamma\zeta\right)}{\pi\Gamma\zeta);
-             */
-            Ref = reflecq (re, QQ, 0, 0, 0);
-          } else {
-            /*interpolate in table*/
-            Ref = reflecq (re, QQ, 0, 0, 0);
-          }
+          /*interpolate in table*/
+          Ref = reflecq (re, QQ, 0, k, fabs (90 - acos (s / k) * RAD2DEG));
 	  /*Reflectivity can't be >1*/
           if (Ref>1) {
             Ref = 1.0;

From ede0ab29e075de500cc6d1ce5045506cac77d3d2 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Tue, 7 Jul 2026 12:00:16 +0200
Subject: [PATCH 15/17] Mark NOACC as an experiment

---
 mcxtrace-comps/optics/Multilayer_elliptic.comp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/mcxtrace-comps/optics/Multilayer_elliptic.comp b/mcxtrace-comps/optics/Multilayer_elliptic.comp
index fe798060e..7a5b400ce 100644
--- a/mcxtrace-comps/optics/Multilayer_elliptic.comp
+++ b/mcxtrace-comps/optics/Multilayer_elliptic.comp
@@ -62,6 +62,7 @@ SETTING PARAMETERS (string coating="Ref_W_B4C.txt",
     theta=1.2,s1=0,s2=0,length=0.5,width=0.2,R0=1, 
     Emin=-1, Emax=-1, Estep=-1, Gamma=0, Lambda=0, rho_AB=0, int N=0,
     xwidth=0, zdepth=0)
+NOACC
 
 /* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */ 
 

From 3e07238a4f55c6e11bd25f853b97b3c66aab7361 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Tue, 7 Jul 2026 18:32:30 +0200
Subject: [PATCH 16/17] Candidate-fix for rectification of Multilayer_elliptic
 on GPU

---
 mcxtrace-comps/optics/Multilayer_elliptic.comp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/mcxtrace-comps/optics/Multilayer_elliptic.comp b/mcxtrace-comps/optics/Multilayer_elliptic.comp
index 7a5b400ce..6b2552d6e 100644
--- a/mcxtrace-comps/optics/Multilayer_elliptic.comp
+++ b/mcxtrace-comps/optics/Multilayer_elliptic.comp
@@ -62,7 +62,6 @@ SETTING PARAMETERS (string coating="Ref_W_B4C.txt",
     theta=1.2,s1=0,s2=0,length=0.5,width=0.2,R0=1, 
     Emin=-1, Emax=-1, Estep=-1, Gamma=0, Lambda=0, rho_AB=0, int N=0,
     xwidth=0, zdepth=0)
-NOACC
 
 /* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */ 
 
@@ -207,6 +206,8 @@ TRACE
         kyo = ky, kzo = kz;
         CoordTransEM (&kx, &ky, &kz, kxn, kyn, kzn, 0, 0, xi);
 
+	NORM (kx, ky, kz);
+
         kx = k * kx;
         ky = k * ky;
         kz = k * kz;

From 82fdc13ae03dafacdd409fc8a66c953b85bcf02a Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Tue, 7 Jul 2026 18:35:50 +0200
Subject: [PATCH 17/17] Formatting + removal of isconst var

---
 .../optics/Multilayer_elliptic.comp           | 22 ++++++++-----------
 1 file changed, 9 insertions(+), 13 deletions(-)

diff --git a/mcxtrace-comps/optics/Multilayer_elliptic.comp b/mcxtrace-comps/optics/Multilayer_elliptic.comp
index 6b2552d6e..20577c4b3 100644
--- a/mcxtrace-comps/optics/Multilayer_elliptic.comp
+++ b/mcxtrace-comps/optics/Multilayer_elliptic.comp
@@ -101,7 +101,6 @@ DECLARE
   double cost0;
   int kinematical;
   t_Reflec re;
-  int isconst;
 %}
 
 INITIALIZE
@@ -109,7 +108,6 @@ INITIALIZE
   /* calculation of the elliptical parameters according to the input mirror parameters:
   ellipse major axis a/2, minor axis b/2, M-magnification factor, Z0&Y0 - position of the mirror centre in the elliptical coordinate system.*/
   double Theta = DEG2RAD * theta;
-  isconst=0;
   if (xwidth)
     width = xwidth;
   if (zdepth)
@@ -133,7 +131,6 @@ INITIALIZE
       status = reflec_Init (&re, COATING_UNDEFINED, coating, NULL);
     } else {
       /*assume a constant reflectivity*/
-      isconst=1;
       status = reflec_Init (&re, CONSTANT, NULL, &(R0));
     }
   } else {
@@ -206,23 +203,22 @@ TRACE
         kyo = ky, kzo = kz;
         CoordTransEM (&kx, &ky, &kz, kxn, kyn, kzn, 0, 0, xi);
 
-	NORM (kx, ky, kz);
+        NORM (kx, ky, kz);
 
         kx = k * kx;
         ky = k * ky;
         kz = k * kz;
 
         double QQ;
-	double Ref=R0;
+        double Ref = R0;
         QQ = sqrt ((kx - kxo) * (kx - kxo) + (ky - kyo) * (ky - kyo) + (kz - kzo) * (kz - kzo));
-	if (!isconst) {
-          /*interpolate in table*/
-          Ref = reflecq (re, QQ, 0, k, fabs (90 - acos (s / k) * RAD2DEG));
-	  /*Reflectivity can't be >1*/
-          if (Ref>1) {
-            Ref = 1.0;
-          }
-	}
+        /* Evaluate reflectivity curve */
+        Ref = reflecq (re, QQ, 0, k, fabs (90 - acos (s / k) * RAD2DEG));
+        /*Reflectivity can't be >1*/
+        if (Ref > 1) {
+          Ref = 1.0;
+        }
+
         /* apply reflectivity */
         p *= Ref;
         SCATTER;