From 633926db85844fc5ce321ecc1167aa5e253346db Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Thu, 22 Oct 2026 17:35:31 +0200 Subject: [PATCH 1/3] cogen: generate _getcomp_type(), _getcomp_rotation() and MCCODE_COMP_TYPE_ COMP_GETPAR3(type, "name", par) lets a component write another instance's SETTING parameters, but needs the target's type at compile time and gives no way to check it at run time. The code generator now emits: - #define MCCODE_COMP_TYPE_ next to each component parameter typedef, so code can guard COMP_GETPAR3 on that type being present; - char* _getcomp_type(char* name): the type of a named instance, or NULL; - int _getcomp_rotation(char* name, Rotation r): its absolute rotation; - #define MCCODE_HAS_GETCOMP_TYPE, so components can detect the feature. No change to generated simulation behaviour. mccode-antlr needs the same additions to support components that use them. Co-Authored-By: Claude Opus 5.5 --- common/lib/share/mccode-r.h.in | 4 ++++ mccode/src/cogen.c.in | 35 ++++++++++++++++++++++++++++++++++ 2 files changed, 39 insertions(+) diff --git a/common/lib/share/mccode-r.h.in b/common/lib/share/mccode-r.h.in index 62613e21de..33465dbf13 100644 --- a/common/lib/share/mccode-r.h.in +++ b/common/lib/share/mccode-r.h.in @@ -535,6 +535,10 @@ int any_set(int n, ...); void* _getvar_parameters(char* compname); int _getcomp_index(char* compname); +/* type name and absolute rotation of a named instance (cogen'd); + * MCCODE_COMP_TYPE_ is defined for every component type present */ +char* _getcomp_type(char* compname); +int _getcomp_rotation(char* compname, Rotation r); /* Note: The two-stage approach to COMP_GETPAR is NOT redundant; without it, * after #define C sample, COMP_GETPAR(C,x) would refer to component C, not to diff --git a/mccode/src/cogen.c.in b/mccode/src/cogen.c.in index ca24eff9fc..c699c52e21 100644 --- a/mccode/src/cogen.c.in +++ b/mccode/src/cogen.c.in @@ -611,6 +611,8 @@ int cogen_comp_declare(struct comp_inst *comp) coutf("}; /* _struct_%s_parameters */", comp->def->name); coutf("typedef struct _struct_%s_parameters _class_%s_parameters;", comp->def->name, comp->def->name); + /* lets other components guard COMP_GETPAR3(type, ...) on the type being present */ + coutf("#define MCCODE_COMP_TYPE_%s 1", comp->def->name); cout(""); coutf("/* Parameters for component type '%s' */", comp->def->name); @@ -1201,6 +1203,37 @@ void cogen_getcompindex_fct(struct instr_def *instr){ cout(""); } +/* + * Generates the functions that return the type name and the absolute rotation + * of a named component. Used to check the type before COMP_GETPAR3 is applied. + */ +void cogen_getcomptype_fct(struct instr_def *instr){ + List_handle iter; + struct comp_inst *comp; + cout("char* _getcomp_type(char* compname)"); + cout("/* Returns the component type of a named instance, or NULL */"); + cout("{"); + iter = list_iterate(instr->complist); + while ((comp = list_next(iter)) != NULL) { + coutf(" if (!strcmp(compname, \"%s\")) return \"%s\";", comp->name, comp->def->name); + } + list_iterate_end(iter); + cout(" return NULL;"); + cout("}"); + cout(""); + cout("int _getcomp_rotation(char* compname, Rotation r)"); + cout("/* Copies the absolute rotation of a named instance into r; returns 0 if not found */"); + cout("{"); + iter = list_iterate(instr->complist); + while ((comp = list_next(iter)) != NULL) { + coutf(" if (!strcmp(compname, \"%s\")) { rot_copy(r, _%s_var._rotation_absolute); return 1; }", comp->name, comp->name); + } + list_iterate_end(iter); + cout(" return 0;"); + cout("}"); + cout(""); +} + void cogen_getdistance_fct(){ cout("double index_getdistance(int first_index, int second_index)"); cout("/* Calculate the distance two components from their indexes*/"); @@ -1314,6 +1347,7 @@ int cogen_decls(struct instr_def *instr) } cout("};"); cout("typedef struct _struct_instrument_parameters _class_instrument_parameters;"); + cout("#define MCCODE_HAS_GETCOMP_TYPE 1 /* _getcomp_type(), _getcomp_rotation(), MCCODE_COMP_TYPE_ */"); cout(""); /* check if we have GROUP, SPLIT and JUMP */ @@ -3303,6 +3337,7 @@ cogen(char *output_name, struct instr_def *instr) cogen_getcompindex_fct(instr); + cogen_getcomptype_fct(instr); embed_file("metadata-r.c"); // functions used to query and display instrument/component-defined metadata strings embed_file("mccode_main.c"); From af7509d06efd7e224d67cbee41ca8dee0056c113 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Thu, 22 Oct 2026 17:35:31 +0200 Subject: [PATCH 2/3] UB_matrix: import a UB matrix and inject it into a crystal component New component for McStas and McXtrace. Placed before a crystal component, it reads a UB matrix (a file in McCode or Mantid ISAW format, 9 numbers, or a cell with two reflections), writes the orientation into the named target's SETTING parameters before the target initialises, and exports the UB to a text file and to NeXus as the run-time value of its 'UB' parameter. Targets: Single_crystal, Single_crystal_process, Single_magnetic_crystal, Single_crystal_inelastic, FluoCrystal (direct lattice vectors ax..cz), and NCrystal_sample, NCrystal_process, Union_NCrystal_material (dir1/dir2 appended to cfg). The target type is looked up from the instance name. Convention: Q = 2*pi*R*UB*(h,k,l) in the target frame (z along the beam, y up), no 2*pi in UB, B in the Busing-Levy frame, R the goniometer rotation. Test_UB_matrix (McStas: Al Single_crystal and NCrystal_sample; McXtrace: Si Single_crystal) compares each input route with a directly oriented reference; all reproduce the reference intensity and Laue pattern. Co-Authored-By: Claude Opus 5.5 --- .../Test_UB_matrix/Test_UB_matrix.instr | 113 +++++ .../Tests_samples/Test_UB_matrix/UB_Al.mat | 8 + mcstas-comps/misc/UB_matrix.comp | 400 ++++++++++++++++++ .../Test_UB_matrix/Test_UB_matrix.instr | 79 ++++ mcxtrace-comps/misc/UB_matrix.comp | 395 +++++++++++++++++ 5 files changed, 995 insertions(+) create mode 100644 mcstas-comps/examples/Tests_samples/Test_UB_matrix/Test_UB_matrix.instr create mode 100644 mcstas-comps/examples/Tests_samples/Test_UB_matrix/UB_Al.mat create mode 100644 mcstas-comps/misc/UB_matrix.comp create mode 100644 mcxtrace-comps/examples/Tests_samples/Test_UB_matrix/Test_UB_matrix.instr create mode 100644 mcxtrace-comps/misc/UB_matrix.comp diff --git a/mcstas-comps/examples/Tests_samples/Test_UB_matrix/Test_UB_matrix.instr b/mcstas-comps/examples/Tests_samples/Test_UB_matrix/Test_UB_matrix.instr new file mode 100644 index 0000000000..69657009ff --- /dev/null +++ b/mcstas-comps/examples/Tests_samples/Test_UB_matrix/Test_UB_matrix.instr @@ -0,0 +1,113 @@ +/******************************************************************************* +* Instrument: Test_UB_matrix +* +* %I +* Written by: Peter Willendrup +* Date: October 2026 +* Origin: DTU +* %INSTRUMENT_SITE: Tests_samples +* +* Test of the UB_matrix component injecting a UB matrix into Single_crystal +* +* %D +* An Al single crystal with (1,1,0) along the beam and (0,0,1) up is oriented +* in several equivalent ways, and illuminated by a white beam: +* SAMPLE=0: Single_crystal reference, direct lattice vectors given +* SAMPLE=1: Single_crystal, UB_matrix with a UB matrix as 9 numbers (UBin) +* SAMPLE=2: Single_crystal, UB_matrix with the cell and two reflections +* SAMPLE=3: Single_crystal, UB_matrix reading a Mantid SaveIsawUB file (UB_Al.mat) +* SAMPLE=4: NCrystal_sample reference, dir1/dir2 given in cfg +* SAMPLE=5: NCrystal_sample, UB_matrix with UBin appending dir1/dir2 to cfg +* SAMPLE=0-3 give the same Laue pattern and intensity, as do SAMPLE=4-5. +* The goniometer angle A3 rotates all samples about the vertical axis. +* +* %Example: SAMPLE=0 Detector: Laue_I=1.2705e-09 +* %Example: SAMPLE=1 Detector: Laue_I=1.2705e-09 +* %Example: SAMPLE=2 Detector: Laue_I=1.2705e-09 +* %Example: SAMPLE=3 Detector: Laue_I=1.2705e-09 +* %Example: SAMPLE=4 Detector: Laue_I=5.6971e-09 +* %Example: SAMPLE=5 Detector: Laue_I=5.6971e-09 +* +* %P +* SAMPLE: [1] Orientation input: 0-3 Single_crystal (reference, UBin, cell, ISAW file), 4-5 NCrystal_sample (reference, UBin) +* A3: [deg] Goniometer rotation about the vertical axis +* +* %L +* A. T. Savici, Introduction to the use of the UB matrix (ORNL, 2011) +* +* %E +*******************************************************************************/ +DEFINE INSTRUMENT Test_UB_matrix(int SAMPLE=1, A3=0) + +TRACE + +COMPONENT Origin = Progress_bar() + AT (0,0,0) ABSOLUTE + +COMPONENT Source = Source_simple(radius=0.001, dist=1, focus_xw=0.002, focus_yh=0.002, lambda0=3, dlambda=2.5) + AT (0,0,0) RELATIVE Origin + +COMPONENT gonio = Arm() + AT (0,0,1) RELATIVE Origin + ROTATED (0,A3,0) RELATIVE Origin + +COMPONENT UB_numbers = UB_matrix(target="Sample_UBin", filename="UB_numbers.txt", + UBin={-0.17459426696, 0.17459426696, 0, 0, 0, 0.246913580247, 0.17459426696, 0.17459426696, 0}) + AT (0,0,0) RELATIVE gonio + +COMPONENT UB_cell = UB_matrix(target="Sample_cell", filename="UB_cell.txt", + lattice_a=4.05, lattice_b=4.05, lattice_c=4.05, + hkl1={1,1,0}, dir1={0,0,1}, hkl2={0,0,1}, dir2={0,1,0}) + AT (0,0,0) RELATIVE gonio + +COMPONENT UB_file = UB_matrix(target="Sample_file", filename="UB_file.txt", file="UB_Al.mat", format="isaw", + lattice_a=4.05, lattice_b=4.05, lattice_c=4.05) + AT (0,0,0) RELATIVE gonio + +COMPONENT UB_nc = UB_matrix(target="Sample_nc", filename="UB_nc.txt", + UBin={-0.17459426696, 0.17459426696, 0, 0, 0, 0.246913580247, 0.17459426696, 0.17459426696, 0}) + AT (0,0,0) RELATIVE gonio + +COMPONENT Sample_ref = Single_crystal(xwidth=0.005, yheight=0.005, zdepth=0.005, reflections="Al.lau", mosaic=5, + ax=-2.863782463806, ay=0, az=2.863782463806, bx=2.863782463806, by=0, bz=2.863782463806, cx=0, cy=4.05, cz=0) + WHEN (SAMPLE==0) AT (0,0,0) RELATIVE gonio +EXTEND %{ + if (!SCATTERED) ABSORB; +%} + +COMPONENT Sample_UBin = Single_crystal(xwidth=0.005, yheight=0.005, zdepth=0.005, reflections="Al.lau", mosaic=5) + WHEN (SAMPLE==1) AT (0,0,0) RELATIVE gonio +EXTEND %{ + if (!SCATTERED) ABSORB; +%} + +COMPONENT Sample_cell = Single_crystal(xwidth=0.005, yheight=0.005, zdepth=0.005, reflections="Al.lau", mosaic=5) + WHEN (SAMPLE==2) AT (0,0,0) RELATIVE gonio +EXTEND %{ + if (!SCATTERED) ABSORB; +%} + +COMPONENT Sample_file = Single_crystal(xwidth=0.005, yheight=0.005, zdepth=0.005, reflections="Al.lau", mosaic=5) + WHEN (SAMPLE==3) AT (0,0,0) RELATIVE gonio +EXTEND %{ + if (!SCATTERED) ABSORB; +%} + +COMPONENT Sample_nc_ref = NCrystal_sample(xwidth=0.005, yheight=0.005, zdepth=0.005, + cfg="Al_sg225.ncmat;mos=5arcmin;incoh_elas=0;inelas=0;dir1=@crys_hkl:1,1,0@lab:0,0,1;dir2=@crys_hkl:0,0,1@lab:0,1,0") + WHEN (SAMPLE==4) AT (0,0,0) RELATIVE gonio +EXTEND %{ + if (!SCATTERED) ABSORB; +%} + +COMPONENT Sample_nc = NCrystal_sample(xwidth=0.005, yheight=0.005, zdepth=0.005, + cfg="Al_sg225.ncmat;mos=5arcmin;incoh_elas=0;inelas=0") + WHEN (SAMPLE==5) AT (0,0,0) RELATIVE gonio +EXTEND %{ + if (!SCATTERED) ABSORB; +%} + +COMPONENT Laue = PSD_monitor_4PI(radius=0.5, nx=180, ny=90, filename="Laue") + AT (0,0,0) RELATIVE gonio + +END diff --git a/mcstas-comps/examples/Tests_samples/Test_UB_matrix/UB_Al.mat b/mcstas-comps/examples/Tests_samples/Test_UB_matrix/UB_Al.mat new file mode 100644 index 0000000000..23c3e83bfe --- /dev/null +++ b/mcstas-comps/examples/Tests_samples/Test_UB_matrix/UB_Al.mat @@ -0,0 +1,8 @@ +0.174594266960 -0.174594266960 0.000000000000 +0.174594266960 0.174594266960 0.000000000000 +0.000000000000 0.000000000000 0.246913580247 +4.050000 4.050000 4.050000 90.000000 90.000000 90.000000 66.430125 +0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 + +Al, (1,1,0) along the beam and (0,0,1) up, in Mantid SaveIsawUB layout: +the first three rows are UB transposed in the IPNS frame (x along beam, z up). diff --git a/mcstas-comps/misc/UB_matrix.comp b/mcstas-comps/misc/UB_matrix.comp new file mode 100644 index 0000000000..3d6a713595 --- /dev/null +++ b/mcstas-comps/misc/UB_matrix.comp @@ -0,0 +1,400 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright (C) 1997-2026, All rights reserved +* DTU Physics, Kgs. Lyngby, Denmark +* +* Component: UB_matrix +* +* %I +* Written by: Peter Willendrup +* Date: October 2026 +* Origin: DTU +* +* Import a UB matrix, inject it into a crystal component, and export it +* +* %D +* Reads a UB matrix and writes the corresponding crystal orientation into the +* SETTING parameters of the crystal component named by 'target', before that +* component initialises. The UB matrix is also written to a text file and is +* stored in the NeXus output as the run-time value of the 'UB' parameter. +* +* Convention: Q = 2*pi*R*UB*(h,k,l), with UB in the frame of the target +* component (z along the beam, y up at zero goniometer angles), no 2*pi in UB +* (|UB*(h,k,l)| = 1/d), B in the Busing-Levy frame (a* along x), and +* R the active rotation of the 'goniometer' component (default: the target). +* +* The UB matrix comes from exactly one of: +* - file: 3 rows of 3 numbers. format="mccode" reads UB in the convention above; +* format="isaw" reads a Mantid SaveIsawUB file (UB transposed, IPNS frame +* with x along the beam and z up). +* - UBin: 9 numbers, row-major, e.g. from Mantid getUB(). +* - lattice_a..lattice_cc with hkl1/dir1 and hkl2/dir2: two reflections whose +* directions in the target frame are given (dir1 is kept exactly, dir2 only +* fixes the plane), as for the NCrystal dir1/dir2 options. +* When a UB is read and lattice_a is also given, the cell of the UB is checked +* against the given cell. +* +* Supported targets (the type is looked up from the instance name): +* Single_crystal, Single_crystal_process, Single_magnetic_crystal, +* Single_crystal_inelastic: the direct lattice vectors ax..cz. +* NCrystal_sample, NCrystal_process, Union_NCrystal_material: dir1/dir2 +* appended to cfg. +* The component must be placed BEFORE its target in the instrument, and the +* target must not already define an orientation. For Union processes, place +* the process AT (0,0,0) ABSOLUTE without rotation and set 'goniometer' to the +* Union volume. +* +* Example: UB_matrix(target="Sample", file="sample.mat", format="isaw") +* +* %P +* INPUT PARAMETERS: +* target: [str] Name of the crystal component to inject into +* file: [str] UB matrix file to read, 0 for none +* format: [str] File format, "mccode" or "isaw" +* UBin: [AA^-1] UB matrix, 9 numbers row-major, used when no file is given +* lattice_a: [AA] Lattice parameter a; with hkl1..dir2 defines the UB, or checks a read UB +* lattice_b: [AA] Lattice parameter b +* lattice_c: [AA] Lattice parameter c +* lattice_aa: [deg] Lattice angle alpha +* lattice_bb: [deg] Lattice angle beta +* lattice_cc: [deg] Lattice angle gamma +* hkl1: [1] Miller indices of the first reflection +* dir1: [1] Direction of hkl1 in the target frame +* hkl2: [1] Miller indices of the second reflection +* dir2: [1] Direction in the target frame fixing the plane of hkl2 +* tol: [rad] Tolerance on the hkl1-hkl2 opening angle +* goniometer: [str] Component whose rotation is R, 0 for the target +* filename: [str] Output file name, 0 for none +* +* OUTPUT PARAMETERS: +* UB: [AA^-1] The UB matrix, 9 numbers row-major (stored in NeXus) +* +* %E +*******************************************************************************/ + +DEFINE COMPONENT UB_matrix + +SETTING PARAMETERS (string target="Sample", string file=0, string format="mccode", + vector UBin={0,0,0,0,0,0,0,0,0}, + lattice_a=0, lattice_b=0, lattice_c=0, lattice_aa=90, lattice_bb=90, lattice_cc=90, + vector hkl1={1,0,0}, vector dir1={0,0,1}, vector hkl2={0,1,0}, vector dir2={1,0,0}, + tol=1e-4, string goniometer=0, string filename="UB.txt", + vector UB={0,0,0,0,0,0,0,0,0}) + +SHARE +%{ +%include "read_table-lib" +#ifndef UB_MATRIX_DECL +#define UB_MATRIX_DECL + /* 3x3 matrices are row-major double[9] */ + void ub_mul(double *A, double *B, double *C) { /* C = A*B */ + double T[9]; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + T[3*i+j] = A[3*i]*B[j] + A[3*i+1]*B[3+j] + A[3*i+2]*B[6+j]; + for (int i = 0; i < 9; i++) C[i] = T[i]; + } + void ub_transpose(double *A, double *T) { + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) T[3*j+i] = A[3*i+j]; + } + double ub_det(double *A) { + return A[0]*(A[4]*A[8]-A[5]*A[7]) - A[1]*(A[3]*A[8]-A[5]*A[6]) + A[2]*(A[3]*A[7]-A[4]*A[6]); + } + int ub_inv(double *A, double *I) { /* returns 0 when singular */ + double d = ub_det(A); + if (fabs(d) < 1e-300) return 0; + I[0] = (A[4]*A[8]-A[5]*A[7])/d; I[1] = (A[2]*A[7]-A[1]*A[8])/d; I[2] = (A[1]*A[5]-A[2]*A[4])/d; + I[3] = (A[5]*A[6]-A[3]*A[8])/d; I[4] = (A[0]*A[8]-A[2]*A[6])/d; I[5] = (A[2]*A[3]-A[0]*A[5])/d; + I[6] = (A[3]*A[7]-A[4]*A[6])/d; I[7] = (A[1]*A[6]-A[0]*A[7])/d; I[8] = (A[0]*A[4]-A[1]*A[3])/d; + return 1; + } + /* B from the cell (Savici eq. 52, Busing-Levy), no 2pi; returns 0 for an impossible cell */ + int ub_B(double a, double b, double c, double al, double be, double ga, double *B) { + double G[9], Gs[9]; + al *= DEG2RAD; be *= DEG2RAD; ga *= DEG2RAD; + G[0] = a*a; G[1] = a*b*cos(ga); G[2] = a*c*cos(be); + G[3] = G[1]; G[4] = b*b; G[5] = b*c*cos(al); + G[6] = G[2]; G[7] = G[5]; G[8] = c*c; + if (a <= 0 || b <= 0 || c <= 0 || ub_det(G) <= 0 || !ub_inv(G, Gs)) return 0; + double as = sqrt(Gs[0]), bs = sqrt(Gs[4]), cs = sqrt(Gs[8]); + double cbs = Gs[2]/(as*cs), cgs = Gs[1]/(as*bs); + B[0] = as; B[1] = bs*cgs; B[2] = cs*cbs; + B[3] = 0; B[4] = bs*sqrt(1 - cgs*cgs); B[5] = -cs*sqrt(1 - cbs*cbs)*cos(al); + B[6] = 0; B[7] = 0; B[8] = 1/c; + return 1; + } + /* orthonormal triad as columns of T: t1 along v1, t3 along v1 x v2; returns 0 if degenerate */ + int ub_triad(double *v1, double *v2, double *T) { + double t1[3], t2[3], t3[3]; + double n1 = sqrt(v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2]); + vec_prod(t3[0], t3[1], t3[2], v1[0], v1[1], v1[2], v2[0], v2[1], v2[2]); + double n3 = sqrt(t3[0]*t3[0] + t3[1]*t3[1] + t3[2]*t3[2]); + if (n1 <= 0 || n3 <= 1e-12*n1*sqrt(v2[0]*v2[0] + v2[1]*v2[1] + v2[2]*v2[2])) return 0; + for (int i = 0; i < 3; i++) { t1[i] = v1[i]/n1; t3[i] /= n3; } + vec_prod(t2[0], t2[1], t2[2], t3[0], t3[1], t3[2], t1[0], t1[1], t1[2]); + for (int i = 0; i < 3; i++) { T[3*i] = t1[i]; T[3*i+1] = t2[i]; T[3*i+2] = t3[i]; } + return 1; + } + double ub_angle(double *u, double *v) { + double d = u[0]*v[0] + u[1]*v[1] + u[2]*v[2]; + double n = sqrt((u[0]*u[0] + u[1]*u[1] + u[2]*u[2])*(v[0]*v[0] + v[1]*v[1] + v[2]*v[2])); + return acos(fmax(-1, fmin(1, d/n))); + } + /* U from two reflections (Savici eqs. 76-81); returns 0 on failure */ + int ub_U2(double *B, double *h1, double *d1, double *h2, double *d2, double tol, double *U) { + double v1[3], v2[3], Tc[9], Tn[9], TcT[9]; + for (int i = 0; i < 3; i++) { + v1[i] = B[3*i]*h1[0] + B[3*i+1]*h1[1] + B[3*i+2]*h1[2]; + v2[i] = B[3*i]*h2[0] + B[3*i+1]*h2[1] + B[3*i+2]*h2[2]; + } + if (!ub_triad(v1, v2, Tc) || !ub_triad(d1, d2, Tn)) return 0; + if (fabs(ub_angle(v1, v2) - ub_angle(d1, d2)) > tol) return 0; + ub_transpose(Tc, TcT); + ub_mul(Tn, TcT, U); + return 1; + } + /* B by Cholesky of (UB)^T(UB) (Savici eq. 88), then U = UB*B^-1; returns 0 if not positive definite */ + int ub_split(double *UB, double *U, double *B) { + double UBT[9], M[9], Bi[9]; + ub_transpose(UB, UBT); + ub_mul(UBT, UB, M); + if (M[0] <= 0) return 0; + B[0] = sqrt(M[0]); B[1] = M[1]/B[0]; B[2] = M[2]/B[0]; + double r = M[4] - B[1]*B[1]; + if (r <= 0) return 0; + B[4] = sqrt(r); B[5] = (M[5] - B[1]*B[2])/B[4]; + r = M[8] - B[2]*B[2] - B[5]*B[5]; + if (r <= 0) return 0; + B[8] = sqrt(r); B[3] = B[6] = B[7] = 0; + if (!ub_inv(B, Bi)) return 0; + ub_mul(UB, Bi, U); + return 1; + } + void ub_fprint(FILE *f, char *title, double *M) { + fprintf(f, "# %s\n", title); + for (int i = 0; i < 3; i++) fprintf(f, "%.12g %.12g %.12g\n", M[3*i], M[3*i+1], M[3*i+2]); + } +#endif +%} + +DECLARE +%{ + double U_out[9]; + double B_out[9]; + double R_out[9]; + double cell[6]; +%} + +INITIALIZE +%{ +#ifndef MCCODE_HAS_GETCOMP_TYPE +#error "UB_matrix needs a code generator that provides _getcomp_type() and MCCODE_COMP_TYPE_ (the C code generator in this McCode release does; mccode-antlr may not yet)." +#endif + double M[9], Bc[9], Uc[9], Ainv[9], A[9]; + int have_file = strlen(file) && strcmp(file, "0") && strcmp(file, "NULL"); + int have_ubin = 0; + for (int i = 0; i < 9; i++) if (UBin[i]) have_ubin = 1; + + if (have_file + have_ubin + (!have_file && !have_ubin && lattice_a > 0) != 1) { + fprintf(stderr, "UB_matrix: %s: ERROR: give exactly one of file, UBin, or lattice_a..lattice_cc with hkl1/dir1/hkl2/dir2.\n", NAME_CURRENT_COMP); + exit(-1); + } + + if (have_file) { + /* the first 3 lines that are not blank or '#' comments; ISAW files carry more lines after them */ + FILE *fp = Open_File(file, "r", NULL); + char line[1024]; + int rows = 0; + while (fp && rows < 3 && fgets(line, sizeof(line), fp)) { + char *s = line + strspn(line, " \t"); + if (*s == '#' || *s == '\n' || *s == '\r' || *s == '\0') continue; + if (sscanf(s, "%lf %lf %lf", &M[3*rows], &M[3*rows+1], &M[3*rows+2]) != 3) break; + rows++; + } + if (fp) fclose(fp); + if (rows < 3) { + fprintf(stderr, "UB_matrix: %s: ERROR: could not read 3 rows of 3 numbers from %s.\n", NAME_CURRENT_COMP, file); + exit(-1); + } + if (!strcmp(format, "isaw")) { + /* file holds UB^T in the IPNS frame: (x,y,z)_McStas = (y,z,x)_IPNS */ + double P[9] = {0,1,0, 0,0,1, 1,0,0}, MT[9]; + ub_transpose(M, MT); + ub_mul(P, MT, UB); + } else if (!strcmp(format, "mccode")) { + for (int i = 0; i < 9; i++) UB[i] = M[i]; + } else { + fprintf(stderr, "UB_matrix: %s: ERROR: unknown format \"%s\" (use \"mccode\" or \"isaw\").\n", NAME_CURRENT_COMP, format); + exit(-1); + } + } else if (have_ubin) { + for (int i = 0; i < 9; i++) UB[i] = UBin[i]; + } else { + if (!ub_B(lattice_a, lattice_b, lattice_c, lattice_aa, lattice_bb, lattice_cc, Bc)) { + fprintf(stderr, "UB_matrix: %s: ERROR: impossible lattice parameters.\n", NAME_CURRENT_COMP); + exit(-1); + } + if (!ub_U2(Bc, hkl1, dir1, hkl2, dir2, tol, Uc)) { + fprintf(stderr, "UB_matrix: %s: ERROR: hkl1/hkl2 or dir1/dir2 are parallel, or their opening angles differ by more than tol=%g rad.\n", + NAME_CURRENT_COMP, tol); + exit(-1); + } + ub_mul(Uc, Bc, UB); + } + + /* every UB passes the same checks */ + if (!ub_split(UB, U_out, B_out)) { + fprintf(stderr, "UB_matrix: %s: ERROR: the UB matrix does not describe a crystal ((UB)^T(UB) not positive definite).\n", NAME_CURRENT_COMP); + exit(-1); + } + if (ub_det(UB) < 0) { + fprintf(stderr, "UB_matrix: %s: ERROR: the UB matrix has a negative determinant (improper U). " + "It was probably written with the opposite Q sign convention; negate it.\n", NAME_CURRENT_COMP); + exit(-1); + } + { /* cell from B^T B = G*, G = (G*)^-1 */ + double BT[9], Gs[9], G[9]; + ub_transpose(B_out, BT); ub_mul(BT, B_out, Gs); ub_inv(Gs, G); + cell[0] = sqrt(G[0]); cell[1] = sqrt(G[4]); cell[2] = sqrt(G[8]); + cell[3] = acos(G[5]/(cell[1]*cell[2]))*RAD2DEG; + cell[4] = acos(G[2]/(cell[0]*cell[2]))*RAD2DEG; + cell[5] = acos(G[1]/(cell[0]*cell[1]))*RAD2DEG; + } + if ((have_file || have_ubin) && lattice_a > 0) { + double given[6] = {lattice_a, lattice_b, lattice_c, lattice_aa, lattice_bb, lattice_cc}; + for (int i = 0; i < 6; i++) + if (fabs(cell[i] - given[i]) > 1e-4*fabs(given[i])) { + fprintf(stderr, "UB_matrix: %s: ERROR: the UB cell (%g %g %g %g %g %g) does not match the given lattice parameters.\n", + NAME_CURRENT_COMP, cell[0], cell[1], cell[2], cell[3], cell[4], cell[5]); + exit(-1); + } + } + + /* goniometer rotation: McStas rotations are passive, R is their transpose */ + char *gonio = (strlen(goniometer) && strcmp(goniometer, "0") && strcmp(goniometer, "NULL")) ? goniometer : target; + Rotation Rp; + if (!_getcomp_rotation(gonio, Rp)) { + fprintf(stderr, "UB_matrix: %s: ERROR: no component named \"%s\".\n", NAME_CURRENT_COMP, gonio); + exit(-1); + } + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) R_out[3*i+j] = Rp[j][i]; + + /* inject into the target */ + char *type = _getcomp_type(target); + if (!type) { + fprintf(stderr, "UB_matrix: %s: ERROR: no component named \"%s\".\n", NAME_CURRENT_COMP, target); + exit(-1); + } + if (_getcomp_index(target) <= _getcomp_index(NAME_CURRENT_COMP)) { + fprintf(stderr, "UB_matrix: %s: ERROR: must be placed before its target \"%s\", which would otherwise already be initialised.\n", + NAME_CURRENT_COMP, target); + exit(-1); + } + /* direct lattice vectors a, b, c are the rows of UB^-1 (A^T UB = I) */ + ub_inv(UB, Ainv); + ub_transpose(Ainv, A); /* columns a, b, c */ + int injected = 0; + + /* direct lattice vectors, as COMP_GETPAR3(type, target, ax) etc. */ + #define UB_INJECT_ABC(T) \ + if (!strcmp(type, #T)) { \ + _class_##T##_parameters *p = (_class_##T##_parameters *) _getvar_parameters(target); \ + if (p->ax || p->ay || p->az || p->bx || p->by || p->bz || p->cx || p->cy || p->cz) { \ + fprintf(stderr, "UB_matrix: %s: ERROR: %s \"%s\" already has lattice vectors ax..cz; remove them.\n", NAME_CURRENT_COMP, type, target); \ + exit(-1); \ + } \ + p->ax = A[0]; p->ay = A[3]; p->az = A[6]; \ + p->bx = A[1]; p->by = A[4]; p->bz = A[7]; \ + p->cx = A[2]; p->cy = A[5]; p->cz = A[8]; \ + p->recip_cell = 0; \ + injected = 1; \ + } + /* NCrystal orientation, appended to cfg: lab directions of (100) and (010) */ + #define UB_INJECT_CFG(T) \ + if (!strcmp(type, #T)) { \ + char *cfgp = ((_class_##T##_parameters *) _getvar_parameters(target))->cfg; \ + if (strstr(cfgp, "dir1") || strstr(cfgp, "dir2")) { \ + fprintf(stderr, "UB_matrix: %s: ERROR: %s \"%s\" cfg already has dir1/dir2; remove them.\n", NAME_CURRENT_COMP, type, target); \ + exit(-1); \ + } \ + char dirs[512]; \ + snprintf(dirs, sizeof(dirs), ";dir1=@crys_hkl:1,0,0@lab:%.17g,%.17g,%.17g;dir2=@crys_hkl:0,1,0@lab:%.17g,%.17g,%.17g", \ + UB[0], UB[3], UB[6], UB[1], UB[4], UB[7]); \ + if (strlen(cfgp) + strlen(dirs) >= 16384) { \ + fprintf(stderr, "UB_matrix: %s: ERROR: cfg of \"%s\" too long to append dir1/dir2.\n", NAME_CURRENT_COMP, target); \ + exit(-1); \ + } \ + strcat(cfgp, dirs); \ + injected = 1; \ + } + + #ifdef MCCODE_COMP_TYPE_Single_crystal + UB_INJECT_ABC(Single_crystal) + #endif + #ifdef MCCODE_COMP_TYPE_Single_crystal_process + UB_INJECT_ABC(Single_crystal_process) + #endif + #ifdef MCCODE_COMP_TYPE_Single_magnetic_crystal + UB_INJECT_ABC(Single_magnetic_crystal) + #endif + #ifdef MCCODE_COMP_TYPE_Single_crystal_inelastic + UB_INJECT_ABC(Single_crystal_inelastic) + #endif + #ifdef MCCODE_COMP_TYPE_FluoCrystal + UB_INJECT_ABC(FluoCrystal) + #endif + #ifdef MCCODE_COMP_TYPE_NCrystal_sample + UB_INJECT_CFG(NCrystal_sample) + #endif + #ifdef MCCODE_COMP_TYPE_NCrystal_process + UB_INJECT_CFG(NCrystal_process) + #endif + #ifdef MCCODE_COMP_TYPE_Union_NCrystal_material + UB_INJECT_CFG(Union_NCrystal_material) + #endif + #undef UB_INJECT_ABC + #undef UB_INJECT_CFG + + if (!injected) { + fprintf(stderr, "UB_matrix: %s: ERROR: target \"%s\" is a %s, which UB_matrix does not support.\n", NAME_CURRENT_COMP, target, type); + exit(-1); + } + + MPI_MASTER( + printf("UB_matrix: %s: injected into %s \"%s\"; cell %g %g %g AA, %g %g %g deg\n", NAME_CURRENT_COMP, type, target, + cell[0], cell[1], cell[2], cell[3], cell[4], cell[5]); + printf("UB_matrix: %s: UB = [%g %g %g; %g %g %g; %g %g %g] 1/AA (no 2pi)\n", NAME_CURRENT_COMP, + UB[0], UB[1], UB[2], UB[3], UB[4], UB[5], UB[6], UB[7], UB[8]); + ); +%} + +SAVE +%{ + if (strlen(filename) && strcmp(filename, "0") && strcmp(filename, "NULL") && !mcdisable_output_files) { + MPI_MASTER( + char *path = mcfull_file(filename, NULL); + FILE *f = fopen(path, "w"); + if (!f) fprintf(stderr, "UB_matrix: %s: WARNING: could not write %s\n", NAME_CURRENT_COMP, path); + else { + double RUB[9]; + ub_mul(R_out, UB, RUB); + fprintf(f, "# UB matrix written by McCode UB_matrix component %s\n", NAME_CURRENT_COMP); + fprintf(f, "# target: %s\n", target); + fprintf(f, "# convention: Q = 2*pi*R*UB*(h,k,l); frame: z along beam, y up; no 2*pi in UB; Q = kf - ki; B Busing-Levy (a* along x)\n"); + fprintf(f, "# cell: %.10g %.10g %.10g %.10g %.10g %.10g\n", cell[0], cell[1], cell[2], cell[3], cell[4], cell[5]); + ub_fprint(f, "UB [1/AA] (read back with format=\"mccode\")", UB); + ub_fprint(f, "U", U_out); + ub_fprint(f, "B [1/AA]", B_out); + ub_fprint(f, "R (goniometer, active)", R_out); + ub_fprint(f, "R*UB [1/AA] (lab frame)", RUB); + fclose(f); + } + free(path); + ); + } +%} + +END diff --git a/mcxtrace-comps/examples/Tests_samples/Test_UB_matrix/Test_UB_matrix.instr b/mcxtrace-comps/examples/Tests_samples/Test_UB_matrix/Test_UB_matrix.instr new file mode 100644 index 0000000000..0a1b6ddf14 --- /dev/null +++ b/mcxtrace-comps/examples/Tests_samples/Test_UB_matrix/Test_UB_matrix.instr @@ -0,0 +1,79 @@ +/******************************************************************************* +* Instrument: Test_UB_matrix +* +* %I +* Written by: Peter Willendrup +* Date: October 2026 +* Origin: DTU +* %INSTRUMENT_SITE: Tests_samples +* +* Test of the UB_matrix component injecting a UB matrix into Single_crystal +* +* %D +* A Si single crystal with (1,1,0) along the beam and (0,0,1) up is oriented +* in three equivalent ways, and illuminated by a white beam: +* SAMPLE=0: reference, direct lattice vectors given to Single_crystal +* SAMPLE=1: UB_matrix with a UB matrix as 9 numbers (UBin) +* SAMPLE=2: UB_matrix with the cell and two reflections +* All three give the same Laue pattern and intensity. The goniometer angle A3 +* rotates all samples about the vertical axis. +* +* %Example: SAMPLE=0 Detector: Laue_I=1.8484e-11 +* %Example: SAMPLE=1 Detector: Laue_I=1.8484e-11 +* %Example: SAMPLE=2 Detector: Laue_I=1.8484e-11 +* +* %P +* SAMPLE: [1] Orientation input: 0 reference, 1 UBin, 2 cell + reflections +* A3: [deg] Goniometer rotation about the vertical axis +* +* %L +* A. T. Savici, Introduction to the use of the UB matrix (ORNL, 2011) +* +* %E +*******************************************************************************/ +DEFINE INSTRUMENT Test_UB_matrix(int SAMPLE=1, A3=0) + +TRACE + +COMPONENT Origin = Progress_bar() + AT (0,0,0) ABSOLUTE + +COMPONENT Source = Source_flat(radius=0.0001, dist=1, focus_xw=0.0002, focus_yh=0.0002, E0=15, dE=10) + AT (0,0,0) RELATIVE Origin + +COMPONENT gonio = Arm() + AT (0,0,1) RELATIVE Origin + ROTATED (0,A3,0) RELATIVE Origin + +COMPONENT UB_numbers = UB_matrix(target="Sample_UBin", filename="UB_numbers.txt", + UBin={-0.130209534095, 0.130209534095, 0, 0, 0, 0.184144089067, 0.130209534095, 0.130209534095, 0}) + AT (0,0,0) RELATIVE gonio + +COMPONENT UB_cell = UB_matrix(target="Sample_cell", filename="UB_cell.txt", + lattice_a=5.43053, lattice_b=5.43053, lattice_c=5.43053, + hkl1={1,1,0}, dir1={0,0,1}, hkl2={0,0,1}, dir2={0,1,0}) + AT (0,0,0) RELATIVE gonio + +COMPONENT Sample_ref = Single_crystal(xwidth=0.0005, yheight=0.0005, zdepth=0.0005, reflections="Si.lau", mosaic=5, + ax=-3.83996459, ay=0, az=3.83996459, bx=3.83996459, by=0, bz=3.83996459, cx=0, cy=5.43053, cz=0) + WHEN (SAMPLE==0) AT (0,0,0) RELATIVE gonio +EXTEND %{ + if (!SCATTERED) ABSORB; +%} + +COMPONENT Sample_UBin = Single_crystal(xwidth=0.0005, yheight=0.0005, zdepth=0.0005, reflections="Si.lau", mosaic=5) + WHEN (SAMPLE==1) AT (0,0,0) RELATIVE gonio +EXTEND %{ + if (!SCATTERED) ABSORB; +%} + +COMPONENT Sample_cell = Single_crystal(xwidth=0.0005, yheight=0.0005, zdepth=0.0005, reflections="Si.lau", mosaic=5) + WHEN (SAMPLE==2) AT (0,0,0) RELATIVE gonio +EXTEND %{ + if (!SCATTERED) ABSORB; +%} + +COMPONENT Laue = PSD_monitor_4PI(radius=0.5, nx=180, ny=90, filename="Laue") + AT (0,0,0) RELATIVE gonio + +END diff --git a/mcxtrace-comps/misc/UB_matrix.comp b/mcxtrace-comps/misc/UB_matrix.comp new file mode 100644 index 0000000000..e51b83de47 --- /dev/null +++ b/mcxtrace-comps/misc/UB_matrix.comp @@ -0,0 +1,395 @@ +/******************************************************************************* +* +* McXtrace, X-ray tracing package +* Copyright (C) 1997-2026, All rights reserved +* DTU Physics, Kgs. Lyngby, Denmark +* +* Component: UB_matrix +* +* %I +* Written by: Peter Willendrup +* Date: October 2026 +* Origin: DTU +* +* Import a UB matrix, inject it into a crystal component, and export it +* +* %D +* Reads a UB matrix and writes the corresponding crystal orientation into the +* SETTING parameters of the crystal component named by 'target', before that +* component initialises. The UB matrix is also written to a text file and is +* stored in the NeXus output as the run-time value of the 'UB' parameter. +* +* Convention: Q = 2*pi*R*UB*(h,k,l), with UB in the frame of the target +* component (z along the beam, y up at zero goniometer angles), no 2*pi in UB +* (|UB*(h,k,l)| = 1/d), B in the Busing-Levy frame (a* along x), and +* R the active rotation of the 'goniometer' component (default: the target). +* +* The UB matrix comes from exactly one of: +* - file: 3 rows of 3 numbers. format="mccode" reads UB in the convention above; +* format="isaw" reads a Mantid SaveIsawUB file (UB transposed, IPNS frame +* with x along the beam and z up). +* - UBin: 9 numbers, row-major, e.g. from Mantid getUB(). +* - lattice_a..lattice_cc with hkl1/dir1 and hkl2/dir2: two reflections whose +* directions in the target frame are given (dir1 is kept exactly, dir2 only +* fixes the plane). +* When a UB is read and lattice_a is also given, the cell of the UB is checked +* against the given cell. +* +* Supported targets (the type is looked up from the instance name): +* Single_crystal, FluoCrystal: the direct lattice vectors ax..cz. +* The component must be placed BEFORE its target in the instrument, and the +* target must not already define an orientation. +* +* Example: UB_matrix(target="Sample", file="sample.mat", format="isaw") +* +* %P +* INPUT PARAMETERS: +* target: [str] Name of the crystal component to inject into +* file: [str] UB matrix file to read, 0 for none +* format: [str] File format, "mccode" or "isaw" +* UBin: [AA^-1] UB matrix, 9 numbers row-major, used when no file is given +* lattice_a: [AA] Lattice parameter a; with hkl1..dir2 defines the UB, or checks a read UB +* lattice_b: [AA] Lattice parameter b +* lattice_c: [AA] Lattice parameter c +* lattice_aa: [deg] Lattice angle alpha +* lattice_bb: [deg] Lattice angle beta +* lattice_cc: [deg] Lattice angle gamma +* hkl1: [1] Miller indices of the first reflection +* dir1: [1] Direction of hkl1 in the target frame +* hkl2: [1] Miller indices of the second reflection +* dir2: [1] Direction in the target frame fixing the plane of hkl2 +* tol: [rad] Tolerance on the hkl1-hkl2 opening angle +* goniometer: [str] Component whose rotation is R, 0 for the target +* filename: [str] Output file name, 0 for none +* +* OUTPUT PARAMETERS: +* UB: [AA^-1] The UB matrix, 9 numbers row-major (stored in NeXus) +* +* %E +*******************************************************************************/ + +DEFINE COMPONENT UB_matrix + +SETTING PARAMETERS (string target="Sample", string file=0, string format="mccode", + vector UBin={0,0,0,0,0,0,0,0,0}, + lattice_a=0, lattice_b=0, lattice_c=0, lattice_aa=90, lattice_bb=90, lattice_cc=90, + vector hkl1={1,0,0}, vector dir1={0,0,1}, vector hkl2={0,1,0}, vector dir2={1,0,0}, + tol=1e-4, string goniometer=0, string filename="UB.txt", + vector UB={0,0,0,0,0,0,0,0,0}) + +SHARE +%{ +%include "read_table-lib" +#ifndef UB_MATRIX_DECL +#define UB_MATRIX_DECL + /* 3x3 matrices are row-major double[9] */ + void ub_mul(double *A, double *B, double *C) { /* C = A*B */ + double T[9]; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + T[3*i+j] = A[3*i]*B[j] + A[3*i+1]*B[3+j] + A[3*i+2]*B[6+j]; + for (int i = 0; i < 9; i++) C[i] = T[i]; + } + void ub_transpose(double *A, double *T) { + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) T[3*j+i] = A[3*i+j]; + } + double ub_det(double *A) { + return A[0]*(A[4]*A[8]-A[5]*A[7]) - A[1]*(A[3]*A[8]-A[5]*A[6]) + A[2]*(A[3]*A[7]-A[4]*A[6]); + } + int ub_inv(double *A, double *I) { /* returns 0 when singular */ + double d = ub_det(A); + if (fabs(d) < 1e-300) return 0; + I[0] = (A[4]*A[8]-A[5]*A[7])/d; I[1] = (A[2]*A[7]-A[1]*A[8])/d; I[2] = (A[1]*A[5]-A[2]*A[4])/d; + I[3] = (A[5]*A[6]-A[3]*A[8])/d; I[4] = (A[0]*A[8]-A[2]*A[6])/d; I[5] = (A[2]*A[3]-A[0]*A[5])/d; + I[6] = (A[3]*A[7]-A[4]*A[6])/d; I[7] = (A[1]*A[6]-A[0]*A[7])/d; I[8] = (A[0]*A[4]-A[1]*A[3])/d; + return 1; + } + /* B from the cell (Savici eq. 52, Busing-Levy), no 2pi; returns 0 for an impossible cell */ + int ub_B(double a, double b, double c, double al, double be, double ga, double *B) { + double G[9], Gs[9]; + al *= DEG2RAD; be *= DEG2RAD; ga *= DEG2RAD; + G[0] = a*a; G[1] = a*b*cos(ga); G[2] = a*c*cos(be); + G[3] = G[1]; G[4] = b*b; G[5] = b*c*cos(al); + G[6] = G[2]; G[7] = G[5]; G[8] = c*c; + if (a <= 0 || b <= 0 || c <= 0 || ub_det(G) <= 0 || !ub_inv(G, Gs)) return 0; + double as = sqrt(Gs[0]), bs = sqrt(Gs[4]), cs = sqrt(Gs[8]); + double cbs = Gs[2]/(as*cs), cgs = Gs[1]/(as*bs); + B[0] = as; B[1] = bs*cgs; B[2] = cs*cbs; + B[3] = 0; B[4] = bs*sqrt(1 - cgs*cgs); B[5] = -cs*sqrt(1 - cbs*cbs)*cos(al); + B[6] = 0; B[7] = 0; B[8] = 1/c; + return 1; + } + /* orthonormal triad as columns of T: t1 along v1, t3 along v1 x v2; returns 0 if degenerate */ + int ub_triad(double *v1, double *v2, double *T) { + double t1[3], t2[3], t3[3]; + double n1 = sqrt(v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2]); + vec_prod(t3[0], t3[1], t3[2], v1[0], v1[1], v1[2], v2[0], v2[1], v2[2]); + double n3 = sqrt(t3[0]*t3[0] + t3[1]*t3[1] + t3[2]*t3[2]); + if (n1 <= 0 || n3 <= 1e-12*n1*sqrt(v2[0]*v2[0] + v2[1]*v2[1] + v2[2]*v2[2])) return 0; + for (int i = 0; i < 3; i++) { t1[i] = v1[i]/n1; t3[i] /= n3; } + vec_prod(t2[0], t2[1], t2[2], t3[0], t3[1], t3[2], t1[0], t1[1], t1[2]); + for (int i = 0; i < 3; i++) { T[3*i] = t1[i]; T[3*i+1] = t2[i]; T[3*i+2] = t3[i]; } + return 1; + } + double ub_angle(double *u, double *v) { + double d = u[0]*v[0] + u[1]*v[1] + u[2]*v[2]; + double n = sqrt((u[0]*u[0] + u[1]*u[1] + u[2]*u[2])*(v[0]*v[0] + v[1]*v[1] + v[2]*v[2])); + return acos(fmax(-1, fmin(1, d/n))); + } + /* U from two reflections (Savici eqs. 76-81); returns 0 on failure */ + int ub_U2(double *B, double *h1, double *d1, double *h2, double *d2, double tol, double *U) { + double v1[3], v2[3], Tc[9], Tn[9], TcT[9]; + for (int i = 0; i < 3; i++) { + v1[i] = B[3*i]*h1[0] + B[3*i+1]*h1[1] + B[3*i+2]*h1[2]; + v2[i] = B[3*i]*h2[0] + B[3*i+1]*h2[1] + B[3*i+2]*h2[2]; + } + if (!ub_triad(v1, v2, Tc) || !ub_triad(d1, d2, Tn)) return 0; + if (fabs(ub_angle(v1, v2) - ub_angle(d1, d2)) > tol) return 0; + ub_transpose(Tc, TcT); + ub_mul(Tn, TcT, U); + return 1; + } + /* B by Cholesky of (UB)^T(UB) (Savici eq. 88), then U = UB*B^-1; returns 0 if not positive definite */ + int ub_split(double *UB, double *U, double *B) { + double UBT[9], M[9], Bi[9]; + ub_transpose(UB, UBT); + ub_mul(UBT, UB, M); + if (M[0] <= 0) return 0; + B[0] = sqrt(M[0]); B[1] = M[1]/B[0]; B[2] = M[2]/B[0]; + double r = M[4] - B[1]*B[1]; + if (r <= 0) return 0; + B[4] = sqrt(r); B[5] = (M[5] - B[1]*B[2])/B[4]; + r = M[8] - B[2]*B[2] - B[5]*B[5]; + if (r <= 0) return 0; + B[8] = sqrt(r); B[3] = B[6] = B[7] = 0; + if (!ub_inv(B, Bi)) return 0; + ub_mul(UB, Bi, U); + return 1; + } + void ub_fprint(FILE *f, char *title, double *M) { + fprintf(f, "# %s\n", title); + for (int i = 0; i < 3; i++) fprintf(f, "%.12g %.12g %.12g\n", M[3*i], M[3*i+1], M[3*i+2]); + } +#endif +%} + +DECLARE +%{ + double U_out[9]; + double B_out[9]; + double R_out[9]; + double cell[6]; +%} + +INITIALIZE +%{ +#ifndef MCCODE_HAS_GETCOMP_TYPE +#error "UB_matrix needs a code generator that provides _getcomp_type() and MCCODE_COMP_TYPE_ (the C code generator in this McCode release does; mccode-antlr may not yet)." +#endif + double M[9], Bc[9], Uc[9], Ainv[9], A[9]; + int have_file = strlen(file) && strcmp(file, "0") && strcmp(file, "NULL"); + int have_ubin = 0; + for (int i = 0; i < 9; i++) if (UBin[i]) have_ubin = 1; + + if (have_file + have_ubin + (!have_file && !have_ubin && lattice_a > 0) != 1) { + fprintf(stderr, "UB_matrix: %s: ERROR: give exactly one of file, UBin, or lattice_a..lattice_cc with hkl1/dir1/hkl2/dir2.\n", NAME_CURRENT_COMP); + exit(-1); + } + + if (have_file) { + /* the first 3 lines that are not blank or '#' comments; ISAW files carry more lines after them */ + FILE *fp = Open_File(file, "r", NULL); + char line[1024]; + int rows = 0; + while (fp && rows < 3 && fgets(line, sizeof(line), fp)) { + char *s = line + strspn(line, " \t"); + if (*s == '#' || *s == '\n' || *s == '\r' || *s == '\0') continue; + if (sscanf(s, "%lf %lf %lf", &M[3*rows], &M[3*rows+1], &M[3*rows+2]) != 3) break; + rows++; + } + if (fp) fclose(fp); + if (rows < 3) { + fprintf(stderr, "UB_matrix: %s: ERROR: could not read 3 rows of 3 numbers from %s.\n", NAME_CURRENT_COMP, file); + exit(-1); + } + if (!strcmp(format, "isaw")) { + /* file holds UB^T in the IPNS frame: (x,y,z)_McStas = (y,z,x)_IPNS */ + double P[9] = {0,1,0, 0,0,1, 1,0,0}, MT[9]; + ub_transpose(M, MT); + ub_mul(P, MT, UB); + } else if (!strcmp(format, "mccode")) { + for (int i = 0; i < 9; i++) UB[i] = M[i]; + } else { + fprintf(stderr, "UB_matrix: %s: ERROR: unknown format \"%s\" (use \"mccode\" or \"isaw\").\n", NAME_CURRENT_COMP, format); + exit(-1); + } + } else if (have_ubin) { + for (int i = 0; i < 9; i++) UB[i] = UBin[i]; + } else { + if (!ub_B(lattice_a, lattice_b, lattice_c, lattice_aa, lattice_bb, lattice_cc, Bc)) { + fprintf(stderr, "UB_matrix: %s: ERROR: impossible lattice parameters.\n", NAME_CURRENT_COMP); + exit(-1); + } + if (!ub_U2(Bc, hkl1, dir1, hkl2, dir2, tol, Uc)) { + fprintf(stderr, "UB_matrix: %s: ERROR: hkl1/hkl2 or dir1/dir2 are parallel, or their opening angles differ by more than tol=%g rad.\n", + NAME_CURRENT_COMP, tol); + exit(-1); + } + ub_mul(Uc, Bc, UB); + } + + /* every UB passes the same checks */ + if (!ub_split(UB, U_out, B_out)) { + fprintf(stderr, "UB_matrix: %s: ERROR: the UB matrix does not describe a crystal ((UB)^T(UB) not positive definite).\n", NAME_CURRENT_COMP); + exit(-1); + } + if (ub_det(UB) < 0) { + fprintf(stderr, "UB_matrix: %s: ERROR: the UB matrix has a negative determinant (improper U). " + "It was probably written with the opposite Q sign convention; negate it.\n", NAME_CURRENT_COMP); + exit(-1); + } + { /* cell from B^T B = G*, G = (G*)^-1 */ + double BT[9], Gs[9], G[9]; + ub_transpose(B_out, BT); ub_mul(BT, B_out, Gs); ub_inv(Gs, G); + cell[0] = sqrt(G[0]); cell[1] = sqrt(G[4]); cell[2] = sqrt(G[8]); + cell[3] = acos(G[5]/(cell[1]*cell[2]))*RAD2DEG; + cell[4] = acos(G[2]/(cell[0]*cell[2]))*RAD2DEG; + cell[5] = acos(G[1]/(cell[0]*cell[1]))*RAD2DEG; + } + if ((have_file || have_ubin) && lattice_a > 0) { + double given[6] = {lattice_a, lattice_b, lattice_c, lattice_aa, lattice_bb, lattice_cc}; + for (int i = 0; i < 6; i++) + if (fabs(cell[i] - given[i]) > 1e-4*fabs(given[i])) { + fprintf(stderr, "UB_matrix: %s: ERROR: the UB cell (%g %g %g %g %g %g) does not match the given lattice parameters.\n", + NAME_CURRENT_COMP, cell[0], cell[1], cell[2], cell[3], cell[4], cell[5]); + exit(-1); + } + } + + /* goniometer rotation: McStas rotations are passive, R is their transpose */ + char *gonio = (strlen(goniometer) && strcmp(goniometer, "0") && strcmp(goniometer, "NULL")) ? goniometer : target; + Rotation Rp; + if (!_getcomp_rotation(gonio, Rp)) { + fprintf(stderr, "UB_matrix: %s: ERROR: no component named \"%s\".\n", NAME_CURRENT_COMP, gonio); + exit(-1); + } + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) R_out[3*i+j] = Rp[j][i]; + + /* inject into the target */ + char *type = _getcomp_type(target); + if (!type) { + fprintf(stderr, "UB_matrix: %s: ERROR: no component named \"%s\".\n", NAME_CURRENT_COMP, target); + exit(-1); + } + if (_getcomp_index(target) <= _getcomp_index(NAME_CURRENT_COMP)) { + fprintf(stderr, "UB_matrix: %s: ERROR: must be placed before its target \"%s\", which would otherwise already be initialised.\n", + NAME_CURRENT_COMP, target); + exit(-1); + } + /* direct lattice vectors a, b, c are the rows of UB^-1 (A^T UB = I) */ + ub_inv(UB, Ainv); + ub_transpose(Ainv, A); /* columns a, b, c */ + int injected = 0; + + /* direct lattice vectors, as COMP_GETPAR3(type, target, ax) etc. */ + #define UB_INJECT_ABC(T) \ + if (!strcmp(type, #T)) { \ + _class_##T##_parameters *p = (_class_##T##_parameters *) _getvar_parameters(target); \ + if (p->ax || p->ay || p->az || p->bx || p->by || p->bz || p->cx || p->cy || p->cz) { \ + fprintf(stderr, "UB_matrix: %s: ERROR: %s \"%s\" already has lattice vectors ax..cz; remove them.\n", NAME_CURRENT_COMP, type, target); \ + exit(-1); \ + } \ + p->ax = A[0]; p->ay = A[3]; p->az = A[6]; \ + p->bx = A[1]; p->by = A[4]; p->bz = A[7]; \ + p->cx = A[2]; p->cy = A[5]; p->cz = A[8]; \ + p->recip_cell = 0; \ + injected = 1; \ + } + /* NCrystal orientation, appended to cfg: lab directions of (100) and (010) */ + #define UB_INJECT_CFG(T) \ + if (!strcmp(type, #T)) { \ + char *cfgp = ((_class_##T##_parameters *) _getvar_parameters(target))->cfg; \ + if (strstr(cfgp, "dir1") || strstr(cfgp, "dir2")) { \ + fprintf(stderr, "UB_matrix: %s: ERROR: %s \"%s\" cfg already has dir1/dir2; remove them.\n", NAME_CURRENT_COMP, type, target); \ + exit(-1); \ + } \ + char dirs[512]; \ + snprintf(dirs, sizeof(dirs), ";dir1=@crys_hkl:1,0,0@lab:%.17g,%.17g,%.17g;dir2=@crys_hkl:0,1,0@lab:%.17g,%.17g,%.17g", \ + UB[0], UB[3], UB[6], UB[1], UB[4], UB[7]); \ + if (strlen(cfgp) + strlen(dirs) >= 16384) { \ + fprintf(stderr, "UB_matrix: %s: ERROR: cfg of \"%s\" too long to append dir1/dir2.\n", NAME_CURRENT_COMP, target); \ + exit(-1); \ + } \ + strcat(cfgp, dirs); \ + injected = 1; \ + } + + #ifdef MCCODE_COMP_TYPE_Single_crystal + UB_INJECT_ABC(Single_crystal) + #endif + #ifdef MCCODE_COMP_TYPE_Single_crystal_process + UB_INJECT_ABC(Single_crystal_process) + #endif + #ifdef MCCODE_COMP_TYPE_Single_magnetic_crystal + UB_INJECT_ABC(Single_magnetic_crystal) + #endif + #ifdef MCCODE_COMP_TYPE_Single_crystal_inelastic + UB_INJECT_ABC(Single_crystal_inelastic) + #endif + #ifdef MCCODE_COMP_TYPE_FluoCrystal + UB_INJECT_ABC(FluoCrystal) + #endif + #ifdef MCCODE_COMP_TYPE_NCrystal_sample + UB_INJECT_CFG(NCrystal_sample) + #endif + #ifdef MCCODE_COMP_TYPE_NCrystal_process + UB_INJECT_CFG(NCrystal_process) + #endif + #ifdef MCCODE_COMP_TYPE_Union_NCrystal_material + UB_INJECT_CFG(Union_NCrystal_material) + #endif + #undef UB_INJECT_ABC + #undef UB_INJECT_CFG + + if (!injected) { + fprintf(stderr, "UB_matrix: %s: ERROR: target \"%s\" is a %s, which UB_matrix does not support.\n", NAME_CURRENT_COMP, target, type); + exit(-1); + } + + MPI_MASTER( + printf("UB_matrix: %s: injected into %s \"%s\"; cell %g %g %g AA, %g %g %g deg\n", NAME_CURRENT_COMP, type, target, + cell[0], cell[1], cell[2], cell[3], cell[4], cell[5]); + printf("UB_matrix: %s: UB = [%g %g %g; %g %g %g; %g %g %g] 1/AA (no 2pi)\n", NAME_CURRENT_COMP, + UB[0], UB[1], UB[2], UB[3], UB[4], UB[5], UB[6], UB[7], UB[8]); + ); +%} + +SAVE +%{ + if (strlen(filename) && strcmp(filename, "0") && strcmp(filename, "NULL") && !mcdisable_output_files) { + MPI_MASTER( + char *path = mcfull_file(filename, NULL); + FILE *f = fopen(path, "w"); + if (!f) fprintf(stderr, "UB_matrix: %s: WARNING: could not write %s\n", NAME_CURRENT_COMP, path); + else { + double RUB[9]; + ub_mul(R_out, UB, RUB); + fprintf(f, "# UB matrix written by McCode UB_matrix component %s\n", NAME_CURRENT_COMP); + fprintf(f, "# target: %s\n", target); + fprintf(f, "# convention: Q = 2*pi*R*UB*(h,k,l); frame: z along beam, y up; no 2*pi in UB; Q = kf - ki; B Busing-Levy (a* along x)\n"); + fprintf(f, "# cell: %.10g %.10g %.10g %.10g %.10g %.10g\n", cell[0], cell[1], cell[2], cell[3], cell[4], cell[5]); + ub_fprint(f, "UB [1/AA] (read back with format=\"mccode\")", UB); + ub_fprint(f, "U", U_out); + ub_fprint(f, "B [1/AA]", B_out); + ub_fprint(f, "R (goniometer, active)", R_out); + ub_fprint(f, "R*UB [1/AA] (lab frame)", RUB); + fclose(f); + } + free(path); + ); + } +%} + +END From 742952f89436ddc99bc3f8164d78337a705820c3 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Sun, 25 Oct 2026 13:29:43 +0100 Subject: [PATCH 3/3] UB_matrix: support contrib Dispersion_relation as a target Dispersion_relation takes direct real-space lattice vectors ax..cz in the component frame, like Single_crystal but without recip_cell, so UB_matrix injects them the same way. Test_UB_matrix gains SAMPLE=6 (Dispersion_relation with ax..cz given) and SAMPLE=7 (injected through UBin) for an orthorhombic 4x5x6 AA cell. The anisotropic dispersion mode is written at compile time by write_ub_dispersion.py (standard library only). Both give 2.26331e-08 at 1e6 rays; the 2D patterns agree to 1.5e-8 relative. Co-Authored-By: Claude Opus 5.5 --- .../Test_UB_matrix/Test_UB_matrix.instr | 33 +++++++++++++++++-- .../Test_UB_matrix/write_ub_dispersion.py | 22 +++++++++++++ mcstas-comps/misc/UB_matrix.comp | 19 ++++++----- mcxtrace-comps/misc/UB_matrix.comp | 17 ++++++---- 4 files changed, 74 insertions(+), 17 deletions(-) create mode 100644 mcstas-comps/examples/Tests_samples/Test_UB_matrix/write_ub_dispersion.py diff --git a/mcstas-comps/examples/Tests_samples/Test_UB_matrix/Test_UB_matrix.instr b/mcstas-comps/examples/Tests_samples/Test_UB_matrix/Test_UB_matrix.instr index 69657009ff..a21fd8b4d9 100644 --- a/mcstas-comps/examples/Tests_samples/Test_UB_matrix/Test_UB_matrix.instr +++ b/mcstas-comps/examples/Tests_samples/Test_UB_matrix/Test_UB_matrix.instr @@ -18,7 +18,12 @@ * SAMPLE=3: Single_crystal, UB_matrix reading a Mantid SaveIsawUB file (UB_Al.mat) * SAMPLE=4: NCrystal_sample reference, dir1/dir2 given in cfg * SAMPLE=5: NCrystal_sample, UB_matrix with UBin appending dir1/dir2 to cfg -* SAMPLE=0-3 give the same Laue pattern and intensity, as do SAMPLE=4-5. +* SAMPLE=6: Dispersion_relation reference, orthorhombic cell (4, 5, 6 AA) with +* direct lattice vectors given, same orientation as above +* SAMPLE=7: Dispersion_relation, UB_matrix with UBin +* SAMPLE=0-3 give the same Laue pattern and intensity, as do SAMPLE=4-5 and +* SAMPLE=6-7. The Dispersion_relation mode is written by the SHELL command +* below (write_ub_dispersion.py, Python standard library only). * The goniometer angle A3 rotates all samples about the vertical axis. * * %Example: SAMPLE=0 Detector: Laue_I=1.2705e-09 @@ -27,9 +32,11 @@ * %Example: SAMPLE=3 Detector: Laue_I=1.2705e-09 * %Example: SAMPLE=4 Detector: Laue_I=5.6971e-09 * %Example: SAMPLE=5 Detector: Laue_I=5.6971e-09 +* %Example: SAMPLE=6 Detector: Laue_I=2.2633e-08 +* %Example: SAMPLE=7 Detector: Laue_I=2.2633e-08 * * %P -* SAMPLE: [1] Orientation input: 0-3 Single_crystal (reference, UBin, cell, ISAW file), 4-5 NCrystal_sample (reference, UBin) +* SAMPLE: [1] Orientation input: 0-3 Single_crystal (reference, UBin, cell, ISAW file), 4-5 NCrystal_sample (reference, UBin), 6-7 Dispersion_relation (reference, UBin) * A3: [deg] Goniometer rotation about the vertical axis * * %L @@ -39,6 +46,9 @@ *******************************************************************************/ DEFINE INSTRUMENT Test_UB_matrix(int SAMPLE=1, A3=0) +// Write the dispersion file of the Dispersion_relation samples, unless it exists. +SHELL "python write_ub_dispersion.py" + TRACE COMPONENT Origin = Progress_bar() @@ -68,6 +78,10 @@ COMPONENT UB_nc = UB_matrix(target="Sample_nc", filename="UB_nc.txt", UBin={-0.17459426696, 0.17459426696, 0, 0, 0, 0.246913580247, 0.17459426696, 0.17459426696, 0}) AT (0,0,0) RELATIVE gonio +COMPONENT UB_disp = UB_matrix(target="Sample_disp", filename="UB_disp.txt", + UBin={-0.176776695296637, 0.14142135623731, 0, 0, 0, 0.166666666666667, 0.176776695296637, 0.14142135623731, 0}) + AT (0,0,0) RELATIVE gonio + COMPONENT Sample_ref = Single_crystal(xwidth=0.005, yheight=0.005, zdepth=0.005, reflections="Al.lau", mosaic=5, ax=-2.863782463806, ay=0, az=2.863782463806, bx=2.863782463806, by=0, bz=2.863782463806, cx=0, cy=4.05, cz=0) WHEN (SAMPLE==0) AT (0,0,0) RELATIVE gonio @@ -107,6 +121,21 @@ EXTEND %{ if (!SCATTERED) ABSORB; %} +COMPONENT Sample_disp_ref = Dispersion_relation(radius=0.0025, yheight=0.005, sigma_abs=0, sigma_inc=0, DW=1, T=10, + folder_path="UB_dispersion", + ax=-2.82842712474619, ay=0, az=2.82842712474619, bx=3.53553390593274, by=0, bz=3.53553390593274, cx=0, cy=6, cz=0) + WHEN (SAMPLE==6) AT (0,0,0) RELATIVE gonio +EXTEND %{ + if (!SCATTERED) ABSORB; +%} + +COMPONENT Sample_disp = Dispersion_relation(radius=0.0025, yheight=0.005, sigma_abs=0, sigma_inc=0, DW=1, T=10, + folder_path="UB_dispersion") + WHEN (SAMPLE==7) AT (0,0,0) RELATIVE gonio +EXTEND %{ + if (!SCATTERED) ABSORB; +%} + COMPONENT Laue = PSD_monitor_4PI(radius=0.5, nx=180, ny=90, filename="Laue") AT (0,0,0) RELATIVE gonio diff --git a/mcstas-comps/examples/Tests_samples/Test_UB_matrix/write_ub_dispersion.py b/mcstas-comps/examples/Tests_samples/Test_UB_matrix/write_ub_dispersion.py new file mode 100644 index 0000000000..61c72fa38b --- /dev/null +++ b/mcstas-comps/examples/Tests_samples/Test_UB_matrix/write_ub_dispersion.py @@ -0,0 +1,22 @@ +"""Write the dispersion file read by the Dispersion_relation samples of Test_UB_matrix.instr. + +Run by the SHELL command of the instrument when it is compiled; the file is only +written if it does not exist yet. The mode is deliberately anisotropic in h, k +and l, so a wrongly oriented lattice changes the scattered intensity. +""" +import math +import os + +OUTPUT = os.path.join("UB_dispersion", "mode_0.dat") +POINTS = 9 + +if not os.path.exists(OUTPUT): + os.makedirs(os.path.dirname(OUTPUT), exist_ok=True) + grid = [i / (POINTS - 1) for i in range(POINTS)] + with open(OUTPUT, "w") as f: + for h in grid: + for k in grid: + for l in grid: + energy = 3 + 4 * math.sin(math.pi * h) ** 2 + 2 * math.sin(math.pi * k) ** 2 + math.sin(math.pi * l) ** 2 + intensity = 1 + h + f.write(f"{h:.6g} {k:.6g} {l:.6g} {energy:.6g} {intensity:.6g}\n") diff --git a/mcstas-comps/misc/UB_matrix.comp b/mcstas-comps/misc/UB_matrix.comp index 3d6a713595..9c58bd5f63 100644 --- a/mcstas-comps/misc/UB_matrix.comp +++ b/mcstas-comps/misc/UB_matrix.comp @@ -37,7 +37,7 @@ * * Supported targets (the type is looked up from the instance name): * Single_crystal, Single_crystal_process, Single_magnetic_crystal, -* Single_crystal_inelastic: the direct lattice vectors ax..cz. +* Single_crystal_inelastic, Dispersion_relation: the direct lattice vectors ax..cz. * NCrystal_sample, NCrystal_process, Union_NCrystal_material: dir1/dir2 * appended to cfg. * The component must be placed BEFORE its target in the instrument, and the @@ -299,7 +299,7 @@ INITIALIZE int injected = 0; /* direct lattice vectors, as COMP_GETPAR3(type, target, ax) etc. */ - #define UB_INJECT_ABC(T) \ + #define UB_INJECT_ABC(T, EXTRA) \ if (!strcmp(type, #T)) { \ _class_##T##_parameters *p = (_class_##T##_parameters *) _getvar_parameters(target); \ if (p->ax || p->ay || p->az || p->bx || p->by || p->bz || p->cx || p->cy || p->cz) { \ @@ -309,7 +309,7 @@ INITIALIZE p->ax = A[0]; p->ay = A[3]; p->az = A[6]; \ p->bx = A[1]; p->by = A[4]; p->bz = A[7]; \ p->cx = A[2]; p->cy = A[5]; p->cz = A[8]; \ - p->recip_cell = 0; \ + EXTRA \ injected = 1; \ } /* NCrystal orientation, appended to cfg: lab directions of (100) and (010) */ @@ -332,19 +332,22 @@ INITIALIZE } #ifdef MCCODE_COMP_TYPE_Single_crystal - UB_INJECT_ABC(Single_crystal) + UB_INJECT_ABC(Single_crystal, p->recip_cell = 0;) #endif #ifdef MCCODE_COMP_TYPE_Single_crystal_process - UB_INJECT_ABC(Single_crystal_process) + UB_INJECT_ABC(Single_crystal_process, p->recip_cell = 0;) #endif #ifdef MCCODE_COMP_TYPE_Single_magnetic_crystal - UB_INJECT_ABC(Single_magnetic_crystal) + UB_INJECT_ABC(Single_magnetic_crystal, p->recip_cell = 0;) #endif #ifdef MCCODE_COMP_TYPE_Single_crystal_inelastic - UB_INJECT_ABC(Single_crystal_inelastic) + UB_INJECT_ABC(Single_crystal_inelastic, p->recip_cell = 0;) #endif #ifdef MCCODE_COMP_TYPE_FluoCrystal - UB_INJECT_ABC(FluoCrystal) + UB_INJECT_ABC(FluoCrystal, p->recip_cell = 0;) + #endif + #ifdef MCCODE_COMP_TYPE_Dispersion_relation + UB_INJECT_ABC(Dispersion_relation, ) #endif #ifdef MCCODE_COMP_TYPE_NCrystal_sample UB_INJECT_CFG(NCrystal_sample) diff --git a/mcxtrace-comps/misc/UB_matrix.comp b/mcxtrace-comps/misc/UB_matrix.comp index e51b83de47..1bc1ae6a38 100644 --- a/mcxtrace-comps/misc/UB_matrix.comp +++ b/mcxtrace-comps/misc/UB_matrix.comp @@ -294,7 +294,7 @@ INITIALIZE int injected = 0; /* direct lattice vectors, as COMP_GETPAR3(type, target, ax) etc. */ - #define UB_INJECT_ABC(T) \ + #define UB_INJECT_ABC(T, EXTRA) \ if (!strcmp(type, #T)) { \ _class_##T##_parameters *p = (_class_##T##_parameters *) _getvar_parameters(target); \ if (p->ax || p->ay || p->az || p->bx || p->by || p->bz || p->cx || p->cy || p->cz) { \ @@ -304,7 +304,7 @@ INITIALIZE p->ax = A[0]; p->ay = A[3]; p->az = A[6]; \ p->bx = A[1]; p->by = A[4]; p->bz = A[7]; \ p->cx = A[2]; p->cy = A[5]; p->cz = A[8]; \ - p->recip_cell = 0; \ + EXTRA \ injected = 1; \ } /* NCrystal orientation, appended to cfg: lab directions of (100) and (010) */ @@ -327,19 +327,22 @@ INITIALIZE } #ifdef MCCODE_COMP_TYPE_Single_crystal - UB_INJECT_ABC(Single_crystal) + UB_INJECT_ABC(Single_crystal, p->recip_cell = 0;) #endif #ifdef MCCODE_COMP_TYPE_Single_crystal_process - UB_INJECT_ABC(Single_crystal_process) + UB_INJECT_ABC(Single_crystal_process, p->recip_cell = 0;) #endif #ifdef MCCODE_COMP_TYPE_Single_magnetic_crystal - UB_INJECT_ABC(Single_magnetic_crystal) + UB_INJECT_ABC(Single_magnetic_crystal, p->recip_cell = 0;) #endif #ifdef MCCODE_COMP_TYPE_Single_crystal_inelastic - UB_INJECT_ABC(Single_crystal_inelastic) + UB_INJECT_ABC(Single_crystal_inelastic, p->recip_cell = 0;) #endif #ifdef MCCODE_COMP_TYPE_FluoCrystal - UB_INJECT_ABC(FluoCrystal) + UB_INJECT_ABC(FluoCrystal, p->recip_cell = 0;) + #endif + #ifdef MCCODE_COMP_TYPE_Dispersion_relation + UB_INJECT_ABC(Dispersion_relation, ) #endif #ifdef MCCODE_COMP_TYPE_NCrystal_sample UB_INJECT_CFG(NCrystal_sample)