/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright (C) 1997-2015, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: PerfectCrystal * * %I * Written by: Markus Appel * Date: 2015-12-21 * Origin: ILL / FAU Erlangen-Nuernberg * Based on a perfect crystal component by: Miguel A. Gonzalez, A. Dianoux June 2013 (ILL) * * Changelog: * Version 1.1 * - BUGFIX: correct neutron energy shift in Doppler mode * - added option 'debyescherrer' to select analyzer geometry * - added option 'facette' to approximate analyzer sphere by small, flat crystals * * Version 1.0 * - inital release * * * %D * Perfect crystal component, primarily for use as monochromator and analyzer in * backscattering spectrometers. Reflection from a single Bragg reflex of a flat or * spherical surface is simulated using a Darwin, Ewald or Gaussian profile. Doppler * movement of the monochromator is supported as well. * * Orientation of the crystal surface is *different* from other monochromator components! * Gravitational energy shift for tall analyzers should work in principle, but it not tested yet. * See the parameter description on how to define the geometry and properties. * * [1] Website for Backscattering Spectroscopy: http://www.ill.eu/sites/BS-review/index.htm * * Examples: * IN16B (ILL) Si111 large angle analyzers (approximate dimensions): * PerfectCrystal(radius=2, lambda=6.2708, sigma=0.24e-3, * ttmin=40, ttmax=165, phimin=-45, phimax=+45, centerfocus=1) * * IN16B (ILL) Si111 Doppler monochromator: * PerfectCrystal(radius=2.165, lambda=6.2708, sigma=0.24e-3, * width = 0.500, height = 0.250, centerfocus=0, * speed = 4.7, amplitude = 0.075, exclusive=1) * * %P * (A) Size, shape and position: * ============================= * radius: [m] Radius of curvature, set to 0 for a flat surface. Default: 0 * centerfocus: [0/1] Component origin is the center of the analyzer sphere if 1, if set to 0 the origin is on the analyzer surface. Default: 0 * * option 1 * ttmin: [deg] * ttmax: [deg] analyzer coverage angle in horizontal xz-plane between -180 and 180 * phimin: [deg] * phimax: [deg] vertical analyzer coverage between -90 and 90 (-180 and 180 if debyescherrer==1) * option 2 * tt0: [deg] * ttwidth: [deg] horizontal coverage as center and full width * phi0: [deg] * phiwidth: [deg] vertical coverage as center and full width * option 3 * tt0: [deg] angular center position in the horizontal plane (only if centerfocus==1) * width: [m] absolute width * phi0: [deg] angular vertical position (only if centerfocus==1) * height: [m] absolute height (only if debyescherrer==0) * * debyescherrer: [-180...180] (0/1) bend analyzer following a Debye-Scherrer ring along scattering angle tt (twotheta) (phi = * facette: [m] width of square crystal facettes arranged on the spherical surface (set 0 to disable). Default: 0 Warning: Facettation will fail at the poles along +-y axis. * facette_xi: [deg] random misalignemt of each facette. Default: 0 * * (B) Neutron optics * =================== * tau: [A^-1] Scattering vector of the reflex (sometimes also called Q...) * lambda: [A] Alternatively to tau: backscattered wavelength * R0: [] Peak reflectivity. Default: 1 * * Ewald/Darwin mode * dtauovertau: [1] Plateau width of the Ewald/Darwin curve (see * dtauovertau_ext: [] Relative variation of tau (randomized for each trajectory, full width) * ewald: [0/1] Use Ewald curve if 1, Darwin curve if 0. Default: 1 (Ewald) * * Gaussian mode * sigma: [meV] Width of the Gaussian reflectivity curve in meV (corresponding to the energy resolution). (The width will be transformed and the Gaussian is actually calculated in k-space.) * * Mirror mode * ismirror: [0/1] Simply reflect all neutrons if 1. Good for debugging/visualization. Default: 0 * * (C) Doppler movement for monochromators (sinusoidal speed profile): * =================================================================== * speed: [m/s] Maximum Doppler speed. The actual monochromator velocity is randomized between -speed and +speed. Default: 0 * amplitude: [m] Amplitude of the Doppler movement. Default: 0 * * smartphase: [0/1] Optimize Doppler phase for better MC efficiency if set to 1. *WARNING:* Experimental option! Always compare to a simulation without smartphase. Better do not use smartwidth with Ewald/Darwin curves due to their endless tails. Default: 0 * smartwidth: [meV] Half width of the possible energy reflection window for smartphase. Default: 5*sigma OR 10*2*dtauovertau*E0 * * (D) Miscellaneous * ================== * exclusive: [0/1] If set to 1, absorb all neutrons that missed the monochromator/analyzer surface. Default: 0 * transmit: [0...1] Monte-Carlo probability of transmitting an event through the monochromator/analyzer surface. (Events with R=1.0 for DarwinE/Ewald curves are always reflected!). Default: 0 * * Output parameters * ================== * tt: [deg] Position where the neutron hit (or missed) the analyzer sphere * phi: [deg] Position where the neutron hit (or missed) the analyzer sphere * xi: [deg] Reflection angle between neutron trajectory and analyzer surface normal * phid: [rad] Actual doppler phase during reflection [0,2*PI] * vd: [m/s] Actual doppler speed during reflection * zd: [m] Actual doppler displacement during reflection (zd>0 => displaced towards inc. beam from -z) * v0: [m/s] Backscattered neutron velocity * E0: [meV] Backscattered neutron energy * R: [0...1] Reflectivity value (on Gauss/Darwin/Ewald curve) * eps: [see 1] 'abs(y)' from the Darwin/Ewald reflectivity formula * * %E *******************************************************************************/ DEFINE COMPONENT PerfectCrystal SETTING PARAMETERS ( ttmin=NAN, ttmax=NAN, tt0=NAN, ttwidth=NAN, width=NAN, phimin=NAN, phimax=NAN, phi0=NAN, phiwidth=NAN, height=NAN, debyescherrer=0, facette=0, facette_xi=0, centerfocus=0, radius=0, tau=NAN, lambda=NAN, R0=1.0, dtauovertau=NAN, dtauovertau_ext=0, ewald=1,sigma=NAN,ismirror=0,speed=0,amplitude=0,smartphase=0, smartwidth=NAN, exclusive=0, transmit=0, verbose=0 ) SHARE %{ // convert between spherical (r,tt,phi) and cartesian (x,y,z) coordinates (angles in deg) // debyescherrer swaps axes void sph2cart (double* x, double* y, double* z, double r, double tt, double phi, int debyescherrer) { tt *= DEG2RAD; phi *= DEG2RAD; if (debyescherrer) { double sintt = sin (tt); *x = -r * cos (phi) * sintt; *y = r * sin (phi) * sintt; *z = r * cos (tt); } else { double cosphi = cos (phi); *x = -r * cosphi * sin (tt); *y = r * sin (phi); *z = r * cosphi * cos (tt); } } /******************************************************************************* * grandvec_target_circle: Choose random direction towards target at (x,y,z) * with given radius and gaussian area distribution. * If radius is zero, choose random direction in full 4PI, no target. ******************************************************************************/ void grandvec_target_circle (double* xo, double* yo, double* zo, double* solid_angle, double xi, double yi, double zi, double radius) { double l2, phi, theta, nx, ny, nz, xt, yt, zt, xu, yu, zu; if (radius == 0.0) { /* No target, choose uniformly a direction in full 4PI solid angle. */ theta = acos (1 - rand0max (2)); phi = rand0max (2 * PI); if (solid_angle) *solid_angle = 4 * PI; nx = 1; ny = 0; nz = 0; yi = sqrt (xi * xi + yi * yi + zi * zi); zi = 0; xi = 0; } else { double costheta0; l2 = xi * xi + yi * yi + zi * zi; /* sqr Distance to target. */ costheta0 = sqrt (l2 / (radius * radius + l2)); if (radius < 0) costheta0 *= -1; if (solid_angle) { /* Compute solid angle of target as seen from origin. */ *solid_angle = 2 * PI * (1 - costheta0); } /* Now choose point uniformly on circle surface within angle theta0 */ double costheta; costheta = (1 - fabs (randnorm () * (1 - costheta0))); if (costheta < -1) costheta = -1; theta = acos (costheta); /* radius on circle */ phi = rand0max (2 * PI); /* rotation on circle at given radius */ /* Now, to obtain the desired vector rotate (xi,yi,zi) angle theta around a perpendicular axis u=i x n and then angle phi around i. */ if (xi == 0 && zi == 0) { nx = 1; ny = 0; nz = 0; } else { nx = -zi; nz = xi; ny = 0; } } /* [xyz]u = [xyz]i x n[xyz] (usually vertical) */ vec_prod (xu, yu, zu, xi, yi, zi, nx, ny, nz); /* [xyz]t = [xyz]i rotated theta around [xyz]u */ rotate (xt, yt, zt, xi, yi, zi, theta, xu, yu, zu); /* [xyz]o = [xyz]t rotated phi around n[xyz] */ rotate (*xo, *yo, *zo, xt, yt, zt, phi, xi, yi, zi); } /* randvec_target_circle */ %} DECLARE %{ // official output variables double tt; double phi; double xi; double phid; double vd; double zd; double v0; double vmin; double vmax; double vperp; double E0; double R; double eps; // internal stuff double sin_facette_xi; %} INITIALIZE %{ // checks and balances ... if (radius < 0) { MPI_MASTER (fprintf (stderr, "PerfectCrystal: %s ERROR: Invalid parameters: negative radius\n", NAME_CURRENT_COMP);); exit (-1); } //****************************************** // position and size if (!radius) { // flat analyzer surface if (isnan (width) || isnan (height)) { MPI_MASTER (fprintf (stderr, "PerfectCrystal: %s ERROR: Invalid parameters: Need width and height for radius==0\n", NAME_CURRENT_COMP);); exit (-1); } } else { // curved analyzer surface //****************************************** // Determine analyzer width // ttmin / ttmax if (!isnan (ttmin) && !isnan (ttmax) && isnan (tt0) && isnan (ttwidth) && isnan (width)) { // all right. } // tt0 / ttwidth else if (isnan (ttmin) && isnan (ttmax) && !isnan (ttwidth) && isnan (width)) { if (isnan (tt0)) { MPI_MASTER (fprintf (stderr, "PerfectCrystal: %s WARNING: Missing parameter: tt0 set to zero\n", NAME_CURRENT_COMP);); tt0 = 0; } ttmin = tt0 - ttwidth / 2.0; ttmax = tt0 + ttwidth / 2.0; } // tt0 / width else if (isnan (ttmin) && isnan (ttmax) && isnan (ttwidth) && !isnan (width)) { if (isnan (tt0)) { MPI_MASTER (fprintf (stderr, "PerfectCrystal: %s WARNING: Missing parameter: tt0 set to zero\n", NAME_CURRENT_COMP);); tt0 = 0; } ttwidth = 2.0 * asin (width / 2.0 / radius) * RAD2DEG; ttmin = tt0 - ttwidth / 2.0; ttmax = tt0 + ttwidth / 2.0; } else { MPI_MASTER (fprintf (stderr, "PerfectCrystal: %s ERROR: Invalid parameters: Cannot determine analyzer width/tt\n", NAME_CURRENT_COMP);); exit (-1); } //****************************************** // Determine analyzer height // phimin / phimax if (!isnan (phimin) && !isnan (phimax) && isnan (phi0) && isnan (phiwidth) && isnan (height)) { // all right, nothing to do. } // phi0 / phiwidth else if (isnan (phimin) && isnan (phimax) && !isnan (phiwidth) && isnan (height)) { if (isnan (phi0)) { MPI_MASTER (fprintf (stderr, "PerfectCrystal: %s WARNING: Missing parameter: phi0 set to zero\n", NAME_CURRENT_COMP);); phi0 = 0; } phimin = phi0 - phiwidth / 2.0; phimax = phi0 + phiwidth / 2.0; } // phi0 / height else if (isnan (phimin) && isnan (phimax) && isnan (phiwidth) && !isnan (height)) { if (isnan (phi0)) { MPI_MASTER (fprintf (stderr, "PerfectCrystal: %s WARNING: Missing parameter: phi0 set to zero\n", NAME_CURRENT_COMP);); phi0 = 0; } phiwidth = 2.0 * asin (height / 2.0 / radius) * RAD2DEG; phimin = phi0 - phiwidth / 2.0; phimax = phi0 + phiwidth / 2.0; } else { MPI_MASTER (fprintf (stderr, "PerfectCrystal: %s ERROR: Invalid parameters: Cannot determine analyzer height/phi\n", NAME_CURRENT_COMP);); exit (-1); } } if (centerfocus && !radius) { MPI_MASTER (fprintf (stderr, "PerfectCrystal: %s ERROR: Invalid parameters: centerfocus doesn't make sense with radius==0\n", NAME_CURRENT_COMP);); exit (-1); } //****************************************** // neutron optics if (!ismirror) { if (isnan (lambda) == isnan (tau)) { MPI_MASTER (fprintf (stderr, "PerfectCrystal: %s ERROR: Invalid parameters: provide either tau or lambda\n", NAME_CURRENT_COMP);); exit (-1); } if (isnan (tau)) tau = 4 * PI / lambda; v0 = tau / 2.0 * K2V; E0 = SQR (v0) * VS2E; if (isnan (dtauovertau) == isnan (sigma)) { MPI_MASTER ( fprintf (stderr, "PerfectCrystal: %s ERROR: Invalid parameters: provide either sigma or dtauovertau, or switch on ismirror\n", NAME_CURRENT_COMP);); exit (-1); } } else { if (!isnan (dtauovertau) || !isnan (sigma)) { MPI_MASTER (fprintf (stderr, "PerfectCrystal: %s WARNING: dtauovertau and/or sigma is ignored with ismirror==1\n", NAME_CURRENT_COMP);); } } if (!(R0 >= 0 && R0 <= 1.0)) { MPI_MASTER (fprintf (stderr, "PerfectCrystal: %s ERROR: Invalid parameter: R0 must be between 0 and 1\n", NAME_CURRENT_COMP);); exit (-1); } if (transmit < 0.0 || transmit > 1.0) { MPI_MASTER (fprintf (stderr, "PerfectCrystal: %s ERROR: Invalid parameter: transmit must be between 0 and 1\n", NAME_CURRENT_COMP);); exit (-1); } // transform Gaussian sigma from energy (meV) to neutron velocity (m/s) if (!isnan (sigma)) sigma /= 3.29106e-3 * tau; // constant is hbar/2 in appropriate units // transform smartwidth from energy (meV) to neutron velocity (m/s) if (!isnan (smartwidth)) smartwidth /= 3.29106e-3 * tau; // constant is hbar/2 in appropriate units // set standard widths for smartphase if (isnan (smartwidth) && !isnan (sigma)) smartwidth = 5.0 * sigma; if (isnan (smartwidth) && !isnan (dtauovertau)) smartwidth = 10.0 * dtauovertau * v0; if (smartphase && !isnan (dtauovertau)) { MPI_MASTER (fprintf (stderr, "PerfectCrystal: %s WARNING: Using smartphase for ewald/darwin curves is probably a bad idea, because these curves have very long wings " "which will be cut off!!!\n", NAME_CURRENT_COMP);); } if (facette_xi) { if (facette_xi < 0) { MPI_MASTER (fprintf (stderr, "PerfectCrystal: %s ERROR: facette_xi must be >= 0\n", NAME_CURRENT_COMP);); exit (-1); } sin_facette_xi = sin (DEG2RAD * facette_xi); } // talk to the user ... if (verbose) { #define PRINTVAR(str,val) fprintf(stderr,"%s --- %s : %g\n",NAME_CURRENT_COMP,str,val); MPI_MASTER (PRINTVAR ("ttmin", ttmin); PRINTVAR ("ttmax", ttmax); PRINTVAR ("ttwidth", ttwidth); PRINTVAR ("width", width); PRINTVAR ("phimin", phimin); PRINTVAR ("phimax", phimax); PRINTVAR ("phi0", phi0); PRINTVAR ("phiwidth", phiwidth); PRINTVAR ("height", height); PRINTVAR ("centerfocus", centerfocus); PRINTVAR ("debyescherrer", debyescherrer); PRINTVAR ("radius", radius); PRINTVAR ("facette", facette); PRINTVAR ("facette_xi", facette_xi); PRINTVAR ("sin_facette_xi", sin_facette_xi); PRINTVAR ("tau", tau); PRINTVAR ("lambda", lambda); PRINTVAR ("v0", v0); PRINTVAR ("E0", E0); PRINTVAR ("dtauovertau", dtauovertau); PRINTVAR ("dtauovertau_ext", dtauovertau_ext); PRINTVAR ("ewald", ewald); PRINTVAR ("R0", R0); PRINTVAR ("sigma[m/s]", sigma); PRINTVAR ("ismirror", ismirror); PRINTVAR ("speed", speed); PRINTVAR ("amplitude", amplitude); PRINTVAR ("smartphase", smartphase); PRINTVAR ("smartwidth[m/s]", smartwidth); PRINTVAR ("exclusive", exclusive);); } %} TRACE %{ double dt1, dt2, q0mod; double nx, ny, nz; int missed = 0; // determine phase of doppler movement if ((speed != 0) && smartphase) { // do something smart vmin = v0 - sqrt (SQR (vx) + SQR (vy) + SQR (vz)) - smartwidth; vmax = vmin + 2.0 * smartwidth; if ((vmin > speed) || (vmax < -speed)) ABSORB; if (vmin < -speed) vmin = -speed; if (vmax > speed) vmax = speed; // fprintf(stderr,"%g %g\n",vmin,vmax); vd = vmin + (vmax - vmin) * rand01 (); phid = acos (vd / speed) + (rand01 () < 0.5 ? 0 : PI); zd = amplitude * sin (phid); p *= (vmax - vmin) / 2.0 / speed; } else if (speed != 0) { // random selection of phase // phid = rand01() * 2.0 * PI; // zd = amplitude * sin(phid); // vd = speed * cos(phid); // random selection of speed phid = acos (randpm1 ()) + (rand01 () < 0.5 ? 0 : PI); vd = speed * cos (phid); zd = amplitude * sin (phid); } else { // no movement zd = 0; vd = 0; } // Propagate to analyzer and determine surface normal if (!radius) { // flat analyzer, use height and width // propoagate to surface if (!vz) ABSORB; dt2 = (-z - zd) / vz; PROP_DT (dt2); // see if the covered area is hit missed = !inside_rectangle (x, y, width, height); // surface normal in this case is simple nx = 0; ny = 0; nz = -1; } else { // spherical analyzer, use spherical coordinates ... // compute neutron path intersection with analyzer sphere if (centerfocus) missed = !sphere_intersect (&dt1, &dt2, x, y, z + zd, vx, vy, vz, radius); else missed = !sphere_intersect (&dt1, &dt2, x, y, z + radius + zd, vx, vy, vz, radius); if (!missed) { // propoagate to surface PROP_DT (dt2); // tt (twotheta) is calculated as in IN16B, positive values downstream rightwards // select coordinate system depending on 'debyescherrer' switch double zprime = z + (centerfocus ? 0 : radius) + zd; if (debyescherrer) { tt = RAD2DEG * atan2 (sqrt (SQR (x) + SQR (y)), zprime); phi = RAD2DEG * atan2 (y, -x); } else { tt = -RAD2DEG * atan2 (x, zprime); phi = RAD2DEG * asin (y / sqrt (SQR (x) + SQR (y) + SQR (zprime))); } missed = (tt < ttmin || tt > ttmax || phi < phimin || phi > phimax); if (!missed) { // analyzer surface normal if (facette) { // calculate center of the facette hit // always use spherical coordinates with y main axis (as in debyescherrer=0), // otherwise the pole will be in the analyzer center! // for the sake of confusion, use radians in this part. double tt1, ttfacette, ttn, phi1, phifacette, phin; if (debyescherrer) { tt1 = -atan2 (x, zprime); phi1 = asin (y / sqrt (SQR (x) + SQR (y) + SQR (zprime))); } else { tt1 = DEG2RAD * tt; phi1 = DEG2RAD * phi; } phifacette = facette / radius; phin = floor (phi1 / phifacette + 0.5) * phifacette; ttfacette = facette / (radius * cos (phin)); ttn = floor (tt1 / ttfacette + 0.5) * ttfacette; nx = cos (phin) * sin (ttn); ny = -sin (phin); nz = -cos (phin) * cos (ttn); if (facette_xi) { grandvec_target_circle (&nx, &ny, &nz, NULL, nx, ny, nz, sin_facette_xi); } } else { nx = -x; ny = -y; nz = -z - (centerfocus ? 0 : radius) - zd; NORM (nx, ny, nz); } } } } // Do the reflection if (!missed) { // velocity vector projected on surface normal // in moving doppler frame vperp = scalar_prod (vx, vy, vz + vd, nx, ny, nz); // angle between surface normal and velocity (only used as output parameter for monitoring) xi = RAD2DEG * acos (-vperp / sqrt (SQR (vx) + SQR (vy) + SQR (vz + vd))); // energy selection if (!ismirror) { if (!isnan (dtauovertau)) { // eps is actually abs(y) // vperp is negative! double this_tau = tau; if (dtauovertau_ext) { this_tau *= 1.0 + 0.5 * dtauovertau_ext * randpm1 (); } eps = fabs (4.0 * vperp * V2K / this_tau + 2.0) / dtauovertau; // Darwin/Ewald curve if (eps > 1) { if (ewald) { // energy selection with Ewald curve R = 1.0 - sqrt (SQR (eps) - 1.0) / eps; } else { // energy selection with Darwin curve R = eps - sqrt (SQR (eps) - 1.0); R *= R; } R *= R0; } else { R = R0; } } else { // Gauss curve (vperp is negative!) eps = fabs (v0 + vperp) / sigma; R = exp (-SQR (eps) / 2.0); R *= R0; } } else { R = 1; eps = NAN; } if (transmit && (R != 1) && (rand01 () < transmit)) { p *= (1.0 - R) / transmit; SCATTER; } else { // reflection: vector kf = ki - qmod0 where qmod0 = 2*n*(n dot ki'), // n is the surface normal vector and ki' the incident wavevector // in the moving frame q0mod = 2.0 * vperp; // q0mod is now directly in velocity units (and contains energy shift due to Doppler effect!) vx -= q0mod * nx; vy -= q0mod * ny; vz -= q0mod * nz; if (R != 1) p *= R / (1.0 - transmit); if (R < 1e-10) ABSORB; else SCATTER; } } else if (!exclusive) { R = NAN; xi = NAN; RESTORE_NEUTRON (INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); } else ABSORB; %} MCDISPLAY %{ int ttsegments; int phisegments; int steps; // const double centerspheresize = 0.1*radius; int i, itt, iphi, istep; double phi, tt, dphi, dtt; double x1, y1, z1, x2, y2, z2; double z0; if (debyescherrer) { ttsegments = 4; phisegments = 20; steps = 50; } else { ttsegments = 10; phisegments = 10; steps = 30; } // mark the center with a sphere // circle("xy",0,0,0,centerspheresize); // circle("xz",0,0,0,centerspheresize); // circle("yz",0,0,z0,centerspheresize); // draw crystal surface if (radius) { z0 = (centerfocus ? 0 : radius); // curved surface dtt = (ttmax - ttmin) / (double)steps; dphi = (phimax - phimin) / (double)steps; for (itt = 0; itt <= ttsegments; itt++) { tt = ttmin + (ttmax - ttmin) * (double)itt / (double)ttsegments; for (istep = 0; istep < steps; istep++) { phi = phimin + dphi * istep; sph2cart (&x1, &y1, &z1, radius, tt, phi, debyescherrer); sph2cart (&x2, &y2, &z2, radius, tt, phi + dphi, debyescherrer); line (x1, y1, z1 - z0, x2, y2, z2 - z0); } } for (iphi = 0; iphi <= phisegments; iphi++) { phi = phimin + (phimax - phimin) * (double)iphi / (double)phisegments; for (istep = 0; istep < steps; istep++) { tt = ttmin + dtt * istep; sph2cart (&x1, &y1, &z1, radius, tt, phi, debyescherrer); sph2cart (&x2, &y2, &z2, radius, tt + dtt, phi, debyescherrer); line (x1, y1, z1 - z0, x2, y2, z2 - z0); } } // additional dashed moving one if (amplitude) { sph2cart (&x1, &y1, &z1, radius, ttmin, phimin, debyescherrer); dashed_line (x1, y1, z1 - z0 - amplitude, x1, y1, z1 - z0 + amplitude, 6); sph2cart (&x1, &y1, &z1, radius, ttmax, phimin, debyescherrer); dashed_line (x1, y1, z1 - z0 - amplitude, x1, y1, z1 - z0 + amplitude, 6); sph2cart (&x1, &y1, &z1, radius, ttmin, phimax, debyescherrer); dashed_line (x1, y1, z1 - z0 - amplitude, x1, y1, z1 - z0 + amplitude, 6); sph2cart (&x1, &y1, &z1, radius, ttmax, phimax, debyescherrer); dashed_line (x1, y1, z1 - z0 - amplitude, x1, y1, z1 - z0 + amplitude, 6); for (i = -1; i <= 1; i += 2) { z0 = (centerfocus ? 0 : radius) + i * amplitude; for (itt = 0; itt <= ttsegments; itt++) { tt = ttmin + (ttmax - ttmin) * (double)itt / (double)ttsegments; for (istep = 0; istep < steps; istep++) { phi = phimin + dphi * istep; sph2cart (&x1, &y1, &z1, radius, tt, phi, debyescherrer); sph2cart (&x2, &y2, &z2, radius, tt, phi + dphi, debyescherrer); dashed_line (x1, y1, z1 - z0, x2, y2, z2 - z0, 1); } } for (iphi = 0; iphi <= phisegments; iphi++) { phi = phimin + (phimax - phimin) * (double)iphi / (double)phisegments; for (istep = 0; istep < steps; istep++) { tt = ttmin + dtt * istep; sph2cart (&x1, &y1, &z1, radius, tt, phi, debyescherrer); sph2cart (&x2, &y2, &z2, radius, tt + dtt, phi, debyescherrer); dashed_line (x1, y1, z1 - z0, x2, y2, z2 - z0, 1); } } } } } else { // flat surface rectangle ("xy", 0.0, 0.0, 0.0, width, height); // additional moving one if (amplitude) { dashed_line (+width / 2.0, +height / 2.0, -amplitude, +width / 2.0, +height / 2.0, +amplitude, 6); dashed_line (+width / 2.0, -height / 2.0, -amplitude, +width / 2.0, -height / 2.0, +amplitude, 6); dashed_line (-width / 2.0, +height / 2.0, -amplitude, -width / 2.0, +height / 2.0, +amplitude, 6); dashed_line (-width / 2.0, -height / 2.0, -amplitude, -width / 2.0, -height / 2.0, +amplitude, 6); for (i = -1; i <= 1; i += 2) rectangle ("xy", 0.0, 0.0, i * amplitude, width, height); } } %} END