From 440bb8cca1dcf87101097c3058dde99149715c87 Mon Sep 17 00:00:00 2001 From: SilasSchack <93326018+SSBNS@users.noreply.github.com> Date: Mon, 10 Nov 2025 13:52:25 +0100 Subject: [PATCH 01/14] Update Phonon_simple.comp Changed focus_aw and focus_ah from deg to rad to comply with randvec_target_angular. --- mcstas-comps/samples/Phonon_simple.comp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mcstas-comps/samples/Phonon_simple.comp b/mcstas-comps/samples/Phonon_simple.comp index 557205eb10..90861b427d 100644 --- a/mcstas-comps/samples/Phonon_simple.comp +++ b/mcstas-comps/samples/Phonon_simple.comp @@ -44,8 +44,8 @@ * focus_r: [m] Radius of sphere containing target. * focus_xw: [m] horiz. dimension of a rectangular area * focus_yh: [m] vert. dimension of a rectangular area -* focus_aw: [deg] horiz. angular dimension of a rectangular area -* focus_ah: [deg] vert. angular dimension of a rectangular area +* focus_aw: [rad] horiz. angular dimension of a rectangular area +* focus_ah: [rad] vert. angular dimension of a rectangular area * target_x: [m] position of target to focus at . Transverse coordinate * target_y: [m] position of target to focus at. Vertical coordinate * target_z: [m] position of target to focus at. Straight ahead. From 8367f33df1bfe9ddeadabfd36f19ad7bb5b46e50 Mon Sep 17 00:00:00 2001 From: SilasSchack Date: Fri, 27 Mar 2026 10:55:05 +0100 Subject: [PATCH 02/14] Added SpinWave_BCO to branch Added both component to contrib and test instrument to examples --- Samples_SpinWave_BCO.instr | 248 +++++++ mcstas-comps/contrib/SpinWave_BCO.comp | 638 ++++++++++++++++++ .../Samples_SpinWave_BCO.instr | 250 +++++++ 3 files changed, 1136 insertions(+) create mode 100644 Samples_SpinWave_BCO.instr create mode 100644 mcstas-comps/contrib/SpinWave_BCO.comp create mode 100644 mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr diff --git a/Samples_SpinWave_BCO.instr b/Samples_SpinWave_BCO.instr new file mode 100644 index 0000000000..8cddd4a72c --- /dev/null +++ b/Samples_SpinWave_BCO.instr @@ -0,0 +1,248 @@ +/***************************************************************************** +* McStas instrument definition URL=http://www.mcstas.org +* +* Instrument: Samples_SpinWave_BCO +* +* %Identification +* Written by: Kim Lefmann +* Date: Feb 2004 +* Origin: RISOE +* Written by: K. Lefmann RISOE, Feb 2004 +* %INSTRUMENT_SITE: Tests_samples +* +* Simple test instrument for the SpinWave_BCO component +* +* +* %Description +* Simple test instrument for the SpinWave_BCO component. +* Refer to the component documentation for further instructions. +* +* Example: h=1 l=0 -n 1e5 Detector: mon1_I=2.86265e-25 +* +* %Parameters +* E: [meV] Mean energy at source +* DE: [meV] Energy spread at source +* HDIV: [deg] Horizontal divergence produced from source +* VDIV: [deg] Vertical divergence produced from source +* TT: [deg] Two-theta detetector-angle +* OM: [deg] Sample rotation angle +* C: [meV/AA^(-1)] Sample velocity of sound +* +* %End +******************************************************************************/ +DEFINE INSTRUMENT Samples_SpinWave_BCO(E=1.06,Ef=14.7,Dlambda=0.1, h=1, l=0,B=0,dA3=-90, Temp=2, width = 0.005, E_steps_low = 50, E_steps_high = 50,mode = 2,verbose=0) + +DECLARE +%{ + double Gqx,Gqy,Gqz; +//double Ef=24.8 +double Ei; +//meV +double thetaM; +double twothetaS; +double thetaA; +double A3; +double QM; +double alpha; +double lambda_i; +double SMALL__NUMBER; +double a; +double c; +double jAFM; +double jFM; +double jextra; +double S; +double D; +%} + +INITIALIZE +%{ +// Set monochromator/analyzer Q-value for PG +QM = 1.8734; +SMALL__NUMBER = 1e-6; + +// MnF2 parameters +// Lattice vectors in AA +a = 4.873; +c = 3.130; + +jAFM = 2*0.152; +jFM = -2*0.028; +jextra = 2*0.004; +S = 2.5; +D = -0.023; + +//MnF2 reciprocal lattice vectors, in 1/AA +double astar = 2*PI/a; +double cstar = 2*PI/c; + +//calculate Ei +Ei=Ef+E; + +//calculate ki, kf, lambda_i, q +double ki = V2K*SE2V*sqrt(Ei); +double kf = V2K*SE2V*sqrt(Ef); +lambda_i=2*PI/ki; +double q = sqrt(h*h*astar*astar+l*l*cstar*cstar); + +//calculate 2thetaM and 2thetaA +thetaM = asin(QM/(2*ki))*RAD2DEG; +thetaA = asin(QM/(2*kf))*RAD2DEG; + +//calculate scattering angle and sample rotation +twothetaS = acos((q*q-ki*ki-kf*kf)/(-2*ki*kf))*RAD2DEG; +alpha = acos((kf*kf-ki*ki-q*q)/(-2*ki*q)); +A3=(asin(l*cstar/(q+SMALL__NUMBER))-alpha)*RAD2DEG; + +%} + +TRACE + +COMPONENT origin = Progress_bar() +AT (0, 0, 0) RELATIVE ABSOLUTE + +COMPONENT source = Source_Maxwell_3( + xwidth=0.01, + yheight=0.01, + Lmin=lambda_i-Dlambda/2, + Lmax=lambda_i+Dlambda/2, + dist=7.5, + focus_yh=width, + focus_xw=width, + T1=300, + T2=300, + T3=300, + I1=1E15, + I2=1E15, + I3=1E15) +AT (0, 0, 0) RELATIVE PREVIOUS + +COMPONENT monochromator_flat = Monochromator_flat( + mosaicv=30, + mosaich=30, + zwidth = width, + yheight = width, + Q=QM) +AT (0, 0, 7.5) RELATIVE source +ROTATED (0, thetaM, 0) RELATIVE source + +COMPONENT arm1 = Arm() +AT (0, 0, 0) RELATIVE PREVIOUS +ROTATED (0, thetaM, 0) RELATIVE PREVIOUS + +COMPONENT collimator_linear1 = Collimator_linear( + xwidth=0.1, + yheight=0.2, + length=0.2, + divergence=40) +AT (0, 0, 0.7) RELATIVE arm1 + +/*COMPONENT e_monitor_before_sample = E_monitor( + nE=100, + filename="E_before_sample", + Emin=Ei-2, + Emax=Ei+2, + xwidth=width, + yheight=2*width, + restore_neutron=1) +AT (0, 0, 0.99) RELATIVE arm1*/ + +COMPONENT arm2 = Arm() +AT (0, 0, 1) RELATIVE arm1 +ROTATED (0, 0, 0) RELATIVE arm1 + +SPLIT 1 + +COMPONENT SAMPLE = SpinWave_BCO( + radius = width/2, + yheight = 2*width, + sigma_abs = 0, + sigma_inc = 0, + T = Temp, + a1 = a, + a2 = a, + a3 = c, + j = jAFM, + jc = jFM, + ja=jextra, + jb=jextra, + S=S, + B=B, + D=D, + FM=0, + mode_input = mode, + e_steps_low = E_steps_low, + e_steps_high = E_steps_high, + target_index=3, + focus_xw=0.005, + focus_yh=0.005, + verbose = verbose) +AT (0, 0, 0) RELATIVE arm2 +ROTATED (0, -A3+dA3, 180) RELATIVE arm2 + +COMPONENT arm3 = Arm() +AT (0, 0, 0) RELATIVE arm2 +ROTATED (0, -twothetaS, 0) RELATIVE arm2 + +COMPONENT collimator_linear2 = Collimator_linear( + xwidth=0.1, + yheight=0.2, + length=0.2, + divergence=40, + xwidth=width, + yheight=width) +AT (0, 0, 0.5) RELATIVE arm3 + +COMPONENT slit = Slit( + xwidth=0.005, + yheight=0.005) +AT (0, 0, 1) RELATIVE arm3 + +/*COMPONENT e_monitor_beforeana = E_monitor( + nE=250, + filename="Ebeforeana", + xwidth=width, + yheight=width, + Emin=0, + Emax=25, + restore_neutron=1) +AT (0, 0, 0.001) RELATIVE PREVIOUS*/ + +COMPONENT analyzer = Monochromator_flat( + zwidth=width, + yheight=width, + mosaicv=30, + mosaich=30, + Q=QM + ) +AT (0, 0, 0.1) RELATIVE PREVIOUS +ROTATED (0, thetaA, 0) RELATIVE PREVIOUS + +COMPONENT arm4 = Arm() +AT (0, 0, 0) RELATIVE PREVIOUS +ROTATED (0, thetaA, 0) RELATIVE PREVIOUS + +COMPONENT collimator_linear3 = Collimator_linear( + xwidth=0.1, + yheight=0.2, + length=0.2, + divergence=40) +AT (0, 0, 0.5) RELATIVE PREVIOUS + +COMPONENT Emon = E_monitor( + nE=250, + filename="Emon", + xwidth=width, + yheight=width, + Emin=0, + Emax=75, + restore_neutron=1) +AT (0, 0, 1) RELATIVE arm4 + + +FINALLY +%{ +%} + +END + diff --git a/mcstas-comps/contrib/SpinWave_BCO.comp b/mcstas-comps/contrib/SpinWave_BCO.comp new file mode 100644 index 0000000000..166148068a --- /dev/null +++ b/mcstas-comps/contrib/SpinWave_BCO.comp @@ -0,0 +1,638 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2000 Risoe National Laboratory. +* +* %I +* Written by: Silas Schack (algoritme based on "phonon_simple" by Kim Lefmann) +* Date: 13.06.25 +* Origin: KU +* +* A sample for FM and AFM magnon scattering from a oI lattice based on cross section expressions from Marshall and Lovesey ch.9. +* +* %D +* Single-cylinder shape. +* Absorption included. +* No multiple scattering. +* No incoherent scattering emitted. +* No attenuation from coherent scattering. No Bragg scattering. +* Bravais lattice only. (i.e. just one spin per sublattice unit cell) +* +* Algorithm: +* 0. Always perform the scattering if possible (otherwise ABSORB) +* 1. Choose direction within a focusing solid angle +* 2. Choose mode seperated by applied B-field +* 3. Calculate the zeros of (E_i-E_f-hbar omega(kappa)) as a function of k_f +* 4. Choose one value of k_f (always at least one is possible!) +* 5. Perform the correct weight transformation +* +* %P +* INPUT PARAMETERS: +* radius: [m] Outer radius of sample in (x,z) plane +* yheight: [m] Height of sample in y direction +* sigma_abs: [barns] Absorption cross section at 2200 m/s per atom +* sigma_inc: [barns] Incoherent scattering cross section per atom +* a1: [AA] lattice constants of orthorhombic body centred lattice +* a2: [AA] lattice constants of orthorhombic body centred lattice +* a3: [AA] lattice constants of orthorhombic body centred lattice +* S: [1] spin +* j: [meV] coupling constant to body centered spins +* ja: [meV] coupling constant to nearest neighbours along a +* jb: [meV] coupling constant to nearest neighbours along b +* jc: [meV] coupling constant to nearest neighbours along c +* j_110: [meV] coupling constant to nearest neighbours +x+y +* j_110_prime: [meV] coupling constant to nearest neighbours +x-y +* D: [meV] uniaxial anisotropy; easy axis is c-axis. Non-positive value assumed +* B: [T] field strength of external magnetic field applied along z-axis +* FM: [1] tag for ferromagnetic system if FM==1 +* DW: [1] Debye-Waller factor +* T: [K] Temperature +* mode_input Index of mode(s) to simulate +* e_steps_low [1] Number of low-energy steps in zrid +* e_steps_high [1] Number of high-energy steps in zrid +* focus_r: [m] Radius of sphere containing target. +* focus_xw: [m] horiz. dimension of a rectangular area +* focus_yh: [m] vert. dimension of a rectangular area +* focus_aw: [deg] horiz. angular dimension of a rectangular area +* focus_ah: [deg] vert. angular dimension of a rectangular area +* target_x: [m] position of target to focus at. Transverse coordinate +* target_y: [m] position of target to focus at. Vertical coordinate +* target_z: [m] position of target to focus at. Straight ahead. +* target_index: [1] relative index of component to focus at, e.g. next is +1 +* verbose [1] Verbosity level +* +* CALCULATED PARAMETERS: +* V_0: [AA] Sublattice unit cell volume +* V_my_s: [m^-1] Attenuation factor due to incoherent scattering +* V_my_a_v: [m^-1] Attenuation factor due to absorbtion +* +******************************************************************************/ + +DEFINE COMPONENT SpinWave_BCO + +SETTING PARAMETERS (radius ,yheight ,sigma_abs =1,sigma_inc=0,a1,a2,a3,j,ja,jb,jc,j_110=0,j_110_prime=0,S,D,B,FM=0,Fq=1,mode_input=2,DW=1,T,e_steps_low,e_steps_high, +target_x=0, target_y=0, target_z=0, int target_index=0,focus_r=0,focus_xw=0,focus_yh=0,focus_aw=0,focus_ah=0,verbose=0) + +/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ +SHARE +%{ +#ifndef MAGNON_oI +#define MAGNON_oI $Revision$ +#define T2E (1/11.605) /* Kelvin to meV */ +#define B2E 0.11590188615547234 /*Tesla to meV; g*Bohr magneton*/ + +#pragma acc routine +double nbose(double omega, double T) /* Other name ?? */ + { + double nb; + + nb= (omega>0) ? 1+1/(exp(omega/(T*T2E))-1) : 1/(exp(-omega/(T*T2E))-1); + return nb; + } +#undef T2E +/* Routine types from Numerical Recipies book */ +#define UNUSED (-1.11e30) +#define MAXRIDD 60 + +void fatalerror_cpu(char *s) + { + fprintf(stderr,"%s \n",s); + exit(1); + } + +#pragma acc routine +void fatalerror(char *s) + { + #ifndef OPENACC + fatalerror_cpu(s); + #endif + } + + #pragma acc routine + double omega_q(double* parms) + { + /* dispersion in units of meV */ + double vi, vf, vv_x, vv_y, vv_z, vi_x, vi_y, vi_z; + double qx, qy, qz, J, J1, J0, J10, omega_magnon, omega_neutron; + double coherence_factor; + double a1, a2, a3; + double S,D; + double j, ja, jb, jc; + double j_110, j_110_prime; + double B; + double coherence_flag; + int branch,FM_TAG; + + vf=parms[0]; + vi=parms[1]; + vv_x=parms[2]; + vv_y=parms[3]; + vv_z=parms[4]; + vi_x=parms[5]; + vi_y=parms[6]; + vi_z=parms[7]; + a1 = parms[8]; + a2 = parms[9]; + a3 = parms[10]; + j = parms[11]; + ja = parms[12]; + jb = parms[13]; + jc = parms[14]; + S = parms[15]; + B = parms[16]; + D = parms[17]; + branch = parms[18]; + coherence_flag = parms[19]; + FM_TAG = parms[20]; + j_110 = parms[21]; + j_110_prime = parms[22]; + + qx=V2K*(vi_x-vf*vv_x); + qy=V2K*(vi_y-vf*vv_y); + qz=V2K*(vi_z-vf*vv_z); + + J = 8*j*cos(a1*qx/2)*cos(a2*qy/2)*cos(a3*qz/2); + J0 = 8*j; + J1 = 2*(ja*cos(qx*a1)+jb*cos(qy*a2)+jc*cos(qz*a3))+2*(j_110+j_110_prime)*cos(qx*a1)*cos(qy*a2); + J10 = 2*(ja+jb+jc)+2*(j_110+j_110_prime); + + if (FM_TAG == 0) { + omega_magnon = sqrt((S*(J0-J10+J1)-2*D*(S-0.5))*(S*(J0-J10+J1)-2*D*(S-0.5))-S*S*(J)*(J)) - (2*branch - 1)*(B*B2E+2*(j_110-j_110_prime)*sin(qx*a1)*sin(qy*a2)); /*AM dispersion*/ + if ((omega_magnon<0) || ((S*(J0-J10+J1)-2*D*(S-0.5))*(S*(J0-J10+J1)-2*D*(S-0.5))-S*S*(J)*(J)<0)) { + fatalerror("Unphysical dispersion: Check parameters"); + } + } else { + omega_magnon = S*(J1+J-(J10+J0))-D*(S-0.5)+B2E*abs(B); /*FM dispersion*/ + if (omega_magnon<0) { + fatalerror("Unphysical dispersion: Check parameters"); + } + } + omega_neutron = fabs(VS2E*(vi*vi-vf*vf)); /*Magnitude of energy transfer*/ + if (coherence_flag == 0) { + return (omega_magnon - omega_neutron); + } else { /*calculate coherence factor*/ + if (FM_TAG == 0) { + coherence_factor = 2*S*((S*(J0-J)-S*(J10-J1)-2*D*(S-0.5)))/(omega_magnon + (2*branch - 1)*(B*B2E+2*(j_110-j_110_prime)*sin(qx*a1)*sin(qy*a2))); + + } else { + coherence_factor = 2*S; + } + if (coherence_factor < 0) { + fatalerror("Negative coherence factor: Check parameters"); /*Stop at negative coherence factor from unphysical parms*/ + } + return coherence_factor; + } + } + + + +double zridd(double (*func)(double*), double x1, double x2, double *parms, double xacc) + { + int j; + double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; + + parms[0]=x1; + fl=(*func)(parms); + parms[0]=x2; + fh=(*func)(parms); + if (fl*fh >= 0) + { + if (fl==0) return x1; + if (fh==0) return x2; + return UNUSED; + } + else + { + xl=x1; + xh=x2; + ans=UNUSED; + for (j=1; j= fh ? 1.0 : -1.0)*fm/s); + if (fabs(xnew-ans) <= xacc) + return ans; + ans=xnew; + parms[0]=ans; + fnew=(*func)(parms); + if (fnew == 0.0) return ans; + if (fabs(fm)*SIGN(fnew) != fm) + { + xl=xm; + fl=fm; + xh=ans; + fh=fnew; + } + else + if (fabs(fl)*SIGN(fnew) != fl) + { + xh=ans; + fh=fnew; + } + else + if(fabs(fh)*SIGN(fnew) != fh) + { + xl=ans; + fl=fnew; + } + else + fatalerror("never get here in zridd"); + if (fabs(xh-xl) <= xacc) + return ans; + } + fatalerror("zridd exceeded maximum iterations"); + } + return 0.0; /* Never get here */ + } + +#pragma acc routine +double zridd_gpu(double x1, double x2, double *parms, double xacc) + { + int j; + double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; + + parms[0]=x1; + fl=omega_q(parms); + parms[0]=x2; + fh=omega_q(parms); + if (fl*fh >= 0) + { + if (fl==0) return x1; + if (fh==0) return x2; + return UNUSED; + } + else + { + xl=x1; + xh=x2; + ans=UNUSED; + for (j=1; j= fh ? 1.0 : -1.0)*fm/s); + if (fabs(xnew-ans) <= xacc) + return ans; + ans=xnew; + parms[0]=ans; + fnew=omega_q(parms); + if (fnew == 0.0) return ans; + if (fabs(fm)*SIGN(fnew) != fm) + { + xl=xm; + fl=fm; + xh=ans; + fh=fnew; + } + else + if (fabs(fl)*SIGN(fnew) != fl) + { + xh=ans; + fh=fnew; + } + else + if(fabs(fh)*SIGN(fnew) != fh) + { + xl=ans; + fl=fnew; + } + else + fatalerror("never get here in zridd"); + if (fabs(xh-xl) <= xacc) + return ans; + } + fatalerror("zridd exceeded maximum iterations"); + } + return 0.0; /* Never get here */ + } + + +#define ROOTACC 1e-8 + + int findroots(double brack_low, double brack_mid, double brack_high, double *list, int* index, double (*f)(double*), double *parms, int steps_low, int steps_high) + { + double root,range; + int i; + range = brack_mid-brack_low; + for (i=0; i<=(steps_low-1); i++) + { + + root = zridd(f, brack_low+range*i/(int)steps_low, + brack_low+range*(i+1)/(int)steps_low, + (double *)parms, ROOTACC); + if (root != UNUSED) + { + list[(*index)++]=root; + } + } + + range = brack_high-brack_mid; + + for (i=0; i<=(steps_high-1); i++) + { + root = zridd(f, brack_mid+range*i/(int)steps_high, + brack_mid+range*(i+1)/(int)steps_high, + (double *)parms, ROOTACC); + if (root != UNUSED) + { + list[(*index)++]=root; + } + } + } + +#pragma acc routine + int findroots_gpu(double brack_low, double brack_mid, double brack_high, double *list, int* index, double *parms, int steps_low, int steps_high) + { + double root,range; + int i; + + range = brack_mid - brack_low; + + for (i=0; i<=(steps_low-1); i++) + { + root = zridd_gpu(brack_low+range*i/(int)steps_low, + brack_low+range*(i+1)/(int)steps_low, + (double *)parms, ROOTACC); + if (root != UNUSED) + { + list[(*index)++]=root; + } + } + + range = brack_high - brack_mid; + + for (i=0; i<=(steps_high-1); i++) + { + root = zridd_gpu(brack_mid+range*i/(int)steps_high, + brack_mid+range*(i+1)/(int)steps_high, + (double *)parms, ROOTACC); + if (root != UNUSED) + { + list[(*index)++]=root; + } + } + } + +#undef UNUSED +#undef MAXRIDD +#endif +%} + +DECLARE +%{ + double V_0; + double V_my_s; + double V_my_a_v; + double DV; + double gamma_n; + double r_0; + double g; +%} +INITIALIZE +%{ + g = 2.0023; /*electron g-factor*/ + gamma_n = -1.913; /*neutron gamma-factor*/ + r_0 = 2.8179; /*electron charge radius in fm*/ + V_0 = a1*a2*a3; /*inverse unit cell volume of sublattice (oP)*/ + V_my_s = (2 * 100 * sigma_inc/V_0); /*Incoherent volume specific cross section.*/ + V_my_a_v = (2 * 100 * sigma_abs/V_0 * 2200); + DV = 0.001; /* Velocity change used for numerical derivative */ + + /* now compute target coords if a component index is supplied */ + if (!target_index && !target_x && !target_y && !target_z) target_index=1; + if (target_index){ + Coords ToTarget; + ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index),POS_A_CURRENT_COMP); + ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget); + coords_get(ToTarget, &target_x, &target_y, &target_z); + } + if (!(target_x || target_y || target_z)) { + printf("Magnon_oI: %s: The target is not defined. Using direct beam (Z-axis).\n", + NAME_CURRENT_COMP); + target_z=1; + } +%} +TRACE +%{ + double t0, t1; /* Entry/exit time for cylinder */ + double v_i, v_f; /* Neutron velocities: initial, final */ + double vx_i, vy_i, vz_i; /* Neutron initial velocity vector */ + double dt0, dt; /* Flight times through sample */ + double l_full; /* Flight path length for non-scattered neutron */ + double l_i, l_o; /* Flight path lenght in/out for scattered neutron */ + double my_a_i; /* Initial attenuation factor */ + double my_a_f; /* Final attenuation factor */ + double solid_angle; /* Solid angle of target as seen from scattering point */ + double aim_x=0, aim_y=0, aim_z=1; /* Position of target relative to scattering point */ + double kappa_x, kappa_y, kappa_z; /* Scattering vector */ + double kappa2; /* Square of the scattering vector */ + double kappa2_z_norm; + double bose_factor; /* Calculated value of the Bose factor */ + double omega; /* energy transfer */ + int nf, index; /* Number of allowed final velocities */ + double vf_list[20]; /* List of allowed final velocities */ + double J_factor; /* Jacobian from delta fnc.s in cross section */ + double coherence_factor; /* Coherence factor in cross section*/ + double f1,f2; /* probed values of omega_q minus omega */ + double w_atten,w_MC,w_magnon,w_Jacobi,w_coherence; /* temporary multipliers */ + int mode; /* index for mode selection */ + double E_max; /* Maximum upper bound of dispersion */ + double parms[23]; + + + if(cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius, yheight)) + { + if(t0 < 0) + ABSORB; /* Neutron came from the sample or begins inside */ + + /* Neutron enters at t=t0. */ + dt0 = t1-t0; /* Time in sample */ + v_i = sqrt(vx*vx + vy*vy + vz*vz); + l_full = v_i * dt0; /* Length of path through sample if not scattered */ + dt = rand01()*dt0; /* Time of scattering (relative to t0) */ + l_i = v_i*dt; /* Penetration in sample at scattering */ + vx_i=vx; + vy_i=vy; + vz_i=vz; + PROP_DT(dt+t0); /* Point of scattering */ + + aim_x = target_x-x; /* Vector pointing at target (e.g. analyzer) */ + aim_y = target_y-y; + aim_z = target_z-z; + + if(focus_aw && focus_ah) { + randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle, + aim_x, aim_y, aim_z, focus_aw*DEG2RAD, focus_ah*DEG2RAD, ROT_A_CURRENT_COMP); + } else if(focus_xw && focus_yh) { + randvec_target_rect(&vx, &vy, &vz, &solid_angle, + aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP); + } else { + randvec_target_circle(&vx,&vy,&vz,&solid_angle,aim_x,aim_y,aim_z, focus_r); + } + + NORM(vx, vy, vz); + nf=0; + + if (FM==0) { + if (B2E*abs(B) > 2*sqrt(D*D*(S-0.5)*(S-0.5)-8*j*D*S*(S-0.5))) { + printf("Magnetic field too strong. Fieldstrength set to zero\n"); /*Error if B induced spin-flop transition*/ + B=0; + } + } else { + if (B<0) { + fatalerror("Field direction flippes spins: Using abs(B)"); + } + } + + mode = 0; + if (FM==0) { + if (mode_input == 2) + mode = round(2*rand01()-0.5); /* Select mode randomly from 2 possibilities */ + else + mode = mode_input; + if ((mode < 0) || (mode > 1)) + { + printf("mode = %d ",mode); + fatalerror("illegal value of mode"); + } + } + + parms[1]=v_i; + parms[2]=vx; + parms[3]=vy; + parms[4]=vz; + parms[5]=vx_i; + parms[6]=vy_i; + parms[7]=vz_i; + parms[8]=a1; + parms[9]=a2; + parms[10]=a3; + parms[11]=j; + parms[12]=ja; + parms[13]=jb; + parms[14]=jc; + parms[15]=S; + parms[16]=B; + parms[17]=D; + parms[18]=mode; + parms[19]=0; + parms[20]=FM; + parms[21]=j_110; + parms[22]=j_110_prime; + + if (FM==0) { + E_max = fabs(S*(8*fabs(j)+4*(fabs(ja)+fabs(jb)+fabs(jc))+2*(fabs(j_110)+fabs(j_110_prime)))-2*D*(S-0.5))+fabs(B2E*B)+2*fabs(j_110-j_110_prime); + } else { + E_max = fabs(S*(2*(8*fabs(j)+2*(fabs(ja)+fabs(jb)+fabs(jc))))-D*(S-0.5))+fabs(B2E*B); + } + + #ifndef OPENACC + findroots(0, v_i, sqrt(v_i*v_i+1/VS2E*E_max)+20, vf_list, &nf, omega_q, parms, e_steps_low,e_steps_high); + #else + findroots_gpu(0, v_i, sqrt(v_i*v_i+1/VS2E*E_max)+20, vf_list, &nf, parms, e_steps_low,e_steps_high); /*+10 to make sure final point is included*/ + #endif + + index=(int)floor(rand01()*nf); + if (nf>0 && index<20) { + v_f = vf_list[index]; + parms[0] = v_f; + + parms[19] = 1; /*coherence flag*/ + coherence_factor = omega_q(parms); + parms[19] = 0; + if (verbose==1) { /* Print functions used for debugging */ + printf("vi=(%g,%g,%g)\n",vx_i,vy_i,vz_i); + printf("vf=(%g,%g,%g)\n",vx,vy,vz); + printf("nf=%d\n",nf); + printf("index=%d\n",index); + printf("omega_q found=%g\n vf=%g\n", omega_q(parms),v_f); + } + parms[0]=v_f-DV; + f1=omega_q(parms); + parms[0]=v_f+DV; + f2=omega_q(parms); + J_factor = fabs(f2-f1)/(2*DV); + omega=VS2E*(v_i*v_i-v_f*v_f); + vx *= v_f; + vy *= v_f; + vz *= v_f; + + kappa_x=V2K*(vx_i-vx); + kappa_y=V2K*(vy_i-vy); + kappa_z=V2K*(vz_i-vz); + + if (verbose==1) { + if (omega>0) { + printf("hkl=(%f,%f,%f)\n",kappa_x/(2*PI/a1),kappa_y/(2*PI/a2),kappa_z/(2*PI/a3)); + }} + + kappa2=kappa_y*kappa_y+kappa_x*kappa_x+kappa_z*kappa_z; + kappa2_z_norm = kappa_z*kappa_z/kappa2; + + if(!cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius, yheight)) + { + /* ??? did not hit cylinder */ + printf("FATAL ERROR: Did not hit cylinder from inside.\n"); + exit(1); + } + dt = t1; + l_o = v_f*dt; + + my_a_i = V_my_a_v/v_i; + my_a_f = V_my_a_v/v_f; + bose_factor=nbose(omega,T); + w_atten = exp(-(V_my_s*(l_i+l_o)+my_a_i*l_i+my_a_f*l_o)); /* Absorption factor */ + w_MC = nf*solid_angle*(l_full*1e10)/V_0; /* Focusing factors; assume random choice of n_f possibilities. Units = AA^-2 */ + + if (FM == 0) { + if (mode_input == 2) { + w_MC *= 2; /*n_m; number of branches/modes in MC choice*/ + } + } else { + w_MC *= 2; /*If FM==1, this accounts for v_0 only being half magnetic unit cell volume*/ + } + + w_magnon = (gamma_n*gamma_n*r_0*r_0/1e10)*(g*g*Fq*Fq/4)*(v_f/v_i)*DW*(1+kappa2_z_norm)*bose_factor/4; /*Constants and magnetic factors in cross section. Units = AA^2*/ + w_Jacobi = 2*VS2E*v_f/J_factor; /* Jacobian of delta functions in cross section. */ + w_coherence = coherence_factor; /* coherence factor */ + p *= w_atten*w_MC*w_magnon*w_Jacobi*w_coherence; + if (verbose ==2) { /* Print functions used for debugging absolute intensity measurements*/ + if (v_i > v_f) { + double delta_E = VS2E*(v_i*v_i-v_f*v_f); + printf("delta_E=%f\n",delta_E); + printf("hkl=(%f,%f,%f)\n",kappa_x/(2*PI/a1),kappa_y/(2*PI/a2),kappa_z/(2*PI/a3)); + printf("vvi=(%f,%f,%f)\n vvf=(%f,%f,%f)\n",vx_i/v_i,vy_i/v_i,vz_i/v_i,vx/v_f,vy/v_f,vz/v_f); + printf("Ei=%f. Ef=%f\n",VS2E*v_i*v_i,VS2E*v_f*v_f); + printf("J_factor=%f\n",J_factor); + printf("vi=%f, vf=%f\n w_Jacobi = %f\n w_coherence = %f\n w_magnon*10^10= %f\n", v_i,v_f,w_Jacobi,w_coherence,w_magnon*1e10); + printf("cross*10^10/N = %f\n",w_magnon*w_Jacobi*w_coherence*1e10); + } + } + + } else { + ABSORB; // findroots returned junk + } + } /* else transmit: Neutron did not hit the sample */ +%} + +MCDISPLAY +%{ + + 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); +%} + +END \ No newline at end of file diff --git a/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr b/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr new file mode 100644 index 0000000000..2800b0e9a5 --- /dev/null +++ b/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr @@ -0,0 +1,250 @@ +/***************************************************************************** +* McStas instrument definition URL=http://www.mcstas.org +* +* Instrument: Samples_SpinWave_BCO +* +* %Identification +* Written by: Kim Lefmann +* Date: Feb 2004 +* Origin: RISOE +* Written by: K. Lefmann RISOE, Feb 2004 +* %INSTRUMENT_SITE: Tests_samples +* +* Simple test instrument for the SpinWave_BCO component +* +* +* %Description +* Simple test instrument for the SpinWave_BCO component. +* Refer to the component documentation for further instructions. +* +* %Example: h=1 l=0 -n 1e6 Detector: Emon_I=0.0517575 +* %Example: h=0 l=0.5 E=6.73 -n 1e6 Detector: Emon_I=0.0108855 +* %Example: h=1 l=1 -n 1e6 Detector: Emon_I=0.000655744 + +* %Parameters +* E: [meV] Mean energy at source +* DE: [meV] Energy spread at source +* HDIV: [deg] Horizontal divergence produced from source +* VDIV: [deg] Vertical divergence produced from source +* TT: [deg] Two-theta detetector-angle +* OM: [deg] Sample rotation angle +* C: [meV/AA^(-1)] Sample velocity of sound +* +* %End +******************************************************************************/ +DEFINE INSTRUMENT Samples_SpinWave_BCO(E=1.06,Ef=14.7,Dlambda=0.1, h=1, l=0,B=0,dA3=-90, Temp=2, width = 0.005, E_steps_low = 50, E_steps_high = 50,mode = 2,verbose=0) + +DECLARE +%{ + double Gqx,Gqy,Gqz; +//double Ef=24.8 +double Ei; +//meV +double thetaM; +double twothetaS; +double thetaA; +double A3; +double QM; +double alpha; +double lambda_i; +double SMALL__NUMBER; +double a; +double c; +double jAFM; +double jFM; +double jextra; +double S; +double D; +%} + +INITIALIZE +%{ +// Set monochromator/analyzer Q-value for PG +QM = 1.8734; +SMALL__NUMBER = 1e-6; + +// MnF2 parameters +// Lattice vectors in AA +a = 4.873; +c = 3.130; + +jAFM = 2*0.152; +jFM = -2*0.028; +jextra = 2*0.004; +S = 2.5; +D = -0.023; + +//MnF2 reciprocal lattice vectors, in 1/AA +double astar = 2*PI/a; +double cstar = 2*PI/c; + +//calculate Ei +Ei=Ef+E; + +//calculate ki, kf, lambda_i, q +double ki = V2K*SE2V*sqrt(Ei); +double kf = V2K*SE2V*sqrt(Ef); +lambda_i=2*PI/ki; +double q = sqrt(h*h*astar*astar+l*l*cstar*cstar); + +//calculate 2thetaM and 2thetaA +thetaM = asin(QM/(2*ki))*RAD2DEG; +thetaA = asin(QM/(2*kf))*RAD2DEG; + +//calculate scattering angle and sample rotation +twothetaS = acos((q*q-ki*ki-kf*kf)/(-2*ki*kf))*RAD2DEG; +alpha = acos((kf*kf-ki*ki-q*q)/(-2*ki*q)); +A3=(asin(l*cstar/(q+SMALL__NUMBER))-alpha)*RAD2DEG; + +%} + +TRACE + +COMPONENT origin = Progress_bar() +AT (0, 0, 0) RELATIVE ABSOLUTE + +COMPONENT source = Source_Maxwell_3( + xwidth=0.01, + yheight=0.01, + Lmin=lambda_i-Dlambda/2, + Lmax=lambda_i+Dlambda/2, + dist=7.5, + focus_yh=width, + focus_xw=width, + T1=300, + T2=300, + T3=300, + I1=1E15, + I2=1E15, + I3=1E15) +AT (0, 0, 0) RELATIVE PREVIOUS + +COMPONENT monochromator_flat = Monochromator_flat( + mosaicv=30, + mosaich=30, + zwidth = width, + yheight = width, + Q=QM) +AT (0, 0, 7.5) RELATIVE source +ROTATED (0, thetaM, 0) RELATIVE source + +COMPONENT arm1 = Arm() +AT (0, 0, 0) RELATIVE PREVIOUS +ROTATED (0, thetaM, 0) RELATIVE PREVIOUS + +COMPONENT collimator_linear1 = Collimator_linear( + xwidth=0.1, + yheight=0.2, + length=0.2, + divergence=40) +AT (0, 0, 0.7) RELATIVE arm1 + +/*COMPONENT e_monitor_before_sample = E_monitor( + nE=100, + filename="E_before_sample", + Emin=Ei-2, + Emax=Ei+2, + xwidth=width, + yheight=2*width, + restore_neutron=1) +AT (0, 0, 0.99) RELATIVE arm1*/ + +COMPONENT arm2 = Arm() +AT (0, 0, 1) RELATIVE arm1 +ROTATED (0, 0, 0) RELATIVE arm1 + +SPLIT 1 + +COMPONENT SAMPLE = SpinWave_BCO( + radius = width/2, + yheight = 2*width, + sigma_abs = 0, + sigma_inc = 0, + T = Temp, + a1 = a, + a2 = a, + a3 = c, + j = jAFM, + jc = jFM, + ja=jextra, + jb=jextra, + S=S, + B=B, + D=D, + FM=0, + mode_input = mode, + e_steps_low = E_steps_low, + e_steps_high = E_steps_high, + target_index=3, + focus_xw=0.005, + focus_yh=0.005, + verbose = verbose) +AT (0, 0, 0) RELATIVE arm2 +ROTATED (0, -A3+dA3, 180) RELATIVE arm2 + +COMPONENT arm3 = Arm() +AT (0, 0, 0) RELATIVE arm2 +ROTATED (0, -twothetaS, 0) RELATIVE arm2 + +COMPONENT collimator_linear2 = Collimator_linear( + xwidth=0.1, + yheight=0.2, + length=0.2, + divergence=40, + xwidth=width, + yheight=width) +AT (0, 0, 0.5) RELATIVE arm3 + +COMPONENT slit = Slit( + xwidth=0.005, + yheight=0.005) +AT (0, 0, 1) RELATIVE arm3 + +/*COMPONENT e_monitor_beforeana = E_monitor( + nE=250, + filename="Ebeforeana", + xwidth=width, + yheight=width, + Emin=0, + Emax=25, + restore_neutron=1) +AT (0, 0, 0.001) RELATIVE PREVIOUS*/ + +COMPONENT analyzer = Monochromator_flat( + zwidth=width, + yheight=width, + mosaicv=30, + mosaich=30, + Q=QM + ) +AT (0, 0, 0.1) RELATIVE PREVIOUS +ROTATED (0, thetaA, 0) RELATIVE PREVIOUS + +COMPONENT arm4 = Arm() +AT (0, 0, 0) RELATIVE PREVIOUS +ROTATED (0, thetaA, 0) RELATIVE PREVIOUS + +COMPONENT collimator_linear3 = Collimator_linear( + xwidth=0.1, + yheight=0.2, + length=0.2, + divergence=40) +AT (0, 0, 0.5) RELATIVE PREVIOUS + +COMPONENT Emon = E_monitor( + nE=250, + filename="Emon", + xwidth=width, + yheight=width, + Emin=0, + Emax=75, + restore_neutron=1) +AT (0, 0, 1) RELATIVE arm4 + + +FINALLY +%{ +%} + +END + From 9929ea73dda9a9b058697f11c2a4b16851a4d5ea Mon Sep 17 00:00:00 2001 From: SilasSchack Date: Fri, 27 Mar 2026 11:55:53 +0100 Subject: [PATCH 03/14] Update Phonon_simple to newest version --- mcstas-comps/contrib/SpinWave_BCO.comp | 7 ++++++- mcstas-comps/samples/Phonon_simple.comp | 10 +++++++--- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/mcstas-comps/contrib/SpinWave_BCO.comp b/mcstas-comps/contrib/SpinWave_BCO.comp index 166148068a..c2beb682f2 100644 --- a/mcstas-comps/contrib/SpinWave_BCO.comp +++ b/mcstas-comps/contrib/SpinWave_BCO.comp @@ -407,6 +407,11 @@ INITIALIZE V_my_a_v = (2 * 100 * sigma_abs/V_0 * 2200); DV = 0.001; /* Velocity change used for numerical derivative */ + if (focus_aw) + focus_aw *= DEG2RAD; + if (focus_ah) + focus_ah *= DEG2RAD; + /* now compute target coords if a component index is supplied */ if (!target_index && !target_x && !target_y && !target_z) target_index=1; if (target_index){ @@ -471,7 +476,7 @@ TRACE if(focus_aw && focus_ah) { randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle, - aim_x, aim_y, aim_z, focus_aw*DEG2RAD, focus_ah*DEG2RAD, ROT_A_CURRENT_COMP); + aim_x, aim_y, aim_z, focus_aw, focus_ah, ROT_A_CURRENT_COMP); } else if(focus_xw && focus_yh) { randvec_target_rect(&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP); diff --git a/mcstas-comps/samples/Phonon_simple.comp b/mcstas-comps/samples/Phonon_simple.comp index 378e3a5612..e7f4e33266 100644 --- a/mcstas-comps/samples/Phonon_simple.comp +++ b/mcstas-comps/samples/Phonon_simple.comp @@ -44,8 +44,8 @@ * focus_r: [m] Radius of sphere containing target. * focus_xw: [m] horiz. dimension of a rectangular area * focus_yh: [m] vert. dimension of a rectangular area -* focus_aw: [rad] horiz. angular dimension of a rectangular area -* focus_ah: [rad] vert. angular dimension of a rectangular area +* focus_aw: [deg] horiz. angular dimension of a rectangular area +* focus_ah: [deg] vert. angular dimension of a rectangular area * target_x: [m] position of target to focus at . Transverse coordinate * target_y: [m] position of target to focus at. Vertical coordinate * target_z: [m] position of target to focus at. Straight ahead. @@ -348,6 +348,10 @@ INITIALIZE V_my_s = (V_rho * 100 * sigma_inc); V_my_a_v = (V_rho * 100 * sigma_abs * 2200); DV = 0.001; /* Velocity change used for numerical derivative */ + if (focus_aw) + focus_aw *= DEG2RAD; + if (focus_ah) + focus_ah *= DEG2RAD; // Set constant parameters for parms object phonon.a_ = a; @@ -506,4 +510,4 @@ MCDISPLAY line (0, -yheight / 2.0, +radius, 0, +yheight / 2.0, +radius); %} -END +END \ No newline at end of file From 75f4f4347ed052c5048398756f99a98e0b36bcdb Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 8 Apr 2026 15:36:24 +0200 Subject: [PATCH 04/14] Rectify mode_input doc line --- mcstas-comps/contrib/SpinWave_BCO.comp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mcstas-comps/contrib/SpinWave_BCO.comp b/mcstas-comps/contrib/SpinWave_BCO.comp index c2beb682f2..bb95714dcf 100644 --- a/mcstas-comps/contrib/SpinWave_BCO.comp +++ b/mcstas-comps/contrib/SpinWave_BCO.comp @@ -47,7 +47,7 @@ * FM: [1] tag for ferromagnetic system if FM==1 * DW: [1] Debye-Waller factor * T: [K] Temperature -* mode_input Index of mode(s) to simulate +* mode_input: [1] Index of mode(s) to simulate * e_steps_low [1] Number of low-energy steps in zrid * e_steps_high [1] Number of high-energy steps in zrid * focus_r: [m] Radius of sphere containing target. @@ -640,4 +640,4 @@ MCDISPLAY line(0, -yheight/2.0, +radius, 0, +yheight/2.0, +radius); %} -END \ No newline at end of file +END From 079699654b60373b559a66037742a35932b73f9d Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 8 Apr 2026 15:38:01 +0200 Subject: [PATCH 05/14] Add template doc line for Fq --- mcstas-comps/contrib/SpinWave_BCO.comp | 1 + 1 file changed, 1 insertion(+) diff --git a/mcstas-comps/contrib/SpinWave_BCO.comp b/mcstas-comps/contrib/SpinWave_BCO.comp index bb95714dcf..71fc47e27e 100644 --- a/mcstas-comps/contrib/SpinWave_BCO.comp +++ b/mcstas-comps/contrib/SpinWave_BCO.comp @@ -45,6 +45,7 @@ * D: [meV] uniaxial anisotropy; easy axis is c-axis. Non-positive value assumed * B: [T] field strength of external magnetic field applied along z-axis * FM: [1] tag for ferromagnetic system if FM==1 +* Fq: [1] -> PLEASE FILL IN <- * DW: [1] Debye-Waller factor * T: [K] Temperature * mode_input: [1] Index of mode(s) to simulate From 1d3bbb7b5d3aa9b5cace43dd789b2fc0909ca150 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 8 Apr 2026 15:39:38 +0200 Subject: [PATCH 06/14] Add hint to elaborate Spin parameter doc / 1st char case changes in doc strings --- mcstas-comps/contrib/SpinWave_BCO.comp | 42 +++++++++++++------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/mcstas-comps/contrib/SpinWave_BCO.comp b/mcstas-comps/contrib/SpinWave_BCO.comp index 71fc47e27e..942e8d4f18 100644 --- a/mcstas-comps/contrib/SpinWave_BCO.comp +++ b/mcstas-comps/contrib/SpinWave_BCO.comp @@ -32,19 +32,19 @@ * yheight: [m] Height of sample in y direction * sigma_abs: [barns] Absorption cross section at 2200 m/s per atom * sigma_inc: [barns] Incoherent scattering cross section per atom -* a1: [AA] lattice constants of orthorhombic body centred lattice -* a2: [AA] lattice constants of orthorhombic body centred lattice -* a3: [AA] lattice constants of orthorhombic body centred lattice -* S: [1] spin -* j: [meV] coupling constant to body centered spins -* ja: [meV] coupling constant to nearest neighbours along a -* jb: [meV] coupling constant to nearest neighbours along b -* jc: [meV] coupling constant to nearest neighbours along c -* j_110: [meV] coupling constant to nearest neighbours +x+y -* j_110_prime: [meV] coupling constant to nearest neighbours +x-y -* D: [meV] uniaxial anisotropy; easy axis is c-axis. Non-positive value assumed -* B: [T] field strength of external magnetic field applied along z-axis -* FM: [1] tag for ferromagnetic system if FM==1 +* a1: [AA] Lattice constants of orthorhombic body centred lattice +* a2: [AA] Lattice constants of orthorhombic body centred lattice +* a3: [AA] Lattice constants of orthorhombic body centred lattice +* S: [1] Spin <-- PLEASE ELABORATE +* j: [meV] Coupling constant to body centered spins +* ja: [meV] Coupling constant to nearest neighbours along a +* jb: [meV] Coupling constant to nearest neighbours along b +* jc: [meV] Coupling constant to nearest neighbours along c +* j_110: [meV] Coupling constant to nearest neighbours +x+y +* j_110_prime: [meV] Coupling constant to nearest neighbours +x-y +* D: [meV] Uniaxial anisotropy; easy axis is c-axis. Non-positive value assumed +* B: [T] Field strength of external magnetic field applied along z-axis +* FM: [1] Tag for ferromagnetic system if FM==1 * Fq: [1] -> PLEASE FILL IN <- * DW: [1] Debye-Waller factor * T: [K] Temperature @@ -52,14 +52,14 @@ * e_steps_low [1] Number of low-energy steps in zrid * e_steps_high [1] Number of high-energy steps in zrid * focus_r: [m] Radius of sphere containing target. -* focus_xw: [m] horiz. dimension of a rectangular area -* focus_yh: [m] vert. dimension of a rectangular area -* focus_aw: [deg] horiz. angular dimension of a rectangular area -* focus_ah: [deg] vert. angular dimension of a rectangular area -* target_x: [m] position of target to focus at. Transverse coordinate -* target_y: [m] position of target to focus at. Vertical coordinate -* target_z: [m] position of target to focus at. Straight ahead. -* target_index: [1] relative index of component to focus at, e.g. next is +1 +* focus_xw: [m] Horiz. dimension of a rectangular area +* focus_yh: [m] Vert. dimension of a rectangular area +* focus_aw: [deg] Horiz. angular dimension of a rectangular area +* focus_ah: [deg] Vert. angular dimension of a rectangular area +* target_x: [m] Position of target to focus at. Transverse coordinate +* target_y: [m] Position of target to focus at. Vertical coordinate +* target_z: [m] Position of target to focus at. Straight ahead. +* target_index: [1] Relative index of component to focus at, e.g. next is +1 * verbose [1] Verbosity level * * CALCULATED PARAMETERS: From 90c8c688dcc803d338a3ad498066e1115d399df2 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 8 Apr 2026 15:41:53 +0200 Subject: [PATCH 07/14] Remove 'top-level' instr file --- Samples_SpinWave_BCO.instr | 248 ------------------------------------- 1 file changed, 248 deletions(-) delete mode 100644 Samples_SpinWave_BCO.instr diff --git a/Samples_SpinWave_BCO.instr b/Samples_SpinWave_BCO.instr deleted file mode 100644 index 8cddd4a72c..0000000000 --- a/Samples_SpinWave_BCO.instr +++ /dev/null @@ -1,248 +0,0 @@ -/***************************************************************************** -* McStas instrument definition URL=http://www.mcstas.org -* -* Instrument: Samples_SpinWave_BCO -* -* %Identification -* Written by: Kim Lefmann -* Date: Feb 2004 -* Origin: RISOE -* Written by: K. Lefmann RISOE, Feb 2004 -* %INSTRUMENT_SITE: Tests_samples -* -* Simple test instrument for the SpinWave_BCO component -* -* -* %Description -* Simple test instrument for the SpinWave_BCO component. -* Refer to the component documentation for further instructions. -* -* Example: h=1 l=0 -n 1e5 Detector: mon1_I=2.86265e-25 -* -* %Parameters -* E: [meV] Mean energy at source -* DE: [meV] Energy spread at source -* HDIV: [deg] Horizontal divergence produced from source -* VDIV: [deg] Vertical divergence produced from source -* TT: [deg] Two-theta detetector-angle -* OM: [deg] Sample rotation angle -* C: [meV/AA^(-1)] Sample velocity of sound -* -* %End -******************************************************************************/ -DEFINE INSTRUMENT Samples_SpinWave_BCO(E=1.06,Ef=14.7,Dlambda=0.1, h=1, l=0,B=0,dA3=-90, Temp=2, width = 0.005, E_steps_low = 50, E_steps_high = 50,mode = 2,verbose=0) - -DECLARE -%{ - double Gqx,Gqy,Gqz; -//double Ef=24.8 -double Ei; -//meV -double thetaM; -double twothetaS; -double thetaA; -double A3; -double QM; -double alpha; -double lambda_i; -double SMALL__NUMBER; -double a; -double c; -double jAFM; -double jFM; -double jextra; -double S; -double D; -%} - -INITIALIZE -%{ -// Set monochromator/analyzer Q-value for PG -QM = 1.8734; -SMALL__NUMBER = 1e-6; - -// MnF2 parameters -// Lattice vectors in AA -a = 4.873; -c = 3.130; - -jAFM = 2*0.152; -jFM = -2*0.028; -jextra = 2*0.004; -S = 2.5; -D = -0.023; - -//MnF2 reciprocal lattice vectors, in 1/AA -double astar = 2*PI/a; -double cstar = 2*PI/c; - -//calculate Ei -Ei=Ef+E; - -//calculate ki, kf, lambda_i, q -double ki = V2K*SE2V*sqrt(Ei); -double kf = V2K*SE2V*sqrt(Ef); -lambda_i=2*PI/ki; -double q = sqrt(h*h*astar*astar+l*l*cstar*cstar); - -//calculate 2thetaM and 2thetaA -thetaM = asin(QM/(2*ki))*RAD2DEG; -thetaA = asin(QM/(2*kf))*RAD2DEG; - -//calculate scattering angle and sample rotation -twothetaS = acos((q*q-ki*ki-kf*kf)/(-2*ki*kf))*RAD2DEG; -alpha = acos((kf*kf-ki*ki-q*q)/(-2*ki*q)); -A3=(asin(l*cstar/(q+SMALL__NUMBER))-alpha)*RAD2DEG; - -%} - -TRACE - -COMPONENT origin = Progress_bar() -AT (0, 0, 0) RELATIVE ABSOLUTE - -COMPONENT source = Source_Maxwell_3( - xwidth=0.01, - yheight=0.01, - Lmin=lambda_i-Dlambda/2, - Lmax=lambda_i+Dlambda/2, - dist=7.5, - focus_yh=width, - focus_xw=width, - T1=300, - T2=300, - T3=300, - I1=1E15, - I2=1E15, - I3=1E15) -AT (0, 0, 0) RELATIVE PREVIOUS - -COMPONENT monochromator_flat = Monochromator_flat( - mosaicv=30, - mosaich=30, - zwidth = width, - yheight = width, - Q=QM) -AT (0, 0, 7.5) RELATIVE source -ROTATED (0, thetaM, 0) RELATIVE source - -COMPONENT arm1 = Arm() -AT (0, 0, 0) RELATIVE PREVIOUS -ROTATED (0, thetaM, 0) RELATIVE PREVIOUS - -COMPONENT collimator_linear1 = Collimator_linear( - xwidth=0.1, - yheight=0.2, - length=0.2, - divergence=40) -AT (0, 0, 0.7) RELATIVE arm1 - -/*COMPONENT e_monitor_before_sample = E_monitor( - nE=100, - filename="E_before_sample", - Emin=Ei-2, - Emax=Ei+2, - xwidth=width, - yheight=2*width, - restore_neutron=1) -AT (0, 0, 0.99) RELATIVE arm1*/ - -COMPONENT arm2 = Arm() -AT (0, 0, 1) RELATIVE arm1 -ROTATED (0, 0, 0) RELATIVE arm1 - -SPLIT 1 - -COMPONENT SAMPLE = SpinWave_BCO( - radius = width/2, - yheight = 2*width, - sigma_abs = 0, - sigma_inc = 0, - T = Temp, - a1 = a, - a2 = a, - a3 = c, - j = jAFM, - jc = jFM, - ja=jextra, - jb=jextra, - S=S, - B=B, - D=D, - FM=0, - mode_input = mode, - e_steps_low = E_steps_low, - e_steps_high = E_steps_high, - target_index=3, - focus_xw=0.005, - focus_yh=0.005, - verbose = verbose) -AT (0, 0, 0) RELATIVE arm2 -ROTATED (0, -A3+dA3, 180) RELATIVE arm2 - -COMPONENT arm3 = Arm() -AT (0, 0, 0) RELATIVE arm2 -ROTATED (0, -twothetaS, 0) RELATIVE arm2 - -COMPONENT collimator_linear2 = Collimator_linear( - xwidth=0.1, - yheight=0.2, - length=0.2, - divergence=40, - xwidth=width, - yheight=width) -AT (0, 0, 0.5) RELATIVE arm3 - -COMPONENT slit = Slit( - xwidth=0.005, - yheight=0.005) -AT (0, 0, 1) RELATIVE arm3 - -/*COMPONENT e_monitor_beforeana = E_monitor( - nE=250, - filename="Ebeforeana", - xwidth=width, - yheight=width, - Emin=0, - Emax=25, - restore_neutron=1) -AT (0, 0, 0.001) RELATIVE PREVIOUS*/ - -COMPONENT analyzer = Monochromator_flat( - zwidth=width, - yheight=width, - mosaicv=30, - mosaich=30, - Q=QM - ) -AT (0, 0, 0.1) RELATIVE PREVIOUS -ROTATED (0, thetaA, 0) RELATIVE PREVIOUS - -COMPONENT arm4 = Arm() -AT (0, 0, 0) RELATIVE PREVIOUS -ROTATED (0, thetaA, 0) RELATIVE PREVIOUS - -COMPONENT collimator_linear3 = Collimator_linear( - xwidth=0.1, - yheight=0.2, - length=0.2, - divergence=40) -AT (0, 0, 0.5) RELATIVE PREVIOUS - -COMPONENT Emon = E_monitor( - nE=250, - filename="Emon", - xwidth=width, - yheight=width, - Emin=0, - Emax=75, - restore_neutron=1) -AT (0, 0, 1) RELATIVE arm4 - - -FINALLY -%{ -%} - -END - From 3b4b2381e2de72074bdb2e5d18cedd15d7e02fbb Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 8 Apr 2026 15:43:43 +0200 Subject: [PATCH 08/14] Correct author info --- .../Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr b/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr index 2800b0e9a5..7c9327551c 100644 --- a/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr +++ b/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr @@ -4,10 +4,9 @@ * Instrument: Samples_SpinWave_BCO * * %Identification -* Written by: Kim Lefmann -* Date: Feb 2004 -* Origin: RISOE -* Written by: K. Lefmann RISOE, Feb 2004 +* Written by: Silas Schack +* Date: Spring 2026 +* Origin: NBI * %INSTRUMENT_SITE: Tests_samples * * Simple test instrument for the SpinWave_BCO component From a9b87e138043026ae4d27715ae76204e8ee57258 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 8 Apr 2026 15:47:38 +0200 Subject: [PATCH 09/14] KU->NBI and case change Phonon_simple --- mcstas-comps/contrib/SpinWave_BCO.comp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mcstas-comps/contrib/SpinWave_BCO.comp b/mcstas-comps/contrib/SpinWave_BCO.comp index 942e8d4f18..f8b8fb8c6f 100644 --- a/mcstas-comps/contrib/SpinWave_BCO.comp +++ b/mcstas-comps/contrib/SpinWave_BCO.comp @@ -4,9 +4,9 @@ * Copyright(C) 2000 Risoe National Laboratory. * * %I -* Written by: Silas Schack (algoritme based on "phonon_simple" by Kim Lefmann) +* Written by: Silas Schack (algorithm based on "Phonon_simple" by Kim Lefmann) * Date: 13.06.25 -* Origin: KU +* Origin: NBI * * A sample for FM and AFM magnon scattering from a oI lattice based on cross section expressions from Marshall and Lovesey ch.9. * From fe6c88959bd2c56c18153e8e37b4abc39f98b277 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 8 Apr 2026 15:57:30 +0200 Subject: [PATCH 10/14] Put version from main back in place --- mcstas-comps/samples/Phonon_simple.comp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mcstas-comps/samples/Phonon_simple.comp b/mcstas-comps/samples/Phonon_simple.comp index e7f4e33266..1931584ced 100644 --- a/mcstas-comps/samples/Phonon_simple.comp +++ b/mcstas-comps/samples/Phonon_simple.comp @@ -510,4 +510,4 @@ MCDISPLAY line (0, -yheight / 2.0, +radius, 0, +yheight / 2.0, +radius); %} -END \ No newline at end of file +END From a1bb59055678980b8620daee8d579fc913bbe1e9 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 8 Apr 2026 15:59:34 +0200 Subject: [PATCH 11/14] Series of edits: * Remove commented monitor-instances * Remove -n specifiers in %Example lines * Add -1 html syntax in header --- .../Samples_SpinWave_BCO.instr | 28 +++---------------- 1 file changed, 4 insertions(+), 24 deletions(-) diff --git a/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr b/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr index 7c9327551c..ad1ff65b85 100644 --- a/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr +++ b/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr @@ -16,9 +16,9 @@ * Simple test instrument for the SpinWave_BCO component. * Refer to the component documentation for further instructions. * -* %Example: h=1 l=0 -n 1e6 Detector: Emon_I=0.0517575 -* %Example: h=0 l=0.5 E=6.73 -n 1e6 Detector: Emon_I=0.0108855 -* %Example: h=1 l=1 -n 1e6 Detector: Emon_I=0.000655744 +* %Example: h=1 l=0 Detector: Emon_I=0.0517575 +* %Example: h=0 l=0.5 E=6.73 Detector: Emon_I=0.0108855 +* %Example: h=1 l=1 Detector: Emon_I=0.000655744 * %Parameters * E: [meV] Mean energy at source @@ -27,7 +27,7 @@ * VDIV: [deg] Vertical divergence produced from source * TT: [deg] Two-theta detetector-angle * OM: [deg] Sample rotation angle -* C: [meV/AA^(-1)] Sample velocity of sound +* C: [meV/AA-1] Sample velocity of sound * * %End ******************************************************************************/ @@ -138,16 +138,6 @@ COMPONENT collimator_linear1 = Collimator_linear( divergence=40) AT (0, 0, 0.7) RELATIVE arm1 -/*COMPONENT e_monitor_before_sample = E_monitor( - nE=100, - filename="E_before_sample", - Emin=Ei-2, - Emax=Ei+2, - xwidth=width, - yheight=2*width, - restore_neutron=1) -AT (0, 0, 0.99) RELATIVE arm1*/ - COMPONENT arm2 = Arm() AT (0, 0, 1) RELATIVE arm1 ROTATED (0, 0, 0) RELATIVE arm1 @@ -199,16 +189,6 @@ COMPONENT slit = Slit( yheight=0.005) AT (0, 0, 1) RELATIVE arm3 -/*COMPONENT e_monitor_beforeana = E_monitor( - nE=250, - filename="Ebeforeana", - xwidth=width, - yheight=width, - Emin=0, - Emax=25, - restore_neutron=1) -AT (0, 0, 0.001) RELATIVE PREVIOUS*/ - COMPONENT analyzer = Monochromator_flat( zwidth=width, yheight=width, From 0ded5b4f19855be087d09b1602f2c7144fa8d5dd Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 8 Apr 2026 16:04:56 +0200 Subject: [PATCH 12/14] Clangformat run --- mcstas-comps/contrib/SpinWave_BCO.comp | 880 ++++++++++++------------- 1 file changed, 422 insertions(+), 458 deletions(-) diff --git a/mcstas-comps/contrib/SpinWave_BCO.comp b/mcstas-comps/contrib/SpinWave_BCO.comp index f8b8fb8c6f..53b223e6f6 100644 --- a/mcstas-comps/contrib/SpinWave_BCO.comp +++ b/mcstas-comps/contrib/SpinWave_BCO.comp @@ -77,315 +77,282 @@ target_x=0, target_y=0, target_z=0, int target_index=0,focus_r=0,focus_xw=0,focu /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE %{ -#ifndef MAGNON_oI -#define MAGNON_oI $Revision$ -#define T2E (1/11.605) /* Kelvin to meV */ -#define B2E 0.11590188615547234 /*Tesla to meV; g*Bohr magneton*/ + #ifndef MAGNON_oI + #define MAGNON_oI $Revision$ + #define T2E (1/11.605) /* Kelvin to meV */ + #define B2E 0.11590188615547234 /*Tesla to meV; g*Bohr magneton*/ -#pragma acc routine -double nbose(double omega, double T) /* Other name ?? */ + #pragma acc routine + double + nbose (double omega, double T) /* Other name ?? */ { double nb; - nb= (omega>0) ? 1+1/(exp(omega/(T*T2E))-1) : 1/(exp(-omega/(T*T2E))-1); + nb = (omega > 0) ? 1 + 1 / (exp (omega / (T * T2E)) - 1) : 1 / (exp (-omega / (T * T2E)) - 1); return nb; } -#undef T2E -/* Routine types from Numerical Recipies book */ -#define UNUSED (-1.11e30) -#define MAXRIDD 60 - -void fatalerror_cpu(char *s) - { - fprintf(stderr,"%s \n",s); - exit(1); + #undef T2E + /* Routine types from Numerical Recipies book */ + #define UNUSED (-1.11e30) + #define MAXRIDD 60 + + void + fatalerror_cpu (char* s) { + fprintf (stderr, "%s \n", s); + exit (1); } - -#pragma acc routine -void fatalerror(char *s) - { - #ifndef OPENACC - fatalerror_cpu(s); - #endif + + #pragma acc routine + void + fatalerror (char* s) { + #ifndef OPENACC + fatalerror_cpu (s); + #endif } #pragma acc routine - double omega_q(double* parms) - { - /* dispersion in units of meV */ - double vi, vf, vv_x, vv_y, vv_z, vi_x, vi_y, vi_z; - double qx, qy, qz, J, J1, J0, J10, omega_magnon, omega_neutron; - double coherence_factor; - double a1, a2, a3; - double S,D; - double j, ja, jb, jc; - double j_110, j_110_prime; - double B; - double coherence_flag; - int branch,FM_TAG; - - vf=parms[0]; - vi=parms[1]; - vv_x=parms[2]; - vv_y=parms[3]; - vv_z=parms[4]; - vi_x=parms[5]; - vi_y=parms[6]; - vi_z=parms[7]; - a1 = parms[8]; - a2 = parms[9]; - a3 = parms[10]; - j = parms[11]; - ja = parms[12]; - jb = parms[13]; - jc = parms[14]; - S = parms[15]; - B = parms[16]; - D = parms[17]; - branch = parms[18]; - coherence_flag = parms[19]; - FM_TAG = parms[20]; - j_110 = parms[21]; - j_110_prime = parms[22]; - - qx=V2K*(vi_x-vf*vv_x); - qy=V2K*(vi_y-vf*vv_y); - qz=V2K*(vi_z-vf*vv_z); - - J = 8*j*cos(a1*qx/2)*cos(a2*qy/2)*cos(a3*qz/2); - J0 = 8*j; - J1 = 2*(ja*cos(qx*a1)+jb*cos(qy*a2)+jc*cos(qz*a3))+2*(j_110+j_110_prime)*cos(qx*a1)*cos(qy*a2); - J10 = 2*(ja+jb+jc)+2*(j_110+j_110_prime); - - if (FM_TAG == 0) { - omega_magnon = sqrt((S*(J0-J10+J1)-2*D*(S-0.5))*(S*(J0-J10+J1)-2*D*(S-0.5))-S*S*(J)*(J)) - (2*branch - 1)*(B*B2E+2*(j_110-j_110_prime)*sin(qx*a1)*sin(qy*a2)); /*AM dispersion*/ - if ((omega_magnon<0) || ((S*(J0-J10+J1)-2*D*(S-0.5))*(S*(J0-J10+J1)-2*D*(S-0.5))-S*S*(J)*(J)<0)) { - fatalerror("Unphysical dispersion: Check parameters"); - } - } else { - omega_magnon = S*(J1+J-(J10+J0))-D*(S-0.5)+B2E*abs(B); /*FM dispersion*/ - if (omega_magnon<0) { - fatalerror("Unphysical dispersion: Check parameters"); + double + omega_q (double* parms) { + /* dispersion in units of meV */ + double vi, vf, vv_x, vv_y, vv_z, vi_x, vi_y, vi_z; + double qx, qy, qz, J, J1, J0, J10, omega_magnon, omega_neutron; + double coherence_factor; + double a1, a2, a3; + double S, D; + double j, ja, jb, jc; + double j_110, j_110_prime; + double B; + double coherence_flag; + int branch, FM_TAG; + + vf = parms[0]; + vi = parms[1]; + vv_x = parms[2]; + vv_y = parms[3]; + vv_z = parms[4]; + vi_x = parms[5]; + vi_y = parms[6]; + vi_z = parms[7]; + a1 = parms[8]; + a2 = parms[9]; + a3 = parms[10]; + j = parms[11]; + ja = parms[12]; + jb = parms[13]; + jc = parms[14]; + S = parms[15]; + B = parms[16]; + D = parms[17]; + branch = parms[18]; + coherence_flag = parms[19]; + FM_TAG = parms[20]; + j_110 = parms[21]; + j_110_prime = parms[22]; + + qx = V2K * (vi_x - vf * vv_x); + qy = V2K * (vi_y - vf * vv_y); + qz = V2K * (vi_z - vf * vv_z); + + J = 8 * j * cos (a1 * qx / 2) * cos (a2 * qy / 2) * cos (a3 * qz / 2); + J0 = 8 * j; + J1 = 2 * (ja * cos (qx * a1) + jb * cos (qy * a2) + jc * cos (qz * a3)) + 2 * (j_110 + j_110_prime) * cos (qx * a1) * cos (qy * a2); + J10 = 2 * (ja + jb + jc) + 2 * (j_110 + j_110_prime); + + if (FM_TAG == 0) { + omega_magnon = sqrt ((S * (J0 - J10 + J1) - 2 * D * (S - 0.5)) * (S * (J0 - J10 + J1) - 2 * D * (S - 0.5)) - S * S * (J) * (J)) + - (2 * branch - 1) * (B * B2E + 2 * (j_110 - j_110_prime) * sin (qx * a1) * sin (qy * a2)); /*AM dispersion*/ + if ((omega_magnon < 0) || ((S * (J0 - J10 + J1) - 2 * D * (S - 0.5)) * (S * (J0 - J10 + J1) - 2 * D * (S - 0.5)) - S * S * (J) * (J) < 0)) { + fatalerror ("Unphysical dispersion: Check parameters"); } + } else { + omega_magnon = S * (J1 + J - (J10 + J0)) - D * (S - 0.5) + B2E * abs (B); /*FM dispersion*/ + if (omega_magnon < 0) { + fatalerror ("Unphysical dispersion: Check parameters"); } - omega_neutron = fabs(VS2E*(vi*vi-vf*vf)); /*Magnitude of energy transfer*/ - if (coherence_flag == 0) { - return (omega_magnon - omega_neutron); - } else { /*calculate coherence factor*/ - if (FM_TAG == 0) { - coherence_factor = 2*S*((S*(J0-J)-S*(J10-J1)-2*D*(S-0.5)))/(omega_magnon + (2*branch - 1)*(B*B2E+2*(j_110-j_110_prime)*sin(qx*a1)*sin(qy*a2))); - - } else { - coherence_factor = 2*S; - } - if (coherence_factor < 0) { - fatalerror("Negative coherence factor: Check parameters"); /*Stop at negative coherence factor from unphysical parms*/ - } - return coherence_factor; - } } - - - -double zridd(double (*func)(double*), double x1, double x2, double *parms, double xacc) - { - int j; - double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; - - parms[0]=x1; - fl=(*func)(parms); - parms[0]=x2; - fh=(*func)(parms); - if (fl*fh >= 0) - { - if (fl==0) return x1; - if (fh==0) return x2; - return UNUSED; + omega_neutron = fabs (VS2E * (vi * vi - vf * vf)); /*Magnitude of energy transfer*/ + if (coherence_flag == 0) { + return (omega_magnon - omega_neutron); + } else { /*calculate coherence factor*/ + if (FM_TAG == 0) { + coherence_factor = 2 * S * ((S * (J0 - J) - S * (J10 - J1) - 2 * D * (S - 0.5))) + / (omega_magnon + (2 * branch - 1) * (B * B2E + 2 * (j_110 - j_110_prime) * sin (qx * a1) * sin (qy * a2))); + + } else { + coherence_factor = 2 * S; } - else - { - xl=x1; - xh=x2; - ans=UNUSED; - for (j=1; j= fh ? 1.0 : -1.0)*fm/s); - if (fabs(xnew-ans) <= xacc) - return ans; - ans=xnew; - parms[0]=ans; - fnew=(*func)(parms); - if (fnew == 0.0) return ans; - if (fabs(fm)*SIGN(fnew) != fm) - { - xl=xm; - fl=fm; - xh=ans; - fh=fnew; - } - else - if (fabs(fl)*SIGN(fnew) != fl) - { - xh=ans; - fh=fnew; - } - else - if(fabs(fh)*SIGN(fnew) != fh) - { - xl=ans; - fl=fnew; - } - else - fatalerror("never get here in zridd"); - if (fabs(xh-xl) <= xacc) - return ans; - } - fatalerror("zridd exceeded maximum iterations"); + if (coherence_factor < 0) { + fatalerror ("Negative coherence factor: Check parameters"); /*Stop at negative coherence factor from unphysical parms*/ } - return 0.0; /* Never get here */ + return coherence_factor; } + } -#pragma acc routine -double zridd_gpu(double x1, double x2, double *parms, double xacc) - { - int j; - double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; - - parms[0]=x1; - fl=omega_q(parms); - parms[0]=x2; - fh=omega_q(parms); - if (fl*fh >= 0) - { - if (fl==0) return x1; - if (fh==0) return x2; - return UNUSED; + double + zridd (double (*func) (double*), double x1, double x2, double* parms, double xacc) { + int j; + double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; + + parms[0] = x1; + fl = (*func) (parms); + parms[0] = x2; + fh = (*func) (parms); + if (fl * fh >= 0) { + if (fl == 0) + return x1; + if (fh == 0) + return x2; + return UNUSED; + } else { + xl = x1; + xh = x2; + ans = UNUSED; + for (j = 1; j < MAXRIDD; j++) { + xm = 0.5 * (xl + xh); + parms[0] = xm; + fm = (*func) (parms); + s = sqrt (fm * fm - fl * fh); + if (s == 0.0) + return ans; + xnew = xm + (xm - xl) * ((fl >= fh ? 1.0 : -1.0) * fm / s); + if (fabs (xnew - ans) <= xacc) + return ans; + ans = xnew; + parms[0] = ans; + fnew = (*func) (parms); + if (fnew == 0.0) + return ans; + if (fabs (fm) * SIGN (fnew) != fm) { + xl = xm; + fl = fm; + xh = ans; + fh = fnew; + } else if (fabs (fl) * SIGN (fnew) != fl) { + xh = ans; + fh = fnew; + } else if (fabs (fh) * SIGN (fnew) != fh) { + xl = ans; + fl = fnew; + } else + fatalerror ("never get here in zridd"); + if (fabs (xh - xl) <= xacc) + return ans; } - else - { - xl=x1; - xh=x2; - ans=UNUSED; - for (j=1; j= fh ? 1.0 : -1.0)*fm/s); - if (fabs(xnew-ans) <= xacc) - return ans; - ans=xnew; - parms[0]=ans; - fnew=omega_q(parms); - if (fnew == 0.0) return ans; - if (fabs(fm)*SIGN(fnew) != fm) - { - xl=xm; - fl=fm; - xh=ans; - fh=fnew; - } - else - if (fabs(fl)*SIGN(fnew) != fl) - { - xh=ans; - fh=fnew; - } - else - if(fabs(fh)*SIGN(fnew) != fh) - { - xl=ans; - fl=fnew; - } - else - fatalerror("never get here in zridd"); - if (fabs(xh-xl) <= xacc) - return ans; - } - fatalerror("zridd exceeded maximum iterations"); + fatalerror ("zridd exceeded maximum iterations"); + } + return 0.0; /* Never get here */ + } + + #pragma acc routine + double + zridd_gpu (double x1, double x2, double* parms, double xacc) { + int j; + double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; + + parms[0] = x1; + fl = omega_q (parms); + parms[0] = x2; + fh = omega_q (parms); + if (fl * fh >= 0) { + if (fl == 0) + return x1; + if (fh == 0) + return x2; + return UNUSED; + } else { + xl = x1; + xh = x2; + ans = UNUSED; + for (j = 1; j < MAXRIDD; j++) { + xm = 0.5 * (xl + xh); + parms[0] = xm; + fm = omega_q (parms); + s = sqrt (fm * fm - fl * fh); + if (s == 0.0) + return ans; + xnew = xm + (xm - xl) * ((fl >= fh ? 1.0 : -1.0) * fm / s); + if (fabs (xnew - ans) <= xacc) + return ans; + ans = xnew; + parms[0] = ans; + fnew = omega_q (parms); + if (fnew == 0.0) + return ans; + if (fabs (fm) * SIGN (fnew) != fm) { + xl = xm; + fl = fm; + xh = ans; + fh = fnew; + } else if (fabs (fl) * SIGN (fnew) != fl) { + xh = ans; + fh = fnew; + } else if (fabs (fh) * SIGN (fnew) != fh) { + xl = ans; + fl = fnew; + } else + fatalerror ("never get here in zridd"); + if (fabs (xh - xl) <= xacc) + return ans; } - return 0.0; /* Never get here */ + fatalerror ("zridd exceeded maximum iterations"); } + return 0.0; /* Never get here */ + } - -#define ROOTACC 1e-8 - - int findroots(double brack_low, double brack_mid, double brack_high, double *list, int* index, double (*f)(double*), double *parms, int steps_low, int steps_high) - { - double root,range; - int i; - range = brack_mid-brack_low; - for (i=0; i<=(steps_low-1); i++) - { - - root = zridd(f, brack_low+range*i/(int)steps_low, - brack_low+range*(i+1)/(int)steps_low, - (double *)parms, ROOTACC); - if (root != UNUSED) - { - list[(*index)++]=root; + #define ROOTACC 1e-8 + + int + findroots (double brack_low, double brack_mid, double brack_high, double* list, int* index, double (*f) (double*), double* parms, int steps_low, + int steps_high) { + double root, range; + int i; + range = brack_mid - brack_low; + for (i = 0; i <= (steps_low - 1); i++) { + + root = zridd (f, brack_low + range * i / (int)steps_low, brack_low + range * (i + 1) / (int)steps_low, (double*)parms, ROOTACC); + if (root != UNUSED) { + list[(*index)++] = root; } - } - - range = brack_high-brack_mid; - - for (i=0; i<=(steps_high-1); i++) - { - root = zridd(f, brack_mid+range*i/(int)steps_high, - brack_mid+range*(i+1)/(int)steps_high, - (double *)parms, ROOTACC); - if (root != UNUSED) - { - list[(*index)++]=root; + } + + range = brack_high - brack_mid; + + for (i = 0; i <= (steps_high - 1); i++) { + root = zridd (f, brack_mid + range * i / (int)steps_high, brack_mid + range * (i + 1) / (int)steps_high, (double*)parms, ROOTACC); + if (root != UNUSED) { + list[(*index)++] = root; } - } } - -#pragma acc routine - int findroots_gpu(double brack_low, double brack_mid, double brack_high, double *list, int* index, double *parms, int steps_low, int steps_high) - { - double root,range; - int i; - - range = brack_mid - brack_low; - - for (i=0; i<=(steps_low-1); i++) - { - root = zridd_gpu(brack_low+range*i/(int)steps_low, - brack_low+range*(i+1)/(int)steps_low, - (double *)parms, ROOTACC); - if (root != UNUSED) - { - list[(*index)++]=root; + } + + #pragma acc routine + int + findroots_gpu (double brack_low, double brack_mid, double brack_high, double* list, int* index, double* parms, int steps_low, int steps_high) { + double root, range; + int i; + + range = brack_mid - brack_low; + + for (i = 0; i <= (steps_low - 1); i++) { + root = zridd_gpu (brack_low + range * i / (int)steps_low, brack_low + range * (i + 1) / (int)steps_low, (double*)parms, ROOTACC); + if (root != UNUSED) { + list[(*index)++] = root; } - } - - range = brack_high - brack_mid; - - for (i=0; i<=(steps_high-1); i++) - { - root = zridd_gpu(brack_mid+range*i/(int)steps_high, - brack_mid+range*(i+1)/(int)steps_high, - (double *)parms, ROOTACC); - if (root != UNUSED) - { - list[(*index)++]=root; + } + + range = brack_high - brack_mid; + + for (i = 0; i <= (steps_high - 1); i++) { + root = zridd_gpu (brack_mid + range * i / (int)steps_high, brack_mid + range * (i + 1) / (int)steps_high, (double*)parms, ROOTACC); + if (root != UNUSED) { + list[(*index)++] = root; } - } } - -#undef UNUSED -#undef MAXRIDD -#endif + } + + #undef UNUSED + #undef MAXRIDD + #endif %} DECLARE @@ -400,13 +367,13 @@ DECLARE %} INITIALIZE %{ - g = 2.0023; /*electron g-factor*/ - gamma_n = -1.913; /*neutron gamma-factor*/ - r_0 = 2.8179; /*electron charge radius in fm*/ - V_0 = a1*a2*a3; /*inverse unit cell volume of sublattice (oP)*/ - V_my_s = (2 * 100 * sigma_inc/V_0); /*Incoherent volume specific cross section.*/ - V_my_a_v = (2 * 100 * sigma_abs/V_0 * 2200); - DV = 0.001; /* Velocity change used for numerical derivative */ + g = 2.0023; /*electron g-factor*/ + gamma_n = -1.913; /*neutron gamma-factor*/ + r_0 = 2.8179; /*electron charge radius in fm*/ + V_0 = a1 * a2 * a3; /*inverse unit cell volume of sublattice (oP)*/ + V_my_s = (2 * 100 * sigma_inc / V_0); /*Incoherent volume specific cross section.*/ + V_my_a_v = (2 * 100 * sigma_abs / V_0 * 2200); + DV = 0.001; /* Velocity change used for numerical derivative */ if (focus_aw) focus_aw *= DEG2RAD; @@ -414,213 +381,210 @@ INITIALIZE focus_ah *= DEG2RAD; /* now compute target coords if a component index is supplied */ - if (!target_index && !target_x && !target_y && !target_z) target_index=1; - if (target_index){ + if (!target_index && !target_x && !target_y && !target_z) + target_index = 1; + if (target_index) { Coords ToTarget; - ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index),POS_A_CURRENT_COMP); - ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget); - coords_get(ToTarget, &target_x, &target_y, &target_z); + ToTarget = coords_sub (POS_A_COMP_INDEX (INDEX_CURRENT_COMP + target_index), POS_A_CURRENT_COMP); + ToTarget = rot_apply (ROT_A_CURRENT_COMP, ToTarget); + coords_get (ToTarget, &target_x, &target_y, &target_z); } if (!(target_x || target_y || target_z)) { - printf("Magnon_oI: %s: The target is not defined. Using direct beam (Z-axis).\n", - NAME_CURRENT_COMP); - target_z=1; + printf ("Magnon_oI: %s: The target is not defined. Using direct beam (Z-axis).\n", NAME_CURRENT_COMP); + target_z = 1; } %} TRACE %{ - double t0, t1; /* Entry/exit time for cylinder */ - double v_i, v_f; /* Neutron velocities: initial, final */ - double vx_i, vy_i, vz_i; /* Neutron initial velocity vector */ - double dt0, dt; /* Flight times through sample */ - double l_full; /* Flight path length for non-scattered neutron */ - double l_i, l_o; /* Flight path lenght in/out for scattered neutron */ - double my_a_i; /* Initial attenuation factor */ - double my_a_f; /* Final attenuation factor */ - double solid_angle; /* Solid angle of target as seen from scattering point */ - double aim_x=0, aim_y=0, aim_z=1; /* Position of target relative to scattering point */ - double kappa_x, kappa_y, kappa_z; /* Scattering vector */ - double kappa2; /* Square of the scattering vector */ + double t0, t1; /* Entry/exit time for cylinder */ + double v_i, v_f; /* Neutron velocities: initial, final */ + double vx_i, vy_i, vz_i; /* Neutron initial velocity vector */ + double dt0, dt; /* Flight times through sample */ + double l_full; /* Flight path length for non-scattered neutron */ + double l_i, l_o; /* Flight path lenght in/out for scattered neutron */ + double my_a_i; /* Initial attenuation factor */ + double my_a_f; /* Final attenuation factor */ + double solid_angle; /* Solid angle of target as seen from scattering point */ + double aim_x = 0, aim_y = 0, aim_z = 1; /* Position of target relative to scattering point */ + double kappa_x, kappa_y, kappa_z; /* Scattering vector */ + double kappa2; /* Square of the scattering vector */ double kappa2_z_norm; - double bose_factor; /* Calculated value of the Bose factor */ - double omega; /* energy transfer */ - int nf, index; /* Number of allowed final velocities */ - double vf_list[20]; /* List of allowed final velocities */ - double J_factor; /* Jacobian from delta fnc.s in cross section */ - double coherence_factor; /* Coherence factor in cross section*/ - double f1,f2; /* probed values of omega_q minus omega */ - double w_atten,w_MC,w_magnon,w_Jacobi,w_coherence; /* temporary multipliers */ - int mode; /* index for mode selection */ - double E_max; /* Maximum upper bound of dispersion */ + double bose_factor; /* Calculated value of the Bose factor */ + double omega; /* energy transfer */ + int nf, index; /* Number of allowed final velocities */ + double vf_list[20]; /* List of allowed final velocities */ + double J_factor; /* Jacobian from delta fnc.s in cross section */ + double coherence_factor; /* Coherence factor in cross section*/ + double f1, f2; /* probed values of omega_q minus omega */ + double w_atten, w_MC, w_magnon, w_Jacobi, w_coherence; /* temporary multipliers */ + int mode; /* index for mode selection */ + double E_max; /* Maximum upper bound of dispersion */ double parms[23]; - - if(cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius, yheight)) - { - if(t0 < 0) + if (cylinder_intersect (&t0, &t1, x, y, z, vx, vy, vz, radius, yheight)) { + if (t0 < 0) ABSORB; /* Neutron came from the sample or begins inside */ /* Neutron enters at t=t0. */ - dt0 = t1-t0; /* Time in sample */ - v_i = sqrt(vx*vx + vy*vy + vz*vz); + dt0 = t1 - t0; /* Time in sample */ + v_i = sqrt (vx * vx + vy * vy + vz * vz); l_full = v_i * dt0; /* Length of path through sample if not scattered */ - dt = rand01()*dt0; /* Time of scattering (relative to t0) */ - l_i = v_i*dt; /* Penetration in sample at scattering */ - vx_i=vx; - vy_i=vy; - vz_i=vz; - PROP_DT(dt+t0); /* Point of scattering */ - - aim_x = target_x-x; /* Vector pointing at target (e.g. analyzer) */ - aim_y = target_y-y; - aim_z = target_z-z; - - if(focus_aw && focus_ah) { - randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle, - aim_x, aim_y, aim_z, focus_aw, focus_ah, ROT_A_CURRENT_COMP); - } else if(focus_xw && focus_yh) { - randvec_target_rect(&vx, &vy, &vz, &solid_angle, - aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP); + dt = rand01 () * dt0; /* Time of scattering (relative to t0) */ + l_i = v_i * dt; /* Penetration in sample at scattering */ + vx_i = vx; + vy_i = vy; + vz_i = vz; + PROP_DT (dt + t0); /* Point of scattering */ + + aim_x = target_x - x; /* Vector pointing at target (e.g. analyzer) */ + aim_y = target_y - y; + aim_z = target_z - z; + + if (focus_aw && focus_ah) { + randvec_target_rect_angular (&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_aw, focus_ah, ROT_A_CURRENT_COMP); + } else if (focus_xw && focus_yh) { + randvec_target_rect (&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP); } else { - randvec_target_circle(&vx,&vy,&vz,&solid_angle,aim_x,aim_y,aim_z, focus_r); + randvec_target_circle (&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_r); } - NORM(vx, vy, vz); - nf=0; + NORM (vx, vy, vz); + nf = 0; - if (FM==0) { - if (B2E*abs(B) > 2*sqrt(D*D*(S-0.5)*(S-0.5)-8*j*D*S*(S-0.5))) { - printf("Magnetic field too strong. Fieldstrength set to zero\n"); /*Error if B induced spin-flop transition*/ - B=0; - } - } else { - if (B<0) { - fatalerror("Field direction flippes spins: Using abs(B)"); + if (FM == 0) { + if (B2E * abs (B) > 2 * sqrt (D * D * (S - 0.5) * (S - 0.5) - 8 * j * D * S * (S - 0.5))) { + printf ("Magnetic field too strong. Fieldstrength set to zero\n"); /*Error if B induced spin-flop transition*/ + B = 0; + } + } else { + if (B < 0) { + fatalerror ("Field direction flippes spins: Using abs(B)"); + } } - } mode = 0; - if (FM==0) { - if (mode_input == 2) - mode = round(2*rand01()-0.5); /* Select mode randomly from 2 possibilities */ - else - mode = mode_input; - if ((mode < 0) || (mode > 1)) - { - printf("mode = %d ",mode); - fatalerror("illegal value of mode"); + if (FM == 0) { + if (mode_input == 2) + mode = round (2 * rand01 () - 0.5); /* Select mode randomly from 2 possibilities */ + else + mode = mode_input; + if ((mode < 0) || (mode > 1)) { + printf ("mode = %d ", mode); + fatalerror ("illegal value of mode"); + } } - } - parms[1]=v_i; - parms[2]=vx; - parms[3]=vy; - parms[4]=vz; - parms[5]=vx_i; - parms[6]=vy_i; - parms[7]=vz_i; - parms[8]=a1; - parms[9]=a2; - parms[10]=a3; - parms[11]=j; - parms[12]=ja; - parms[13]=jb; - parms[14]=jc; - parms[15]=S; - parms[16]=B; - parms[17]=D; - parms[18]=mode; - parms[19]=0; - parms[20]=FM; - parms[21]=j_110; - parms[22]=j_110_prime; - - if (FM==0) { - E_max = fabs(S*(8*fabs(j)+4*(fabs(ja)+fabs(jb)+fabs(jc))+2*(fabs(j_110)+fabs(j_110_prime)))-2*D*(S-0.5))+fabs(B2E*B)+2*fabs(j_110-j_110_prime); + parms[1] = v_i; + parms[2] = vx; + parms[3] = vy; + parms[4] = vz; + parms[5] = vx_i; + parms[6] = vy_i; + parms[7] = vz_i; + parms[8] = a1; + parms[9] = a2; + parms[10] = a3; + parms[11] = j; + parms[12] = ja; + parms[13] = jb; + parms[14] = jc; + parms[15] = S; + parms[16] = B; + parms[17] = D; + parms[18] = mode; + parms[19] = 0; + parms[20] = FM; + parms[21] = j_110; + parms[22] = j_110_prime; + + if (FM == 0) { + E_max = fabs (S * (8 * fabs (j) + 4 * (fabs (ja) + fabs (jb) + fabs (jc)) + 2 * (fabs (j_110) + fabs (j_110_prime))) - 2 * D * (S - 0.5)) + fabs (B2E * B) + + 2 * fabs (j_110 - j_110_prime); } else { - E_max = fabs(S*(2*(8*fabs(j)+2*(fabs(ja)+fabs(jb)+fabs(jc))))-D*(S-0.5))+fabs(B2E*B); + E_max = fabs (S * (2 * (8 * fabs (j) + 2 * (fabs (ja) + fabs (jb) + fabs (jc)))) - D * (S - 0.5)) + fabs (B2E * B); } #ifndef OPENACC - findroots(0, v_i, sqrt(v_i*v_i+1/VS2E*E_max)+20, vf_list, &nf, omega_q, parms, e_steps_low,e_steps_high); + findroots (0, v_i, sqrt (v_i * v_i + 1 / VS2E * E_max) + 20, vf_list, &nf, omega_q, parms, e_steps_low, e_steps_high); #else - findroots_gpu(0, v_i, sqrt(v_i*v_i+1/VS2E*E_max)+20, vf_list, &nf, parms, e_steps_low,e_steps_high); /*+10 to make sure final point is included*/ + findroots_gpu (0, v_i, sqrt (v_i * v_i + 1 / VS2E * E_max) + 20, vf_list, &nf, parms, e_steps_low, e_steps_high); /*+10 to make sure final point is included*/ #endif - index=(int)floor(rand01()*nf); - if (nf>0 && index<20) { + index = (int)floor (rand01 () * nf); + if (nf > 0 && index < 20) { v_f = vf_list[index]; parms[0] = v_f; - parms[19] = 1; /*coherence flag*/ - coherence_factor = omega_q(parms); + parms[19] = 1; /*coherence flag*/ + coherence_factor = omega_q (parms); parms[19] = 0; - if (verbose==1) { /* Print functions used for debugging */ - printf("vi=(%g,%g,%g)\n",vx_i,vy_i,vz_i); - printf("vf=(%g,%g,%g)\n",vx,vy,vz); - printf("nf=%d\n",nf); - printf("index=%d\n",index); - printf("omega_q found=%g\n vf=%g\n", omega_q(parms),v_f); + if (verbose == 1) { /* Print functions used for debugging */ + printf ("vi=(%g,%g,%g)\n", vx_i, vy_i, vz_i); + printf ("vf=(%g,%g,%g)\n", vx, vy, vz); + printf ("nf=%d\n", nf); + printf ("index=%d\n", index); + printf ("omega_q found=%g\n vf=%g\n", omega_q (parms), v_f); } - parms[0]=v_f-DV; - f1=omega_q(parms); - parms[0]=v_f+DV; - f2=omega_q(parms); - J_factor = fabs(f2-f1)/(2*DV); - omega=VS2E*(v_i*v_i-v_f*v_f); + parms[0] = v_f - DV; + f1 = omega_q (parms); + parms[0] = v_f + DV; + f2 = omega_q (parms); + J_factor = fabs (f2 - f1) / (2 * DV); + omega = VS2E * (v_i * v_i - v_f * v_f); vx *= v_f; vy *= v_f; vz *= v_f; - kappa_x=V2K*(vx_i-vx); - kappa_y=V2K*(vy_i-vy); - kappa_z=V2K*(vz_i-vz); + kappa_x = V2K * (vx_i - vx); + kappa_y = V2K * (vy_i - vy); + kappa_z = V2K * (vz_i - vz); - if (verbose==1) { - if (omega>0) { - printf("hkl=(%f,%f,%f)\n",kappa_x/(2*PI/a1),kappa_y/(2*PI/a2),kappa_z/(2*PI/a3)); - }} + if (verbose == 1) { + if (omega > 0) { + printf ("hkl=(%f,%f,%f)\n", kappa_x / (2 * PI / a1), kappa_y / (2 * PI / a2), kappa_z / (2 * PI / a3)); + } + } - kappa2=kappa_y*kappa_y+kappa_x*kappa_x+kappa_z*kappa_z; - kappa2_z_norm = kappa_z*kappa_z/kappa2; + kappa2 = kappa_y * kappa_y + kappa_x * kappa_x + kappa_z * kappa_z; + kappa2_z_norm = kappa_z * kappa_z / kappa2; - if(!cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius, yheight)) - { - /* ??? did not hit cylinder */ - printf("FATAL ERROR: Did not hit cylinder from inside.\n"); - exit(1); - } + if (!cylinder_intersect (&t0, &t1, x, y, z, vx, vy, vz, radius, yheight)) { + /* ??? did not hit cylinder */ + printf ("FATAL ERROR: Did not hit cylinder from inside.\n"); + exit (1); + } dt = t1; - l_o = v_f*dt; - - my_a_i = V_my_a_v/v_i; - my_a_f = V_my_a_v/v_f; - bose_factor=nbose(omega,T); - w_atten = exp(-(V_my_s*(l_i+l_o)+my_a_i*l_i+my_a_f*l_o)); /* Absorption factor */ - w_MC = nf*solid_angle*(l_full*1e10)/V_0; /* Focusing factors; assume random choice of n_f possibilities. Units = AA^-2 */ - - if (FM == 0) { - if (mode_input == 2) { - w_MC *= 2; /*n_m; number of branches/modes in MC choice*/ - } - } else { - w_MC *= 2; /*If FM==1, this accounts for v_0 only being half magnetic unit cell volume*/ + l_o = v_f * dt; + + my_a_i = V_my_a_v / v_i; + my_a_f = V_my_a_v / v_f; + bose_factor = nbose (omega, T); + w_atten = exp (-(V_my_s * (l_i + l_o) + my_a_i * l_i + my_a_f * l_o)); /* Absorption factor */ + w_MC = nf * solid_angle * (l_full * 1e10) / V_0; /* Focusing factors; assume random choice of n_f possibilities. Units = AA^-2 */ + + if (FM == 0) { + if (mode_input == 2) { + w_MC *= 2; /*n_m; number of branches/modes in MC choice*/ } + } else { + w_MC *= 2; /*If FM==1, this accounts for v_0 only being half magnetic unit cell volume*/ + } - w_magnon = (gamma_n*gamma_n*r_0*r_0/1e10)*(g*g*Fq*Fq/4)*(v_f/v_i)*DW*(1+kappa2_z_norm)*bose_factor/4; /*Constants and magnetic factors in cross section. Units = AA^2*/ - w_Jacobi = 2*VS2E*v_f/J_factor; /* Jacobian of delta functions in cross section. */ - w_coherence = coherence_factor; /* coherence factor */ - p *= w_atten*w_MC*w_magnon*w_Jacobi*w_coherence; - if (verbose ==2) { /* Print functions used for debugging absolute intensity measurements*/ - if (v_i > v_f) { - double delta_E = VS2E*(v_i*v_i-v_f*v_f); - printf("delta_E=%f\n",delta_E); - printf("hkl=(%f,%f,%f)\n",kappa_x/(2*PI/a1),kappa_y/(2*PI/a2),kappa_z/(2*PI/a3)); - printf("vvi=(%f,%f,%f)\n vvf=(%f,%f,%f)\n",vx_i/v_i,vy_i/v_i,vz_i/v_i,vx/v_f,vy/v_f,vz/v_f); - printf("Ei=%f. Ef=%f\n",VS2E*v_i*v_i,VS2E*v_f*v_f); - printf("J_factor=%f\n",J_factor); - printf("vi=%f, vf=%f\n w_Jacobi = %f\n w_coherence = %f\n w_magnon*10^10= %f\n", v_i,v_f,w_Jacobi,w_coherence,w_magnon*1e10); - printf("cross*10^10/N = %f\n",w_magnon*w_Jacobi*w_coherence*1e10); + w_magnon = (gamma_n * gamma_n * r_0 * r_0 / 1e10) * (g * g * Fq * Fq / 4) * (v_f / v_i) * DW * (1 + kappa2_z_norm) * bose_factor + / 4; /*Constants and magnetic factors in cross section. Units = AA^2*/ + w_Jacobi = 2 * VS2E * v_f / J_factor; /* Jacobian of delta functions in cross section. */ + w_coherence = coherence_factor; /* coherence factor */ + p *= w_atten * w_MC * w_magnon * w_Jacobi * w_coherence; + if (verbose == 2) { /* Print functions used for debugging absolute intensity measurements*/ + if (v_i > v_f) { + double delta_E = VS2E * (v_i * v_i - v_f * v_f); + printf ("delta_E=%f\n", delta_E); + printf ("hkl=(%f,%f,%f)\n", kappa_x / (2 * PI / a1), kappa_y / (2 * PI / a2), kappa_z / (2 * PI / a3)); + printf ("vvi=(%f,%f,%f)\n vvf=(%f,%f,%f)\n", vx_i / v_i, vy_i / v_i, vz_i / v_i, vx / v_f, vy / v_f, vz / v_f); + printf ("Ei=%f. Ef=%f\n", VS2E * v_i * v_i, VS2E * v_f * v_f); + printf ("J_factor=%f\n", J_factor); + printf ("vi=%f, vf=%f\n w_Jacobi = %f\n w_coherence = %f\n w_magnon*10^10= %f\n", v_i, v_f, w_Jacobi, w_coherence, w_magnon * 1e10); + printf ("cross*10^10/N = %f\n", w_magnon * w_Jacobi * w_coherence * 1e10); } } @@ -632,13 +596,13 @@ TRACE MCDISPLAY %{ - - 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); %} END From ce90938ad5d39e7913ae9084ded8505fb2e5a8ec Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 8 Apr 2026 16:20:20 +0200 Subject: [PATCH 13/14] Add back missing * in header (should rectify instr table) --- .../Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr b/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr index ad1ff65b85..4362b11a05 100644 --- a/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr +++ b/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr @@ -19,7 +19,7 @@ * %Example: h=1 l=0 Detector: Emon_I=0.0517575 * %Example: h=0 l=0.5 E=6.73 Detector: Emon_I=0.0108855 * %Example: h=1 l=1 Detector: Emon_I=0.000655744 - +* * %Parameters * E: [meV] Mean energy at source * DE: [meV] Energy spread at source From fb4c36787e4012abff58c24620ab09350db300e0 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 8 Apr 2026 16:31:21 +0200 Subject: [PATCH 14/14] Doc header edits --- .../Samples_SpinWave_BCO.instr | 22 ++++++++++++------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr b/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr index 4362b11a05..6b90afebcc 100644 --- a/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr +++ b/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr @@ -21,17 +21,23 @@ * %Example: h=1 l=1 Detector: Emon_I=0.000655744 * * %Parameters -* E: [meV] Mean energy at source -* DE: [meV] Energy spread at source -* HDIV: [deg] Horizontal divergence produced from source -* VDIV: [deg] Vertical divergence produced from source -* TT: [deg] Two-theta detetector-angle -* OM: [deg] Sample rotation angle -* C: [meV/AA-1] Sample velocity of sound +* E: [meV] Mean energy at source +* Ef: [meV] Final energy +* Dlambda: [AA] Wavelength spread at source +* h: [1] 1st sample Miller index +* l: [1] 3rd sample Miller index +* B: [T] Field strength of external magnetic field applied along z-axis +* dA3: [deg] Offset from A3 nominal value +* Temp: [K] Sample temperature +* width: [deg] Width of sample +* E_steps_low: [ ] -> PLEASE FILL IN <- +* E_steps_high: [ ] -> PLEASE FILL IN <- +* mode: [1] -> PLEASE FILL IN <- +* verbose: [1] Verbosity flag * * %End ******************************************************************************/ -DEFINE INSTRUMENT Samples_SpinWave_BCO(E=1.06,Ef=14.7,Dlambda=0.1, h=1, l=0,B=0,dA3=-90, Temp=2, width = 0.005, E_steps_low = 50, E_steps_high = 50,mode = 2,verbose=0) +DEFINE INSTRUMENT Samples_SpinWave_BCO(E=1.06,Ef=14.7,Dlambda=0.1, h=1, l=0,B=0,dA3=-90, Temp=2, width = 0.005, E_steps_low = 50, E_steps_high = 50, int mode = 2,verbose=0) DECLARE %{