Skip to content

Commit

Permalink
Merge pull request #1738 from farhi/main
Browse files Browse the repository at this point in the history 
McStas: Fix #1732: Single_crystal powder mode was broken by #1521
  • Loading branch information
@willend
willend authored Oct 11, 2024
2 parents 1d63187 + 7219d96 commit bbf7e8b
Show file tree
Hide file tree
Showing 3 changed files with 47 additions and 67 deletions.
1 change: 1 addition & 0 deletions mcstas-comps/examples/Tests_optics/Test_Monochromators/Test_Monochromators.instr
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
111 changes: 45 additions & 66 deletions mcstas-comps/samples/Single_crystal.comp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand All @@ -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);

Expand Down Expand Up @@ -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;
}

Expand Down Expand Up @@ -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 */

Expand Down Expand Up @@ -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)
{
Expand All @@ -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);
}


Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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. */
Expand All @@ -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 */
Expand All @@ -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 */
Expand Down Expand Up @@ -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);
Expand All @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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; }
Expand Down Expand Up @@ -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;
Expand All @@ -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);
}
%}
Expand Down
2 changes: 1 addition & 1 deletion mcxtrace-comps/optics/Capillary.comp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit bbf7e8b

Please sign in to comment.