/******************************************************************************* * * 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