/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright 1997-2002, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: Guide * * %I * Written by: Kristian Nielsen * Date: September 2 1998 * Origin: Risoe * * Neutron guide. * * %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. * The reflectivity profile may either use an analytical mode (see Component * Manual) or a 2-columns reflectivity free text file with format * [q(Angs-1) R(0-1)]. * * Example: Guide(w1=0.1, h1=0.1, w2=0.1, h2=0.1, l=2.0, R0=0.99, Qc=0.021, alpha=6.07, m=2, W=0.003 * * %VALIDATION * May 2005: extensive internal test, no bugs found * Validated by: K. Lieutenant * * %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. glass/SiO2 Si Ni Ni58 supermirror Be Diamond m= 0.65 0.47 1 1.18 2-6 1.01 1.12 * W: [AA-1] Width of supermirror cut-off * * %D * Example values: m=4 Qc=0.0219 W=1/300 alpha=6.49 R0=1 * * %E *******************************************************************************/ DEFINE COMPONENT Guide_simple DEFINITION PARAMETERS () SETTING PARAMETERS (w1, h1, w2=0, h2=0, l, R0=0.99, Qc=0.0219, alpha=6.07, m=2, W=0.003) OUTPUT PARAMETERS () /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE %{ //%include "read_table-lib" %include "ref-lib" %} DECLARE %{ //t_Table pTable; %} INITIALIZE %{ if (mcgravitation) fprintf(stderr,"WARNING: Guide: %s: " "This component produces wrong results with gravitation !\n" "Use Guide_gravity.\n", NAME_CURRENT_COMP); if (!w2) w2=w1; if (!h2) h2=h1; // if (reflect && strlen(reflect) && strcmp(reflect,"NULL") && strcmp(reflect,"0")) { // if (Table_Read(&pTable, reflect, 1) <= 0) /* read 1st block data from file into pTable */ // exit(fprintf(stderr,"Guide: %s: can not read file %s\n", NAME_CURRENT_COMP, reflect)); //} else { if (W < 0 || R0 < 0 || Qc < 0 || m < 0) { fprintf(stderr,"Guide: %s: W R0 Qc must be >0.\n", NAME_CURRENT_COMP); exit(-1); } //} %} TRACE %{ double t1,t2; /* Intersection times. */ double av,ah,bv,bh,cv1,cv2,ch1,ch2,d; /* Intermediate values */ double weight; /* Internal probability weight */ double vdotn_v1,vdotn_v2,vdotn_h1,vdotn_h2; /* Dot products. */ int i; /* Which mirror hit? */ double q; /* Q [1/AA] of reflection */ double nlen2; /* Vector lengths squared */ /* ToDo: These could be precalculated. */ double ww = .5*(w2 - w1), hh = .5*(h2 - h1); double whalf = .5*w1, hhalf = .5*h1; /* 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 = ww*vz; ah = l*vy; bh = hh*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*ww; cv2 = x*l; ch1 = -hhalf*l - z*hh; ch2 = y*l; /* Compute intersection times. */ t1 = (l - z)/vz; i = 0; if(vdotn_v1 < 0 && (t2 = (cv1 - cv2)/vdotn_v1) < t1) { t1 = t2; i = 1; } if(vdotn_v2 < 0 && (t2 = (cv1 + cv2)/vdotn_v2) < t1) { t1 = t2; i = 2; } if(vdotn_h1 < 0 && (t2 = (ch1 - ch2)/vdotn_h1) < t1) { t1 = t2; i = 3; } if(vdotn_h2 < 0 && (t2 = (ch1 + ch2)/vdotn_h2) < t1) { t1 = t2; i = 4; } if(i == 0) break; /* Neutron left guide. */ PROP_DT(t1); switch(i) { case 1: /* Left vertical mirror */ nlen2 = l*l + ww*ww; q = V2Q*(-2)*vdotn_v1/sqrt(nlen2); d = 2*vdotn_v1/nlen2; vx = vx - d*l; vz = vz - d*ww; break; case 2: /* Right vertical mirror */ nlen2 = l*l + ww*ww; q = V2Q*(-2)*vdotn_v2/sqrt(nlen2); d = 2*vdotn_v2/nlen2; vx = vx + d*l; vz = vz - d*ww; break; case 3: /* Lower horizontal mirror */ nlen2 = l*l + hh*hh; q = V2Q*(-2)*vdotn_h1/sqrt(nlen2); d = 2*vdotn_h1/nlen2; vy = vy - d*l; vz = vz - d*hh; break; case 4: /* Upper horizontal mirror */ nlen2 = l*l + hh*hh; q = V2Q*(-2)*vdotn_h2/sqrt(nlen2); d = 2*vdotn_h2/nlen2; vy = vy + d*l; vz = vz - d*hh; break; } /* Now compute reflectivity. */ weight = 1.0; /* Initial internal weight factor */ if(m == 0) ABSORB; // if (reflect && strlen(reflect) && strcmp(reflect,"NULL") && strcmp(reflect,"0")) // TableReflecFunc(q, &pTable, &weight); //else { double par[] = {R0, Qc, alpha, m, W}; StdReflecFunc(q, par, &weight); //} if (weight > 0) p *= weight; else ABSORB; SCATTER; } %} MCDISPLAY %{ 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