/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright (C) 1997-2011, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: Virtual_mcnp_ss_Guide * * %I * Written by: Esben klinkby and Peter Willendrup * Date: Marts 2012 * Origin: Risoe-DTU * * Neutron guide initiated using Virtual_mcnp_ss_input.comp, and replacing Virtual_mcnp_ss_output.comp - see examples//Test_SSR_SSW_Guide.instr * * %D * Based on Kristian Nielsens Guide.comp * Models a rectangular guide tube centered on the Z axis. The entrance lies * in the X-Y plane. * The component must be initiated after Virtual_mcnp_ss_input.comp, * and replaces Virtual_mcnp_ss_output.comp - see examples/Test_SSR_SSW_Guide.instr. * The basic idea is, that rather than discarding unreflected (i.e. absorbed) * neutrons at the guide mirrors, these neutron states are stored on disk. * Thus, after the McStas simulation a MCNP simulation can be performed based * on the un-reflected neutrons - intended for shielding studies. * (details: we don't deal with actual neutrons, so what is transferred between * simulations suites is neutron state parameters: pos,mom,time,weight. The latter is * whatever remains after reflection. * * 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: Virtual_mcnp_ss_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 * Upcomming in 2012 based on ESS shielding * Validated by: D. Ene & E. Klinkby * * %BUGS * This component does not work with gravitation on. Use component Guide_gravity then * (doesn't work with SSR/SSW unfortunately) * * %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 * reflect: [str] Reflectivity file name. Format * * %D * Example values: m=4 Qc=0.0219 W=1/300 alpha=6.49 R0=1 * * %E *******************************************************************************/ DEFINE COMPONENT Virtual_mcnp_ss_Guide SETTING PARAMETERS (string reflect=0, w1, h1, w2, h2, l, R0=0.99, Qc=0.0219, alpha=6.07, m=2, W=0.003) NOACC SHARE %{ %include "read_table-lib" %} DECLARE %{ t_Table pTable; double from_mcstas[8]; unsigned int ntrk; unsigned int nhis; %} INITIALIZE %{ writeheader_(&ntrk,&nhis); //Important that ntrk & nhis does not exceed actuall rssa file content if (mcgravitation) fprintf(stderr,"WARNING: Virtual_mcnp_ss_Guide: %s: " "This component produces wrong results with gravitation !\n" "Use Guide_gravity (doesn't work with SSR/SSW).\n", NAME_CURRENT_COMP); if (reflect && strlen(reflect)) { if (Table_Read(&pTable, reflect, 1) <= 0) /* read 1st block data from file into pTable */ exit(fprintf(stderr,"Virtual_mcnp_ss_Guide: %s: can not read file %s\n", NAME_CURRENT_COMP, reflect)); } else { if (W < 0 || R0 < 0 || Qc < 0 || m < 0) { fprintf(stderr,"Virtual_mcnp_ss_Guide: %s: W R0 Qc must be >0.\n", NAME_CURRENT_COMP); exit(-1); } if (m < 1 && m != 0) fprintf(stderr,"WARNING: Virtual_mcnp_ss_Guide: %s: m < 1 behaves as if m=1.\n", NAME_CURRENT_COMP); } %} 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; double vx_org=0.,vy_org=0.,vz_org=0.; /* 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(;;) { vx_org=vx; vy_org=vy; vz_org=vz; /* 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)) weight = Table_Value(pTable, q, 1); else if(q > Qc) { double arg = (q-m*Qc)/W; // printf(" %e %e %e %e %e \n ", arg, q, m, Qc, W ); if(arg < 10) weight = .5*(1-tanh(arg))*(1-alpha*(q-Qc)); else ABSORB; /* Cutoff ~ 1E-10 */ weight *= R0; } else { /* q <= Qc */ weight *= R0; } p *= weight; SCATTER; from_mcstas[0]=x; from_mcstas[1]=y; from_mcstas[2]=z; from_mcstas[3]=vx_org; from_mcstas[4]=vy_org; from_mcstas[5]=vz_org; from_mcstas[6]=p*(1.-weight); //Hmm... assuming that weigths has the same meening in MCNP and McStas. from_mcstas[7]=t; writeneutron_(&from_mcstas); // double v2=vx_org*vx_org+vy_org*vy_org+vz_org*vz_org; // printf("%e %e %e %e \n ", weight, 1-weight, from_mcstas[6], v2); } %} 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