Skip to content

Add comp-formatter to McStas samples #2317

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
willend merged 1 commit into from
Feb 19, 2026
+6,462 −6,546
Merged

Add comp-formatter to McStas samples #2317

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
422 changes: 216 additions & 206 deletions mcstas-comps/samples/Incoherent.comp
View file
Open in desktop

Large diffs are not rendered by default.

4,578 changes: 2,295 additions & 2,283 deletions mcstas-comps/samples/Isotropic_Sqw.comp
View file
Open in desktop

Large diffs are not rendered by default.

806 changes: 380 additions & 426 deletions mcstas-comps/samples/Magnon_bcc.comp
View file
Open in desktop

Large diffs are not rendered by default.

290 changes: 140 additions & 150 deletions mcstas-comps/samples/NCrystal_sample.comp
View file
Open in desktop

Large diffs are not rendered by default.

635 changes: 300 additions & 335 deletions mcstas-comps/samples/Phonon_simple.comp
View file
Open in desktop

Large diffs are not rendered by default.

181 changes: 92 additions & 89 deletions mcstas-comps/samples/Powder1.comp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,143 +73,146 @@ pack= 1, j= 6, DW= 1, F2= 56.8, Vc= 85.0054, sigma_abs= 0.463)

DECLARE
%{
double my_s_v2;
double my_a_v;
double q_v;
char isrect;
double my_s_v2;
double my_a_v;
double q_v;
char isrect;
%}

INITIALIZE
%{
isrect=0;
isrect = 0;

if (yheight) yheight=yheight;
if (yheight)
yheight = yheight;
if (!radius || !yheight) {
if (!xwidth || !yheight || !zdepth) exit(fprintf(stderr,"Powder1: %s: sample has no volume (zero dimensions)\n", NAME_CURRENT_COMP));
else isrect=1; }
if (!xwidth || !yheight || !zdepth)
exit (fprintf (stderr, "Powder1: %s: sample has no volume (zero dimensions)\n", NAME_CURRENT_COMP));
else
isrect = 1;
}

my_a_v = pack*sigma_abs/Vc*2200*100; /* Is not yet divided by v */
my_s_v2 = 4*PI*PI*PI*pack*j*F2*DW/(Vc*Vc*V2K*V2K*q)*100;
my_a_v = pack * sigma_abs / Vc * 2200 * 100; /* Is not yet divided by v */
my_s_v2 = 4 * PI * PI * PI * pack * j * F2 * DW / (Vc * Vc * V2K * V2K * q) * 100;
/* Is not yet divided by v^2. 100: convert from barns to fm^2 */
/* Squires [3.103] */
if (d) q=2*PI/d;
q_v = q*K2V;
if (d)
q = 2 * PI / d;
q_v = q * K2V;
%}
TRACE
%{
double t0, t1, v, l_full, l, l_1, dt, dphi_in,d_phi0, theta, my_s;
double t0, t1, v, l_full, l, l_1, dt, dphi_in, d_phi0, theta, my_s;
double arg, tmp_vx, tmp_vy, tmp_vz, vout_x, vout_y, vout_z;
char intersect=0;
char intersect = 0;

dphi_in = d_phi;
if (isrect)
intersect = box_intersect(&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);
intersect = box_intersect (&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);
else
intersect = cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius, yheight);
if(intersect)
{
if(t0 < 0)
intersect = cylinder_intersect (&t0, &t1, x, y, z, vx, vy, vz, radius, yheight);
if (intersect) {
if (t0 < 0)
ABSORB;
/* Neutron enters at t=t0. */
v = sqrt(vx*vx + vy*vy + vz*vz);
l_full = v * (t1 - t0); /* Length of full path through sample */
dt = rand01()*(t1 - t0); /* Time of scattering */
PROP_DT(dt+t0); /* Point of scattering */
l = v*dt; /* Penetration in sample */
v = sqrt (vx * vx + vy * vy + vz * vz);
l_full = v * (t1 - t0); /* Length of full path through sample */
dt = rand01 () * (t1 - t0); /* Time of scattering */
PROP_DT (dt + t0); /* Point of scattering */
l = v * dt; /* Penetration in sample */

/* choose line theta */
arg = q_v/(2.0*v);
if(arg > 1)
ABSORB; /* No bragg scattering possible*/
theta = asin(arg); /* Bragg scattering law */

/* Choose point on Debye-Scherrer cone */
if (dphi_in)
{ /* relate height of detector to the height on DS cone */
arg = sin(dphi_in*DEG2RAD/2)/sin(2*theta);
if (arg < -1 || arg > 1) dphi_in = 0;
else dphi_in = 2*asin(arg);
}
if (dphi_in) {
dphi_in = fabs(dphi_in);
d_phi0= 2*rand01()*dphi_in;
if (d_phi0 > dphi_in) arg = 1; else arg = 0;
if (arg) {
d_phi0=PI+(d_phi0-1.5*dphi_in);
} else {
d_phi0=d_phi0-0.5*dphi_in;
}
p *= dphi_in/PI;
}
arg = q_v / (2.0 * v);
if (arg > 1)
ABSORB; /* No bragg scattering possible*/
theta = asin (arg); /* Bragg scattering law */

/* Choose point on Debye-Scherrer cone */
if (dphi_in) { /* relate height of detector to the height on DS cone */
arg = sin (dphi_in * DEG2RAD / 2) / sin (2 * theta);
if (arg < -1 || arg > 1)
dphi_in = 0;
else
dphi_in = 2 * asin (arg);
}
if (dphi_in) {
dphi_in = fabs (dphi_in);
d_phi0 = 2 * rand01 () * dphi_in;
if (d_phi0 > dphi_in)
arg = 1;
else
d_phi0 = PI*randpm1();
arg = 0;
if (arg) {
d_phi0 = PI + (d_phi0 - 1.5 * dphi_in);
} else {
d_phi0 = d_phi0 - 0.5 * dphi_in;
}
p *= dphi_in / PI;
} else
d_phi0 = PI * randpm1 ();

/* now find a nearly vertical rotation axis:
* (v along Z) x (X axis) -> nearly Y axis
*/
vec_prod(tmp_vx,tmp_vy,tmp_vz, vx,vy,vz, 1,0,0);
* (v along Z) x (X axis) -> nearly Y axis
*/
vec_prod (tmp_vx, tmp_vy, tmp_vz, vx, vy, vz, 1, 0, 0);

/* handle case where v and aim are parallel */
if (!tmp_vx && !tmp_vy && !tmp_vz) { tmp_vx=tmp_vz=0; tmp_vy=1; }
if (!tmp_vx && !tmp_vy && !tmp_vz) {
tmp_vx = tmp_vz = 0;
tmp_vy = 1;
}

/* v_out = rotate 'v' by 2*theta around tmp_v: Bragg angle */
rotate(vout_x,vout_y,vout_z, vx,vy,vz, 2*theta, tmp_vx,tmp_vy,tmp_vz);
rotate (vout_x, vout_y, vout_z, vx, vy, vz, 2 * theta, tmp_vx, tmp_vy, tmp_vz);

/* tmp_v = rotate v_out by d_phi0 around 'v' (Debye-Scherrer cone) */
rotate(tmp_vx,tmp_vy,tmp_vz, vout_x,vout_y,vout_z, d_phi0, vx, vy, vz);
rotate (tmp_vx, tmp_vy, tmp_vz, vout_x, vout_y, vout_z, d_phi0, vx, vy, vz);
vx = tmp_vx;
vy = tmp_vy;
vz = tmp_vz;

arg=0;
if (isrect && !box_intersect(&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zdepth)) arg=1;
else if(!isrect && !cylinder_intersect(&t0, &t1, x, y, z,
vx, vy, vz, radius, yheight)) arg=1;
arg = 0;
if (isrect && !box_intersect (&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zdepth))
arg = 1;
else if (!isrect && !cylinder_intersect (&t0, &t1, x, y, z, vx, vy, vz, radius, yheight))
arg = 1;

if (arg) {
/* Strange error: did not hit cylinder */
fprintf(stderr, "PowderN: FATAL ERROR: Did not hit sample from inside.\n");
fprintf (stderr, "PowderN: FATAL ERROR: Did not hit sample from inside.\n");
ABSORB;
}
l_1 = v*t1; /* go to exit */
l_1 = v * t1; /* go to exit */

my_s = my_s_v2/(v*v);
p *= l_full*my_s*exp(-(my_a_v/v+my_s)*(l+l_1));
my_s = my_s_v2 / (v * v);
p *= l_full * my_s * exp (-(my_a_v / v + my_s) * (l + l_1));
SCATTER;
}
%}

MCDISPLAY
%{

if (!isrect) {
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);
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 {
double xmin = -0.5*xwidth;
double xmax = 0.5*xwidth;
double ymin = -0.5*yheight;
double ymax = 0.5*yheight;
double zmin = -0.5*zdepth;
double zmax = 0.5*zdepth;
multiline(5, xmin, ymin, zmin,
xmax, ymin, zmin,
xmax, ymax, zmin,
xmin, ymax, zmin,
xmin, ymin, zmin);
multiline(5, xmin, ymin, zmax,
xmax, ymin, zmax,
xmax, ymax, zmax,
xmin, ymax, zmax,
xmin, ymin, zmax);
line(xmin, ymin, zmin, xmin, ymin, zmax);
line(xmax, ymin, zmin, xmax, ymin, zmax);
line(xmin, ymax, zmin, xmin, ymax, zmax);
line(xmax, ymax, zmin, xmax, ymax, zmax);
double xmin = -0.5 * xwidth;
double xmax = 0.5 * xwidth;
double ymin = -0.5 * yheight;
double ymax = 0.5 * yheight;
double zmin = -0.5 * zdepth;
double zmax = 0.5 * zdepth;
multiline (5, xmin, ymin, zmin, xmax, ymin, zmin, xmax, ymax, zmin, xmin, ymax, zmin, xmin, ymin, zmin);
multiline (5, xmin, ymin, zmax, xmax, ymin, zmax, xmax, ymax, zmax, xmin, ymax, zmax, xmin, ymin, zmax);
line (xmin, ymin, zmin, xmin, ymin, zmax);
line (xmax, ymin, zmin, xmax, ymin, zmax);
line (xmin, ymax, zmin, xmin, ymax, zmax);
line (xmax, ymax, zmin, xmax, ymax, zmax);
}
%}
END
Loading