About component version and origin
Emmanuel Farhi
farhi at ill.fr
Fri Mar 16 10:14:54 CET 2001
Hello all of you,
About the Gravity_guide (attached file):
The 'G' constant should be a protected variable, to avoid conflicts with other macros:
OUTPUT PARAMETERS (ww,hh,whalf,hhalf,lwhalf,lhhalf)
...
double G;
...
G=9.82;
About the Bender:
The Bender doc header is not valid. I here attach a full version.
About version control for comp:
I do not think the proposition of Per-Olof will solve the problems of indentifying what are the
component versions and origin. At the end you would have only 'official updated' and
'unofficial' components, and this info would not be very useful !
The best way would be to have a 'version' column in the McDoc automatic page (list of comps).
Also, the column 'Help' is not very informative (except when you click on 'More...'. perhaps
this 'More...' item could be moved to 'More (version number)'. If McStas is to be opened and
shared by users, I think the real origin of the components should be displayed clearly, even
when these comps are included in the official release. I like to know that the bender comes from
Erlangen.
When you put 'official release', just add also the McStas version it is associated to (modify
headers to have this displayed).
Comp update
Also, perhaps it is the time, for version 1.4.1 to come, to update the POLARISATION PARAMETERS
(sx, sy, sz) in every comp ? (I have not done that yet)
Version 1.4.1 to come
I have put in the ILL McStas page (http://www.ill.fr/tas/mcstas/) an updated version of McStas
1.4 called mcstas-1.4a. It includes all the updates that were discussed these two last weeks:
* new lib/mcstas-r.h and .c to remove the 2e9 limit (passed run_num counter from int to
double)
* signal handler for Unix/Linux (to see the simulation status with kill -USR1 <pid>, and also
to stop it saving results with other signals/interupts) in lib/mcstas-r.h and .c
* a new internal mcstas function 'double get_run_num()' to get the number of iterations
performed until now in a simulation.
* new Monitor_nD (0.13.9) with 3He gas handling, intermediate results saving during long
simulations (works with the new mcstas-r lib only). a PreMonitor_nD to make correlations
between two positions.
* bug corrected in Source_Optimizer (0.09), and asociated Monitor_Optimizer (0.08) that
caused a SIGFPE error (div per 0 during optimisation of the spin divergence)
* new Gravity_Guide (1.2) protecting internal variables, and having G as a protected comp
symbol (not a macro)
* new Bender (1.1)
* an update of the 'cogen.c' file for correct handling of string parametrs in instruments
* an update of lib/mcfrontlib.pl to plot results for simulations using instrument string
parameters
This is to wait until 1.4.1 appears !!
Cheers. Emmanuel.
Per-Olof Åstrand wrote:
> Hi,
>
> Thanks for pointing this out. I have put the new version in the McDoc source and both the
> Risø and ILL web-pages will be updated automatically by the McDoc program.
>
> Since Bender.comp is available at the web-pages but is not a part of the official McStas
> release, it gives me the opportunity to elaborate on this. The web-pages at Risø and ILL
> generated by the McDoc program contain components from the official release, unofficial
> components provided by users, and updated versions of official components. Just by inspection
> on the web-page, it may be very difficult to be able to decide which of these categories a
> component belongs to. For the official components developed at Risø, the "Origin" entry
> always say "McStas release", even if it is an updated component different from the version in
> the official release. For components kindly provided by users, the "Origin" entry is mostly
> used for giving a geographic origin regardless of if it is an official component or not.
>
> For me (and I guess also for others), some time is spent for each component (e.g. the Bender
> component) to figure out where its "origin" is. I therefore suggest that the "Origin"
> definition in a component can only contain three different entries: "McStas release", "McStas
> release, updated version" and "Unofficial", and that the "Authors" definition is used for
> specifying the geographic origin. I would like to have some response on this, otherwise I
> will just implement it in a week or so.
>
> Finally, the McDoc program is a very efficient tool maintained by ILL and Risø for sharing
> components and keep them up-to-date. Users are therefore encouraged to make their components
> public either by sending them to us or to establish an additional web-server handling McDoc
> (which is not too difficult) in addition to ILL and Risø. If a component is of general
> interest, it may eventually be, with the permission of the author, included in the official
> McStas release and documented in the manual.
>
> Best regards,
>
> Per-Olof Åstrand
>
> Stuart Rycroft wrote:
>
> > Hello
> >
> > I just wanted to mention a problem with the bender component, which can be
> > downloaded from the Riso or ILL website. They both have an early version
> > of the component (Feb. 7th 1999) but from the email archive, a later
> > version (Feb. 22nd 1999) was sent by the author (Philipp Bernhardt) which
> > corrects a small bug in the calculation of the reflectivity for the top
> > and bottom mirrors. Without this correction the output from my simulated
> > guides is obviously not correct when looking at the divergence, but it is
> > very hard to see a problem when just looking at an intensity since the
> > only effect I have seen is a small loss, maybe around five percent.
> >
> > The updated component and the problem are given by the author in the
> > mailing archive:
> >
> > http://neutron.risoe.dk/cgi-bin/wilma_hiliter/neutron-mc/9902/msg00023.html?line=24#hilite
> >
> > Maybe the component list on the Mcstas webpages should have this version
> > instead of the old one.
> >
> > Regards,
> > Stuart Rycroft
> > IRI, Delft
> >
> >
--
What's up Doc ?
--------------------------------------------
Emmanuel FARHI, http://www.ill.fr/tas/people/Farhi.html \|/ ____ \|/
CS-Group ILL4/156, Institut Laue-Langevin (ILL) Grenoble ~@-/ oO \-@~
6 rue J. Horowitz, BP 156, 38042 Grenoble Cedex 9,France /_( \__/ )_\
Work :Tel (33/0) 4 76 20 71 35. Fax (33/0) 4 76 48 39 06 \__U_/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20010316/dfb6e851/attachment.html>
-------------- next part --------------
/*******************************************************************************
*
* Component: Bender
*
* %Identification
* Written by: Philipp Bernhardt
* Date: Februar 7 1999
* Origin: Uni. Erlangen (Germany)
* Version: 1.1
* Modified by: PB, Februar 22 1999
*
* Models a curved neutron guide.
*
* %Description
* Models a curved neutron guide.The entrance lies in the X-Y plane, centered
* on the Z axis. The neutrons will also leave the bender in the X-Y plane at
* the z-value r*Win, so you do not have to calculate the real exit coordinates
* and you do not need a new arm. The bender is bent to the negative X axis;
* it behaves like a parallel guide in the Y axis.
* There is no warning, if the input parameters are wrong, so please make sure
* that
* w,h,r,d,Win,k are positive numbers and
* k*d is smaller than the width w.
*
* Example values: w=0.05 h=0.12 r=250 d=0.001 Win=0.04 k=3 ma=3 mi=2 ms=2
* See also <a href="http://neutron.risoe.dk/mcstas/support/misc/Bender.html">Additional note</a> from <a href="mailto:philipp.bernhardt at krist.uni-erlangen.de">Philipp Bernhardt</a>.
*
*
* %Parameters
* INPUT PARAMETERS:
*
* w: (m) Width at the bender entry and exit
* h: (m) Height at the bender entry and exit
* r: (m) Radius of the bender
* d: (m) Thickness of the partition, which separates the channels
* Win: (rad) Angle of the deflection of the whole bender
* k: (1) Number of channels inside the bender
* R0a: (1) Low-angle reflectivity at the bender <b>concave</b> side
* Qca: (AA-1) Critical scattering vector
* alphaa: (AA) Slope of reflectivity
* ma: (1) m-value of material
* Wa: (AA-1) Width of supermirror cut-off
* R0i: (1) Low-angle reflectivity at the bender <b>convex</b> side
* Qci: (AA-1) Critical scattering vector
* alphai: (AA) Slope of reflectivity
* mi: (1) m-value of material
* Wi: (AA-1) Width of supermirror cut-off
* R0s: (1) Low-angle reflectivity at bender <b>top and bottom</b> side
* Qcs: (AA-1) Critical scattering vector
* alphas: (AA) Slope of reflectivity
* ms: (1) m-value of material
* Ws: (AA-1) Width of supermirror cut-off
*
* %End
*******************************************************************************/
DEFINE COMPONENT Bender
DEFINITION PARAMETERS (w,h,r,d,Win,k,R0a,Qca,alphaa,ma,Wa,
R0i,Qci,alphai,mi,Wi,R0s,Qcs,alphas,ms,Ws)
SETTING PARAMETERS ()
OUTPUT PARAMETERS (bk)
STATE PARAMETERS (x, y, z, vx, vy, vz, t, s1, s2, p)
DECLARE
%{
double bk;
#ifndef BENDER_DECLARE
#define BENDER_DECLARE
double sgn(double x) {
if (x>=0)
return 1.0;
else
return -1.0; }
#endif
%}
INITIALIZE
%{
/* width of one channel + thickness d of partition */
bk=(w+d)/k;
%}
TRACE
%{
int i,num,numa,numi;
double dru,ab,dab,R,Q,arg,arga,argi,Ta,vpl;
double einwin,auswin,zykwin,aeuwin,innwin,ref,innref,aeuref;
double einzei,auszei,zykzei;
/* does the neutron hit the bender at the entrance? */
PROP_Z0;
if ((fabs(x)>w/2) || (fabs(y)>h/2))
ABSORB;
/*** reflections in the XZ-plane ***/
/* distance between neutron and concave side of the channel at the
entrance */
dru=floor((w/2-x)/bk)*bk;
ab=w/2.0-x-dru;
/* radius of the channel */
R=r-dru;
/* does the neutron hit the partition at the entrance? */
if (ab>bk-d)
ABSORB;
/* velocity in the XZ-plane */
vpl=sqrt(vx*vx+vz*vz);
/* divergence of the neutron at the entrance */
einwin=atan(vx/vz);
/* maximal distance between neutron and concave side of the channel */
dab=R-cos(einwin)*(R-ab);
/* reflection angle at the concave side */
aeuwin=acos((R-dab)/R);
/* reflection coefficient at the concave side */
arga=0.0;
Q=2.0*V2K*vpl*sin(aeuwin);
if (Q<=Qca)
aeuref=R0a;
else {
arga=(Q-ma*Qca)/Wa;
aeuref=0.5*R0a*(1.0-tanh(arga))*(1.0-alphaa*(Q-Qca)); }
/* does the neutron hit the convex side of the channel? */
innwin=0.0;
innref=1.0;
argi=0.0;
if (dab>bk-d) {
/* reflection coefficient at the convex side */
innwin=acos((R-dab)/(R-bk+d));
Q=2.0*V2K*vpl*sin(innwin);
if (Q<=Qci)
innref=R0i;
else {
argi=(Q-mi*Qci)/Wi;
innref=0.5*R0i*(1.0-tanh(argi))*(1.0-alphai*(Q-Qci)); } }
/* divergence of the neutron at the exit */
zykwin=2.0*(aeuwin-innwin);
auswin=fmod(Win+einwin+aeuwin-innwin*(1.0+sgn(einwin)),zykwin)-zykwin/2.0;
auswin+=innwin*sgn(auswin);
/* number of reflections at the concave side */
numa=(Win+einwin+aeuwin-innwin*(1.0+sgn(einwin)))/zykwin;
/* number of reflections at the convex side */
numi=numa;
if (auswin*einwin<0) {
if (auswin-einwin>0)
numi++;
else
numi--; }
/* is the reflection coefficient too small? */
if (((numa>0) && (arga>10.0)) || ((numi>0) && (argi>10.0)))
ABSORB;
/* calculation of the neutron probability weight p */
for (i=1;i<=numa;i++)
p*=aeuref;
for (i=1;i<=numi;i++)
p*=innref;
/* time to cross the bender */
Ta=(2*numa*(tan(aeuwin)-tan(innwin))+tan(auswin)-tan(einwin)-tan(innwin)*
(sgn(auswin)-sgn(einwin)))*(R-dab)/vpl;
t+=Ta;
/* distance between neutron and concave side of channel at the exit */
ab=R-(R-dab)/cos(auswin);
/* calculation of the exit coordinates in the XZ-plane */
x=w/2.0-ab-dru;
z=r*Win;
vx=sin(auswin)*vpl;
vz=cos(auswin)*vpl;
/*** reflections at top and bottom side (Y axis) ***/
if (vy!=0.0) {
/* reflection coefficent at the top and bottom side */
arg=0.0;
Q=2.0*V2K*fabs(vy);
if (Q<=Qcs)
ref=R0s;
else {
arg=(Q-ms*Qcs)/Ws;
ref=0.5*R0s*(1.0-tanh(arg))*(1.0-alphas*(Q-Qcs)); }
/* number of reflections at top and bottom */
einzei=h/2.0/fabs(vy)+y/vy;
zykzei=h/fabs(vy);
num=(Ta+einzei)/zykzei;
/* time between the last reflection and the exit */
auszei=fmod(Ta+einzei,zykzei);
/* is the reflection coefficient too small? */
if ((num>0) && (arg>10.0))
ABSORB;
/* calculation of the probability weight p */
for (i=1;i<=num;i++) {
p*=ref;
vy*=-1.0; }
/* calculation of the exit coordinate */
y=auszei*vy-vy*h/fabs(vy)/2.0; }
%}
END
-------------- next part --------------
/*******************************************************************************
*
* McStas, the neutron ray-tracing package
* Maintained by Per-Olof Astrand and Kim Lefmann,
* Copyright 1997-2001 Risoe National Laboratory, Roskilde, Denmark
*
* %I
* Written by: KL
* Date: November 22 2000
* Version: $Revision: 1.2 $
* Origin: McStas release
* Modified by: KL, EF, March 15, 2001 protect comp variable in OUTPUT macro
* Modified by: POA, March 15, 2001; small changes of documentation
*
* Neutron guide with gravity.
* To do: information about gravity should be included in the kernel.
*
* %D
* Models a rectangular guide tube centered on the Z axis. The entrance lies
* in the X-Y plane.
* A constant gravity is applied in the Y direction of the LOCAL coordinate
* system. No gravity is applied for the flight path before the guide.
* For details on the geometry calculation see the description in the McStas
* reference manual.
*
* %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
*
* %D
* Example values: m=4 Qc=0.02 W=1/300 alpha=6.49 R0=1
*
* %E
*******************************************************************************/
DEFINE COMPONENT Gravity_guide
DEFINITION PARAMETERS ()
SETTING PARAMETERS (w1, h1, w2, h2, l, R0, Qc, alpha, m, W)
OUTPUT PARAMETERS (ww,hh,whalf,hhalf,lwhalf,lhhalf, G)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
double G;
double ww,hh,whalf,hhalf,lwhalf,lhhalf;
%}
INITIALIZE
%{
ww = .5*(w2 - w1), hh = .5*(h2 - h1);
whalf = .5*w1, hhalf = .5*h1;
lwhalf = l*whalf, lhhalf = l*hhalf;
G = 9.82;
%}
TRACE
%{
double t1,tv1,tv2,th1,th2; /* 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 D_1,b_1,D_2,b_2; /* coeefficients for 2nd order equations */
int i; /* Which mirror hit? */
double q; /* Q [1/AA] of reflection */
double vlen2,nlen2; /* Vector lengths squared */
/* Propagate neutron to guide entrance. */
/* NO GRAVITY HERE! */
PROP_Z0;
if(x <= -whalf || x >= whalf || y <= -hhalf || y >= hhalf)
ABSORB;
for(;;)
{
/* Compute the dot products of v and n for the two vertical mirrors. */
av = l*vx; bv = ww*vz;
vdotn_v1 = bv + av; /* Left vertical */
vdotn_v2 = bv - av; /* Right vertical */
/* Compute the dot products of (O - r) and n as c1+c2 and c1-c2 */
cv1 = -whalf*l - z*ww; cv2 = x*l;
/* Compute intersection times. WARNING, there may be division by zero */
t1 = (l - z)/vz;
tv1= (cv1 - cv2)/vdotn_v1;
tv2 = (cv1 + cv2)/vdotn_v2;
/* Now: solve the 2nd order equation for the "horizontal" surfaces*/
b_1 = -vy+hh*vz/l;
D_1 = b_1*b_1+2*G*(y-hhalf);
if(D_1 < 0)
th1= -1;
else /* Choose first intersection with upper surface */
th1= (-b_1-sqrt(D_1))/G;
b_2 = -vy-hh*vz/l;
D_2 = b_2*b_2+2*G*(y+hhalf);
if(D_2 < 0)
{
printf("Fatal error, neutron cannot intersect lower guide surface \n");
exit(1);
}
else /* Choose second intersection with lower surface */
th2= (-b_2+sqrt(D_2))/G;
i = 0;
if(vdotn_v1 < 0 && tv1 < t1)
{
t1 = tv1;
i = 1;
}
if(vdotn_v2 < 0 && tv2 < t1)
{
t1 = tv2;
i = 2;
}
if(th1>0 && th1 < t1)
{
t1 = th1;
i = 3;
}
if(th2>0 && th2 < t1)
{
t1 = th2;
i = 4;
}
if(i == 0)
break; /* Neutron left guide. */
ah = l*vy; bh = hh*vz;
vdotn_h1 = bh + ah; /* Lower "horizontal" */
vdotn_h2 = bh - ah; /* Upper "horizontal" */
ch1 = -hhalf*l - z*hh; ch2 = y*l;
x+=vx*t1;
z+=vz*t1;
y+=vy*t1-G*t1*t1/2;
vy-= G*t1; /* TO BE REPLACED BY A MACRO: PROP_GRAV(dt) */
t+=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. */
if(m == 0)
ABSORB;
if(q > Qc)
{
double arg = (q-m*Qc)/W;
if(arg < 10)
p *= .5*(1-tanh(arg))*(1-alpha*(q-Qc));
else
ABSORB; /* Cutoff ~ 1E-10 */
}
p *= R0;
SCATTER;
}
x+=vx*t1;
z+=vz*t1;
y+=vy*t1-G*t1*t1/2;
vy-= G*t1; /* TO BE REPLACED WITH A MACRO: PROP_GRAV(dt) */
t+=t1;
%}
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
More information about the mcstas-users
mailing list