/******************************************************************************* * * 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 of magnetic ions * 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] Magnetic form factor * 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, Jq, J1q, 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, verbose; 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]; verbose = parms[23]; qx = V2K * (vi_x - vf * vv_x); qy = V2K * (vi_y - vf * vv_y); qz = V2K * (vi_z - vf * vv_z); Jq = 8 * j * cos (a1 * qx / 2) * cos (a2 * qy / 2) * cos (a3 * qz / 2); J0 = 8 * j; J1q = 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 + J1q) - 2 * D * (S - 0.5)) * (S * (J0 - J10 + J1q) - 2 * D * (S - 0.5)) - S * S * (Jq) * (Jq)) - (2 * branch - 1) * (B * B2E + 2 * S * (j_110 - j_110_prime) * sin (qx * a1) * sin (qy * a2)); /*AM dispersion*/ if ((omega_magnon < 0) || ((S * (J0 - J10 + J1q) - 2 * D * (S - 0.5)) * (S * (J0 - J10 + J1q) - 2 * D * (S - 0.5)) - S * S * (Jq) * (Jq) < 0)) { fatalerror ("Unphysical dispersion: Check parameters"); } } else { omega_magnon = S * (J1q + Jq - (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) { if ((verbose == 2) && fabs (omega_magnon - omega_neutron) < 1e-3 && (vi > vf)) { printf ("omega_q called with parameters vf= %g, vi=%g (%g %g %g) vv=(%g, %g, %g) q=(%g %g %g)\n", vf, vi, vi_x, vi_y, vi_z, vv_x, vv_y, vv_z, qx, qy, qz); printf ("omega_q gives: J0 = %g , Jq = %g, J10 = %g, J1q = %g, D = %g\n", J0, Jq, J10, J1q, D); printf ("in omega_q: q=(%g %g %g) omega_magnon=%g, omega_neutron=%g\n", qx, qy, qz, omega_magnon, omega_neutron); } return (omega_magnon - omega_neutron); } else { /*calculate coherence factor*/ if (FM_TAG == 0) { coherence_factor = 2 * S * ((S * (J0 - Jq) - S * (J10 - J1q) - 2 * D * (S - 0.5))) / (omega_magnon + (2 * branch - 1) * (B * B2E + 2 * S * (j_110 - j_110_prime) * sin (qx * a1) * sin (qy * a2))); if (verbose = 3) { printf ("Coherence factor = %g\n", coherence_factor); } } 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[24]; 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; parms[23] = verbose; 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 * S * 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; 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); 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; SCATTER; if (verbose == 1) { printf ("w factors : %g %g %g %g %g Omega: %g \n", w_atten, w_MC, w_magnon, w_Jacobi, w_coherence, omega); printf ("J_factor %g l_full %g, v_f/v_i %g, DW %g, kappa2 %g, bose_factor%g, fabs(omega) %g, coherence %g \n", J_factor, l_full, v_f / v_i, DW, kappa2, bose_factor, fabs (omega), coherence_factor); } } else { ABSORB; // findroots returned junk // Simply absorb if we can not hit, no fatal errors. (Typically indication of v_f==0) } } /* 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