/******************************************************************************* * * McStas, neutron ray-tracing pacxkage * Copyright 1997-2002, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: Guide_honeycomb * * %I * Written by: G. Venturi * Date: Apr 2003 * Origin: ILL (France). * Modified by: G. Venturi, Made Honeycomb from Guide_gravity (Apr 2003) * Modified by: E. Farhi, corrected gravitation bug (Feb 2005) * Modified by: E. Farhi, uniformize parameter names (Jul 2008) * * Neutron guide with gravity and honeycomb geometry. Can be channeled and focusing. * * %D * Models a honeycomb guide tube centered on the Z axis. The entrance lies * in the X-Y plane. Gravitation applies also when reaching the guide input * window. The guide can be channeled (hexagonal channel section; nslit,d parameters). * The guide coating specifications may be entered via different ways (global, * or for each wall m-value). * For details on the geometry calculation see the description in the McStas * reference manual. * * Example: Guide_honeycomb(w1=0.1, w2=0.1, l=12, * R0=0.99, Qc=0.0219, alpha=6.07, m=1.0, W=0.003, nslit=1, d=0.0005) * * %P * INPUT PARAMETERS: * * w1: [m] Width at the guide entry * w2: [m] Width at the guide exit. If zero, sets w2=w1. * 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 * d: [m] Thickness of subdividing walls [0] * nslit: [1] Number of horizontal channels in the guide (>= 1) [1] * (nslit-1 vertical dividing walls) * * Optional input parameters: (different ways for m-specifications) * * mright: [1] m-value of material for right. vertical mirror * mleft: [1] m-value of material for left. vertical mirror * mleftup: [1] m-value of material for leftup. oblique mirror * mrightup: [1] m-value of material for rightdown. oblique mirror * mrightdown: [1] m-value of material for leftup. oblique mirror * mleftdown: [1] m-value of material for rightdown. oblique mirror * G: [m/s2] Gravitation norm. 0 value disables G effects. * * OUTPUT PARAMETERS * * GVars: [1] internal variables * GVars.N_reflection: (1) Array of the cumulated Number of reflections * N_reflection[0] total nb of reflections * N_reflection[1,2,3,4.5.6] reflections * N_reflection[7] total nb neutrons exiting guide * N_reflection[8] total nb neutrons entering guide * * %D * Example values: m=4 Qc=0.02 W=1/300 alpha=6.49 R0=1 * * %E *******************************************************************************/ DEFINE COMPONENT Guide_honeycomb DEFINITION PARAMETERS () SETTING PARAMETERS (w1, w2=0, l, R0=0.995, Qc=0.0218, alpha=4.38, m=1.0, W=0.003, nslit=1, d=0.0005, mleftup=-1, mrightup=-1, mleftdown=-1, mrightdown=-1,mleft=-1, mright=-1, G=0) OUTPUT PARAMETERS () /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE %{ %include "ref-lib" #ifndef Honeycomb_guide_Version #define Honeycomb_guide_Version "$Revision$" #ifndef PROP_GRAV_DT #error McStas : You need PROP_GRAV_DT (McStas >= 1.4.3) to run this component #endif /* * G: (m/s^2) Gravitation acceleration along y axis [-9.81] * Gx: (m/s^2) Gravitation acceleration along x axis [0] * Gy: (m/s^2) Gravitation acceleration along y axis [-9.81] * Gz: (m/s^2) Gravitation acceleration along z axis [0] * mh: (1) m-value of material for left/right vert. mirrors * mv: (1) m-value of material for top/bottom horz. mirrors * mx: (1) m-value of material for left/right vert. mirrors * my: (1) m-value of material for top/bottom horz. mirrors */ typedef struct Honeycomb_guide_Vars { double gx; double gy; double gz; double nx[8], ny[8], nz[8]; double wx[8], wy[8], wz[8]; double A[8], norm_n2[8], norm_n[8]; long N_reflection[9]; double w1c, w2c; double M[7]; double nzC[7], norm_n2xy[7], Axy[7]; char compcurname[256]; double warnings; } Honeycomb_guide_Vars_type; void Honeycomb_guide_Init(Honeycomb_guide_Vars_type *aVars, MCNUM a_w1, MCNUM a_w2, MCNUM a_l, MCNUM a_R0, MCNUM a_Qc, MCNUM a_alpha, MCNUM a_m, MCNUM a_W, MCNUM a_nslit, MCNUM a_d, MCNUM a_Gx, MCNUM a_Gy, MCNUM a_Gz, MCNUM a_mright, MCNUM a_mleft, MCNUM a_mleftup, MCNUM a_mrightdown, MCNUM a_mrightup, MCNUM a_mleftdown) { int i; for (i=0; i<8; aVars->N_reflection[i++] = 0); for (i=0; i<7; aVars->M[i++] = 0); aVars->gx = a_Gx; /* The gravitation vector in the current component axis system */ aVars->gy = a_Gy; aVars->gz = a_Gz; aVars->warnings=0; if (a_nslit <= 0) { fprintf(stderr,"%s: Fatal: no channel in this guide (kh or nslit=0).\n", aVars->compcurname); exit(-1); } if (a_d < 0) { fprintf(stderr,"%s: Fatal: subdividing walls have negative thickness in this guide (d<0).\n", aVars->compcurname); exit(-1); } aVars->w1c = 0.5*(a_w1 - (a_nslit-1) *2*a_d)/(double)a_nslit; /*INPUT APOTHEM*/ aVars->w2c = 0.5*(a_w2 - (a_nslit-1) *2*a_d)/(double)a_nslit; /*OUTPUT APOTHEM*/ for (i=0; i <= 6; aVars->M[i++]=a_m); if (a_mright >= 0) aVars->M[1] =a_mright; if (a_mleft >= 0) aVars->M[2] =a_mleft; if (a_mleftup >= 0) aVars->M[3] =a_mleftup; if (a_mrightdown >= 0) aVars->M[4] =a_mrightdown; if (a_mrightup >= 0) aVars->M[5] =a_mrightup; if (a_mleftdown >= 0) aVars->M[6] =a_mleftdown; /* n: normal vectors to surfaces */ aVars->nx[1] = -a_l; aVars->ny[1] = 0; aVars->nz[1] = (aVars->w2c-aVars->w1c); /* 1:+X right */ aVars->nx[2] = +a_l; aVars->ny[2] = 0; aVars->nz[2] = -aVars->nz[1]; /* 2:-X left */ aVars->nx[3] = +a_l*0.5; aVars->ny[3] = -0.866025*a_l; aVars->nz[3] = (aVars->w2c-aVars->w1c); /* 3:+Y leftup*/ aVars->nx[4] = -a_l*0.5; aVars->ny[4] = +0.866025*a_l; aVars->nz[4] = -(aVars->w2c-aVars->w1c); /* 4:+Y rightdown*/ aVars->nx[5] = -a_l*0.5; aVars->ny[5] = -0.866025*a_l; aVars->nz[5] = (aVars->w2c-aVars->w1c); /* 5:+Y rightup */ aVars->nx[6] = +a_l*0.5; aVars->ny[6] = +0.866025*a_l; aVars->nz[6] = -(aVars->w2c-aVars->w1c); /* 6:+Y leftdown */ aVars->nx[7] = 0; aVars->ny[7] = 0; aVars->nz[7] = a_l; aVars->nx[0] = 0; aVars->ny[0] = 0; aVars->nz[0] = -a_l; /* w: a point on these surfaces */ aVars->wx[1] = (aVars->w1c); aVars->wy[1] = 0; aVars->wz[1] = 0; /* 1: right */ aVars->wx[2] = -(aVars->w1c); aVars->wy[2] = 0; aVars->wz[2] = 0; /* 2: left */ aVars->wx[3] = -0.5*(aVars->w1c); aVars->wy[3] = +0.866025*(aVars->w1c); aVars->wz[3] = 0; /* 3: leftup */ aVars->wx[4] = +0.5*(aVars->w1c); aVars->wy[4] = -0.866025*(aVars->w1c); aVars->wz[4] = 0; /* 4: rightdown */ aVars->wx[5] = +0.5*(aVars->w1c); aVars->wy[5] = +0.866025*(aVars->w1c); aVars->wz[5] = 0; /* 5: rightup */ aVars->wx[6] = -0.5*(aVars->w1c); aVars->wy[6] = -0.866025*(aVars->w1c); aVars->wz[6] = 0; /* 6: leftdown */ aVars->wx[7] = 0; aVars->wy[7] = 0; aVars->wz[7] = a_l; /* 7:+Z exit */ aVars->wx[0] = 0; aVars->wy[0] = 0; aVars->wz[0] = 0; /* 0:Z0 input */ for (i=0; i <= 7; i++) { aVars->A[i] = scalar_prod(aVars->nx[i], aVars->ny[i], aVars->nz[i], aVars->gx, aVars->gy, aVars->gz)/2; aVars->norm_n2[i] = aVars->nx[i]*aVars->nx[i] + aVars->ny[i]*aVars->ny[i] + aVars->nz[i]*aVars->nz[i]; if (aVars->norm_n2[i] <= 0) { fprintf(stderr,"%s: Fatal: normal vector norm %i is null/negative ! check guide dimensions.\n", aVars->compcurname, i); exit(-1); } /* should never occur */ else aVars->norm_n[i] = sqrt(aVars->norm_n2[i]); } /* partial computations for sides, to save computing time */ for (i=1; i <= 6; i++) { aVars->nzC[i] = aVars->nz[i]; aVars->norm_n2xy[i]= aVars->nx[i]*aVars->nx[i] + aVars->ny[i]*aVars->ny[i]; aVars->Axy[i] =(aVars->nx[i]*aVars->gx + aVars->ny[i]*aVars->gy)/2; } } #pragma acc routine seq int Honeycomb_guide_Trace(double *dt, Honeycomb_guide_Vars_type *aVars, double cx, double cy, double cz, double cvx, double cvy, double cvz, double cxnum, int nslit, double cynum) { double B, C; int ret=0; int side=0; double n1; double dt0, dt_min=0; int i; double loc_num; int i_slope=3; /* look if there is a previous intersection with guide sides */ /* A = 0.5 n.g; B = n.v; C = n.(r-W); */ /* 5=+Z side: n=(0, 0, -l) ; W = (0, 0, l) (at z=l, guide exit)*/ B = aVars->nz[7]*cvz; C = aVars->nz[7]*(cz - aVars->wz[7]); ret = solve_2nd_order(&dt0, NULL, aVars->A[7], B, C); if (ret && dt0>10e-10) { dt_min = dt0; side=7; } loc_num = (3*cynum+cxnum)/2; for (i=6; i>0; i--) { if (i == 4) { i_slope=1; loc_num =(3*cynum-cxnum)/2; } if (i == 2) { i_slope=1; loc_num = cxnum;} if (aVars->nzC[i_slope] != 0) { n1=loc_num-1; loc_num++; loc_num++; aVars->nz[i] = aVars->nzC[i]*n1; aVars->A[i] = aVars->Axy[i] + aVars->nz[i]*aVars->gz/2; } B = aVars->nx[i]*cvx + aVars->ny[i]*cvy + aVars->nz[i]*cvz; /* n.v */ C = aVars->nx[i]*(cx-aVars->wx[i]) + aVars->ny[i]*(cy-aVars->wy[i]) + aVars->nz[i]*cz; /* n.(r-W) */ ret = solve_2nd_order(&dt0, NULL, aVars->A[i], B, C); if (ret && dt0>10e-10 && (dt0nzC[i] != 0) { aVars->norm_n2[i] = aVars->norm_n2xy[i] + aVars->nz[i]*aVars->nz[i]; aVars->norm_n[i] = sqrt(aVars->norm_n2[i]); } } } *dt = dt_min; return (side); } #endif %} DECLARE %{ Honeycomb_guide_Vars_type GVars; //#prama acc routine declare create(GVars) %} INITIALIZE %{ double Gx=0, Gy=9.81, Gz=0; Coords mcLocG; int i; if (!w2) w2=w1; if (W < 0 || nslit <= 0 || R0 < 0 || Qc < 0) { fprintf(stderr,"%s:Guide_gravity: W nslit R0 Qc must be >0.\n", NAME_CURRENT_COMP); exit(-1); } if (mcgravitation) G=-9.81; mcLocG = rot_apply(ROT_A_CURRENT_COMP, coords_set(Gx,G,Gz)); coords_get(mcLocG, &Gx, &Gy, &Gz); strcpy(GVars.compcurname, NAME_CURRENT_COMP); Honeycomb_guide_Init(&GVars, w1, w2, l, R0, Qc, alpha, m, W, nslit, d, Gx, Gy, Gz,mright, mleft, mleftup, mrightdown, mrightup, mleftdown); if (!G) for (i=0; i<7; GVars.A[i++] = 0); //#pragma acc update device(GVars) %} TRACE %{ double B, C, dt; int ret, bounces = 0; float n,m1,nv,mv; int nup=-1, mright1=-1; double xhole,yhole, y_min; double cn; double inside; double w_edge, w_adj; /* Channel displacement */ double h_edge, h_adj; double w_chnum,h_chnum,w_c,h_c; /* channel indexes */ /* propagate to box input (with gravitation) in comp local coords */ /* 0=Z0 side: n=(0, 0, l) ; W = (0, 0, 0) (at z=0, guide input)*/ B = -l*vz; C = -l*z; ret = solve_2nd_order(&dt, NULL, GVars.A[0], B, C); if (ret && dt>0) { PROP_GRAV_DT(dt, GVars.gx, GVars.gy, GVars.gz); GVars.N_reflection[8]++; } /* check if we are in the box input, else absorb */ if(dt > 0 && fabs(x) <= w1/2 && fabs(y) <= w1/2) { for(n=0; x>((2*n+1)*(GVars.w1c+d)) ; n++); nv=n; if (n==0) { for(n=0; x<((2*n-1)*(GVars.w1c+d)); n--); nv=n; } xhole=2*n*(GVars.w1c+d); if(x-xhole>0) nup=1; y_min=1.732*(GVars.w1c+d); for(m1=0;y>((2*m1+1)*y_min); m1++); mv=m1; if (m1==0) { for (m1=0; y<((2*m1-1)*(y_min)); m1--); mv=m1; } yhole=2*m1*1.732*(GVars.w1c+d); if(y-yhole>0) mright1=1; cn=1.1547*(GVars.w1c+d); inside=(fabs(y-yhole)+0.577351*fabs(x-xhole)-cn); if(inside>0) { xhole +=nup*(GVars.w1c+d); yhole +=mright1*(1.732*(GVars.w1c+d)); } /* double w_chnum,h_chnum;*/ /* channel indexes */ /* Shift origin to center of channel hit (absorb if hit dividing walls) */ w_c=xhole/(GVars.w1c+d); h_c=yhole/(1.732*(GVars.w1c+d)); w_chnum=rint(w_c); h_chnum=rint(h_c); x -=xhole; y -=yhole; w_adj = xhole; h_adj = yhole; if(fabs(x) > GVars.w1c) { x += xhole; /* Re-adjust origin */ y += yhole; ABSORB; } if(fabs(x*0.5+y*0.866025) > GVars.w1c) { x += xhole; /* Re-adjust origin */ y += yhole; ABSORB; } if(fabs(-x*0.5+y*0.866025) > GVars.w1c) { x += xhole; /* Re-adjust origin */ y += yhole; ABSORB; } /* neutron is now in the input window of the guide */ /* do loops on reflections in the box */ for(;;) { /* get intersections for all box sides */ double q; int side=0; bounces++; /* now look for intersection with guide sides and exit */ side = Honeycomb_guide_Trace(&dt, &GVars, x, y, z, vx, vy, vz, w_chnum, nslit, h_chnum); /* only positive dt are valid */ /* exit reflection loops if no intersection (neutron is after box) */ if (side == 0 || dt <= 0) { if (GVars.warnings < 100) fprintf(stderr,"%s: warning: neutron has entered guide, but can not exit !\n", GVars.compcurname); GVars.warnings++; x += w_adj; y += h_adj; ABSORB; } /* should never occur */ /* propagate to dt */ PROP_GRAV_DT(dt, GVars.gx, GVars.gy, GVars.gz); /* do reflection on speed for l/r/u/d sides */ if (side == 7) /* neutron reaches end of guide: end loop and exit comp */ { GVars.N_reflection[side]++; x += w_adj; y += h_adj; SCATTER; x -= w_adj; y -= h_adj; break; } /* else reflection on a guide wall */ if(GVars.M[side] == 0 || Qc == 0) /* walls are absorbing */ { x += w_adj; y += h_adj; ABSORB; } /* change/mirror velocity: h_f = v - n.2*n.v/|n|^2 */ B = GVars.nx[side]*vx + GVars.ny[side]*vy + GVars.nz[side]*vz; /* n.v */ GVars.N_reflection[side]++; /* GVars.norm_n2 > 0 was checked at INIT */ dt = 2*B/GVars.norm_n2[side]; /* 2*n.v/|n|^2 */ vx -= GVars.nx[side]*dt; vy -= GVars.ny[side]*dt; vz -= GVars.nz[side]*dt; /* compute q and modify neutron weight */ /* scattering q=|nslit_i-nslit_f| = V2Q*|vf - v| = V2Q*2*n.v/|n| */ q = 2*V2Q*fabs(B)/GVars.norm_n[side]; { double par[] = {R0, Qc, alpha, GVars.M[side], W}; StdReflecFunc(q, par, &B); } if (B <= 0) { x += w_adj; y += h_adj; ABSORB; } else p *= B; x += w_adj; y += h_adj; SCATTER; x -= w_adj; y -= h_adj; GVars.N_reflection[0]++; /* go to the next reflection */ if (bounces > 1000) ABSORB; } /* end for */ x += w_adj; y += h_adj; /* Re-adjust origin after SCATTER */ } else ABSORB; %} FINALLY %{ if (GVars.warnings > 100) { fprintf(stderr,"%s: warning: neutron has entered guide, but can not exit !\n", GVars.compcurname); fprintf(stderr,"%s: warning: This message has been repeated %g times\n", GVars.compcurname, GVars.warnings); } %} MCDISPLAY %{ int i,j; double a,b,c; double x0,x01,x1,x11,x2,x21; double y0,y01,y1,y11,y2,y21,y3,y31,y4,y41; for(j = -nslit/2; j <= nslit/2; j++) { y0 = j*(GVars.w1c+d)*1.732; y01= j*(GVars.w2c+d)*1.732; y1 =y0 +(GVars.w1c+d)/1.732; y11=y01+(GVars.w2c+d)/1.732; y2 =y0 -(GVars.w1c+d)/1.732; y21=y01-(GVars.w2c+d)/1.732; y3 =y0 +(GVars.w1c+d)*1.1547; y31=y01+(GVars.w2c+d)*1.1547; y4 =y0 -(GVars.w1c+d)*1.1547; y41=y01-(GVars.w2c+d)*1.1547; for(i = -nslit; i <= nslit; i++) { a=i+j; b=a/2+0.1; c=rint(b); if(fabs(c-b)<0.3) { x0 = i*(GVars.w1c+d); x01= i*(GVars.w2c+d); x1 =x0 +(GVars.w1c+d); x11=x01+(GVars.w2c+d); x2 =x0 -(GVars.w1c+d); x21=x01-(GVars.w2c+d); multiline(5, x1, y1, 0.0, x11, y11, (double)l, x21, y11, (double)l, x2, y1, 0.0, x1, y1, 0.0); multiline(5, x1, y1, 0.0, x11, y11, (double)l, x01, y31, (double)l, x0, y3, 0.0, x1, y1, 0.0); multiline(5, x0, y3, 0.0, x01, y31, (double)l, x21, y11, (double)l, x2, y1, 0.0, x0, y3, 0.0); multiline(5, x2, y1, 0.0, x21, y11, (double)l, x21, y21, (double)l, x2, y2, 0.0, x2, y1, 0.0); multiline(5, x2, y2, 0.0, x21, y21, (double)l, x01, y41, (double)l, x0, y4, 0.0, x2, y2, 0.0); multiline(5, x0, y4, 0.0, x01, y41, (double)l, x11, y21, (double)l, x1, y2, 0.0, x0, y4, 0.0); } } } %} END