[mcstas-users] Best way to simulate a non-symmetric guide
Peter Link
Peter.Link at frm2.tum.de
Tue May 23 13:21:07 CEST 2017
Hi Inigo,
as the geometry is still rather simple, I would use the guide_anyshape
component. You have to define a .off geometry file to describe your
guide shape. I guess the m-values of your different mirrors aren't
identical? The original component doesn't provide that - I have done a
workaround, not fully tested and no warranty that it work in all cases...
Best regards,
Peter
Am 23.05.2017 um 12:46 schrieb Iñigo Herranz:
>
> Dear all. I am Iñigo Herranz, working for the instrument MIRACLES.
> We have to simulate a guide that is non-symetrically in the vertical
> plane, as figure below shows.
> Could anyone suggest us the best way to do it?. We are not sure
> whether there is a component to do that, or you have to define a group
> of mirrors, or with rotations...
> We look forward to hearing from you.
> In behalf of the MIRACLES team,
> thank you very much.
>
>
> Imágenes integradas 2
>
>
>
>
> _______________________________________________
> mcstas-users mailing list
> mcstas-users at mcstas.org
> http://mailman.mcstas.org/cgi-bin/mailman/listinfo/mcstas-users
--
Dr. Peter Link
Leiter Neutronenoptik
Technische Universität München
Heinz Maier-Leibnitz Zentrum (MLZ)
Lichtenbergstr. 1
85747 Garching
Tel.: +49 (0)89 289 14622
Fax: +49 (0) 89 289 11694
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20170523/42919118/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: nonsymetricalguide.png
Type: image/png
Size: 11364 bytes
Desc: not available
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20170523/42919118/attachment.png>
-------------- next part --------------
/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright (C) 1997-2010, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Guide
*
* %I
* Written by: Emmanuel Farhi
* Date: August 4th 2010
* Version: $Revision: 1.0$
* Origin: ILL
* Release: McStas 1.12b
*
* Reflecting surface (guide and mirror) with any shape, defined from an OFF file.
*
* %D
*
* This is a reflecting object component. Its shape is defined from an OFF file,
* given with its file name. The object size is set as given from the file, where
* dimensions should be in meters. The bounding box may be re-scaled by specifying
* xwidth,yheight,zdepth parameters. The object may optionally be centered when
* using 'center=1'.
*
* The reflectivity is specified either from the usual parametric description
* R0,Qc,alpha,W,m, or from a reflectivity file 'reflect' with a 2 column
* format [q(Angs-1) R(0-1)].
*
* The complex OFF/PLY geometry option handles any closed non-convex polyhedra.
* It supports the OFF and NOFF file format but not COFF (colored faces).
* Such files may be generated from XYZ data using:
* qhull < coordinates.xyz Qx Qv Tv o > geomview.off
* or
* powercrust coordinates.xyz
* and viewed with geomview or java -jar jroff.jar (see below).
* The default size of the object depends of the OFF file data, but its
* bounding box may be resized using xwidth,yheight and zdepth.
* PLY geometry files are also supported.
*
* %BUGS
* This component does take into account gravitation accurately.
*
* %P
* INPUT PARAMETERS:
*
* geometry:(str) Name of the OFF/PLY file describing the guide geometry
* xwidth: (m) Redimension the object bounding box on X axis is non-zero
* yheight: (m) Redimension the object bounding box on Y axis is non-zero
* zdepth: (m) Redimension the object bounding box on Z axis is non-zero
* center: (1) When set to 1, the object will be centered w.r.t. the local coordinate frame
* 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
* m0: (1) m-value of material for color index 0 mirror
* m1: (1) m-value of material for color index 1 mirror
* m2: (1) m-value of material for color index 2 mirror
* m3: (1) m-value of material for color index 3 mirror
* PL: this is not at all elegant, should be replaced by some kind of list or table
* reflect: (str) Reflectivity file name. Format [q(Angs-1) R(0-1)]
* transmit:(1) When true, non reflected neutrons are transmitted through
* the surfaces, instead of being absorbed. No material absorption
* is taken into account though
*
* OUTPUT PARAMETERS:
* SCATTERED: number of reflected events
*
* %D
* Example values: m=4 Qc=0.0219 W=1/300 alpha=6.49 R0=1
*
* %L
* <a href="http://www.geomview.org">Geomview and Object File Format (OFF)</a>
* %L
* Java version of Geomview (display only) <a href="http://www.holmes3d.net/graphics/roffview/">jroff.jar</a>
* %L
* <a href="http://qhull.org">qhull</a>
* %L
* <a href="http://www.cs.ucdavis.edu/~amenta/powercrust.html">powercrust</a>
*
* %E
*******************************************************************************/
DEFINE COMPONENT Guide_anyshape_pl
DEFINITION PARAMETERS (string reflect=0, string geometry=0)
SETTING PARAMETERS (xwidth=0, yheight=0, zdepth=0, center=0, transmit=0,
R0=0.99, Qc=0.0219, alpha=3, m=0, W=0.0033,
m0=-1, m1=-1, m2=-1, m3=-1)
OUTPUT PARAMETERS (pTable, offdata)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
%include "read_table-lib"
%include "interoff-lib"
%include "ref-lib"
%}
DECLARE
%{
off_struct offdata;
t_Table pTable;
%}
INITIALIZE
%{
/* initialize OFF object from the file(s) */
if (!off_init( geometry, xwidth, yheight, zdepth, !center, &offdata )) exit(-1);
/* optionally initialize reflectivity curves */
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_anyshape: %s: can not read file %s. Aborting.\n", NAME_CURRENT_COMP, reflect));
else
Table_Rebin(&pTable); /* rebin as evenly, increasing array */
}
%}
TRACE
%{
int intersect=0;
int counter=0;
/* main loop for multiple reflections */
do {
double t0=0, t3=0, dt=0;
unsigned long faceindex0=0, faceindex3=0, fi=0;
Coords n0, n3, n={0,0,0};
/* determine intersections with object; PL: get index of the intersected face */
intersect = off_intersect(&t0, &t3, &n0, &n3, &faceindex0, &faceindex3,
x,y,z, vx, vy, vz, offdata );
/* get the smallest positive */
if (t0 > 0) { dt = t0; n=n0; fi=faceindex0;}
if (intersect > 1 && dt <= 0 && t3 > dt) { dt = t3; n=n3; fi=faceindex3;}
/* exit loop when no intersection forward */
if (dt <= 0 || !intersect) break;
double nx,ny,nz;
coords_get(n, &nx, &ny, &nz);
/* test if the angle is large in case the object has an internal coating */
double n2 = nx*nx+ny*ny+nz*nz;
double n_dot_v = scalar_prod(vx,vy,vz,nx,ny,nz);
double q = 2*fabs(n_dot_v)*V2K/sqrt(n2);
/* propagate neutron to reflection point */
PROP_DT(dt);
/* handle surface intersection */
double R=0;
/* Reflectivity (see component Guide). */
if (reflect && strlen(reflect) && strcmp(reflect,"NULL") && strcmp(reflect,"0"))
TableReflecFunc(q, &pTable, &R);
else {
/* very crude first attempt to code up to four different m-values */
double m_value=m;
if (offdata.facepropsArray[fi]==0) m_value=m0;
if (offdata.facepropsArray[fi]==1) m_value=m1;
if (offdata.facepropsArray[fi]==2) m_value=m2;
if (offdata.facepropsArray[fi]==3) m_value=m3;
double par[] = {R0, Qc, alpha, m_value, W};
StdReflecFunc(q, par, &R);
}
if (R > 1) {
fprintf(stderr,"Guide_anyshape: %s: Warning: Reflectivity R=%g > 1 lowered to R=1.\n", NAME_CURRENT_COMP, R);
R=1;
}
/* now handle either probability when transmit or reflect */
if (R > 0) {
/* when allowing transmission, we should check if indeed we reflect */
if (!transmit || (transmit && rand01() < R)) {
/* reflect velocity: -q -> -2*n*n.v/|n|^2 */
if (!transmit) p *= R;
n_dot_v *= 2/n2;
vx -= nx*n_dot_v;
vy -= ny*n_dot_v;
vz -= nz*n_dot_v;
SCATTER;
} else {
if (transmit) {
p *= (1-R); /* transmitted beam has non reflected weight */
} else ABSORB;
}
} else {
/* R=0: no reflection: absorb or transmit through when allowed */
if (!transmit) ABSORB;
}
/* leave surface by a small amount so that next intersection is not the same one */
PROP_DT(1e-9);
} while (intersect && counter++<CHAR_BUF_LENGTH);
/* end of main loop */
%}
MCDISPLAY
%{
/* display the object */
off_display(offdata);
/* and its bounding box */
/*
mcdis_box(( offdata.xmin+offdata.xmax)/2,
( offdata.ymin+offdata.ymax)/2,
( offdata.zmin+offdata.zmax)/2,
(-offdata.xmin+offdata.xmax),
(-offdata.ymin+offdata.ymax),
(-offdata.zmin+offdata.zmax));
*/
%}
END
-------------- next part --------------
OFF
8 4 0
-0.071 0.0 0
0.0 0.071 0
0.071 0.0 0
0.0 -0.071 0
-0.071 0.0 20.0
0.0 0.071 20.0
0.071 0.0 20.0
0.0 -0.071 20.0
4 0 1 5 4 0
4 1 2 6 5 1
4 2 3 7 6 0
4 3 0 4 7 1
-------------- next part --------------
/*******************************************************************************
* McStas instrument definition URL=http://www.mcstas.org
*
* Instrument: TEST
*
* %Identification
* Written by: Peter Link (peter.link at frm2.tum.de)
* Date: 08.03.2017
* Origin: FRM II
* Release: McStas
* Version: 2.3
* %INSTRUMENT_SITE: FRMII
*
* Instrument short description
*
* %Description
* Simulation to test Guide_anyshape
*
* Example: mcrun test.instr <parameters=values>
*
* %Parameters
* Par1: [unit] Parameter1 description
*
* %Link
* A reference/HTML link for more information
*
* %End
*******************************************************************************/
/* Change name of instrument and input parameters with default values */
DEFINE INSTRUMENT Test()
/* The DECLARE section allows us to declare variables or small */
/* functions in C syntax. These may be used in the whole instrument. */
DECLARE
%{
%}
/* The INITIALIZE section is executed when the simulation starts */
/* (C code). You may use them as component parameter values. */
INITIALIZE
%{
%}
TRACE
COMPONENT origin = Progress_bar()
AT (0,0,0) ABSOLUTE
/*####################################################################*/
/* Source */
/*####################################################################*/
COMPONENT source = Source_gen(
// T1=35.0, I1=9.38e12,T2=547.5,I2=2.23e12,T3=195.4,I3=1.26e13, // "FRM2 cold" from mcstas.org
Lmin = 2, Lmax = 9, // wavelength range can be adapted to optimise stats
xwidth = 0.14, yheight= 0.21, // SR1 neutron emitting area
dist = 1.9, // focus on 1.9m
focus_xw = 0.11, focus_yh = 0.11) // entry of test guide + 1cm overillumination
AT (0,0,0) RELATIVE origin
/*
COMPONENT testguide = Guide_gravity(
R0 = .99, W = 0.003 ,Qc=0.02174, alpha = 6.07, mtop=2, mbottom=2, mleft = 1.2, mright = 1.2,
w1 = 0.1, h1 = 0.1, w2 = 0.1, h2 = 0.1, l = 20.0)
AT (0, 0, 1.9) RELATIVE origin
ROTATED (0,0,45) RELATIVE origin
*/
/* using the original Guide_anyshape component
COMPONENT testany = Guide_anyshape (
geometry= "testguide.off",
R0=0.99,Qc=0.02174,alpha=6.07,m=2.0,W=0.0033)
AT (0,0,1.9) RELATIVE origin
*/
COMPONENT testany = Guide_anyshape_pl (
geometry= "testguide_45deg.off",
R0=0.99,Qc=0.02174,alpha=6.07,W=0.0033,m0=1.2,m1=2.0)
AT (0,0,1.9) RELATIVE origin
/******************************************************************************/
/* Monitor section at the exit */
/******************************************************************************/
COMPONENT det_arm = Arm()
AT (0,0,21.901) RELATIVE origin
COMPONENT PSD_Det1 = PSD_monitor(
yheight = 0.15, xwidth = 0.15,
filename = "PSD_det1",
nx= 150, ny = 150,
restore_neutron = 1)
AT (0,0, 0.0001) RELATIVE PREVIOUS
COMPONENT Hdiv_Det1 = Monitor_nD(xmin=-0.06, xmax=0.06, ymin=-0.060, ymax=0.060,
filename="Hdivx_Det1.dat",
options= "x bins=120 limits[-0.06 0.06] ; hdiv bins=40 limits[-2.0 2.0]")
AT (0,0,0.0001) RELATIVE PREVIOUS
ROTATED (0,0,0) RELATIVE PREVIOUS
COMPONENT Vdiv_Det1 = Monitor_nD(xmin=-0.06, xmax=0.06, ymin=-0.06, ymax=0.06,
filename="Vdivy_Det1.dat",
options= "vdiv bins=40 limits[-2.0 2.0] ; y bins=120 limits[-0.06 0.06]")
AT (0,0,0.0001) RELATIVE PREVIOUS
ROTATED (0,0,0) RELATIVE PREVIOUS
COMPONENT Vdivlam_Det1 = Monitor_nD(xmin=-0.06, xmax=0.06, ymin=-0.06, ymax=0.06,
filename="Vdiv_det1.dat",
options= "lamda wavelength bins=70 limits[2.0 9.0] ; vdiv bins=80 limits[-2.0 2.0]")
AT (0,0,0.0001) RELATIVE PREVIOUS
ROTATED (0,0,0) RELATIVE PREVIOUS
COMPONENT Hdivlam_Det1 = Monitor_nD(xmin=-0.06, xmax=0.06, ymin=-0.06, ymax=0.06,
filename="Hdiv_det1.dat",
options= "lamda wavelength bins=70 limits[2.0 9.0] ; hdiv bins=80 limits[-2.0 2.0]")
AT (0,0,0.0001) RELATIVE PREVIOUS
ROTATED (0,0,0) RELATIVE PREVIOUS
END
-------------- next part --------------
OFF
8 4 0
-0.05 -0.05 0
-0.05 0.05 0
0.05 0.05 0
0.05 -0.05 0
-0.05 -0.05 20.0
-0.05 0.05 20.0
0.05 0.05 20.0
0.05 -0.05 20.0
4 0 1 5 4 0
4 1 2 6 5 1
4 2 3 7 6 0
4 3 0 4 7 1
More information about the mcstas-users
mailing list