Skip to content

Commit

Permalink
Merge pull request #1521 from McStasMcXtrace/Single_crystal_PG_powder…
Browse files Browse the repository at this point in the history 
…_fix

Single crystal PG powder fix
  • Loading branch information
@willend
willend committed Nov 21, 2023
2 parents 5782a94 + 71d88d2 commit 239b2a2
Show file tree
Hide file tree
Showing 2 changed files with 51 additions and 19 deletions.
38 changes: 33 additions & 5 deletions mcstas-comps/samples/Single_crystal.comp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1237,6 +1237,11 @@ 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 @@ -1326,10 +1331,16 @@ 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);
}

Expand Down Expand Up @@ -1447,9 +1458,9 @@ TRACE

if (force_transmit) {
/* Exit due to truncated order, weight with relevant cross-sections to distance l_full */
p*=exp(-abs_xlen*l_full);
p*=exp(-abs_xlen*l_full);
intersect=0;
break;
break;
}

/* (5). Transmission */
Expand All @@ -1471,10 +1482,10 @@ TRACE
PGderotate(&vx, &vy, &vz, Alpha, hkl_info.csx, hkl_info.csy, hkl_info.csz);
}

type = 't';
if (!itype) itype = 1;
type = 't';
if (!itype) itype = 1;
#ifndef OPENACC
hkl_info.type = type;
hkl_info.type = type;
#endif

break;
Expand Down Expand Up @@ -1510,9 +1521,26 @@ TRACE
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);
event_counter++;

if (PG || powder) {
/* In case of PG or powder return to crystalite frame after propagation */
vx = _vx;
vy = _vy;
vz = _vz;
}

/* (4). Account for the probability of sigma_abs */
p *= (coh_xlen + inc_xlen)/tot_xlen;
Expand Down
32 changes: 18 additions & 14 deletions mcstas-comps/union/Single_crystal_process.comp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -719,16 +719,18 @@ SHARE

/* Functions for "reorientation", powder and PG modes */
/* Powder, forward */
void randrotate_union(double *nx, double *ny, double *nz, double a, double b) {
double nvx, nvy, nvz;
rotate(nvx,nvy,nvz, *nx, *ny, *nz, a, 1, 0, 0);
rotate(*nx,*ny,*nz, nvx,nvy,nvz, b, 0, 1, 0);
void randrotate_union(double *nx, double *ny, double *nz, double a, double b, double c) {
double x1, y1, z1, x2, y2, z2;
rotate(x1, y1, z1, *nx,*ny,*nz, a, 1, 0, 0); /* <1> = rot(<n>,a) */
rotate(x2, y2, z2, x1, y1, z1, b, 0, 1, 0); /* <2> = rot(<1>,b) */
rotate(*nx,*ny,*nz, x2, y2, z2, c, 0, 0, 1); /* <n> = rot(<2>,c) */
}
/* Powder, back */
void randderotate_union(double *nx, double *ny, double *nz, double a, double b) {
double nvx, nvy, nvz;
rotate(nvx,nvy,nvz,*nx,*ny,*nz, -b, 0, 1, 0);
rotate(*nx, *ny, *nz, nvx,nvy,nvz, -a, 1, 0, 0);
void randderotate_union(double *nx, double *ny, double *nz, double a, double b, double c) {
double x1, y1, z1, x2, y2, z2;
rotate(x1, y1, z1, *nx,*ny,*nz, -c, 0,0,1);
rotate(x2, y2, z2, x1, y1, z1, -b, 0,1,0);
rotate(*nx,*ny,*nz, x2, y2, z2, -a, 1,0,0);
}
/* PG, forward */
void PGrotate_union(double *nx, double *ny, double *nz, double a, double csx, double csy, double csz) {
Expand Down Expand Up @@ -760,6 +762,7 @@ struct Single_crystal_physics_storage_struct{
double powder_setting; // 0 if powder mode is disabled, 1 if enabled. Values between approximates texture?
double Alpha; // random angle between 0 and 2*Pi*powder
double Beta; // random angle between -Pi/2 and Pi/2
double Gamma; // random angle between -Pi and Pi

struct hkl_info_struct_union *hkl_info_storage; // struct containing all necessary info for SC
double pack; // packing factor
Expand All @@ -775,19 +778,19 @@ int Single_crystal_physics_my(double *my, double *k_initial, union data_transfer

struct hkl_info_struct_union *hkl_info = data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->hkl_info_storage;

// Need the Powder/PG/curvature mode rotate added here

// Taken from Single_crystal and changed hkl_info to a pointer.
// The split optimization is less useful here than normally

if (data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->powder_setting) {
//orientation of crystallite is random
data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->Alpha = randpm1()*PI*data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->powder_setting;
data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->Beta = randpm1()*PI/2;
data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->Gamma = randpm1()*PI;

randrotate_union(&kix, &kiy, &kiz,
data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->Alpha,
data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->Beta);
data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->Beta,
data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->Gamma);
}
if (data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->PG_setting) {
// orientation of crystallite is random
Expand Down Expand Up @@ -831,8 +834,6 @@ int Single_crystal_physics_my(double *my, double *k_initial, union data_transfer
hkl_info->coh_xsect = coh_xsect;
hkl_info->nb_refl += hkl_info->tau_count;
hkl_info->nb_refl_count++;


}

/* (3). Probabilities of the different possible interactions. */
Expand Down Expand Up @@ -925,7 +926,10 @@ int Single_crystal_physics_scattering(double *k_final, double *k_initial, double

if (data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->powder_setting) {
// orientation of crystallite is no longer random
randderotate_union(&k_final[0], &k_final[1], &k_final[2], data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->Alpha, data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->Beta);
randderotate_union(&k_final[0], &k_final[1], &k_final[2],
data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->Alpha,
data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->Beta,
data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->Gamma);
}
if (data_transfer.pointer_to_a_Single_crystal_physics_storage_struct->PG_setting) {
// orientation of crystallite is longer random
Expand Down

0 comments on commit 239b2a2

Please sign in to comment.