diff --git a/mcstas-comps/contrib/SpinWave_BCO.comp b/mcstas-comps/contrib/SpinWave_BCO.comp new file mode 100644 index 000000000..53b223e6f --- /dev/null +++ b/mcstas-comps/contrib/SpinWave_BCO.comp @@ -0,0 +1,608 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2000 Risoe National Laboratory. +* +* %I +* Written by: Silas Schack (algorithm based on "Phonon_simple" by Kim Lefmann) +* Date: 13.06.25 +* 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. +* +* %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 <-- 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 +* 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. +* 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 < 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; + } + 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; + } + 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 */ + + 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) { + 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, 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); + } + + 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 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 000000000..6b90afebc --- /dev/null +++ b/mcstas-comps/examples/Tests_samples/Samples_SpinWave_BCO/Samples_SpinWave_BCO.instr @@ -0,0 +1,235 @@ +/***************************************************************************** +* McStas instrument definition URL=http://www.mcstas.org +* +* Instrument: Samples_SpinWave_BCO +* +* %Identification +* Written by: Silas Schack +* Date: Spring 2026 +* Origin: NBI +* %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 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 +* 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, int 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 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 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/samples/Phonon_simple.comp b/mcstas-comps/samples/Phonon_simple.comp index bf45c8c98..1931584ce 100644 --- a/mcstas-comps/samples/Phonon_simple.comp +++ b/mcstas-comps/samples/Phonon_simple.comp @@ -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;