/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright (C) 1997-2011, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: curved_mirror * * %I * Written by: Henrik Jacobsen * Date: April 2012 * Origin: RNBI * * Single mirror plate that is curved and fits into an elliptic guide. * * %D * * * %P * INPUT PARAMETERS: * * focus_start_w: [m] Horizontal position of first focal point of ellipse in ABSOLUTE COORDINATES * focus_start_h: [m] Vertical position of first focal point of ellipse in ABSOLUTE COORDINATES * focus_end_w: [m] Horizontal position of second focal point of ellipse in ABSOLUTE COORDINATES * focus_end_h: [m] VErtical position of second focal point of ellipse in ABSOLUTE COORDINATES * mirror_start: [m] start of mirror in ABSOLUTE COORDINATES * guide_start: [m] start of guide in ABSOLUTE COORDINATES * smallaxis_w: [m] Small-axis dimension of horizontal ellipse * smallaxis_h: [m] Small-axis dimension of vertical ellipse * yheight: [m] the height of the mirror when not in the guide * smallaxis: [m] the smallaxis of the guide * length: [m] the length of the mirror * R0=0.99: [] low angle reflectivity * Qc: [AA-1] critical scattering vector * alpha: [AA] slope of reflectivity * m: [] m-value of material * W: [AA-1] width of supermirror cutoff * transmit: [0/1] 0: non reflected neutrons are absorbed. 1: non reflected neutrons pass through * substrate_thickness: [m] thickness of substrate (for absorption) * coating_thickness: [m] thickness of coating (for absorption) * theta_1: [deg] angle of mirror wrt propagation direction at start of mirror * theta_2: [deg] angle of mirror wrt propagation direction at center of mirror * theta_3: [deg] angle of mirror wrt propagation direction at end of mirror * reflect: [q(Angs-1) R(0-1)] (str) Name of relfectivity file. Format * * * * %D * * %E *******************************************************************************/ DEFINE COMPONENT Mirror_Elliptic_Bispectral SETTING PARAMETERS (string reflect=0,focus_start_w,focus_end_w, focus_start_h, focus_end_h, mirror_start, m, smallaxis_w, smallaxis_h, length, transmit=0, substrate_thickness=0.0005, coating_thickness=10e-6) SHARE %{ %include "read_table-lib" %} DECLARE %{ t_Table pTable; %} INITIALIZE %{ if (reflect && strlen (reflect) && strcmp (reflect, "NULL") && strcmp (reflect, "0")) { if (Table_Read (&pTable, reflect, 1) <= 0) /* read 1st block data from file into pTable */ exit (fprintf (stderr, "Mirror: %s: can not read file %s\n", NAME_CURRENT_COMP, reflect)); } %} TRACE %{ double f; // half of distance between focal points double asquared; double a; // half of ellipse length double b; // half of ellipse width double xprime; // x in coordinates with center of ellipse at (xprime,zprime)=(0,0) double ymirror; // height of the mirror // Defining the mirror double a1; double b1; double c1; // solving the time the neutron hits the sample double A, B, C, D, E, P, Q, R, U, V, I, J, K; // finding rotation of mirror double alpha1, beta1, gamma1; double theta_m; double xhi, zeta; double b_w; double f_w; double asquared_w; double A_w; double B_w; double C_w; double determinant_w; double b_h; double f_h; double asquared_h; double xprime_h; double xprime_w; double zprime_w; double asquare_z; double bsquare_x; double v_n; // speed of neutron perpendicular to surface double Ref; // reflectivity double dt; double q; int intersect; double discriminant; double xprime_start_w; double zprime_start_w; double z_test; double x_test; double z_prime_test; int hit_back_flag; int prop_case; double x_2; double y_2; double z_2; double t_2; int x_hit; int x_hit_2; double xprime_h_2; double ymirror_2; int intersect_2; intersect = 0; x_2 = 0; y_2 = 0; z_2 = 0; t_2 = -1; prop_case = 0; double old_x = x, old_y = y, old_z = z, old_t = t, old_vx = vx, old_vz = vz, old_vy = vy; // Check if neutron hits mirror. First find which z,x coordinates it hits. // define the ellipse b_w = smallaxis_w / 2; f_w = (focus_end_w - focus_start_w) * 0.5; asquared_w = f_w * f_w + b_w * b_w; // in coordinate system of mirror: xprime_w(t)=xprime_w+vx*t. z-value is zprime(t)=b*sqrt(1-x'^2/a^2)+old_z+vz*t. This gives equation for t: // A_w*t^2+B_w*t+C_w=0; xprime_start_w = old_x - f_w - focus_start_w + mirror_start; zprime_start_w = z; A_w = b_w * b_w * vx * vx + asquared_w * vz * vz; B_w = 2 * b_w * b_w * xprime_start_w * vx + 2 * asquared_w * old_z * vz; C_w = b_w * b_w * xprime_start_w * xprime_start_w + asquared_w * old_z * old_z - asquared_w * b_w * b_w; // this equation must now be solved for t; if (A_w != 0) { determinant_w = B_w * B_w - 4.0 * A_w * C_w; if (determinant_w >= 0) { I = (-B_w - sqrt (determinant_w)) / (2.0 * A_w); J = (-B_w + sqrt (determinant_w)) / (2.0 * A_w); if (I <= 0 && J <= 0) { dt = -1.0; } // both times are negative. if (I <= 0 && J > 0) { dt = J; } // set dt to only positive value. if (I > 0 && J <= 0) { dt = I; } // set dt to only positive value. if (I > 0 && J > 0) { prop_case = 1; if (I > J) { dt = J; } else { dt = I; } } // set dt to smallest positive value } else { dt = -1.0; } // only complex solutions: set dt negative so no scattering } else { // end if (A!=0) if (B_w != 0) { dt = -C / B_w; } else { // end if (B!)=0 printf ("warning: A_w=B_w=C_w=0. Neutron is ignored\n"); } } // end if (A!=0){}else{ // now intersection time has been found. // printf("dt=%f\n",dt); if (dt > 0) { // if time is positive, propagate neutron to where it hits mirror. This is done without gravity. x += vx * dt; y += vy * dt; z += vz * dt; t += dt; if (prop_case > 0) // also check if neutron can hit mirror at second solution - it might not be in y-range for first solution { x_2 = x + vx * fabs (J - I); y_2 = y + vy * fabs (J - I); z_2 = z + vz * fabs (J - I); t_2 = t + fabs (J - I); } else { x_2 = x; y_2 = y; z_2 = z; t_2 = t; } x_hit = (x >= 0 && x <= length); x_hit_2 = (x_2 >= 0 && x_2 <= length); // printf("x=%f,y=%f,z=%f\n",x,y,z); // if (x >=0 && x<=length){ //check if neutron is within x limits of the mirror. If so, check if it is within y limits. if (x_hit || x_hit_2) { // define the ellipse b_h = smallaxis_h / 2; f_h = (focus_end_h - focus_start_h) * 0.5; asquared_h = f_h * f_h + b_h * b_h; xprime_h = -f_h - focus_start_h + mirror_start + x; // xprime is the x-coordinate in a coordinate system centered at the center of the ellipse ymirror = b_h * sqrt (1 - xprime_h * xprime_h / asquared_h); xprime_h_2 = -f_h - focus_start_h + mirror_start + x_2; // xprime is the x-coordinate in a coordinate system centered at the center of the ellipse ymirror_2 = b_h * sqrt (1 - xprime_h_2 * xprime_h_2 / asquared_h); intersect = (y >= -ymirror && y <= ymirror && x >= 0 && x <= length && zprime_start_w + vz * dt >= 0); intersect_2 = (y_2 >= -ymirror_2 && y_2 <= ymirror_2 && x_2 >= 0 && x_2 <= length && zprime_start_w + vz * (dt + fabs (J - I)) >= 0); if (!intersect && intersect_2) { // if neutron doesn't hit mirror with smallest t, but hits with largest t, propagte to largest t intersect = intersect_2; x = x_2; y = y_2; z = z_2; t = t_2; dt = t_2 - old_t; // printf("x=%f,y=%f,z=%f\n",x,y,z); } if (intersect) { // if neutron is within ylimits of the mirror handle reflection/transmission // now perform reflection. // First find out if neutron hits front or back of mirror: propagate backwards a bit and check if neutron is outside ellipse: if so, it hits back of // mirror b_w = smallaxis_w / 2.0; f_w = (focus_end_w - focus_start_w) * 0.5; asquared_w = f_w * f_w + b_w * b_w; z_test = zprime_start_w + vz * (dt - 1e-6); x_test = xprime_start_w + vx * (dt - 1e-6); z_prime_test = b_w * sqrt (1 - x_test * x_test / asquared_w); // find velocity in q direction. xprime_w = xprime_start_w + vx * dt; zprime_w = zprime_start_w + vz * dt; asquare_z = asquared_w * zprime_w; bsquare_x = b_w * b_w * xprime_w; zeta = (asquare_z) / (sqrt (asquare_z * asquare_z + bsquare_x * bsquare_x)); xhi = -(bsquare_x) / (sqrt (asquare_z * asquare_z + bsquare_x * bsquare_x)); // printf("z_test=%f, z_prime_test=%f\n",z_test,z_prime_test); if (z_test > z_prime_test) { hit_back_flag = 1; } // printf("vx=%f, vz=%f, vy=%f, xhi=%f, zeta=%f, prop_case=%d\n",vx,vz,vy,xhi,zeta,prop_case); v_n = -xhi * vx + zeta * vz; q = fabs (2.0 * v_n * V2Q); // Reflectivity parameters calculated from SWISS neutronics data. double R0 = 0.99; double Qc = 0.0217; double m_value = m * 0.9853 + 0.1978; double W = -0.0002 * m_value + 0.0022; double alpha = 0.1204 * m_value + 5.0944; double beta = -7.6251 * m_value + 68.1137; if (m_value <= 3) { alpha = m_value; beta = 0; } /* Reflectivity (see component Guide). */ if (m == 0) ABSORB; if (reflect && strlen (reflect) && strcmp (reflect, "NULL") && strcmp (reflect, "0")) Ref = Table_Value (pTable, q, 1); else { Ref = R0; if (q > Qc) { double arg = (q - m_value * Qc) / W; if (arg < 10) Ref *= .5 * (1 - tanh (arg)) * (1 - alpha * (q - Qc) + beta * (q - Qc) * (q - Qc)); // matches data from Swiss Neutronics else Ref = 0; } } if (Ref < 0) Ref = 0; else if (Ref > 1) Ref = 1; // Now comes actual reflection/transmission if (!transmit) { // all neutrons are reflected // printf("v_n=%f,q=%f, Ref=%f, lambda=%f, theta=%f, 1p_before=%f",v_n,q, Ref, (2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz),asin((2*PI/V2K)/sqrt(vx*vx + vy*vy // + vz*vz)*q/(4*PI))*RAD2DEG,p); if (!Ref) ABSORB; p *= Ref; // printf("p_after=%f\n",p); // handle reflection: change v_n -->-v_n vx = old_vx * (zeta * zeta - xhi * xhi) + old_vz * (2 * zeta * xhi); vz = +old_vx * (2 * zeta * xhi) + old_vz * (xhi * xhi - zeta * zeta); SCATTER; // after transmission or reflection } else { // if neutrons can be transmitted // calculate absorption. // substrate double lambda = (2 * PI / V2K) / sqrt (vx * vx + vy * vy + vz * vz); double sin_theta = lambda * q / (4 * PI); double substrate_path_length = substrate_thickness / sin_theta; double mu_substrate = 0.0318 / lambda + 0.0055 * lambda - 0.0050; // unit: cm^-1 mu_substrate = mu_substrate * 100; // unit: m^-1; // For nickel and titanium coating, the following formular is used: // mu = rho/m(atom)*sigma_a,thermal*lambda/lambda_thermal // lambda_thermal=1.798 Å // rho_nickel=8.908g/cm^3 // m(atom)_nickel=58.6934*1.661*10^-27 kg // sigma_a,thermal_nickel=4.49*10^-28 m^2 // rho_titanium=4.506g/cm^3 // m(atom)_titanium=47.867*1.661*10^-27 kg // sigma_a,thermal_titanium=6.09*10^-28 m^2 double coating_path_length = coating_thickness / sin_theta; double Ni_coefficient = 22.8180; double Ti_coefficient = 19.1961; double mu_coating = (0.5 * Ni_coefficient + 0.5 * Ti_coefficient) * lambda; // it is roughly 50% nickel and 50% titanium // transmit when rand > R if (Ref == 0 || rand01 () >= Ref) { // transmit if (substrate_thickness > 0) { p = p * exp (-mu_substrate * substrate_path_length - mu_coating * coating_path_length); } // reduce weight of neutrons due to attenuation in the mirror } else { // neutron is reflected if (hit_back_flag == 1 && substrate_thickness > 0) { // if neutron comes from behind the mirror p = p * exp (-2 * mu_substrate * substrate_path_length - 2 * mu_coating * coating_path_length); } // reduce weight of neutrons due to attenuation in the mirror // handle reflection: change v_n -->-v_n vx = old_vx * (zeta * zeta - xhi * xhi) - old_vz * (2 * zeta * xhi); vz = -old_vx * (2 * zeta * xhi) + old_vz * (xhi * xhi - zeta * zeta); } // printf("p_before=%f, q=%f",p,q); SCATTER; // after transmission or reflection // printf("p_after=%f\n",p); } // end } else { after if (!transmit) { } } // end if (x >=-length/2 && x<=length/2) if (!intersect) { /* No intersection: restore neutron position. */ x = old_x; y = old_y; z = old_z; t = old_t; } } // end if (dt>0) %} MCDISPLAY %{ double xstart; double xend; double xprime_start; double ystart; double xprime_end; double yend; double focus_2_h; double focus_1_h; double b_h; double f_h; double asquared_h; double focus_2_w; double focus_1_w; double b_w; double f_w; double asquared_w; int n_lines; int j = 0; double xprimepos[51]; double ypos[51]; double x_plot[51]; double zpos[51]; double xstep; double xprime_w[51]; focus_2_h = focus_end_h; // focus in local coordinates focus_1_h = focus_start_h; b_h = smallaxis_h / 2.0; f_h = (focus_2_h - focus_1_h) * 0.5; asquared_h = f_h * f_h + b_h * b_h; focus_2_w = focus_end_w; // focus in local coordinates focus_1_w = focus_start_w; b_w = smallaxis_w / 2.0; f_w = (focus_2_w - focus_1_w) * 0.5; asquared_w = f_w * f_w + b_w * b_w; xstart = 0; xend = length; n_lines = 50; xstep = length / n_lines; for (j = 0; j < n_lines + 1; j++) { xprimepos[j] = -f_h - focus_start_h + mirror_start + xstart + xstep * j; ypos[j] = b_h * sqrt (1 - xprimepos[j] * xprimepos[j] / asquared_h); // correct x_plot[j] = xstart + xstep * j; xprime_w[j] = -f_w - focus_start_w + mirror_start + xstart + xstep * j; // xprime is the x-coordinate in a coordinate system centered at the center of the ellipse zpos[j] = b_w * sqrt (1 - xprime_w[j] * xprime_w[j] / asquared_w); // printf("xprime=%f, zpos=%f\n",xprime_w[j], zpos[j]); } for (j = 0; j < n_lines; j++) { line (x_plot[j], -ypos[j], zpos[j], x_plot[j + 1], -ypos[j + 1], zpos[j + 1]); line (x_plot[j], ypos[j], zpos[j], x_plot[j + 1], ypos[j + 1], zpos[j + 1]); } line (x_plot[0], -ypos[0], zpos[0], x_plot[0], ypos[0], zpos[0]); line (x_plot[50], -ypos[50], zpos[50], x_plot[50], ypos[50], zpos[50]); // printf("ypos0=%f xpos0=%f ypos50=%f xpos50=%f",ypos[0], x_plot[0], ypos[50], x_plot[50]); /* double xmax, xmin, ymax, ymin; if (center == 0) { xmax= x1; xmin=0; ymax= yheight; ymin=0; } else { xmax= x1/2; xmin=-xmax; ymax= yheight/2; ymin=-ymax; } multiline(5, (double)xmin, (double)ymin, 0.0, (double)xmax, (double)ymin, 0.0, (double)xmax, (double)ymax, 0.0, (double)xmin, (double)ymax, 0.0, (double)xmin, (double)ymin, 0.0); */ %} END