/******************************************************************************* * * 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 * Version: $Revision: 1.11 $ * Origin: Risoe * Release: McStas 1.6 * * 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. * * %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 * alpha1: (AA) Slope of reflectivity * m1: (1) m-value of material, left * W1: (AA-1) Width of supermirror cut-off * alpha2: (AA) Slope of reflectivity * m2: (1) m-value of material, right. * W2: (AA-1) Width of supermirror cut-off * alpha3: (AA) Slope of reflectivity * m3: (1) m-value of material, top. * W3: (AA-1) Width of supermirror cut-off * alpha4: (AA) Slope of reflectivity * m4: (1) m-value of material, bottom. * W4: (AA-1) Width of supermirror cut-off * wavy_z: (deg) Waviness in the z-(flight-)direction * wavy_xy: (deg) Waviness in the transverse direction * * %D * Example values: m=2 Qc=0.0218 (nat. Ni) W=1/300 alpha=4.38 R0=0.995 * (given by Daniel Clemens, PSI) * %E ******************************************************************************/ DEFINE COMPONENT Guide_wavy DEFINITION PARAMETERS () SETTING PARAMETERS (w1, h1, w2, h2, l, R0, Qc, alpha1, m1, W1, alpha2, m2, W2, alpha3, m3, W3, alpha4, m4, W4, wavy_z, wavy_xy) OUTPUT PARAMETERS (whalf,hhalf, lwhalf, lhhalf, norm_nv, norm_nh, f_h, f_v, eta_z, eta_xy) STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p) DECLARE %{ double whalf,hhalf, lwhalf, lhhalf, norm_nv, norm_nh, f_h, f_v, eta_z, eta_xy; %} INITIALIZE %{ 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 m, alpha, 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; 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 */ m=m1; W=W1; alpha=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 */ m=m2; W=W2; alpha=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 */ m=m3; W=W3; alpha=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 */ m=m4; W=W4; alpha=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. */ if(m == 0) ABSORB; if(q > Qc) { double arg = (q-m*Qc)/W; if(arg < 10) { R = .5*(1-tanh(arg))*(1-alpha*(q-Qc)); p *= R; /* printf("q,Qc,alpha,arg,R: %g %g %g %g %g \n",q,Qc,alpha,arg,R); */ } else ABSORB; /* Cutoff ~ 1E-10 */ } p *= R0; /* Correct for max. reflectivity for both cases */ vx += dvx; vy += dvy; vz += dvz; SCATTER; } %} MCDISPLAY %{ double x; int i; magnify("xy"); 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