/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2003, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Guide_curved
*
* %I
* Written by: Ross Stewart
* Date: November 18 2003
* Origin: ILL (France).
* Modified by: E. Farhi, uniformize parameter names (Jul 2008)
*
* Non-focusing curved neutron guide.
*
* %D
* Models a rectangular curved guide tube with entrance centered on the Z axis.
* The entrance lies in the X-Y plane. Draws a true depiction
* of the guide, and trajectories. Guide is not focusing.
*
* Example: Guide_curved(w1=0.1, h1=0.1, l=2.0, R0=0.99, Qc=0.021,
* alpha=6.07, m=2, W=0.003, curvature=2700)
*
* %BUGS
* This component does not work with gravitation on. Use component Guide_gravity then.
* Systematic error on transmitted flux is found to be about 10%.
*
* %P
* INPUT PARAMETERS:
*
* w1: [m] Width at the guide entry
* h1: [m] Height at the guide entry
* l: [m] length of guide
* R0: [1] Low-angle reflectivity
* Qc: [AA-1] Critical scattering vector
* alpha: [AA] Slope of reflectivity
* m: [1] m-value of material. Zero means completely absorbing.
* W: [AA-1] Width of supermirror cut-off
* curvature: [m] Radius of curvature of the guide
*
* %L
* Bender
*
* %E
*******************************************************************************/
DEFINE COMPONENT Guide_curved
DEFINITION PARAMETERS ()
SETTING PARAMETERS (w1, h1, l, R0=0.995, Qc=0.0218, alpha=4.38, m=2, W=0.003, curvature=2700)
OUTPUT PARAMETERS ()
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
%include "ref-lib"
%}
INITIALIZE
%{
if (mcgravitation) fprintf(stderr,"WARNING: Guide_curved: %s: "
"This component produces wrong results with gravitation !\n"
"Use Guide_gravity.\n",
NAME_CURRENT_COMP);
%}
TRACE
%{
double t11, t12, t21, t22, theta, alphaAng, endtime, phi;
double time, time1, time2, q, R;
int ii, i_bounce;
double whalf = 0.5*w1, hhalf = 0.5*h1; /* half width and height of guide */
double z_off = curvature*sin(l/curvature); /* z-component of total guide length */
double R1 = curvature - whalf; /* radius of curvature of inside mirror */
double R2 = curvature + whalf; /* radius of curvature of outside mirror */
double vel = sqrt(vx*vx + vy*vy + vz*vz); /* neutron velocity */
double vel_xz = sqrt(vx*vx + vz*vz); /* in plane velocity */
double K = V2K*vel; /* neutron wavevector */
double lambda = 2.0*PI/K; /* neutron wavelength */
/* Propagate neutron to guide entrance. */
PROP_Z0;
if(x <= -whalf || x >= whalf || y <= -hhalf || y >= hhalf)
ABSORB;
SCATTER;
for(;;)
{
double par[]={R0, Qc, alpha, m, W};
/* Find itersection points of neutron with inside and outside guide walls */
ii = cylinder_intersect(&t11, &t12 ,x - curvature, y, z, vx, vy, vz, R1, h1);
ii = cylinder_intersect(&t21, &t22 ,x - curvature, y, z, vx, vy, vz, R2, h1);
/* Choose appropriate reflection time */
time1 = (t11 < 1e-7) ? t12 : t11;
time2 = (t21 < 1e-7) ? t22 : t21;
time = (time1 < 1e-7 || time2 < time1) ? time2 : time1;
/* Has neutron left the guide? */
endtime = (z_off - z)/vz;
if (time > endtime || time <= 1e-7) break;
PROP_DT(time);
/* Find reflection surface */
R = (time == time1) ? R1 : R2;
i_bounce = (fabs(y - hhalf) < 1e-7 || fabs(y + hhalf) < 1e-7) ? 2 : 1;
switch(i_bounce) {
case 1: /* Inside or Outside wall */
phi = atan(vx/vz); /* angle of neutron trajectory */
alphaAng = asin(z/R); /* angle of guide wall */
theta = fabs(phi - alphaAng); /* angle of reflection */
vz = vel_xz*cos(2.0*alphaAng - phi);
vx = vel_xz*sin(2.0*alphaAng - phi);
break;
case 2: /* Top or Bottom wall */
theta = fabs(atan(vy/vz));
vy = -vy;
break;
}
/* Now compute reflectivity. */
if (m == 0 || !R0) ABSORB;
q = 4.0*PI*sin(theta)/lambda;
StdReflecFunc(q, par, &R);
if (R >= 0) p *= R; else ABSORB;
SCATTER;
}
%}
MCDISPLAY
%{
double x1, x2, z1, z2;
double xplot1[100], xplot2[100], zplot1[100], zplot2[100];
int n = 100;
int j = 1;
double R1 = (curvature - 0.5*w1); /* radius of inside arc */
double R2 = (curvature + 0.5*w1); /* radius of outside arc */
for(j=0;j