/******************************************************************************* * * McXtrace, x-ray tracing package * Copyright, All rights reserved * DTU Physics, Kgs. Lyngby, Denmark * Synchrotron SOLEIL, Saint-Aubin, France * * Component: Mirror_toroid * * %Identification * * Written by: Erik B Knudsen * Date: Jul 2016 * Version: 1.0 * Origin: DTU Physics * * Toroidal shape mirror (in XZ) * * %Description * This is an implementation of a toroidal mirror which may be curved in two dimensions. * To avoid solving quartic equations, the intersection is compited as a combination of * two intersections. First, the ray is intersected with a cylinder to catch (almost) the small * radius curvature. Secondly, the ray is the intersected with an ellipsoid, with the curvatures * matching that of the torus. * * In the first incarnation we assume the mirror to be curving outwards (a bump). * * Example: Mirror_toroid(zdepth=0.340,xwidth=0.020,radius=246.9254,radius_o=246.9254,R0=1) * * %Parameters * Input parameters: * R0: [1] Reflectivity of mirror. * xwidth: [m] Width of mirror. * zdepth: [m] Length of mirror. * coating: [str] Datafile containing either mirror material constants or reflectivity numbers. * radius: [m] Curvature radius * radius_o:[m] Curvature radius, outwards * * %End *******************************************************************************/ DEFINE COMPONENT Mirror_toroid SETTING PARAMETERS (string coating="", zdepth=0.1, xwidth=0.01, radius, radius_o,R0=1) /* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */ SHARE %{ %include "read_table-lib" %include "reflectivity-lib" %} DECLARE %{ t_Reflec re; %} INITIALIZE %{ int status; if (coating && strlen (coating) && strcmp (coating, "NULL")) { status = reflec_Init (&re, COATING_UNDEFINED, coating, NULL); } else { status = reflec_Init_const (&re, R0); } if (status != 0) { fprintf (stderr, "ERROR (%s): Could not interpret reflectivity. Aborting.\n", NAME_CURRENT_COMP); exit (-1); } %} TRACE %{ int status; double k, l0, l1, l2, l3, tx, ty, tz; do { k = sqrt (scalar_prod (kx, ky, kz, kx, ky, kz)); /*TODO check the permutation of coordinates and the actual position of the cylinder.*/ status = cylinder_intersect (&l0, &l1, x, z, y + radius, kx, kz, ky, radius, zdepth); if (!status) break; /* no interaction with the cylinder*/ if (status & (02 | 04)) break; /*we exit the top/bottom of the cylinder*/ /*if the first intersection is behind the particle - this means the ray is on the wrong side of the mirror*/ if (l0 < 0) break; tx = x + l0 * kx / k; ty = y + l0 * ky / k; tz = z + l0 * kz / k; /*check mirror limits of intersection point.*/ if (tz < -zdepth / 2.0 || tz > zdepth / 2.0) break; /*check if the width is OK*/ double xmax = sin(xwidth / (2.0 * radius)) * radius; if (tx < -xmax || tx > xmax) break; status = ellipsoid_intersect(&l2, &l3, x, y+radius, z, kx, ky, kz, radius, radius, sqrt(radius*(radius_o+radius)), NULL); if (!status) break; if (l2 < 0) break; /*This shouldn't be possible*/ /*the mirror is indeed hit*/ PROP_DL (l2); SCATTER; double nx, ny, nz; nx = 2 * x / (radius * radius); ny = 2 * (y + radius) / (radius * radius); /*ellipsoid is displaced to put Origin on the surface*/ nz = 2 * z / (radius*(radius+radius_o)); NORM (nx, ny, nz); /*By definition the normal vector points out from the ellipsoid*/ double s = scalar_prod (kx, ky, kz, nx, ny, nz); kx = kx - 2 * s * nx; ky = ky - 2 * s * ny; kz = kz - 2 * s * nz; /*the get ref. by call to reflectivity lib*/ double R, q; q = 2.0 * s; R = reflecq (re, q, 0, k, fabs (90 - acos (s / k) * RAD2DEG)); p *= R; phi += M_PI; } while (0); %} MCDISPLAY %{ // See parametric equation used for the 3D view here: https://en.wikipedia.org/wiki/Torus#Geometry or here:https://web.maths.unsw.edu.au/~rsw/Torus/index.php // Because the code above doesn't actually reproduce a torus exactly but uses a cylinder and an ellipsoid to construct the toroidal mirror (to avoid solving // quartic equations), there's a slight mismatch between the 3D view (exact torus) and the actual toroidal mirror coded above. double theta, phi, phi_total, theta_total, theta_beginning, phi_beginning, theta_end, phi_end; // int N=50; too many points, makes the 3D matlab viewer bug int N = 10; phi_total = RAD2DEG * (zdepth / (radius + radius_o)); // degrees theta_total = RAD2DEG * (xwidth / radius); // fix theta=0 and phi=180 degrees respectively theta = 0; phi = 180; // then start at the beginning of the theta or phi range in order to map the flat surface (xwidth, zdepth) on the torus. theta_beginning = theta - (theta_total / 2); phi_beginning = phi - (phi_total / 2); theta_end = theta + (theta_total / 2); phi_end = phi + (phi_total / 2); int i, j; for (i = 0; i < N; i++) { for (j = 0; j < N; j++) { double t0 = DEG2RAD * (theta_beginning + i * theta_total / N), t1 = DEG2RAD * (theta_beginning + (i + 1) * theta_total / N); double p0 = DEG2RAD * (phi_beginning + j * phi_total / N), p1 = DEG2RAD * (phi_beginning + (j + 1) * phi_total / N); polygon (4, radius * sin (t0), (radius_o + radius * cos (t0)) * cos (p0) + (radius_o + radius), (radius_o + radius * cos (t0)) * sin (p0), radius * sin (t1), (radius_o + radius * cos (t1)) * cos (p0) + (radius_o + radius), (radius_o + radius * cos (t1)) * sin (p0), radius * sin (t1), (radius_o + radius * cos (t1)) * cos (p1) + (radius_o + radius), (radius_o + radius * cos (t1)) * sin (p1), radius * sin (t0), (radius_o + radius * cos (t0)) * cos (p1) + (radius_o + radius), (radius_o + radius * cos (t0)) * sin (p1)); } } line (-xwidth / 2.0, 0, -zdepth / 2.0, xwidth / 2.0, 0, -zdepth / 2.0); line (-xwidth / 2.0, 0, zdepth / 2.0, xwidth / 2.0, 0, zdepth / 2.0); line (-xwidth / 2.0, 0, -zdepth / 2.0, -xwidth / 2.0, 0, zdepth / 2.0); line (xwidth / 2.0, 0, -zdepth / 2.0, xwidth / 2.0, 0, zdepth / 2.0); %} END