From b9a3c9c660f6d5b4f5aa478487bfd2636ec43384 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Tue, 13 Oct 2026 07:47:30 +0200 Subject: [PATCH 1/6] conic.h: fix radial normal in degenerate branch of getNormConic The det <= 0 branch (surface parallel to the z-axis, e.g. at the waist of an ellipsoid or for a cylinder) divided the y-component by r^2 instead of r, giving a non-unit, non-radial normal. Co-Authored-By: Claude Opus 5.5 --- mcstas-comps/share/conic.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mcstas-comps/share/conic.h b/mcstas-comps/share/conic.h index 8aa22c5fb..abe4841be 100755 --- a/mcstas-comps/share/conic.h +++ b/mcstas-comps/share/conic.h @@ -1321,8 +1321,8 @@ MAKE SURE p IS ON SURFACE Vec getNormConic(Point p, ConicSurf s) { double det = s.k2*s.k2+4*s.k3*(p.x*p.x+p.y*p.y-s.k1); if (det <= 0.){ - - return makeVec(-p.x/sqrt(p.x*p.x + p.y*p.y),-p.y/(p.x*p.x + p.y*p.y),0); + double r = sqrt(p.x*p.x + p.y*p.y); + return makeVec(-p.x/r,-p.y/r,0); } double den = sqrt(det); double nx = -2*p.x/den; From 2125398924e0a0304018dce97bd20515571bb67c Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Tue, 13 Oct 2026 07:47:30 +0200 Subject: [PATCH 2/6] conic.h: make the figure error of ConicSurf mirrors configurable reflectNeutronConic() applied a hard-coded Gaussian slope error of 1 urad (meridional) and 250 urad (sagittal). These are now ConicSurf fields (figerr_m, figerr_s), set by the make*() functions to the previous values (CONIC_SIGMA_M/S) and exposed as sigma_m/sigma_s on Conics_EH/HE/PH/PP with unchanged defaults, i.e. results are unchanged by default. With both set to 0 the rotation and its two randnorm() draws are skipped. Dead debug output is removed. Mirror_conic and the test instruments pass the parameters through; Mirror_conic defaults to no figure error. Co-Authored-By: Claude Opus 5.5 --- mcstas-comps/contrib/Conics_EH.comp | 8 ++- mcstas-comps/contrib/Conics_HE.comp | 8 ++- mcstas-comps/contrib/Conics_PH.comp | 8 ++- mcstas-comps/contrib/Conics_PP.comp | 8 ++- mcstas-comps/contrib/Mirror_conic.comp | 14 ++-- .../Test_Conics_pairs/Test_Conics_pairs.instr | 20 +++--- .../Test_Mirror_conic/Test_Mirror_conic.instr | 6 +- mcstas-comps/share/conic.h | 71 ++++++++----------- 8 files changed, 82 insertions(+), 61 deletions(-) diff --git a/mcstas-comps/contrib/Conics_EH.comp b/mcstas-comps/contrib/Conics_EH.comp index 4d5fedbc1..7fe5ab38f 100644 --- a/mcstas-comps/contrib/Conics_EH.comp +++ b/mcstas-comps/contrib/Conics_EH.comp @@ -41,6 +41,8 @@ * alpha: [AA] Slope of reflectivity for reflectivity curve approximation * W: [AA^-1] Width of supermirror cut-off * transmit: [1] Fraction of statistics to assign to transmitted beam - NOT YET IMPLEMENTED +* sigma_m: [rad] RMS meridional slope (figure) error of the mirror surfaces +* sigma_s: [rad] RMS sagittal slope (figure) error of the mirror surfaces * mirr_thick: [m] Thickness of mirror shell surfaces - NOT YET IMPLEMENTED * %End * @@ -51,7 +53,7 @@ DEFINE COMPONENT Conics_EH SETTING PARAMETERS ( rmin=0.0031416, rmax=0.05236,focal_length_u=10, focal_length_d=10, le=0.25, lh=0.25, int nshells=4, m=1, mirr_thick=0, int disk=1, vector radii=NULL, - R0 = 0.99, Qc = 0.021, W = 0.003, alpha = 6.07, transmit = 0) + R0 = 0.99, Qc = 0.021, W = 0.003, alpha = 6.07, transmit = 0, sigma_m=1e-6, sigma_s=250e-6) SHARE %{ @@ -106,6 +108,10 @@ INITIALIZE addHyperboloid (focal_length_d, focal_length_d + 2.0 * cH, pi, 0, lh, m, R0, Qc, alpha, W, &s); } } + for (i = 0; i < s.num_c; i++) { + s.c[i].figerr_m = sigma_m; + s.c[i].figerr_s = sigma_s; + } if (disk) { addDisk (pm->zs, 0.0, rConic (pm->ze, *pm), &s); } diff --git a/mcstas-comps/contrib/Conics_HE.comp b/mcstas-comps/contrib/Conics_HE.comp index ef685423d..13f5a63de 100644 --- a/mcstas-comps/contrib/Conics_HE.comp +++ b/mcstas-comps/contrib/Conics_HE.comp @@ -41,6 +41,8 @@ * alpha: [AA] Slope of reflectivity for reflectivity curve approximation * W: [AA^-1] Width of supermirror cut-off * transmit: [1] Fraction of statistics to assign to transmitted beam - NOT YET IMPLEMENTED +* sigma_m: [rad] RMS meridional slope (figure) error of the mirror surfaces +* sigma_s: [rad] RMS sagittal slope (figure) error of the mirror surfaces * mirr_thick: [m] Thickness of mirror shell surfaces - NOT YET IMPLEMENTED * %End * @@ -51,7 +53,7 @@ DEFINE COMPONENT Conics_HE SETTING PARAMETERS ( rmin=0.0031416, rmax=0.05236,focal_length_u=10, focal_length_d=10, le=0.25, lh=0.25, int nshells=4, m=1, mirr_thick=0, int disk=1, vector radii=NULL, - R0 = 0.99, Qc = 0.021, W = 0.003, alpha = 6.07, transmit = 0) + R0 = 0.99, Qc = 0.021, W = 0.003, alpha = 6.07, transmit = 0, sigma_m=1e-6, sigma_s=250e-6) SHARE %{ @@ -106,6 +108,10 @@ INITIALIZE addEllipsoid (-focal_length_u - 2.0 * cH, focal_length_d, pi, 0, le, m, R0, Qc, alpha, W, &s); } } + for (i = 0; i < s.num_c; i++) { + s.c[i].figerr_m = sigma_m; + s.c[i].figerr_s = sigma_s; + } if (disk) { addDisk (pm->zs, 0.0, rConic (pm->ze, *pm), &s); } diff --git a/mcstas-comps/contrib/Conics_PH.comp b/mcstas-comps/contrib/Conics_PH.comp index a54d07d95..e6ce00269 100644 --- a/mcstas-comps/contrib/Conics_PH.comp +++ b/mcstas-comps/contrib/Conics_PH.comp @@ -39,6 +39,8 @@ * alpha: [AA] Slope of reflectivity for reflectivity curve approximation * W: [AA^-1] Width of supermirror cut-off * transmit: [1] Fraction of statistics to assign to transmitted beam - NOT YET IMPLEMENTED +* sigma_m: [rad] RMS meridional slope (figure) error of the mirror surfaces +* sigma_s: [rad] RMS sagittal slope (figure) error of the mirror surfaces * mirr_thick: [m] Thickness of mirror shell surfaces - NOT YET IMPLEMENTED * * %L @@ -52,7 +54,7 @@ DEFINE COMPONENT Conics_PH SETTING PARAMETERS ( rmin=0.325,rmax=0.615, focal_length=10.070, lp=0.84, lh=0.84, int nshells=4, m=1, mirr_thick=0, int disk=1, vector radii=NULL, - R0=0.99, Qc=0.021, W=0.003, alpha=6.07, transmit=0) + R0=0.99, Qc=0.021, W=0.003, alpha=6.07, transmit=0, sigma_m=1e-6, sigma_s=250e-6) SHARE %{ @@ -99,6 +101,10 @@ INITIALIZE pm = addParaboloid (focal_length + 2.0 * cH, pi, -lp, 0, m, R0, Qc, alpha, W, &s); addHyperboloid (focal_length, focal_length + 2.0 * cH, pi, 0, lh, m, R0, Qc, alpha, W, &s); } + for (i = 0; i < s.num_c; i++) { + s.c[i].figerr_m = sigma_m; + s.c[i].figerr_s = sigma_s; + } if (disk) { addDisk (pm->zs, 0.0, rConic (pm->ze, *pm), &s); } diff --git a/mcstas-comps/contrib/Conics_PP.comp b/mcstas-comps/contrib/Conics_PP.comp index b56885cdb..8fe81cfb4 100644 --- a/mcstas-comps/contrib/Conics_PP.comp +++ b/mcstas-comps/contrib/Conics_PP.comp @@ -41,6 +41,8 @@ * alpha: [AA] Slope of reflectivity for reflectivity curve approximation * W: [AA^-1] Width of supermirror cut-off * transmit: [1] Fraction of statistics to assign to transmitted beam - NOT YET IMPLEMENTED +* sigma_m: [rad] RMS meridional slope (figure) error of the mirror surfaces +* sigma_s: [rad] RMS sagittal slope (figure) error of the mirror surfaces * mirr_thick: [m] Thickness of mirror shell surfaces - NOT YET IMPLEMENTED * %End * @@ -51,7 +53,7 @@ DEFINE COMPONENT Conics_PP SETTING PARAMETERS ( rmin=0.0031416, rmax=0.05236, focal_length_u=10.070, focal_length_d=10.070, lp1=0.84, lp2=0.84, int nshells=4, m=1, mirr_thick=0, int disk=1, vector radii=NULL, - R0 = 0.99, Qc = 0.021, W = 0.003, alpha = 6.07, transmit = 0) + R0 = 0.99, Qc = 0.021, W = 0.003, alpha = 6.07, transmit = 0, sigma_m=1e-6, sigma_s=250e-6) SHARE %{ @@ -95,6 +97,10 @@ INITIALIZE pm = addParaboloid (-focal_length_u, pi, -lp1, 0, m, R0, Qc, alpha, W, &s); addParaboloid (focal_length_d, pi, 0, lp2, m, R0, Qc, alpha, W, &s); } + for (i = 0; i < s.num_c; i++) { + s.c[i].figerr_m = sigma_m; + s.c[i].figerr_s = sigma_s; + } if (disk) { addDisk (pm->zs, 0.0, rConic (pm->ze, *pm), &s); } diff --git a/mcstas-comps/contrib/Mirror_conic.comp b/mcstas-comps/contrib/Mirror_conic.comp index ab8fac0c6..0c3ad09be 100644 --- a/mcstas-comps/contrib/Mirror_conic.comp +++ b/mcstas-comps/contrib/Mirror_conic.comp @@ -29,9 +29,9 @@ * unaffected. Compare with Mirror(transmit=1). * conic=1: the mirror is a conic.h ConicSurf cylinder of (very large) radius 1 km, * handled by reflectNeutronConic(), which always reflects and weights by R. -* Compare with Mirror(transmit=0). Note that reflectNeutronConic() absorbs -* neutrons beyond m times the critical angle of Ni (hard cut-off) and -* applies a small hard-coded figure error to the surface normal. +* Compare with Mirror(transmit=0). A Gaussian figure (slope) error can be +* applied via sigma_m and sigma_s; the default is 0 here, whereas the +* Conics_* components default to the conic.h values of 1 and 250 urad. * * Example: Mirror_conic(length=1, yheight=0.1, m=2) * @@ -44,13 +44,15 @@ * m: [1] m-value of material * W: [AA^-1] Width of supermirror cut-off * conic: [1] 0: use FlatSurf/reflectNeutronFlat(), 1: use ConicSurf/reflectNeutronConic() +* sigma_m: [rad] RMS meridional slope (figure) error, conic=1 only +* sigma_s: [rad] RMS sagittal slope (figure) error, conic=1 only * * %E *******************************************************************************/ DEFINE COMPONENT Mirror_conic -SETTING PARAMETERS (length=1, yheight=0.1, R0=0.99, Qc=0.021, alpha=6.07, m=2, W=0.003, int conic=0) +SETTING PARAMETERS (length=1, yheight=0.1, R0=0.99, Qc=0.021, alpha=6.07, m=2, W=0.003, int conic=0, sigma_m=0, sigma_s=0) SHARE %{ @@ -78,6 +80,8 @@ SHARE c->alpha = par[2]; c->m = par[3]; c->W = par[4]; + c->figerr_m = par[5]; + c->figerr_s = par[6]; c->type = PARA; } %} @@ -95,7 +99,7 @@ INITIALIZE coordinates; the local mirror plane x=0 is mapped to x=xshift. */ xshift = 1e3; - double par[5] = { R0, Qc, alpha, m, W }; + double par[7] = { R0, Qc, alpha, m, W, sigma_m, sigma_s }; mirror_conic_init (&fs, &cs, xshift, length, yheight, par); %} diff --git a/mcstas-comps/examples/Tests_optics/Test_Conics_pairs/Test_Conics_pairs.instr b/mcstas-comps/examples/Tests_optics/Test_Conics_pairs/Test_Conics_pairs.instr index 7231de614..739374510 100644 --- a/mcstas-comps/examples/Tests_optics/Test_Conics_pairs/Test_Conics_pairs.instr +++ b/mcstas-comps/examples/Tests_optics/Test_Conics_pairs/Test_Conics_pairs.instr @@ -38,13 +38,15 @@ * rmin: [m] Radius of the innermost shell * rmax: [m] Radius of the outermost shell * quadratic: [ ] Use the quadratic shell-radius mode instead of the radii vector +* sigma_m: [rad] RMS meridional slope (figure) error of the mirrors +* sigma_s: [rad] RMS sagittal slope (figure) error of the mirrors * * %L * * * %E *******************************************************************************/ -DEFINE INSTRUMENT Test_Conics_pairs(int OPTIC=1, ssize=1e-6, fs=10, fi=10, R0=0.99, m=3, W=0.003, alpha=6.07, int nshells=4, rmin=0.0031416, rmax=0.05236, int quadratic=0) +DEFINE INSTRUMENT Test_Conics_pairs(int OPTIC=1, ssize=1e-6, fs=10, fi=10, R0=0.99, m=3, W=0.003, alpha=6.07, int nshells=4, rmin=0.0031416, rmax=0.05236, int quadratic=0, sigma_m=1e-6, sigma_s=250e-6) DECLARE %{ @@ -102,7 +104,7 @@ COMPONENT psd_o0 = PSD_monitor(xwidth=0.12, yheight=0.12, restore_neutron=1, filename="psd_o0") AT(0,0,fs) RELATIVE PREVIOUS -COMPONENT cEH = Conics_EH(R0=R0,alpha=alpha,W=W,m=m, +COMPONENT cEH = Conics_EH(R0=R0,alpha=alpha,W=W,m=m, sigma_m=sigma_m, sigma_s=sigma_s, nshells=nshells, focal_length_u=fs, focal_length_d=fi, radii=radii, le=0.25, lh=0.25, disk=1 ) WHEN(OPTIC==1 && !quadratic) AT(0,0,fs) RELATIVE ss @@ -111,7 +113,7 @@ EXTEND number_of_reflections = SCATTERED - 1; %} -COMPONENT cEH_quadratic = Conics_EH(R0=R0,alpha=alpha,W=W,m=m, +COMPONENT cEH_quadratic = Conics_EH(R0=R0,alpha=alpha,W=W,m=m, sigma_m=sigma_m, sigma_s=sigma_s, nshells=nshells, focal_length_u=fs, focal_length_d=fi, rmin=rmin, rmax=rmax, le=0.25, lh=0.25, disk=1 ) WHEN(OPTIC==1 && quadratic) AT(0,0,fs) RELATIVE ss @@ -120,7 +122,7 @@ EXTEND number_of_reflections = SCATTERED - 1; %} -COMPONENT cPH = Conics_PH(R0=R0,alpha=alpha,W=W,m=m, +COMPONENT cPH = Conics_PH(R0=R0,alpha=alpha,W=W,m=m, sigma_m=sigma_m, sigma_s=sigma_s, nshells=nshells, focal_length=fi, radii=radii, lp=0.25, lh=0.25, disk=0 ) WHEN(OPTIC==2 && !quadratic) AT(0,0,fs) RELATIVE ss @@ -129,7 +131,7 @@ EXTEND number_of_reflections = SCATTERED - 1; %} -COMPONENT cPH_quadratic = Conics_PH(R0=R0,alpha=alpha,W=W,m=m, +COMPONENT cPH_quadratic = Conics_PH(R0=R0,alpha=alpha,W=W,m=m, sigma_m=sigma_m, sigma_s=sigma_s, nshells=nshells, focal_length=fi, rmin=rmin, rmax=rmax, lp=0.25, lh=0.25, disk=0 ) WHEN(OPTIC==2 && quadratic) AT(0,0,fs) RELATIVE ss @@ -138,7 +140,7 @@ EXTEND number_of_reflections = SCATTERED - 1; %} -COMPONENT cHE = Conics_HE(R0=R0,alpha=alpha,W=W,m=m, +COMPONENT cHE = Conics_HE(R0=R0,alpha=alpha,W=W,m=m, sigma_m=sigma_m, sigma_s=sigma_s, nshells=nshells, focal_length_u=fs, focal_length_d=fi, radii=radii, le=0.25, lh=0.25, disk=1 ) WHEN(OPTIC==3 && !quadratic) AT(0,0,fs) RELATIVE ss @@ -147,7 +149,7 @@ EXTEND number_of_reflections = SCATTERED - 1; %} -COMPONENT cHE_quadratic = Conics_HE(R0=R0,alpha=alpha,W=W,m=m, +COMPONENT cHE_quadratic = Conics_HE(R0=R0,alpha=alpha,W=W,m=m, sigma_m=sigma_m, sigma_s=sigma_s, nshells=nshells, focal_length_u=fs, focal_length_d=fi, rmin=rmin, rmax=rmax, le=0.25, lh=0.25, disk=1 ) WHEN(OPTIC==3 && quadratic) AT(0,0,fs) RELATIVE ss @@ -156,7 +158,7 @@ EXTEND number_of_reflections = SCATTERED - 1; %} -COMPONENT cPP = Conics_PP(R0=R0,alpha=alpha,W=W,m=m, +COMPONENT cPP = Conics_PP(R0=R0,alpha=alpha,W=W,m=m, sigma_m=sigma_m, sigma_s=sigma_s, nshells=nshells, focal_length_u=fs, focal_length_d=fi, radii=radii, lp1=0.25, lp2=0.25, disk=0 ) WHEN(OPTIC==4 && !quadratic) AT(0,0,fs) RELATIVE ss @@ -165,7 +167,7 @@ EXTEND number_of_reflections = SCATTERED - 1; %} -COMPONENT cPP_quadratic = Conics_PP(R0=R0,alpha=alpha,W=W,m=m, +COMPONENT cPP_quadratic = Conics_PP(R0=R0,alpha=alpha,W=W,m=m, sigma_m=sigma_m, sigma_s=sigma_s, nshells=nshells, focal_length_u=fs, focal_length_d=fi, rmin=rmin, rmax=rmax, lp1=0.25, lp2=0.25, disk=0 ) WHEN(OPTIC==4 && quadratic) AT(0,0,fs) RELATIVE ss diff --git a/mcstas-comps/examples/Tests_optics/Test_Mirror_conic/Test_Mirror_conic.instr b/mcstas-comps/examples/Tests_optics/Test_Mirror_conic/Test_Mirror_conic.instr index 2e930f2f2..3d9aa77c5 100644 --- a/mcstas-comps/examples/Tests_optics/Test_Mirror_conic/Test_Mirror_conic.instr +++ b/mcstas-comps/examples/Tests_optics/Test_Mirror_conic/Test_Mirror_conic.instr @@ -43,10 +43,12 @@ * W: [AA^-1] Width of supermirror cut-off * Lmin: [AA] Minimum source wavelength * Lmax: [AA] Maximum source wavelength +* sigma_m: [rad] RMS meridional slope (figure) error, COMP=3 only +* sigma_s: [rad] RMS sagittal slope (figure) error, COMP=3 only * * %E *******************************************************************************/ -DEFINE INSTRUMENT Test_Mirror_conic(int COMP=0, theta=1, length=1, m=2, R0=0.99, Qc=0.0217, alpha=6.07, W=0.003, Lmin=1, Lmax=11) +DEFINE INSTRUMENT Test_Mirror_conic(int COMP=0, theta=1, length=1, m=2, R0=0.99, Qc=0.0217, alpha=6.07, W=0.003, Lmin=1, Lmax=11, sigma_m=0, sigma_s=0) DECLARE %{ @@ -102,7 +104,7 @@ AT (0, 0, length/2) RELATIVE mirror_arm ROTATED (0, 90, 0) RELATIVE mirror_arm COMPONENT mirror_conic_conic = Mirror_conic(length=length, yheight=0.1, conic=1, - R0=R0, Qc=Qc, alpha=alpha, m=m, W=W) + R0=R0, Qc=Qc, alpha=alpha, m=m, W=W, sigma_m=sigma_m, sigma_s=sigma_s) WHEN (COMP==3) AT (0, 0, 0) RELATIVE mirror_arm diff --git a/mcstas-comps/share/conic.h b/mcstas-comps/share/conic.h index abe4841be..2df0dc3b2 100755 --- a/mcstas-comps/share/conic.h +++ b/mcstas-comps/share/conic.h @@ -177,6 +177,14 @@ int main() { #endif //! Max number of ConicSurf allowed in a Scene +/* Default RMS figure (slope) errors of ConicSurf mirrors, see reflectNeutronConic() */ +#ifndef CONIC_SIGMA_M +#define CONIC_SIGMA_M 1e-6 +#endif +#ifndef CONIC_SIGMA_S +#define CONIC_SIGMA_S 250e-6 +#endif + #ifndef MAX_CONICSURF #define MAX_CONICSURF 100 #endif @@ -544,6 +552,8 @@ typedef struct { double R0; //!< R0 reflectivity for mirror (0.99 for Nickel substrate) double alpha; //!< alpha, slope of reflectivity (6.07 AA for m=2 mirror) double W; //!< Width of supermirror cut-off, 0.003 AA-1 for m=2 mirror + double figerr_m; //!< RMS meridional slope error [rad] (figure error), default CONIC_SIGMA_M + double figerr_s; //!< RMS sagittal slope error [rad] (figure error), default CONIC_SIGMA_S //Only for reference double f1; //!< z-axis position of first focus @@ -945,6 +955,8 @@ ConicSurf makeHyperboloid(double f1, double f2, Point p, s.Qc = Qc; s.alpha = alpha; s.W = W; + s.figerr_m = CONIC_SIGMA_M; + s.figerr_s = CONIC_SIGMA_S; s.f1 = f1; s.f2 = f2; s.a = a; @@ -992,6 +1004,8 @@ ConicSurf makeEllipsoid(double f1, double f2, Point p, s.Qc = Qc; s.alpha = alpha; s.W = W; + s.figerr_m = CONIC_SIGMA_M; + s.figerr_s = CONIC_SIGMA_S; s.f1 = f1; s.f2 = f2; s.a = a; @@ -1091,6 +1105,8 @@ ConicSurf makeParaboloid(double f, Point p, double zstart, s.Qc = Qc; s.alpha = alpha; s.W = W; + s.figerr_m = CONIC_SIGMA_M; + s.figerr_s = CONIC_SIGMA_S; s.f1 = f; s.a = a; @@ -1517,44 +1533,20 @@ before computing reflection double reflectNeutronConic(_class_particle* _particle, ConicSurf s) { Vec n = getNormConic(get_class_particlePos(*_particle),s); Vec pv = get_class_particleVel(*_particle); - - int disp = 0; - if (disp>0) { - printf("\n n = %f, %f, %f",n.x,n.y,n.z); - printf("\n pv = %f, %f, %f",pv.x,pv.y,pv.z); - } - //add figure error by D. Liu - Vec n_p = makeVec(-n.y/sqrt(n.x*n.x+n.y*n.y),n.x/sqrt(n.x*n.x+n.y*n.y),0); //define a vector perpendicular to norm - //Vec n_p = makeVec(0,-n.z/sqrt(n.z*n.z+n.y*n.y),n.y/sqrt(n.z*n.z+n.y*n.y)); //define a vector perpendicular to norm - //Vec n_p = makeVec(-n.z/sqrt(n.z*n.z+n.x*n.x),0,n.x/sqrt(n.z*n.z+n.x*n.x)); //define a vector perpendicular to norm - if (disp>0) { - printf("\n n_p %f, %f, %f", n_p.x, n_p.y, n_p.z); - } - double len_p2 = sqrt((n.x*n.x+n.y*n.y)*(n.x*n.x+n.y*n.y+n.z*n.z)); - Vec n_p2 = makeVec(-n.x*n.z/len_p2, -n.y*n.z/len_p2, (n.x*n.x+n.y*n.y)/len_p2); //define a vector perpendicular to both n and n_p - - //double phi, theta = rand01()*2*M_PI; - double phi1, theta1, sigma1; - //sigma = 37.7e-6; - sigma1 = 1e-6; // sigma1 is half-sigma (unit in rad): sigma1=5e-6 means the FWHM is 23.6 urad - phi1 = sigma1*randnorm();//randgaussian(0.0,sigma1); - theta1 = 250*sigma1*randnorm();//randgaussian(0.0,sigma1); - if (disp>0) { - printf("\n phi and theta %f, %f", phi1, theta1); - } - - Vec R1 = rotateM(n,phi1,n_p); //rotate n by phi according to n_p - if (disp>0) { - printf("\n R1 %f, %f, %f", R1.x, R1.y, R1.z); - } - Vec R2 = rotateM(R1,theta1,n_p2); //rotate R1 by theta according to n_p2 - //Vec R2 = rotateM(R1,theta1,n); //rotate R1 by theta according to n_p2 - n = R2; - // - //n.y = n.y+theta1; - //n.z = n.z; - - // end of adding figure error + + //add figure error by D. Liu: Gaussian slope errors with RMS s.figerr_m (meridional) + //and s.figerr_s (sagittal), i.e. reflected rays deviate by twice these values + double nr = sqrt(n.x*n.x+n.y*n.y); + if ((s.figerr_m > 0 || s.figerr_s > 0) && nr > 0) { + Vec n_p = makeVec(-n.y/nr,n.x/nr,0); //define a vector perpendicular to norm + double len_p2 = sqrt((n.x*n.x+n.y*n.y)*(n.x*n.x+n.y*n.y+n.z*n.z)); + Vec n_p2 = makeVec(-n.x*n.z/len_p2, -n.y*n.z/len_p2, (n.x*n.x+n.y*n.y)/len_p2); //define a vector perpendicular to both n and n_p + double phi1 = s.figerr_m*randnorm(); + double theta1 = s.figerr_s*randnorm(); + Vec R1 = rotateM(n,phi1,n_p); //rotate n by phi according to n_p + n = rotateM(R1,theta1,n_p2); //rotate R1 by theta according to n_p2 + } + // end of adding figure error double v = getMagVec(pv); double vn = dotVec(pv,n); @@ -1580,9 +1572,6 @@ double reflectNeutronConic(_class_particle* _particle, ConicSurf s) { double par[5] = {s.R0, s.Qc, s.alpha, s.m, s.W}; StdReflecFunc(q, par, &ref); _particle->p = _particle->p * ref; - if (disp>0) { - printf("\n pv final = %f, %f, %f",pv.x,pv.y,pv.z); - } } return ga; } From 84ea859ddf144535676aec4b659137acb70b0e7a Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Tue, 13 Oct 2026 07:47:30 +0200 Subject: [PATCH 3/6] conic.h: use StdReflecFunc only, drop hard grazing-angle cut-off reflectNeutronConic() absorbed neutrons with grazing angle above m*theta_c(Ni), ignoring the Qc parameter and cutting the tanh tail of StdReflecFunc. The reflectivity is now given by StdReflecFunc alone, and neutrons with R=0 are absorbed (as in Mirror.comp). The unused cut-off in reflectNeutronFlat() is removed as well. Small numerical change: slightly more intensity just above m*Qc. Mirror_conic conic=1 now agrees exactly with Mirror(transmit=0). Co-Authored-By: Claude Opus 5.5 --- .../Test_Mirror_conic/Test_Mirror_conic.instr | 7 +++--- mcstas-comps/share/conic.h | 24 ++++++++----------- 2 files changed, 13 insertions(+), 18 deletions(-) diff --git a/mcstas-comps/examples/Tests_optics/Test_Mirror_conic/Test_Mirror_conic.instr b/mcstas-comps/examples/Tests_optics/Test_Mirror_conic/Test_Mirror_conic.instr index 3d9aa77c5..d7676a9bf 100644 --- a/mcstas-comps/examples/Tests_optics/Test_Mirror_conic/Test_Mirror_conic.instr +++ b/mcstas-comps/examples/Tests_optics/Test_Mirror_conic/Test_Mirror_conic.instr @@ -21,16 +21,15 @@ * COMP=2: Mirror, transmit=0 (reference for COMP=3) * COMP=3: Mirror_conic, conic=1 (ConicSurf, reflectNeutronConic) * -* COMP=0 and COMP=1 should agree exactly. COMP=3 agrees with COMP=2 except above -* the cut-off: reflectNeutronConic() absorbs neutrons with grazing angles beyond -* m times the critical angle of Ni, whereas Mirror keeps the tail of width W. +* COMP=0 and COMP=1 should agree exactly, and so should COMP=2 and COMP=3 when +* no figure error is applied (sigma_m=sigma_s=0). * * %Example: COMP=0 Detector: L_refl_I=8.98208e-07 * %Example: COMP=1 Detector: L_refl_I=8.98208e-07 * %Example: COMP=0 Detector: L_trans_I=7.01789e-07 * %Example: COMP=1 Detector: L_trans_I=7.01789e-07 * %Example: COMP=2 Detector: L_refl_I=8.9763e-07 -* %Example: COMP=3 Detector: L_refl_I=8.83094e-07 +* %Example: COMP=3 Detector: L_refl_I=8.9763e-07 * * %P * COMP: [1] Mirror implementation, see description diff --git a/mcstas-comps/share/conic.h b/mcstas-comps/share/conic.h index 2df0dc3b2..6944deea3 100755 --- a/mcstas-comps/share/conic.h +++ b/mcstas-comps/share/conic.h @@ -1518,7 +1518,7 @@ double calcSupermirrorReflectivity(double q, double m, double R_0, double Q_c){ /*! \brief Function to handle reflection of neutron for a ConicSurf. -@note Uses step function for reflectivity +@note Uses StdReflecFunc for reflectivity, applied as a weight factor @warning Make sure particle has been moved to surface of mirror before computing reflection @@ -1558,21 +1558,18 @@ double reflectNeutronConic(_class_particle* _particle, ConicSurf s) { } double ga = fabs(acos(vn/v)) - M_PI/2; - double gc = 6.84459399932*s.m/v; double ref=1.0; - if (ga > gc) { + double q = V2Q*(-2)*vn; + double par[5] = {s.R0, s.Qc, s.alpha, s.m, s.W}; + StdReflecFunc(q, par, &ref); + if (ref <= 0) { absorb_class_particle(_particle); return -1; } - else { - _particle->vx = _particle->vx-2*vn*n.x; - _particle->vy = _particle->vy-2*vn*n.y; - _particle->vz = _particle->vz-2*vn*n.z; - double q = V2Q*(-2)*vn; - double par[5] = {s.R0, s.Qc, s.alpha, s.m, s.W}; - StdReflecFunc(q, par, &ref); - _particle->p = _particle->p * ref; - } + _particle->vx = _particle->vx-2*vn*n.x; + _particle->vy = _particle->vy-2*vn*n.y; + _particle->vz = _particle->vz-2*vn*n.z; + _particle->p = _particle->p * ref; return ga; } @@ -1621,7 +1618,7 @@ void refractNeutronFlat(_class_particle* p, Vec n, double m1, double m2) { /*! \brief Function to handle reflection of neutron for a FlatSurf. -@note Uses step function for reflectivity +@note Uses StdReflecFunc for reflectivity, applied as the probability for reflection @warning Make sure particle has been moved to surface of mirror before computing reflection @@ -1643,7 +1640,6 @@ double reflectNeutronFlat(_class_particle* _particle, FlatSurf s) { //Hitting shell from outside For FlatSurface this has to be checked // make it able to reflect from the outside double ga = fabs(acos(vn/v)) - M_PI/2; - double gc = 6.84459399932*s.m/v; double ref=1.0; double q = V2Q*(-2)*vn; From 33ea7b96dec78714466ed655af8c5941f58bd23d Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Tue, 13 Oct 2026 07:47:30 +0200 Subject: [PATCH 4/6] conic.h: remove dead negative-reflectivity branch in reflectNeutronFlat Co-Authored-By: Claude Opus 5.5 --- mcstas-comps/share/conic.h | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/mcstas-comps/share/conic.h b/mcstas-comps/share/conic.h index 6944deea3..ce0543031 100755 --- a/mcstas-comps/share/conic.h +++ b/mcstas-comps/share/conic.h @@ -1646,11 +1646,7 @@ double reflectNeutronFlat(_class_particle* _particle, FlatSurf s) { double par[5] = {s.R0, s.Qc, s.alpha, s.m, s.W}; StdReflecFunc(q, par, &ref); - if (ref < 0) { - printf("this happens?"); - absorb_class_particle(_particle); - return -1; - } else {//here we need to implement the refraction + {//here we need to implement the refraction if (rand01() <= ref){//to be updated to use the mcstas random function or quasoi deterministic model //printf("oh a reflections\n"); //printf("this total reflection?"); From a9d899e7b0792b115b738e3cc5d74ffbfe229bfc Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Tue, 13 Oct 2026 07:47:30 +0200 Subject: [PATCH 5/6] FlatEllipse_finite_mirror: apply sourceDist default before it is used With sourceDist unset (0), the mirror spacing was calculated by get_r_at_z0() for a source at the mirror centre, since the default sourceDist=LStart was applied only afterwards. Also allocate with sizeof(double) instead of sizeof(double_t). Affects NMO as well (inherits INITIALIZE). Bundled instruments set sourceDist explicitly and are unchanged. Co-Authored-By: Claude Opus 5.5 --- mcstas-comps/contrib/FlatEllipse_finite_mirror.comp | 10 +++++----- mcstas-comps/contrib/NMO.comp | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/mcstas-comps/contrib/FlatEllipse_finite_mirror.comp b/mcstas-comps/contrib/FlatEllipse_finite_mirror.comp index 7356b3e73..bbeda4463 100644 --- a/mcstas-comps/contrib/FlatEllipse_finite_mirror.comp +++ b/mcstas-comps/contrib/FlatEllipse_finite_mirror.comp @@ -17,7 +17,7 @@ * The component relies on an updated version of conic.h from MIT. * * %P -* sourceDist: [m] Distance used for calculating the spacing of the mirrors +* sourceDist: [m] Distance used for calculating the spacing of the mirrors, 0: use LStart * LStart: [m] Left focal point * LEnd: [m] Right focal point * lStart: [m] z-Value of the mirror start @@ -86,7 +86,7 @@ SHARE double* get_r_at_z0 (int number, double z_0, double r_0, double z_extract, double LStart, double LEnd, double lStart, double lEnd) { int n = number; - double* r_zExtracts = malloc (n * sizeof (double_t)); /* n is an array of 10 integers */ + double* r_zExtracts = malloc (n * sizeof (double)); /* n is an array of 10 integers */ if (!r_zExtracts) { fprintf (stderr, "NMO comp function get_r_at_z0: malloc() failed. Exit! \n"); exit (-1); @@ -143,6 +143,9 @@ DECLARE INITIALIZE %{ + if (sourceDist == 0) { // default: mirror spacing is calculated for a source at the left focal point + sourceDist = LStart; + } if (rfront_inner_file && strlen (rfront_inner_file) && strcmp (rfront_inner_file, "NULL") && strcmp (rfront_inner_file, "0")) { if (Table_Read (&rsTable, rfront_inner_file, 1) <= 0) { /* read 1st block data from file into pTable */ exit (fprintf (stderr, "FlatEllipse_finite_mirror: %s: can not read file %s\n", NAME_CURRENT_COMP, rfront_inner_file)); @@ -163,9 +166,6 @@ INITIALIZE rfront_inner = get_r_at_z0 (nummirror, 0, r_0, lStart, sourceDist, LEnd, lStart, lEnd); // calculate the r-distances of all mirrors at the entry of the NMO, we will need this later } - if (sourceDist == 0) { // obsolete? - sourceDist = LStart; - } silicon = (mirror_width == 0) ? 0 : -1; // neutron starts in air by default // Load Reflectivity Data File TODO diff --git a/mcstas-comps/contrib/NMO.comp b/mcstas-comps/contrib/NMO.comp index fa989b0eb..b8dcad7d1 100644 --- a/mcstas-comps/contrib/NMO.comp +++ b/mcstas-comps/contrib/NMO.comp @@ -19,7 +19,7 @@ * NMO.comp contains no actual code, uses INHERIT of all sections from FlatEllipse_finite_mirror * * %P -* sourceDist: [m] Distance used for calculating the spacing of the mirrors +* sourceDist: [m] Distance used for calculating the spacing of the mirrors, 0: use LStart * LStart: [m] Left focal point * LEnd: [m] Right focal point * lStart: [m] z-Value of the mirror start From 70973d07d68b91ad15c584c403f382d84ad93f7c Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Tue, 13 Oct 2026 07:47:30 +0200 Subject: [PATCH 6/6] Test_Conics_pairs: update reference values Re-recorded (1e6 rays, seed 1) after the conic.h reflectivity fixes of #2789 and this branch. Co-Authored-By: Claude Opus 5.5 --- .../Test_Conics_pairs/Test_Conics_pairs.instr | 24 +++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/mcstas-comps/examples/Tests_optics/Test_Conics_pairs/Test_Conics_pairs.instr b/mcstas-comps/examples/Tests_optics/Test_Conics_pairs/Test_Conics_pairs.instr index 739374510..fe92c5028 100644 --- a/mcstas-comps/examples/Tests_optics/Test_Conics_pairs/Test_Conics_pairs.instr +++ b/mcstas-comps/examples/Tests_optics/Test_Conics_pairs/Test_Conics_pairs.instr @@ -12,18 +12,18 @@ * %D * Example instrument that shows some ways of using the Conics pair components * - * %Example: Test_Conics_pairs OPTIC=1 Detector: psd_i_I=6.91402e-13 - * %Example: Test_Conics_pairs OPTIC=1 Detector: psd_i_reflections_2_I=6.91402e-13 - * %Example: Test_Conics_pairs OPTIC=2 fs=100000 Detector: psd_i_I=2.85546e-21 - * %Example: Test_Conics_pairs OPTIC=2 fs=100000 Detector: psd_i_reflections_2_I=2.85515e-21 - * %Example: Test_Conics_pairs OPTIC=3 Detector: psd_i_I=6.91402e-13 - * %Example: Test_Conics_pairs OPTIC=4 ssize=1e-4 Detector: psd_i_I=8.66267e-09 - * %Example: Test_Conics_pairs quadratic=1 Detector: psd_i_I=5.03026e-13 - * %Example: Test_Conics_pairs quadratic=1 fs=8 fi=6 Detector: psd_i_I=1.26085e-12 - * %Example: Test_Conics_pairs OPTIC=2 quadratic=1 fs=100000 Detector: psd_i_I=2.09716e-21 - * %Example: Test_Conics_pairs OPTIC=3 quadratic=1 Detector: psd_i_I=5.03026e-13 - * %Example: Test_Conics_pairs OPTIC=3 quadratic=1 fs=8 fi=6 Detector: psd_i_I=9.00287e-13 - * %Example: Test_Conics_pairs OPTIC=4 quadratic=1 ssize=1e-4 Detector: psd_i_I=6.56678e-09 + * %Example: Test_Conics_pairs OPTIC=1 Detector: psd_i_I=6.71864e-13 + * %Example: Test_Conics_pairs OPTIC=1 Detector: psd_i_reflections_2_I=6.71864e-13 + * %Example: Test_Conics_pairs OPTIC=2 fs=100000 Detector: psd_i_I=2.75574e-21 + * %Example: Test_Conics_pairs OPTIC=2 fs=100000 Detector: psd_i_reflections_2_I=2.75574e-21 + * %Example: Test_Conics_pairs OPTIC=3 Detector: psd_i_I=6.71864e-13 + * %Example: Test_Conics_pairs OPTIC=4 ssize=1e-4 Detector: psd_i_I=8.42691e-09 + * %Example: Test_Conics_pairs quadratic=1 Detector: psd_i_I=4.93027e-13 + * %Example: Test_Conics_pairs quadratic=1 fs=8 fi=6 Detector: psd_i_I=1.2236e-12 + * %Example: Test_Conics_pairs OPTIC=2 quadratic=1 fs=100000 Detector: psd_i_I=2.03618e-21 + * %Example: Test_Conics_pairs OPTIC=3 quadratic=1 Detector: psd_i_I=4.93027e-13 + * %Example: Test_Conics_pairs OPTIC=3 quadratic=1 fs=8 fi=6 Detector: psd_i_I=8.80918e-13 + * %Example: Test_Conics_pairs OPTIC=4 quadratic=1 ssize=1e-4 Detector: psd_i_I=6.41135e-09 * * %P * OPTIC: [ ] Flag to choose between 1: EH pair, 2: PH pair, 3: HE pair, 4: PP pair