Merge pull request #1738 from farhi/main

McStas: Fix #1732: Single_crystal powder mode was broken by #1521

Author: willend

mcstas-comps/examples/Tests_optics/Test_Monochromators/Test_Monochromators.instr

@@ -24,6 +24,7 @@
 * %Example: Mono=5 Detector: psd1_I=8.9e-05
 * %Example: Mono=6 Detector: psd1_I=1.0e-04
 * %Example: Mono=6 PG=0.15 Detector: psd1_I=2.8e-05
+* %Example: Mono=6 powder=0.15 Detector: psd1_I=3.3e-07
 * %Example: Mono=6 ay=2.13389 az=-1.232 bz=2.464 cx=6.711 Detector: psd1_I=1.0e-04
 * %Example: Mono=6 ay=2.94447 by=1.47224 bz=2.54999 cx=0.936252 recip=1 Detector: psd1_I=1.0e-04
 * %Example: Mono=7 ay=2.94447 by=1.47224 bz=2.54999 cx=0.936252 recip=1 Detector: psd1_I=9.6e-05

mcstas-comps/samples/Single_crystal.comp

@@ -598,24 +598,24 @@ double tau;                 /* Length of (tau_x, tau_y, tau_z) */
         (flag ? "INC" : SC_file), as, bs, cs, info->m_aa, info->m_bb, info->m_cc);
     } else {
       if (!info->recip) {
-	printf("Mode: Direct mode lattice\n");
+        printf("Mode: Direct mode lattice\n");
         printf("Single_crystal: %s structure a=[%g,%g,%g] b=[%g,%g,%g] c=[%g,%g,%g] ",
-	       (flag ? "INC" : SC_file), info->m_ax ,info->m_ay ,info->m_az,
-	       info->m_bx ,info->m_by ,info->m_bz,
-	       info->m_cx ,info->m_cy ,info->m_cz);
+               (flag ? "INC" : SC_file), info->m_ax ,info->m_ay ,info->m_az,
+               info->m_bx ,info->m_by ,info->m_bz,
+               info->m_cx ,info->m_cy ,info->m_cz);
       } else {
-	printf("Mode: Reciprocal mode lattice\n");
+        printf("Mode: Reciprocal mode lattice\n");
         printf("Single_crystal: %s structure a*=[%g,%g,%g] b*=[%g,%g,%g] c*=[%g,%g,%g] ",
-	       (flag ? "INC" : SC_file), info->m_ax ,info->m_ay ,info->m_az,
-	       info->m_bx ,info->m_by ,info->m_bz,
-	       info->m_cx ,info->m_cy ,info->m_cz);
+               (flag ? "INC" : SC_file), info->m_ax ,info->m_ay ,info->m_az,
+               info->m_bx ,info->m_by ,info->m_bz,
+               info->m_cx ,info->m_cy ,info->m_cz);
       }
     }
     /* Compute reciprocal or direct lattice vectors. */
     if (!info->recip) {
       vec_prod(tmp_x, tmp_y, tmp_z,
-	       info->m_bx, info->m_by, info->m_bz,
-	       info->m_cx, info->m_cy, info->m_cz);
+               info->m_bx, info->m_by, info->m_bz,
+               info->m_cx, info->m_cy, info->m_cz);
       info->V0 = fabs(scalar_prod(info->m_ax, info->m_ay, info->m_az, tmp_x, tmp_y, tmp_z));
       printf("V0=%g\n", info->V0);
 
@@ -642,8 +642,8 @@ double tau;                 /* Length of (tau_x, tau_y, tau_z) */
       info->csz = info->m_cz;
 
       vec_prod(tmp_x, tmp_y, tmp_z,
-	       info->bsx/(2*PI), info->bsy/(2*PI), info->bsz/(2*PI),
-	       info->csx/(2*PI), info->csy/(2*PI), info->csz/(2*PI));
+               info->bsx/(2*PI), info->bsy/(2*PI), info->bsz/(2*PI),
+               info->csx/(2*PI), info->csy/(2*PI), info->csz/(2*PI));
       info->V0 = 1/fabs(scalar_prod(info->asx/(2*PI), info->asy/(2*PI), info->asz/(2*PI), tmp_x, tmp_y, tmp_z));
       printf("V0=%g\n", info->V0);
 
@@ -1174,7 +1174,7 @@ INITIALIZE
 
   if (PG && !(order==1)) {
     fprintf(stderr,"Single_crystal: %s: PG mode means implicit choice of no multiple scattering!\n"
-	    "WARNING setting order=1\n", NAME_CURRENT_COMP);
+            "WARNING setting order=1\n", NAME_CURRENT_COMP);
     order=1;
   }
 
@@ -1241,11 +1241,6 @@ TRACE
   double _vy;
   double _vz;
 
-  double component_vx;          /* velocities not rotated by Powder / PG for non-rotated propagation*/
-  double component_vy;
-  double component_vz;
-    
-
   char   type;      /* type of last event: t=transmit,c=coherent or i=incoherent */
   int    itype;     /* type of last event: t=1,c=2 or i=3 */
 
@@ -1312,7 +1307,7 @@ TRACE
         intersect = sphere_intersect(&t1, &t2, x, y, z, vx, vy, vz, radius);
       #ifdef USE_OFF
       else if (hkl_info.shape == 3)
-	      intersect = off_intersect(&t1, &t2, NULL, NULL, x, y, z, vx, vy, vz, 0, 0, 0, thread_offdata );
+              intersect = off_intersect(&t1, &t2, NULL, NULL, x, y, z, vx, vy, vz, 0, 0, 0, thread_offdata );
       #endif
       if(!intersect || t2*v < -1e-9 || t1*v > 1e-9)
       {
@@ -1335,17 +1330,11 @@ TRACE
         Alpha = randpm1()*PI*powder;
         Beta  = randpm1()*PI/2;
         Gamma = randpm1()*PI;
-        component_vx = vx;
-        component_vy = vy;
-        component_vz = vz;
         randrotate(&vx, &vy, &vz, Alpha, Beta, Gamma);
       }
       if (PG) { /* orientation of crystallite is random along <c> axis */
-	Alpha = randpm1()*PI*PG;
-        component_vx = vx;
-        component_vy = vy;
-        component_vz = vz;
-	PGrotate(&vx, &vy, &vz, Alpha, hkl_info.csx, hkl_info.csy, hkl_info.csz);
+        Alpha = randpm1()*PI*PG;
+        PGrotate(&vx, &vy, &vz, Alpha, hkl_info.csx, hkl_info.csy, hkl_info.csz);
       }
 
 
@@ -1400,10 +1389,10 @@ TRACE
         && fabs(kiz - hkl_info.kiz) < deltak) {
         hkl_info.nb_reuses++;
 
-	/* Restore in case of matching event (e.g. SPLIT) */
-	coh_refl = hkl_info.coh_refl;
+        /* Restore in case of matching event (e.g. SPLIT) */
+        coh_refl = hkl_info.coh_refl;
         coh_xsect = hkl_info.coh_xsect;
-	tau_count = hkl_info.tau_count;
+        tau_count = hkl_info.tau_count;
 
       } else {
 #endif
@@ -1417,8 +1406,8 @@ TRACE
               kix, kiy, kiz, tau_max,
               &coh_refl, &coh_xsect, oclContext_SX,
               d_L, d_T, d_tau_count, d_coh_refl, d_coh_xsect);
-	  if (tau_count != 0)
-	    printf("\nGPU tau_count:%i\n",tau_count);
+          if (tau_count != 0)
+            printf("\nGPU tau_count:%i\n",tau_count);
         }
         else 
         #endif
@@ -1437,13 +1426,13 @@ TRACE
           hkl_info.kiy = kiy;
           hkl_info.kiz = kiz;
 
-	  /* Store for potential re-use (e.g. SPLIT) */
-	  hkl_info.coh_refl  = coh_refl;
-	  hkl_info.coh_xsect = coh_xsect;
-	  hkl_info.tau_count = tau_count;
-	  hkl_info.nb_refl += tau_count;
-	  hkl_info.nb_refl_count++;
-	}
+          /* Store for potential re-use (e.g. SPLIT) */
+          hkl_info.coh_refl  = coh_refl;
+          hkl_info.coh_xsect = coh_xsect;
+          hkl_info.tau_count = tau_count;
+          hkl_info.nb_refl += tau_count;
+          hkl_info.nb_refl_count++;
+        }
       }
 #endif
       /* (3). Probabilities of the different possible interactions. */
@@ -1461,10 +1450,10 @@ TRACE
       }
 
       if (force_transmit) {
-	/* Exit due to truncated order, weight with relevant cross-sections to distance l_full */
-	p*=exp(-abs_xlen*l_full);
+        /* Exit due to truncated order, weight with relevant cross-sections to distance l_full */
+        p*=exp(-abs_xlen*l_full);
         intersect=0; 
-	break;
+        break;
       }
 
       /* (5). Transmission */
@@ -1487,14 +1476,14 @@ TRACE
         }
 
         type = 't';
-	    if (!itype) itype = 1;
+            if (!itype) itype = 1;
         #ifndef OPENACC
-	      hkl_info.type = type;
+              hkl_info.type = type;
         #endif
 
         break; 
         /* This break means that we are leaving the while-loop, exiting the
-	   crystal by "tunneling". */
+           crystal by "tunneling". */
       }
 
       /* Scattering "proper", i.e. coh or incoh */
@@ -1524,16 +1513,6 @@ TRACE
         l = rand0max(l_full);
       else
         l = -log(1 - rand0max((1 - exp(-tot_xlen*l_full))))/tot_xlen;
- 
-      if (PG || powder) {
-        /* If PG or powder mode, we need to propagate the neutron in the component frame instead of the crystalite frame*/
-        _vx = vx;
-        _vy = vy;
-        _vz = vz;
-        vx = component_vx;
-        vy = component_vy;
-        vz = component_vz;
-      }
 
       /* Propagate to scattering point */
       PROP_DT(l/v);
@@ -1556,9 +1535,9 @@ TRACE
         vx = kix; /* ki vector is used as tmp var with norm v */
         vy = kiy;
         vz = kiz; /* Go for next scattering event */
-	
-	type = 'i';
-	if (!itype) itype = 2;
+        
+        type = 'i';
+        if (!itype) itype = 2;
 #ifndef OPENACC
         hkl_info.type = type;
 #endif
@@ -1601,9 +1580,9 @@ TRACE
         vx = K2V*(L[i].u1x*kfx + L[i].u2x*kfy + L[i].u3x*kfz);
         vy = K2V*(L[i].u1y*kfx + L[i].u2y*kfy + L[i].u3y*kfz);
         vz = K2V*(L[i].u1z*kfx + L[i].u2z*kfy + L[i].u3z*kfz);
-	
-	type = 'c';
-	if (!itype) itype = 3;
+        
+        type = 'c';
+        if (!itype) itype = 3;
 #ifndef OPENACC
         hkl_info.type = type;
         hkl_info.h    = L[i].h;
@@ -1636,7 +1615,7 @@ TRACE
         randderotate(&vx, &vy, &vz, Alpha, Beta, Gamma);
       }
       if (PG) { /* orientation of crystallite is longer random */
-	PGderotate(&vx, &vy, &vz, Alpha, hkl_info.csx, hkl_info.csy, hkl_info.csz);
+        PGderotate(&vx, &vy, &vz, Alpha, hkl_info.csx, hkl_info.csy, hkl_info.csz);
       }
       /* exit if multiple scattering order has been reached */
       if (order && event_counter >= order) { force_transmit=1; }
@@ -1695,15 +1674,15 @@ FINALLY
 
 MCDISPLAY
 %{
-  if (hkl_info.shape == 0) {	/* cylinder */
+  if (hkl_info.shape == 0) {        /* cylinder */
     circle("xz", 0,  yheight/2.0, 0, radius);
     circle("xz", 0, -yheight/2.0, 0, radius);
     line(-radius, -yheight/2.0, 0, -radius, +yheight/2.0, 0);
     line(+radius, -yheight/2.0, 0, +radius, +yheight/2.0, 0);
     line(0, -yheight/2.0, -radius, 0, +yheight/2.0, -radius);
     line(0, -yheight/2.0, +radius, 0, +yheight/2.0, +radius);
   }
-  else if (hkl_info.shape == 1) { 	/* box */
+  else if (hkl_info.shape == 1) {         /* box */
     double xmin = -0.5*xwidth;
     double xmax =  0.5*xwidth;
     double ymin = -0.5*yheight;
@@ -1725,12 +1704,12 @@ MCDISPLAY
     line(xmin, ymax, zmin, xmin, ymax, zmax);
     line(xmax, ymax, zmin, xmax, ymax, zmax);
   }
-  else if (hkl_info.shape == 2) {	/* sphere */
+  else if (hkl_info.shape == 2) {        /* sphere */
     circle("xy", 0,  0.0, 0, radius);
     circle("xz", 0,  0.0, 0, radius);
     circle("yz", 0,  0.0, 0, radius);
   }
-  else if (hkl_info.shape == 3) {	/* OFF file */
+  else if (hkl_info.shape == 3) {        /* OFF file */
     off_display(offdata);
   }
 %}

mcxtrace-comps/optics/Capillary.comp

@@ -32,7 +32,7 @@
 * length:   [m]   Length of the unbent mirror.
 * coating:  [str] Name of file containing the material data (i.e. f1 and f2) for the coating
 * R0:       [0-1] Fixed constant reflectivity
-* rtable:   [0/1] If nonzero, the coating file contains an E,theta paramterized matrix of raw reflectivities.
+* rtable:   [0/1] If nonzero, the coating file contains an E,theta parameterized matrix of raw reflectivities.
 * waviness: [rad] The momentaneous waviness is uniformly distributed in the range [-waviness,waviness].
 * longw:    [0/1] If non-zero, waviness is purely longitudinal in its nature.
 * %End