/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright 1997-2002, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: Guide_wavy * * %I * Written by: Kim Lefmann * Origin: Risoe * * Neutron guide with gaussian waviness. * * %D * Models a rectangular guide tube centered on the Z axis. The entrance lies * in the X-Y plane. * For details on the geometry calculation see the description in the McStas * reference manual. * * Example: m=2 Qc=0.0218 (nat. Ni) W=1/300 alpha=4.38 R0=0.995 (given by Daniel Clemens, PSI) * * %BUGS * This component does not work with gravitation on. Use component Guide_gravity then. * * %P * INPUT PARAMETERS: * * w1: [m] Width at the guide entry * h1: [m] Height at the guide entry * w2: [m] Width at the guide exit * h2: [m] Height at the guide exit * 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 for all mirrors * * alpha1: [AA] Slope of reflectivity left * m1: [1] m-value of material, left * W1: [AA-1] Width of supermirror cut-off left * alpha2: [AA] Slope of reflectivity right * m2: [1] m-value of material, right. * W2: [AA-1] Width of supermirror cut-off right * alpha3: [AA] Slope of reflectivity top * m3: [1] m-value of material, top. * W3: [AA-1] Width of supermirror cut-off top * alpha4: [AA] Slope of reflectivity bottom * m4: [1] m-value of material, bottom. * W4: [AA-1] Width of supermirror cut-off bottom * * wavy_z: [deg] Waviness in the z-(flight-)direction * wavy_xy: [deg] Waviness in the transverse direction * * %E ******************************************************************************/ DEFINE COMPONENT Guide_wavy DEFINITION PARAMETERS () SETTING PARAMETERS (w1, h1, w2=0, h2=0, l, R0=0.995, Qc=0.0218, alpha=0, m=0, W=0, alpha1=4.38, m1=2, W1=0.003, alpha2=4.38, m2=2, W2=0.003, alpha3=4.38, m3=2, W3=0.003, alpha4=4.38, m4=2, W4=0.003, wavy_z=0, wavy_xy=0) OUTPUT PARAMETERS (whalf,hhalf, lwhalf, lhhalf, norm_nv, norm_nh, f_h, f_v, eta_z, eta_xy) /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE %{ %include "ref-lib" %} DECLARE %{ double whalf,hhalf, lwhalf, lhhalf, norm_nv, norm_nh, f_h, f_v, eta_z, eta_xy; %} INITIALIZE %{ if (m) { m1=m2=m3=m4=m; } if (alpha) { alpha1=alpha2=alpha3=alpha4=alpha; } if (W) { W1=W2=W3=W4=W; } if (!w2) w2=w1; if (!h2) h2=h1; f_h = 0.5*(w2 - w1), f_v = 0.5*(h2 - h1); whalf = 0.5*w1, hhalf = 0.5*h1; lwhalf = l*whalf, lhhalf = l*hhalf; norm_nv = sqrt(l*l+f_v*f_v); norm_nh = sqrt(l*l+f_h*f_h); eta_z = wavy_z/(sqrt(8*log(2))); /* Convert from FWHM to Gaussian sigma */ eta_xy = wavy_xy/(sqrt(8*log(2))); if (mcgravitation) fprintf(stderr,"WARNING: Guide_wavy: %s: " "This component produces wrong results with gravitation !\n" "Use Guide_gravity.\n", NAME_CURRENT_COMP); %} TRACE %{ double t_min,t_tmp; /* Intersection times. */ double av,ah,bv,bh,cv1,cv2,ch1,ch2,d; /* Intermediate values */ double vdotn_v1,vdotn_v2,vdotn_h1,vdotn_h2; /* Dot products. */ double vdotn; double d_xy, d_z; /* Random angles */ double m0, alpha0, W; int i; /* Which mirror hit? */ double q; /* Q [1/AA] of reflection */ double dvx, dvy, dvz; /* Velocity change */ double vlen2,nlen2, norm_n2; /* Vector lengths squared */ double nperp,pz,nxy; /* for dot products */ double R; /* Reflectivity */ /* Propagate neutron to guide entrance. */ PROP_Z0; /* Scatter here to ensure that fully transmitted neutrons will not be absorbed in a GROUP construction, e.g. all neutrons - even the later absorbed ones are scattered at the guide entry. */ SCATTER; if(x <= -whalf || x >= whalf || y <= -hhalf || y >= hhalf) ABSORB; for(;;) { /* Compute the dot products of v and n for the four mirrors. */ av = -l*vx ; bv = f_h*vz; ah = -l*vy ; bh = f_v*vz; vdotn_v1 = bv + av; /* Left vertical */ vdotn_v2 = bv - av; /* Right vertical */ vdotn_h1 = bh + ah; /* Lower horizontal */ vdotn_h2 = bh - ah; /* Upper horizontal */ /* Compute the dot products of (O - r) and n as c1+c2 and c1-c2 */ cv1 = -whalf*l - z*f_h; cv2 = x*l; ch1 = -hhalf*l - z*f_v; ch2 = y*l; /* Compute intersection times. */ t_min = (l - z)/vz; /* printf(" (x,y,z)=(%g %g %g) (vx,vy,vz)=(%g %g %g) Exit time : %g \n", x,y,z,vx,vy,vz,t_min); */ i = 0; if(vdotn_v1 < 0 && (t_tmp = (cv1 + cv2)/vdotn_v1) < t_min) { t_min = t_tmp; i = 1; /* printf("Left vertical: t=%g \n",t_min); */ } if(vdotn_v2 < 0 && (t_tmp = (cv1 - cv2)/vdotn_v2) < t_min) { t_min = t_tmp; i = 2; /* printf("Right vertical: t=%g \n",t_min); */ } if(vdotn_h1 < 0 && (t_tmp = (ch1 + ch2)/vdotn_h1) < t_min) { t_min = t_tmp; i = 3; /* printf("Lower horizontal: t=%g \n",t_min); */ } if(vdotn_h2 < 0 && (t_tmp = (ch1 - ch2)/vdotn_h2) < t_min) { t_min = t_tmp; i = 4; /* printf("Upper horizontal: t=%g \n",t_min); */ } if(i == 0) break; /* Neutron left guide. */ PROP_DT(t_min); /* ******* Recalculate dot products ********* */ d_z=DEG2RAD*eta_z*randnorm(); d_xy=DEG2RAD*eta_xy*randnorm(); /* Now the normal vector is rotated. To 1st order in waviness and 2nd order in f=(w2-w1)/l and (f times waviness) the rotation matrix for the left vertical mirror is: { 1 d_xy d_z } { -d_xy 1 d_xy f_h } (for the right vertical mirror the ) { d_z -d_xy f_h 1 } (terms d_xy f_h changes sign ) the left vertical normal vector is { -l, 0, f_h } giving the rotated normal vector { -l + f_h d_z, l d_xy, f_h - l d_z } for the right vertical mirror the normal vector is { l, 0, f_h } and the rotated right normal vector is { l + f_h d_z, -l d_xy, f_h + l d_z } The top horizontal mirror must be (something like) { 1 -d_xy d_xy f_h } { d_xy 1 d_z } (for the bottom horizontal mirror the ) { -d_xy f_h d_z 1 } (terms d_xy f_h changes sign ) The top horizontal normal vector is { 0, -l, f_v} giving the rotated normal vector { l d_xy, -l + f_v d_z, f_v - l d_z } for the bottom mirror the normal vector is { 0, l, f_h } and the rotated bottom normal vector is { -l d_xy, l + f_v d_z, f_v + l d_z } */ switch(i) { case 1: /* Left vertical mirror */ m0=m1; W=W1; alpha0=alpha1; norm_n2 = (-l+d_z*f_h)*(-l+d_z*f_h)+(d_xy*l)*(d_xy*l) +(f_h-d_z*l)*(f_h-d_z*l); /* Square of length of n vector */ vdotn = (vx*(-l+f_h*d_z)+ vy*(l*d_xy)+ vz*(f_h-l*d_z) ); q = 2 * V2Q * fabs(vdotn) / sqrt(norm_n2); dvx = -2*(-l+f_h*d_z)*vdotn/norm_n2; dvy = -2*(l*d_xy)*vdotn/norm_n2; dvz = -2*(f_h-l*d_z)*vdotn/norm_n2; break; case 2: /* Right vertical mirror */ m0=m2; W=W2; alpha0=alpha2; norm_n2 = (l+d_z*f_h)*(l+d_z*f_h)+(-d_xy*l)*(-d_xy*l) +(f_h+d_z*l)*(f_h+d_z*l); /* Square of length of n vector */ vdotn = (vx*(l+f_h*d_z)+ vy*(-l*d_xy)+ vz*(f_h+l*d_z) ); q = 2 * V2Q * fabs(vdotn) / sqrt(norm_n2); dvx = -2*(l+f_h*d_z)*vdotn/norm_n2; dvy = -2*(-l*d_xy)*vdotn/norm_n2; dvz = -2*(f_h+l*d_z)*vdotn/norm_n2; break; case 3: /* Lower horizontal mirror */ m0=m3; W=W3; alpha0=alpha3; norm_n2 = (d_xy*l)*(d_xy*l)+(-l+d_z*f_v)*(-l+d_z*f_v) +(f_v-d_z*l)*(f_v-d_z*l); /* Square of length of n vector */ vdotn = (vx*(l*d_xy)+ vy*(-l+f_v*d_z)+ vz*(f_v-l*d_z) ); q = 2 * V2Q * fabs(vdotn) / sqrt(norm_n2); dvx = -2*(l*d_xy)*vdotn/norm_n2; dvy = -2*(-l+f_v*d_z)*vdotn/norm_n2; dvz = -2*(f_v-l*d_z)*vdotn/norm_n2; break; case 4: /* Upper horizontal mirror */ m0=m4; W=W4; alpha0=alpha4; norm_n2 = (-d_xy*l)*(-d_xy*l)+(l+d_z*f_v)*(l+d_z*f_v) +(f_v+d_z*l)*(f_v+d_z*l); /* Square of length of n vector */ vdotn = (vx*(-l*d_xy)+ vy*(l+f_v*d_z)+ vz*(f_v+l*d_z) ); q = 2 * V2Q * fabs(vdotn) / sqrt(norm_n2); dvx = -2*(-l*d_xy)*vdotn/norm_n2; dvy = -2*(l+f_v*d_z)*vdotn/norm_n2; dvz = -2*(f_v+l*d_z)*vdotn/norm_n2; break; default: printf("Fatal error: No guide wall hit"); exit(1); } /* Now compute reflectivity. */ { double par[] = {R0, Qc, alpha0, m0, W}; StdReflecFunc(q, par, &R); if (R > 0) p *= R; else ABSORB; } vx += dvx; vy += dvy; vz += dvz; SCATTER; } %} MCDISPLAY %{ double x; int i; multiline(5, -w1/2.0, -h1/2.0, 0.0, w1/2.0, -h1/2.0, 0.0, w1/2.0, h1/2.0, 0.0, -w1/2.0, h1/2.0, 0.0, -w1/2.0, -h1/2.0, 0.0); multiline(5, -w2/2.0, -h2/2.0, (double)l, w2/2.0, -h2/2.0, (double)l, w2/2.0, h2/2.0, (double)l, -w2/2.0, h2/2.0, (double)l, -w2/2.0, -h2/2.0, (double)l); line(-w1/2.0, -h1/2.0, 0, -w2/2.0, -h2/2.0, (double)l); line( w1/2.0, -h1/2.0, 0, w2/2.0, -h2/2.0, (double)l); line( w1/2.0, h1/2.0, 0, w2/2.0, h2/2.0, (double)l); line(-w1/2.0, h1/2.0, 0, -w2/2.0, h2/2.0, (double)l); %} END