diff --git a/mcstas-comps/contrib/Monochromator_bent.comp b/mcstas-comps/contrib/Monochromator_bent.comp new file mode 100644 index 0000000000..6e84ad8289 --- /dev/null +++ b/mcstas-comps/contrib/Monochromator_bent.comp @@ -0,0 +1,905 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright 1997-2002, All rights reserved +* Risoe National Laboratory, Roskilde, Denmark +* Institut Laue Langevin, Grenoble, France +* +* Component: Monochromator +* +* %I +* Written by: Daniel Lomholt Christensen +* Based on the model implemented by Jan Šaroun in the paper published in +* Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment +* Volume 529, Issues 1 through 3, 21 August 2004, Pages 162-165 +* +* +* Date: 24 August 2023 +* Origin: ILL / Niels Bohr Institute, University of Copenhagen. +* +* A perfect bent crystal monochromator. +* +* %D +* This monochromator is a number of lamella of a perfect crystal. +* The lamella are placed in the horizontal plane, behind each other, centered at 0. +* The crystal is bent, so that it follows a curve on a cylinder of radius_x. +* The monochromator lies along the z plane, so when a diffraction angle of theta +* is desired, it should just be inserted in the ROTATED parameter around +* the y-axis. +* Instruments that showcase the use of this component is the +* "Test_monochromator_bent.instr", and the "ILL_SALSA.instr" under the examples folder. +* SALSA showcases its complex use in a real instrument, while Test_monochromator_bent +* makes a simple show of its capabilities. +* +* +* %P +* INPUT PARAMETERS: +* +* zwidth: [m] Width of each lamella without bending. +* yheight: [m] Height of each lamella without bending. +* xthickness: [m] Thickness of each lamella without bending. +* radius_x: [m] Radius of the circle the monochromator bends on in the plane. +* lamella_slabs = 1 Amount of horizontal lamella in you monochromator. +* lamella_gap_size = 0 Gap between said horizontal lamella. +* plane_of_reflection: ["Si400"] The plane of reflection from the material. The list of possible reflections can +* can be seen in the source code. +* angle_to_cut_horizontal [degrees] Angle between cut and normal of crystal slab, horizontally +* angle_to_cut_vertical [degrees] Angle between cut and normal of crystal slab, vertically +* domainthickness [micro meter] Thickness of the crystal domains. +* temperature: [K] Temperature of the monochromator in Kelvin. +* verbose: [0] Verbosity of the monochromator. Used for debugging. +* +* %E +*******************************************************************************/ +DEFINE COMPONENT Monochromator_bent +DEFINITION PARAMETERS () +SETTING PARAMETERS (zwidth = 0.2, + yheight = 0.1, + xthickness = 0.0005, + radius_x = 2, + lamella_slabs = 1, + lamella_gap_size = 0, + string plane_of_reflection = "Si400", + angle_to_cut_horizontal = 0, + angle_to_cut_vertical = 0, + domainthickness = 10, + temperature = 300, + int verbose = 0) +OUTPUT PARAMETERS () +/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ + +SHARE +%{ + %include "read_table-lib" + #include + /////////////////////////////////////////////////////////////////////////// + /////////////// Structs and simple functions specific for this component + /////////////////////////////////////////////////////////////////////////// + + struct Monochromator_values{ + double length, height, thickness; + int type; + double radius_horizontal; + double Debye_Waller_factor; + double lattice_spacing; + double Maier_Leibnitz_reflectivity; + double poisson_ratio; + double bound_atom_scattering_cross_section; + double absorption_for_1AA_Neutrons; + double incoherent_scattering_cross_section; + double volume; + double Constant_from_Freund_paper; + double debye_temperature; + double atomic_number; + double temperature_mono; + double B0; + double BT; + double single_phonon_absorption; + double multiple_phonon_absorption; + double nuclear_capture_absorption; + double total_absorption; + double tau[3]; + double perp_to_tau[3]; + double lattice_spacing_gradient_field[3][3]; + double gradient_of_bragg_angle; + int lamellas; + double gap_size; + double domain_thickness; + double max_angle; + double min_angle; + }; + + struct neutron_values { + /* Statically allocate vectors that are always 3 */ + double* ki; // Incoming wavevector + double* kf; // outgoig wavevector + double* r; + double* v; // velocity of neutron + double* tau; //Reciprocal lattice vector + double ki_size; // size of incoming wavevector + double v_size; // speed + double tau_size; // size of reciprocal lattice vector + double kf_size; // size of outgoing wavevector + double* Bragg_angle_gradient; // Deviation from bragg condition + double absorption; // Absorption factor + double path; // Length of the path the neutron follows + double wavelength; // De Broglie wavelength of neutron + double kinematic_reflectivity; // The Q value from the paper this code is based on. + double* path_length; /* The time spent in crystals, to add to path for attenuation */ + double* entry_time; /* Time from start of crystal, to entrance of each lamella */ + double* exit_time; /* Time from start of crystal, to exit of each lamella */ + double* probabilities; /* Probability of reflection in each lamella */ + double* accumulating_probabilities; /* Accumulating probability in each lamella */ + double* time_of_reflection; /* The time in s from crystal edge to reflection */ + }; + + enum crystal_type {flat, bent, mosaic, bent_mosaic}; + /////////////// Define sign function + + double sign(double x){ + if (x > 0) return 1; + if (x < 0) return -1; + } + + double square(double x){ + return x*x; + } + enum crystal_plane {Cu111, Cu200, Cu220, Cu311, Cu400, Cu331, Cu420, Cu440, Ge111, Ge220, Ge311, + Ge400, Ge331, Ge422, Ge511, Ge533, Ge711, Ge551, Si111, Si220, Si311, Si400, Si331, + Si422, Si333, Si511, Si440, Si711, Si551, Be10, Be100, Be102, Be103, Be110, Be112, Be200, + Be00_2, Be10_1, PG00_2,PG00_4,PG00_6, Fe110, HS111,HS222,HS111star,Di111,Di220, Di311, Di400, + Di331, Di422, Di333, Di511, Di440}; + + /* An array containing all the possible strings that will be accepted if given as an + argument to the parameter plane_of_reflection */ + const char* crystal_planeStrings[] = { + "Cu111", "Cu200", "Cu220", "Cu311", "Cu400", "Cu331", "Cu420", "Cu440", "Ge111", "Ge220", "Ge311", + "Ge400", "Ge331", "Ge422", "Ge511", "Ge533", "Ge711", "Ge551", "Si111", "Si220", "Si311", "Si400", "Si331", + "Si422", "Si333", "Si511", "Si440", "Si711", "Si551"," Be10", "Be100", "Be102", "Be103", "Be110", "Be112", "Be200", + "Be00_2", "Be10_1", "PG00_2","PG00_4","PG00_6", "Fe110", "HS111","HS222","HS111star","Di111","Di220", "Di311", "Di400", + "Di331", "Di422", "Di333", "Di511", "Di440"}; + + // Function to convert a string to an enum value + enum crystal_plane stringToEnum(const char* plane) { + for (int i = 0; i < sizeof(crystal_planeStrings) / sizeof(crystal_planeStrings[0]); ++i) { + if (strcmp(plane, crystal_planeStrings[i]) == 0) { + return (enum crystal_plane)i; + } + } + } + /* TITLE Crystal table for perfect crystal bent monochromator + Table copied from SIMRES, current url: https://github.com/saroun/simres + Contents: dhkl, QML,sigmab,sigmaa,V0,A,thetaD,C2,poi + dhkl ... Lattice spacing of crystal plane. + QML = 4*PI*(F*dhkl/V0)**2 [ A^-1 cm^-1] + sigmab ... bound-atom scattering cross-section [barn] + sigmaa ... absorption for 1A neutrons [barn*A^-1] + sigmai ... incoherent scattering cross-section [barn] + V0 .... volume [A^3]/atom + A .... atomic number + thetaD .... Debye temperature (K) + C2 .... constant from the Freund's paper [A^-2 eV^-1] + poi .... Poisson elastic constant */ + + + double crystal_table[56][10] = {{ 2.087063, 0.23391E+00 ,7.485, 2.094, 0.55, 11.81, 63.54, 315, 12.00, 0.30000E+00}, + { 1.80745 , 0.17544E+00 ,7.485, 2.094, 0.55, 11.81, 63.54, 315, 12.00, 0.30000E+00}, + { 1.27806 , 0.87718E-01 ,7.485, 2.094, 0.55, 11.81, 63.54, 315, 12.00, 0.30000E+00}, + { 1.089933, 0.63795E-01 ,7.485, 2.094, 0.55, 11.81, 63.54, 315, 12.00, 0.30000E+00}, + { 0.903725, 0.43859E-01 ,7.485, 2.094, 0.55, 11.81, 63.54, 315, 12.00, 0.30000E+00}, + { 0.829315, 0.36934E-01 ,7.485, 2.094, 0.55, 11.81, 63.54, 315, 12.00, 0.30000E+00}, + { 0.808316, 0.35087E-01 ,7.485, 2.094, 0.55, 11.81, 63.54, 315, 12.00, 0.30000E+00}, + { 0.63903 , 0.21930E-01 ,7.485, 2.094, 0.55, 11.81, 63.54, 315, 12.00, 0.30000E+00}, + { 3.26665 , 0.87700E-01 ,8.42 , 1.216, 0.18, 22.63, 72.6, 290, 9.0, 0.15450E+00}, + { 2.00041 , 0.65760E-01 ,8.42 , 1.216, 0.18, 22.63, 72.6, 290, 9.0, 0.30000E+00}, + { 1.70595 , 0.23920E-01 ,8.42 , 1.216, 0.18, 22.63, 72.6, 290, 9.0, 0.15430E+00}, + { 1.41450 , 0.32880E-01 ,8.42 , 1.216, 0.18, 22.63, 72.6, 290, 9.0, 0.27300E+00}, + { 1.29803 , 0.13850E-01 ,8.42 , 1.216, 0.18, 22.63, 72.6, 290, 9.0, 0.15430E+00}, + { 1.15493 , 0.21925E-01 ,8.42 , 1.216, 0.18, 22.63, 72.6, 290, 9.0, 0.27270E+00}, + { 1.08888 , 0.97400E-02 ,8.42 , 1.216, 0.18, 22.63, 72.6, 290, 9.0, 0.27270E+00}, + { 0.86284 , 0.61200E-02 ,8.42 , 1.216, 0.18, 22.63, 72.6, 290, 9.0, 0.27270E+00}, + { 0.79228 , 0.51588E-02 ,8.42 , 1.216, 0.18, 22.63, 72.6, 290, 9.0, 0.27270E+00}, + { 0.79228 , 0.51600E-02 ,8.42 , 1.216, 0.18, 22.63, 72.6, 290, 9.0, 0.27270E+00}, + { 3.13536 , 0.25970E-01 ,2.18 , 0.0889, 0.0, 20.02, 28.09, 420, 6.36, 0.18080E+00}, + { 1.92001 , 0.19480E-01 ,2.18 , 0.0889, 0.0, 20.02, 28.09, 420, 6.36, 0.30000E+00}, + { 1.63739 , 0.70800E-02 ,2.18 , 0.0889, 0.0, 20.02, 28.09, 420, 6.36, 0.28000E+00}, + { 1.35765 , 0.97400E-02 ,2.18 , 0.0889, 0.0, 20.02, 28.09, 420, 6.36, 0.28000E+00}, + { 1.24587 , 0.41000E-02 ,2.18 , 0.0889, 0.0, 20.02, 28.09, 420, 6.36, 0.18080E+00}, + { 1.10852 , 0.64930E-02 ,2.18 , 0.0889, 0.0, 20.02, 28.09, 420, 6.36, 0.28000E+00}, + { 1.04512 , 0.28900E-02 ,2.18 , 0.0889, 0.0, 20.02, 28.09, 420, 6.36, 0.28000E+00}, + { 1.04512 , 0.28900E-02 ,2.18 , 0.0889, 0.0, 20.02, 28.09, 420, 6.36, 0.28000E+00}, + { 0.96000 , 0.48700E-02 ,2.18 , 0.0889, 0.0, 20.02, 28.09, 420, 6.36, 0.28000E+00}, + { 0.76044 , 0.15277E-02 ,2.18 , 0.0889, 0.0, 20.02, 28.09, 420, 6.36, 0.28000E+00}, + { 0.76044 , 0.15277E-02 ,2.18 , 0.0889, 0.0, 20.02, 28.09, 420, 6.36, 0.28000E+00}, + { 1.97956 , 0.11361 ,7.62579, 0.00422655, 0.002, 8.10926, 9.012, 1100, 7.62, 0.30000E+00}, + { 1.97956 , 0.11361 ,7.62579, 0.00422655, 0.002, 8.10926, 9.012, 1100, 7.62, 0.28000E+00}, + { 1.32857 , 0.05117 ,7.62579, 0.00422655, 0.002, 8.10926, 9.012, 1100, 7.62, 0.28000E+00}, + { 1.02290 , 0.091 ,7.62579, 0.00422655, 0.002, 8.10926, 9.012, 1100, 7.62, 0.28000E+00}, + { 1.14290 , 0.15147 ,7.62579, 0.00422655, 0.002, 8.10926, 9.012, 1100, 7.62, 0.28000E+00}, + { 0.96363 , 0.10768 ,7.62579, 0.00422655, 0.002, 8.10926, 9.012, 1100, 7.62, 0.28000E+00}, + { 0.98978 , 0.0284 ,7.62579, 0.00422655, 0.002, 8.10926, 9.012, 1100, 7.62, 0.28000E+00}, + { 1.79215 , 0.37245 ,7.62579, 0.00422655, 0.002, 8.10926, 9.012, 1100, 7.62, 0.30000E+00}, + { 1.73285 , 0.26116 ,7.62579, 0.00422655, 0.002, 8.10926, 9.012, 1100, 7.62, 0.30000E+00}, + { 3.35500 , 0.79500E+00 ,5.555, 0.0019, 0.0, 8.80, 12.01, 1050, 20.00, 0.30000E+00}, + { 1.67750 , 0.18000E+00 ,5.555, 0.0019, 0.0, 8.80, 12.01, 1050, 20.00, 0.30000E+00}, + { 1.11830 , 0.08833E+00 ,5.555, 0.0019, 0.0, 8.80, 12.01, 1050, 20.00, 0.30000E+00}, + { 2.02660 , 0.34031E+00 ,11.43, 2.53, 0.4 , 11.75 , 55.85, 411, 10.67 , 0.30000E+00}, + { 3.43500 , 0.11020E+00 ,1.79, 2.88, 0.55, 13.16, 48.0, 300, 12.00 , 0.30000E+00}, + { 1.71750 , 0.13130E+00 ,1.79, 2.88, 0.55, 13.16, 48.0, 300, 12.00 , 0.30000E+00}, + { 3.43500 , 0.55100E-01 ,1.79, 2.88, 0.55, 13.16, 48.0, 300, 12.00 , 0.30000E+00}, + { 2.05929 , 0.36606 ,5.55449 ,0.00194444, 0.0, 5.67213, 12.01, 1860, 3.00, 0.30000E+00}, + { 1.26105 , 0.27455 ,5.55449 ,0.00194444, 0.0, 5.67213, 12.01, 1860, 3.00, 0.30000E+00}, + { 1.07543 , 0.09984 ,5.55449 ,0.00194444, 0.0, 5.67213, 12.01, 1860, 3.00, 0.30000E+00}, + { 0.89170 , 0.13727 ,5.55449 ,0.00194444, 0.0, 5.67213, 12.01, 1860, 3.00, 0.30000E+00}, + { 0.81828 , 0.0578 ,5.55449 ,0.00194444, 0.0, 5.67213, 12.01, 1860, 3.00, 0.30000E+00}, + { 0.72807 , 0.09152 ,5.55449 ,0.00194444, 0.0, 5.67213, 12.01, 1860, 3.00, 0.30000E+00}, + { 0.68643 , 0.04067 ,5.55449 ,0.00194444, 0.0, 5.67213, 12.01, 1860, 3.00, 0.30000E+00}, + { 0.68643 , 0.04067 ,5.55449 ,0.00194444, 0.0, 5.67213, 12.01, 1860, 3.00, 0.30000E+00}, + { 0.63053 , 0.06864 ,5.55449 ,0.00194444, 0.0, 5.67213, 12.01, 1860, 3.00, 0.30000E+00}, + { 0.63053 , 0.06864 ,5.55449 ,0.00194444, 0.0, 5.67213, 12.01, 1860, 3.00, 0.30000E+00} + }; + + /////////////////////////////////////////////////////////////////////////// + /////////////// Testing function + /////////////////////////////////////////////////////////////////////////// + void print_neutron_state(struct neutron_values* neutron){ + printf("Neutron state:\nki %g, %g, %g\ntau %g, %g, %g\nkf %g, %g, %g\nv %g, %g, %g\nr %g, %g, %g\nki size %g, tau size %g, kf size %g, v size %g\n\n", + neutron->ki[0], neutron->ki[1], neutron->ki[2], + neutron->tau[0], neutron->tau[1], neutron->tau[2], + neutron->kf[0], neutron->kf[1], neutron->kf[2], + neutron->v[0], neutron->v[1], neutron->v[2], + neutron->r[0], neutron->r[1], neutron->r[2], + neutron->ki_size, neutron->tau_size, neutron->kf_size, neutron->v_size + ); + } + /////////////////////////////////////////////////////////////////////////// + /////////////// Calculations for absorption factor + /////////////// Based on the cross sections from + /////////////// A. K. Freund in Nuclear Instruments and Methods 213 (1983) 495-501 + /////////////////////////////////////////////////////////////////////////// + + // Integral needed for debye factor + + double calculate_phi_integral(double x){ + // Asymptotic approximation + if (x > 5) return PI * PI / 6 - exp(-x)/(x+1); + // Integate with Simpson/3. I dont know what this means + double z = 1 + x/(exp(x)-1); + double dx = x/100; + double ksi; + for (int i = 2; i <= 100; i++) { + ksi = (i-1)*dx; + switch (i%2){ + case 1: + z = z + 4 * ksi/(exp(ksi)-1); + break; + case 0: + z = z + 2 * ksi/(exp(ksi)-1); + break; + } + } + return z*dx/3; + } + + + int neutron_is_inside_monochromator(struct Monochromator_values* monochromator, + struct neutron_values* neutron, int lamella){ + double t1, t2; + if (neutron->time_of_reflection[lamella] <= 0) return 0; + /* Check that the neutron is not beyond */ + cylinder_intersect(&t1, &t2, neutron->r[0]-monochromator->radius_horizontal + monochromator->thickness/2, + neutron->r[1],neutron->r[2],neutron->v[0],neutron->v[1],neutron->v[2], + monochromator->radius_horizontal, monochromator->height); + if (t2 <= 0) return 0; + /* Check that the neutron is with the angle ranges*/ + double angle_on_inner_cylinder = PI - asin(neutron->r[2]/monochromator->radius_horizontal); + if (monochromator-> max_angle + 0.0001 <= angle_on_inner_cylinder || angle_on_inner_cylinder <= monochromator->min_angle - 0.0001){ + return 0; + } + return 1; + } + + /* B0 and BT are values used for the Debye factor */ + void calculate_B0_and_BT(struct Monochromator_values *monochromator){ + double x; + monochromator->B0 = 2872.556/monochromator->atomic_number + /monochromator->debye_temperature; + + if (monochromator->temperature_mono>0.1) x = monochromator->debye_temperature/monochromator->temperature_mono; + else x =monochromator->debye_temperature/0.1; + double phi = calculate_phi_integral(x); + + monochromator->BT = 4 * monochromator->B0 * phi / square(x); + } + + double calculate_kinematic_reflectivity(struct Monochromator_values* monochromator, + struct neutron_values* neutron){ + double sine_of_bragg_angle = neutron->wavelength/2/monochromator->lattice_spacing; + double cosine_of_bragg_angle = sqrt(1-square(sine_of_bragg_angle)); + double extinction_length = monochromator->lattice_spacing/neutron->wavelength*sqrt(4*PI/monochromator->Maier_Leibnitz_reflectivity); + + // Kinenatic reflectivity = QML*DHKL*sin(theta_B)**2/PI/cos(theta_B) + double kinematic_reflectivity = monochromator->Maier_Leibnitz_reflectivity; + kinematic_reflectivity *= monochromator->lattice_spacing; + kinematic_reflectivity *= square(sine_of_bragg_angle); + kinematic_reflectivity *= 1/PI/cosine_of_bragg_angle; + kinematic_reflectivity *= monochromator->Debye_Waller_factor; + // Primary extinction factor, using the approximation in G.E Bacon and R.D. Lowde, Acta Cryst. (1948). 1, 303 + kinematic_reflectivity *= tanh(monochromator->domain_thickness/extinction_length)/monochromator->domain_thickness*extinction_length; + return kinematic_reflectivity; + } + + double calculate_attenuation_coefficient(struct Monochromator_values* monochromator, + struct neutron_values* neutron){ + double E = square(neutron->v_size)*VS2E; // Neutron energy in meV + // Get factor for single phonon cross section + + double Bernoulli_sequence[31] = {1,-0.5,0.166667,0,-0.033333,0,0.0238095,0,-0.033333, + 0,0.0757576,0,-0.253114,0,1.16667,0,-7.09216,0,54.9712, + 0,-529.124,0,6192.12,0,-86580.3,0,1.42551717e6,0,-2.7298231e7, + 0,6.01580874e8}; + double x; + if (monochromator->temperature_mono - 0.1 <= 0){ + x = monochromator->debye_temperature/0.1; + } + else{ + x = monochromator->debye_temperature/monochromator->temperature_mono; + } + double R, Ifact, Xn; + if (x<6){ + R = 0; + Ifact = 1; + Xn = 1/x; + for (int i=0; i<30; i++){ + R += Bernoulli_sequence[i]*Xn/Ifact*(i + 2.5); + Xn *= x; + Ifact *= i + 1; + } + } + else R = 3.3/sqrt(x*x*x*x*x*x*x); + + // Define boltzmann_constant in units of (meV/K) + double boltzmann_constant = 0.08617333262; + double DWMF = 1-exp(-(monochromator->B0+monochromator->BT)*monochromator->Constant_from_Freund_paper*E/1000); + // Set the cross sections, as written in freunds paper + monochromator->nuclear_capture_absorption = monochromator->incoherent_scattering_cross_section + +monochromator->absorption_for_1AA_Neutrons*neutron->wavelength; + + monochromator->multiple_phonon_absorption = monochromator->bound_atom_scattering_cross_section + *square(monochromator->atomic_number/(monochromator->atomic_number + 1)) + *DWMF; + + monochromator->single_phonon_absorption = 3*monochromator->bound_atom_scattering_cross_section/monochromator->atomic_number + * sqrt(boltzmann_constant * monochromator->debye_temperature/E) * R; + + double attenuation_coefficient = (monochromator->nuclear_capture_absorption + + monochromator->single_phonon_absorption + + monochromator->multiple_phonon_absorption) + /monochromator->volume * 100; // *100 to change to per mm? + return attenuation_coefficient; + } + /////////////////////////////////////////////////////////////////////////// + /////////////// Function that retrieves local scattering vector G or tau. + /////////////// TODO: + /////////////// Add mosaic stuff to this function. + /////////////////////////////////////////////////////////////////////////// + void calculate_local_scattering_vector(struct Monochromator_values* monochromator, + struct neutron_values* neutron, int direction){ + neutron->tau[0] = monochromator->tau[0]; + neutron->tau[1] = monochromator->tau[1]; + neutron->tau[2] = monochromator->tau[2]; + for (int i=0 ; i<3; i++) { + neutron->tau[i] += monochromator->lattice_spacing_gradient_field[i][0]*neutron->r[0] + +monochromator->lattice_spacing_gradient_field[i][1]*neutron->r[1] + +monochromator->lattice_spacing_gradient_field[i][2]*neutron->r[2]; + } + /* Renormalize local scat vect */ + double normalization_factor = sqrt(square(monochromator->tau[0]) + square(monochromator->tau[1]) + square(monochromator->tau[2])) + /sqrt(square(neutron->tau[0]) + square(neutron->tau[1]) + square(neutron->tau[2])); + + neutron->tau[0] *= direction*normalization_factor; + neutron->tau[1] *= direction*normalization_factor; + neutron->tau[2] *= direction*normalization_factor; + } + /////////////////////////////////////////////////////////////////////////// + /////////////// Function that sets the neutron structs values at a point + /////////////////////////////////////////////////////////////////////////// + void set_neutron_values( + struct neutron_values* neutron, + double x, double y, double z, + double vx, double vy, double vz){ + neutron->r[0] = x; + neutron->r[1] = y; + neutron->r[2] = z; + neutron->v[0] = vx; + neutron->v[1] = vy; + neutron->v[2] = vz; + neutron->v_size = 0; + neutron->ki_size = 0; + neutron->tau_size = 0; + neutron->kf_size = 0; + for (int i =0; i<3; i++){ + neutron->ki[i] = neutron->v[i]*V2K; + neutron->ki_size += square(neutron->ki[i]); + neutron->v_size += square(neutron->v[i]); + } + neutron->v_size = sqrt(neutron->v_size); + neutron->ki_size = sqrt(neutron->ki_size); + neutron->wavelength = 3956/neutron->v_size;// Wavelength in Angstrom. + } + /////////////////////////////////////////////////////////////////////////// + /////////////// Function that solves the Laue condition + /////////////// Solves |k + tau|^2 - k^2 = 0 after expanding with + /////////////// bending terms + /////////////////////////////////////////////////////////////////////////// + + void solve_Bragg_condition(struct neutron_values* neutron, + struct Monochromator_values* monochromator, + int direction, int lamella) { + calculate_local_scattering_vector(monochromator, neutron, direction); + neutron->tau_size = 0; + neutron->kf_size = 0; + for (int i=0; i<3; i++){ + neutron->kf[i] = neutron->ki[i] + neutron->tau[i]; + neutron->tau_size += square(neutron->tau[i]); + neutron->kf_size += square(neutron->kf[i]); + } + neutron->tau_size = sqrt(neutron->tau_size); + neutron->kf_size = sqrt(neutron->kf_size); + double a = 0; + double b = 0; + double c = square(neutron->kf_size) - square(neutron->ki_size); + double ksi = 0; + + for (int i=0; i<3; i++){ + ksi = (monochromator->lattice_spacing_gradient_field[i][0]*neutron->ki[0] + + monochromator->lattice_spacing_gradient_field[i][1]*neutron->ki[1] + + monochromator->lattice_spacing_gradient_field[i][2]*neutron->ki[2]); + a += ksi*ksi; + b += neutron->kf[i] * ksi; + } + neutron->Bragg_angle_gradient[lamella] = b/neutron->tau_size/( + neutron->kf[0]*monochromator->perp_to_tau[0] + + neutron->kf[1]*monochromator->perp_to_tau[1] + + neutron->kf[2]*monochromator->perp_to_tau[2] + ); + b *= direction; + //printf("a %g, b %g, c %g \n", a, b, c); + double det = b*b - a*c; + if (a> 1e-20 && det > 1e-20){ + neutron->time_of_reflection[lamella] = sign(b)*((-fabs(b)+sqrt(det))/a); + } + else {neutron->time_of_reflection[lamella] = -0.000001;} + neutron->time_of_reflection[lamella] *= neutron->ki_size/neutron->v_size; + } + /////////////////////////////////////////////////////////////////////////// + /////////////// Simple function to reflect a neutron + /////////////////////////////////////////////////////////////////////////// + void reflect_neutron(struct neutron_values* neutron, double* vx, double* vy, double* vz){ + *vx = (neutron->ki[0] + neutron->tau[0])*K2V; + *vy = (neutron->ki[1] + neutron->tau[1])*K2V; + *vz = (neutron->ki[2] + neutron->tau[2])*K2V; + } + + /////////////////////////////////////////////////////////////////////////// + /////////////// Function that scans the neutron path + /////////////// and finds the probability of reflection in each lamella, + /////////////// the accumulative probability, the path length in each lamella, + /////////////// as well as the point of reflection in each lamella. + /////////////////////////////////////////////////////////////////////////// + void get_crossing_times_of_lamellas(struct Monochromator_values* monochromator, struct neutron_values* neutron, + int lamella, int direction){ + double transposed_x; + double inner_t0; + double outer_t0; + double inner_t1; + double outer_t1; + transposed_x = neutron->r[0]-monochromator->radius_horizontal + -(monochromator->lamellas)*monochromator->thickness/2 + - (monochromator->lamellas-1)*monochromator->gap_size/2 + + monochromator->thickness * lamella + monochromator->gap_size * lamella; + cylinder_intersect(&inner_t0,&inner_t1, + transposed_x,neutron->r[1],neutron->r[2], + neutron->v[0],neutron->v[1],neutron->v[2], + monochromator->radius_horizontal, + monochromator->height); + cylinder_intersect(&outer_t0,&outer_t1, + transposed_x + monochromator->thickness,neutron->r[1],neutron->r[2], + neutron->v[0],neutron->v[1],neutron->v[2], + monochromator->radius_horizontal, + monochromator->height); + if (direction == 1){ + neutron->entry_time[lamella] = inner_t1; + neutron->exit_time[lamella] = outer_t1; + } else { + neutron->entry_time[lamella] = outer_t0; + neutron->exit_time[lamella] = inner_t0; + } + neutron->path_length[lamella] = neutron->exit_time[lamella] - neutron->entry_time[lamella]; + } + void transport_neutron_to_lamella_coordinates(struct Monochromator_values* monochromator, struct neutron_values* neutron, + int lamella, int direction){ + neutron->r[0] += neutron->v[0]*neutron->entry_time[lamella]; + neutron->r[1] += neutron->v[1]*neutron->entry_time[lamella]; + neutron->r[2] += neutron->v[2]*neutron->entry_time[lamella]; + neutron->r[0] += monochromator->thickness * lamella + monochromator->gap_size * lamella + -(monochromator->lamellas-1)*monochromator->thickness/2 + - (monochromator->lamellas-1)*monochromator->gap_size/2; + } + void transport_neutron_back_to_entry(struct Monochromator_values* monochromator, struct neutron_values* neutron, + int lamella, int direction){ + neutron->r[0] -= neutron->v[0]*neutron->entry_time[lamella]; + neutron->r[1] -= neutron->v[1]*neutron->entry_time[lamella]; + neutron->r[2] -= neutron->v[2]*neutron->entry_time[lamella]; + neutron->r[0] -= monochromator->thickness * lamella + monochromator->gap_size * lamella + -(monochromator->lamellas-1)*monochromator->thickness/2 + - (monochromator->lamellas-1)*monochromator->gap_size/2; + } + void propagate_neutrons_to_point_of_reflection(struct neutron_values* neutron, int lamella){ + neutron->r[0] += neutron->v[0]*neutron->time_of_reflection[lamella]; + neutron->r[1] += neutron->v[1]*neutron->time_of_reflection[lamella]; + neutron->r[2] += neutron->v[2]*neutron->time_of_reflection[lamella]; + } + void propagate_neutrons_to_lamella_entry(struct neutron_values* neutron, int lamella){ + neutron->r[0] -= neutron->v[0]*neutron->time_of_reflection[lamella]; + neutron->r[1] -= neutron->v[1]*neutron->time_of_reflection[lamella]; + neutron->r[2] -= neutron->v[2]*neutron->time_of_reflection[lamella]; + } + + void scan_lamellas(struct Monochromator_values* monochromator, struct neutron_values* neutron, + int current_lamella, int direction, int neutron_just_reflected){ + double t0, inner_t1, outer_t1; + double kinematic_reflectivity; + double r[3] = {neutron->r[0], neutron->r[1], neutron->r[2]}; + int intersected; + double transposed_x; + for (int i = current_lamella; i < monochromator->lamellas && i >= 0;){ + get_crossing_times_of_lamellas(monochromator, neutron, i, direction); + transport_neutron_to_lamella_coordinates(monochromator, neutron, i, direction); + solve_Bragg_condition(neutron, monochromator, direction, i); + propagate_neutrons_to_point_of_reflection(neutron, i); + /* Assign probabilities in arrays */ + if (!neutron_is_inside_monochromator(monochromator, neutron, i)) { + neutron->probabilities[i] = 0; + } else{ + kinematic_reflectivity = calculate_kinematic_reflectivity(monochromator, neutron); + neutron->probabilities[i] = 1 - exp(-neutron->ki_size*kinematic_reflectivity/neutron->Bragg_angle_gradient[i]); + } + + /* Assigning accumulative probabilities. + * If lamella being checked is the same as the starting lamella, set accumulating probability to the probability of that lamella. + * TODO: Ask Kim about this. It is done by SIMRES, but it makes me uncertain. If reflection has already occured within this lamella, set probability to 0. */ + if (i == current_lamella){ + if (neutron_just_reflected) neutron->probabilities[i] = 0; + neutron->accumulating_probabilities[i] = neutron->probabilities[i]; + } else { + neutron->accumulating_probabilities[i] = 1 - (1-neutron->accumulating_probabilities[i-direction])*(1-neutron->probabilities[i]); + } + /* Move neutron back to starting position */ + propagate_neutrons_to_lamella_entry(neutron, i); + transport_neutron_back_to_entry(monochromator, neutron, i, direction); + i += direction; + } + + } +%} + +DECLARE +%{ + double angle_range; + double neutron_counter; + double changed_neutrons_counter; + double curvature; + double tau_size_zero; + + struct neutron_values; + struct neutron_values neutron; + struct Monochromator_values monochromator; +%} + +INITIALIZE +%{ + /////////////////////////////////////////////////////////////////////////// + /////////////// ERROR FUNCTIONS + /////////////////////////////////////////////////////////////////////////// + if (radius_x <= 0) + exit(printf("tau_approach_bent_perfect_crystal: %s: incorrect radius_x=%g\n", NAME_CURRENT_COMP, radius_x)); + if (xthickness <= 0) + exit(printf("tau_approach_bent_perfect_crystal: %s: invalid monochromator xthickness=%g\n", NAME_CURRENT_COMP, xthickness)); + if (zwidth <= 0) + exit(printf("tau_approach_bent_perfect_crystal: %s: invalid monochromator zwidth=%g\n", NAME_CURRENT_COMP, zwidth)); + + /////////////////////////////////////////////////////////////////////////// + /////////////// INITIALIZING PARAMETERS + /////////////////////////////////////////////////////////////////////////// + neutron_counter = 0; + changed_neutrons_counter = 0; + /* Initialize angles of the monochromator */ + angle_range = zwidth/radius_x; + monochromator.max_angle = angle_range/2 + PI; + monochromator.min_angle = -angle_range/2 + PI; + /* Read the designated plane of reflection, for use in the monochromator. */ + /* Set Monochromator values */ + + enum crystal_plane plane = stringToEnum(&plane_of_reflection); + + monochromator.length = zwidth; + monochromator.type = bent; + monochromator.height = yheight; + monochromator.thickness = xthickness; + monochromator.radius_horizontal = radius_x; + monochromator.lamellas = lamella_slabs; + monochromator.gap_size = lamella_gap_size; + monochromator.domain_thickness = domainthickness; + monochromator.temperature_mono = temperature; + monochromator.lattice_spacing = crystal_table[plane][0]; + monochromator.Maier_Leibnitz_reflectivity = crystal_table[plane][1]*100; /* Convert to SI and Angstrom */ + monochromator.bound_atom_scattering_cross_section = crystal_table[plane][2]; + monochromator.absorption_for_1AA_Neutrons = crystal_table[plane][3]; + monochromator.incoherent_scattering_cross_section = crystal_table[plane][4]; + monochromator.volume = crystal_table[plane][5]; + monochromator.atomic_number = crystal_table[plane][6]; + monochromator.debye_temperature = crystal_table[plane][7]; + monochromator.Constant_from_Freund_paper = crystal_table[plane][9]; + monochromator.poisson_ratio = crystal_table[plane][9]; + calculate_B0_and_BT(&monochromator); + monochromator.Debye_Waller_factor = exp(-(monochromator.B0 + monochromator.BT)/2/square(monochromator.lattice_spacing)); + + /* Initialize reciprocal lattice vector G or tau in some texts, and perp_to_tau. */ + + angle_to_cut_horizontal *= DEG2RAD; + angle_to_cut_vertical *= DEG2RAD; + + tau_size_zero = 2*PI/monochromator.lattice_spacing; + + monochromator.tau[0] = tau_size_zero*cos(angle_to_cut_horizontal)*cos(angle_to_cut_vertical); + monochromator.tau[1] = tau_size_zero*sin(angle_to_cut_vertical); + monochromator.tau[2] = tau_size_zero*sin(angle_to_cut_horizontal)*cos(angle_to_cut_vertical); + + monochromator.perp_to_tau[0] = sin(angle_to_cut_horizontal)*cos(angle_to_cut_vertical); + monochromator.perp_to_tau[1] = sin(angle_to_cut_vertical); + monochromator.perp_to_tau[2] = -cos(angle_to_cut_horizontal)*cos(angle_to_cut_vertical); + + /* Initialize lattice_spacing_gradient_field */ + + curvature = 1/radius_x; + monochromator.lattice_spacing_gradient_field[0][0] = -monochromator.poisson_ratio*cos(angle_to_cut_horizontal)*tau_size_zero*curvature; + monochromator.lattice_spacing_gradient_field[0][1] = 0; + monochromator.lattice_spacing_gradient_field[0][2] = sin(angle_to_cut_horizontal)*tau_size_zero*curvature; + monochromator.lattice_spacing_gradient_field[1][0] = 0; + monochromator.lattice_spacing_gradient_field[1][1] = 0; + monochromator.lattice_spacing_gradient_field[1][2] = 0; + monochromator.lattice_spacing_gradient_field[2][0] = sin(angle_to_cut_horizontal)*tau_size_zero*curvature; + monochromator.lattice_spacing_gradient_field[2][1] = 0; + monochromator.lattice_spacing_gradient_field[2][2] = -cos(angle_to_cut_horizontal)*tau_size_zero*curvature; + + /* Initialize neutron structs values */ + + + neutron.ki = (double*) calloc (3, sizeof(double)); + neutron.r = (double*) calloc (3, sizeof(double)); + neutron.v = (double*) calloc (3, sizeof(double)); + neutron.tau = (double*) calloc (3, sizeof(double)); + neutron.kf = (double*) calloc (3, sizeof(double)); + neutron.Bragg_angle_gradient = (double*) calloc (lamella_slabs, sizeof(double)); + neutron.path_length = (double*) calloc (lamella_slabs, sizeof(double)); + neutron.entry_time = (double*) calloc (lamella_slabs, sizeof(double)); + neutron.exit_time = (double*) calloc (lamella_slabs, sizeof(double)); + neutron.probabilities = (double*) calloc (lamella_slabs, sizeof(double)); + neutron.accumulating_probabilities = (double*) calloc (lamella_slabs, sizeof(double)); + neutron.time_of_reflection = (double*) calloc (lamella_slabs, sizeof(double)); + if (verbose)printf("PLane of reflection %d, lattice vector %g\n", plane, monochromator.lattice_spacing); +%} + +TRACE +%{ + /* Initialize variables for use in TRACE */ + int current_lamella = 0; + double weight_init = p; + if (weight_init == 0.0) { + changed_neutrons_counter ++; + ABSORB; + } + double init_v_size; + double cutoff_minimum_relative_weight = 1e-10; + int neutron_is_in_crystal = 1; + int neutron_just_reflected = 0; + double reflect_condition; + int start_lamella; + double intersected; + double intersect_time_1 = 0; + double intersect_time_2 = 0; + int direction = 1; + double attenuation_coefficient; + intersected = cylinder_intersect(&intersect_time_1, + &intersect_time_2, + x-radius_x - lamella_slabs*xthickness/2 - (lamella_slabs-1)*lamella_gap_size/2, + y,z, + vx,vy,vz, + radius_x, yheight); + + if (verbose) printf("Neutron %i entered comp (from %i)\n",_particle->_uid,INDEX_CURRENT_COMP); + /* If the neutron has either intersected not at all, the top, or the bottom, ABSORB the neutron. */ + if (intersected == 0 || ((int)intersected)%8 == 2 || ((int)intersected)%8 == 4 ) { + if (verbose) printf("Neutron %i should proceed to next comp (from %i, did not intersect)\n",_particle->_uid,INDEX_CURRENT_COMP); + ABSORB; + } + /* If the neutron is not within the angle width of the monochromator now, ABSORB it. */ + PROP_DT(intersect_time_2); + double angle_on_inner_cylinder = PI - asin(z/monochromator.radius_horizontal); + if (monochromator.max_angle + 0.0001 <= angle_on_inner_cylinder || angle_on_inner_cylinder <= monochromator.min_angle - 0.0001){ + if (verbose) printf("Neutron %i should proceed to next comp (from %i, not within angles)\n",_particle->_uid,INDEX_CURRENT_COMP); + ABSORB; + } + if (verbose) print_neutron_state(&neutron); + /* Choose monochromator type. TODO: Implement more types of monochromators. */ + switch (monochromator.type){ + case flat: + break; + case bent: + neutron.path = 0; + /* This is an infinite loop that the neutrons should break out of. */ + while (neutron_is_in_crystal){ + + start_lamella = current_lamella; + /* We set the neutron values for the neutron, to reduce overhead in the share functions */ + set_neutron_values(&neutron, x,y,z,vx,vy,vz); + if (!neutron_just_reflected) init_v_size = neutron.v_size; + /* scan_lamellas calculates time of reflection, entry/exit times, probability of reflection, + * accumulating_propability of reflection, as well as the Bragg angle gradient for each lamella */ + scan_lamellas(&monochromator, &neutron, current_lamella, direction, neutron_just_reflected); + /* If the neutron has reflected in a lamella, the neutron cannot reflect in the same segment. Therefore, add + * the path length left in the lamella to the neutron path.*/ + if (neutron_just_reflected){ + neutron.path += neutron.exit_time[current_lamella]; + } + /* Choose reflection condition depending on the whether the neutron is incoming or outgoing. */ + switch(direction){ + case (1): + reflect_condition = neutron.accumulating_probabilities[monochromator.lamellas - 1]*rand01(); + break; + case(-1): + reflect_condition = 1*rand01(); + break; + } + /* Find the lamella the neutron reflects from, or the final lamella the neutron is in. */ + while(neutron.accumulating_probabilities[current_lamella]< reflect_condition && + current_lamella < monochromator.lamellas && current_lamella >= 0){ + current_lamella += direction; + } + /* Add the neutrons path through the lamellas it has not reflected from, to the total neutron path. */ + for (int i = start_lamella + direction; direction*(current_lamella - i) > 0; i += direction){ + neutron.path += neutron.path_length[i]; + } + /* Allow neutron to break from the component if the relative weight is too low.*/ + if (p*neutron.accumulating_probabilities[current_lamella]/weight_init < cutoff_minimum_relative_weight){ + break; + } + /* If the neutron did not reflect on it's way out, allow it to go to the next component.*/ + if (current_lamella == -1 && direction == -1) { + break; + } + /* Add the distance to the point of reflection in the lamella, to the neutron path, propagate to that point, + * reflect and modify the weight of the neutron. */ + neutron.path += neutron.time_of_reflection[current_lamella]; + PROP_DT(neutron.time_of_reflection[current_lamella] + neutron.entry_time[current_lamella]); + SCATTER; + if (verbose) printf("Neutron %i SCATTERED (from %i)\n\nGROUP BREAK!!\n",_particle->_uid,INDEX_CURRENT_COMP); + p *= neutron.accumulating_probabilities[current_lamella]; + + neutron.r[0] = x; + neutron.r[1] = y; + neutron.r[2] = z; + calculate_local_scattering_vector(&monochromator, &neutron, direction); + reflect_neutron(&neutron, &vx, &vy, &vz); + direction *= -1; + neutron_just_reflected++; + } + /* Attenuate the neutron with the remaining length through the crystal. */ + switch (direction){ + case (1): + neutron.path += neutron.exit_time[monochromator.lamellas - 1]; + case (-1): + neutron.path += neutron.exit_time[0]; + break; + } + attenuation_coefficient = calculate_attenuation_coefficient(&monochromator, &neutron); + p *= exp(-attenuation_coefficient*neutron.path*neutron.v_size); + /* Temporary measure */ + if (p > weight_init + 0.00000001*p) ABSORB; + /* End of algorithm */ + break; + case mosaic: + break; + + case bent_mosaic: + break; + } +// if (!SCATTERED) {RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);} + neutron_counter += 1; + + if (fabs(neutron.v_size - init_v_size)/neutron.v_size>0.003 && verbose){ + printf("init speed %g \nFinal speed %g\nInit weight %f final weight %g\n", + init_v_size, neutron.v_size, weight_init, p); + print_neutron_state(&neutron); + // ABSORB; + } + //if (lambda_init != neutron.wavelength) print_neutron_state(&neutron); +// printf("neutron wavelength incoming %g outgoing %g\n\n", lambda_init, neutron.wavelength); + +%} + +FINALLY +%{ + if (verbose) printf("\nThe amount of neutrons reflected in the monochromator are %g.\nThe amount of neutrons that changed %g\n", neutron_counter, changed_neutrons_counter); +%} + +MCDISPLAY +%{ + /* Draw the monochromator as two curving lines, between 10 points. */ + + double x_inner [2]; + double x_outer [2]; + double y_top; + double y_bottom; + double z_inner [2]; + double z_outer [2]; + double i = 0; + double inner_radii = radius_x + xthickness*lamella_slabs/2 + (lamella_slabs -1)*lamella_gap_size/2; + double outer_radii = radius_x - xthickness*lamella_slabs/2 - (lamella_slabs -1)*lamella_gap_size/2; + x_inner[1] = radius_x - cos(angle_range/2)*inner_radii; + y_top = yheight/2; + z_inner[1] = -sin(angle_range/2)*inner_radii; + + x_outer[1] = radius_x - cos(angle_range/2)*outer_radii; + y_bottom = -yheight/2; + z_outer[1] = -sin(angle_range/2)*outer_radii; + + + for (i = 0; i < 2; i = i + 0.2) { + + x_inner[0] = x_inner[1]; + + z_inner[0] = z_inner[1]; + x_inner[1] = radius_x - cos(fabs(i-1)*angle_range/2)*inner_radii; + z_inner[1] = -sin(angle_range/2)*radius_x + sin(angle_range/2)*inner_radii*i; + + x_outer[0] = x_outer[1]; + z_outer[0] = z_outer[1]; + + x_outer[1] = radius_x - cos(fabs(i-1)*angle_range/2)*outer_radii; + z_outer[1] = -sin(angle_range/2)*outer_radii + sin(angle_range/2)*outer_radii*i; + + multiline(5, + x_inner[0], y_top, z_inner[0], + x_outer[0], y_top, z_outer[0], + x_outer[0], y_bottom, z_outer[0], + x_inner[0], y_bottom, z_inner[0], + x_inner[0], y_top, z_inner[0]); + + + line(x_inner[0], y_top, z_inner[0], + x_inner[1], y_top, z_inner[1]); + line(x_outer[0], y_top, z_outer[0], + x_outer[1], y_top, z_outer[1]); + line(x_inner[0], y_bottom, z_inner[0], + x_inner[1], y_bottom, z_inner[1]); + line(x_outer[0], y_bottom, z_outer[0], + x_outer[1], y_bottom, z_outer[1]); + + } + multiline(5, + x_inner[1], y_top, z_inner[1], + x_outer[1], y_top, z_outer[1], + x_outer[1], y_bottom, z_outer[1], + x_inner[1], y_bottom, z_inner[1], + x_inner[1], y_top, z_inner[1]); + /* line(0,0,0, + monochromator.tau[0],monochromator.tau[1],monochromator.tau[2]); + line(0,0,0, + monochromator.perp_to_tau[0], monochromator.perp_to_tau[1], monochromator.perp_to_tau[2]); */ +%} + +END + + diff --git a/mcstas-comps/examples/Tests/Test_Monochromator_bent.instr b/mcstas-comps/examples/Tests/Test_Monochromator_bent.instr new file mode 100644 index 0000000000..7f07b9a31d --- /dev/null +++ b/mcstas-comps/examples/Tests/Test_Monochromator_bent.instr @@ -0,0 +1,93 @@ +/******************************************************************************** +* +* McStas, neutron ray-tracing package +* Copyright (C) 1997-2008, All rights reserved +* Risoe National Laboratory, Roskilde, Denmark +* Institut Laue Langevin, Grenoble, France +* +* This file was written by McStasScript, which is a +* python based McStas instrument generator written by +* Mads Bertelsen in 2019 while employed at the +* European Spallation Source Data Management and +* Software Centre +* +* Instrument Test_Monochromator_bent +* +* %Identification +* Written by: Daniel Lomholt Christensen +* Date: 14:13:59 on January 23, 2024 +* Origin: ILL +* %INSTRUMENT_SITE: Generated_instruments +* +* +* %Parameters +* +* %End +********************************************************************************/ + +DEFINE INSTRUMENT Test_Monochromator_bent ( +double mono_rotation = 43.5337 // Rotation of the monochromator, relative to the z axis. +) + +DECLARE +%{ +int reflections; +%} + +INITIALIZE +%{ +// Start of initialize for generated Test_Monochromator_bent +%} + +TRACE +COMPONENT Source = Source_gen( + dist = 10, focus_xw = 0.001, + focus_yh = 0.01, lambda0 = 1.5, + dlambda = 7.5E-06, I1 = 7.95775E+12, + yheight = 0.01, xwidth = 0.001, + T1 = 315) +AT (0,0,0) ABSOLUTE + +COMPONENT Slit = Slit( + xwidth = 0.001, yheight = 0.01) +AT (0,0,10) RELATIVE Source + +COMPONENT entry_monitor = PSD_monitor( + nx = 200, ny = 200, + filename = "entry_monitor", xwidth = 0.02, + yheight = 0.02, restore_neutron = 1) +AT (0,0,0) RELATIVE Slit + +COMPONENT monochromator_arm = Arm() +AT (0,0,10.15) RELATIVE Source +ROTATED (0,mono_rotation,0) RELATIVE Source + +COMPONENT Monochromator = Monochromator_bent( + zwidth = 0.07, yheight = 0.012, + xthickness = 0.008, radius_x = 10, + lamella_slabs = 1, lamella_gap_size = 0.0001, + plane_of_reflection = "Ge511", angle_to_cut_horizontal = -19.47, + angle_to_cut_vertical = 0) +AT (0,0,0) RELATIVE monochromator_arm +ROTATED (0,-19.47,0) RELATIVE monochromator_arm + +COMPONENT arm_after_monochromator = Arm() +AT (0,0,0) RELATIVE monochromator_arm +ROTATED (0,43.533702255651555,0) RELATIVE monochromator_arm + +COMPONENT monitor_arm = Arm() +AT (0,0,0) RELATIVE arm_after_monochromator +ROTATED (0,0,0) RELATIVE arm_after_monochromator + +COMPONENT PSD_exit_monitor = PSD_monitor( + nx = 200, ny = 200, + filename = "det2d.dat", xwidth = 0.1, + yheight = 0.02, restore_neutron = 1) +AT (0,0,0.5) RELATIVE monitor_arm + +FINALLY +%{ +// Start of finally for generated Test_Monochromator_bent +%} + +END