[mcstas-users] resonant spin echo flipper in McStas 3.4
Lukas Vogl
Lukas.Vogl at frm2.tum.de
Wed Oct 18 22:44:52 CEST 2023
Hi McStas Team,
At the neutron resonance spin echo instrument RESEDA (FRM II), we are
working on implementing a virtual clone of the instrument, which
combines our instrument control software with a McStas simulation in the
background.
At the heart of the instrument are radio frequency spin flippers in a
longitudinal field geometry. The static field is parallel for the
neutron fight path and the rotating RF-field perpendicular to it.
The build-in McStas functionality for neutron spin evolution seems to
work up to circular frequencies of 2pi*500 1/s. We want to perform
simulations with frequencies in the MHz range.
The particular problem we see, is that the spin is no longer flipped by
Pi for higher frequencies, and only a partial polarisation in the
opposite direction after the flipper is seen.
I have attached the necessary files to reproduce this problem.
pol-lib.c, pol-lib.h and Pol_Bfield.comp have to be placed in the
corresponding folders (share, resp. optics) of the McStas installation
version 3.4.
pol-lib.* and Pol_Bfield.comp are modified to implement a rotating
B-field, which is not yet present in the current McStas version.
I already have tried to change the threshold condition for the time step
calculation to allow smaller time steps but this did not help.
Changing the length of the flipper does seem to influence the
performance which should also not happen.
Thus our question is, if the current McStas spin propagation code is at
all capable to handle RF-fields with such frequencies and if so what we
have to change in our implementation.
Best Regards
Lukas Vogl
Working student at
RESEDA - FRM II
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pol-lib.h
Type: text/x-chdr
Size: 4558 bytes
Desc: not available
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20231018/8fc2a4b3/attachment-0002.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pol-lib.c
Type: text/x-csrc
Size: 17898 bytes
Desc: not available
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20231018/8fc2a4b3/attachment-0003.bin>
-------------- next part --------------
/**************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2006, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Pol_Bfield
*
* %I
* Written by: Erik B Knudsen, Peter Christiansen and Peter Willendrup
* Date: July 2011
* Origin: RISOE
*
* Magnetic field component.
*
* %D
*
* Region with a definable magnetic field.
*
* The component is nestable if the concentric flag is set (the default). This means that it may have a
* construction like the following:
* // START MAGNETIC FIELD
* COMPONENT msf =
* Pol_Bfield(xwidth=0.08, yheight=0.08, zdepth=0.2, Bx=0, By=-0.678332e-4, Bz=0)
* AT (0, 0, 0) RELATIVE armMSF
*
* // OTHER COMPONENTS INSIDE THE MAGNETIC FIELD CAN BE PLACED HERE
*
* // STOP MAGNETIC FIELD
* COMPONENT msf_stop = Pol_simpleBfield_stop(magnet_comp_stop=msf)
* AT (0, 0, 0) RELATIVE armMSF
*
* Note that if you have objects within the magnetic field that extend outside it you may get
* wrong results, because propagation within the field will then possibly extend outside, e.g.
* when using a tabled field. The evaluated field will then use the nearest defined field point
* _also_ outside the defintion area. If these outside-points have a non-zero field precession will
* continue - even after the neutron has left the field region.
*
* In between the two component instances the propagation routine
* PROP_DT also handles spin propagation.
* The current algorithm used for spin propagation is:
* SimpleNumMagnetPrecession
* in pol-lib.
*
* Example: Pol_Bfield(xwidth=0.1, yheight=0.1, zdepth=0.2, Bx=0, By=1, Bz=0)
* Pol_Bfield(xwidth=0.1, yheight=0.1, zdepth=0.2, field_type=1)
*
* Functions supplied by the system are (the number is the ID of the function to be given as the field_type parameter:
* 1. Constant magnetic field: Constant field (Bx,By,Bz) within the region
* 2. Rotating magnetic_field: Field is initially (0,By,0) but after a length of zdepth
* has rotated to (By,0,0)
* 3. Majorana magnetic_field: Field is initially (Bx,By,0) with By<<Bx, then linearly transforms to
* (-Bx,By,0) at z=zdepth.
* 4. MSF field:
* 5. RF field: A oscillating B-field in the x plane with the amplitude given by Bx, freqency given by By and phase given by Bz.
* 6. Gradient field: Similar to Majorana by without an x-component to the field. I.e. it (0,By,0) at z=0, and
* becomes (0,-By,0) AT z=zdepth.
* 7. Rotating RF field: A rotating B-field in the x-y plane with the amplitude given by Bx, freqency given by By and phase given by Bz.
*
* If the concentric parameter is set to 1 the magnetic field is turned on by an
* instance of this component, and turned off by an instance of Pol_simpleBfield_stop.
* Anything in between is considered inside the field.
* If concentric is zero the field is considered a closed component - neutrons are propagated
* through the field region, and no other components are permitted inside the field.
*
* %P
* INPUT PARAMETERS:
*
* xwidth: [m] Width of opening.
* yheight: [m] Height of opening.
* zdepth: [m] Length of field.
* radius: [m] Radius of field if it is cylindrical or spherical.
* Bx: [T] Parameter used for x composant of field.
* By: [T] Parameter used for y composant of field.
* Bz: [T] Parameter used for z composant of field.
* field_type: [] ID of function describing the magnetic field. 1:constant field, 2: rotating fiel, 3: majorana type field, 4: MSF, 5: RF-field, 6: gradient field, 7: rotating RF-field.
* concentric: [ ] Allow components and/or other fields inside field. Requires a subsequent Pol_simpleBfield_stop component.
*
* %E
****************************************************************************/
DEFINE COMPONENT Pol_Bfield
DEFINITION PARAMETERS ()
SETTING PARAMETERS (int field_type=0, xwidth=0, yheight=0,zdepth=0, radius=0, Bx=0, By=0, Bz=0, concentric=1)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
%include "pol-lib"
double fmax(double, double);
double fmin(double, double);
%}
DECLARE
%{
/* Larmor frequency and scalar product threshold*/
double Bprms[8];
int shape;
%}
INITIALIZE
%{
enum shapes {NONE=0, BOX, WINDOW, CYLINDER, SPHERE, ANY};
enum field_functions ff=field_type;
/* these default field functions are part of pol-lib */
if (ff==constant){
double *t=Bprms;
t[0]=Bx;t[1]=By;t[2]=Bz;
} else if (ff==rotating){
double *t=Bprms;
t[0]=By;t[1]=zdepth;
} else if (ff==majorana){
double *t=Bprms;
t[0]=Bx; t[1]=By; t[2]=zdepth;
} else if (ff==RF){
double *t=Bprms;
t[0]=Bx; t[1]=By; t[2]=Bz;
} else if (ff==rotatingRF){
double *t=Bprms;
t[0]=Bx; t[1]=By; t[2]=Bz;
}
if(xwidth && yheight && zdepth){
shape=BOX;
}else if (xwidth && yheight && !zdepth){
shape=WINDOW;
}else if(radius && yheight){
shape=CYLINDER;
}else if (radius) {
shape=SPHERE;
}else{
shape=NONE;
}
if(shape == WINDOW && !concentric){
printf("Warning (%s): Magnetic field has window shape and concentric==0 => Field volume == 0.\n",NAME_CURRENT_COMP);
}
if(shape == NONE) {
printf("Warning (%s): Magnetic field has no geometry. Full Z=0-plane is used as boundary.\n",NAME_CURRENT_COMP);
}
%}
TRACE
%{
double t0,t1;
int hit;
enum shapes {NONE=0, BOX, WINDOW, CYLINDER, SPHERE, ANY};
/*enter through whatever object we are*/
switch (shape){
case BOX:
hit=box_intersect(&t0,&t1,x,y,z,vx,vy,vz,xwidth,yheight,zdepth);
/*terminate neutrons which miss the component*/
if(!hit) ABSORB;
/*If we do hit - propagate to the start of the field unless the neutron is already there.*/
if(t0>FLT_EPSILON) {
PROP_DT(t0);
t1-=t0;
}else if (t0<-FLT_EPSILON && concentric){
PROP_DT(t1);
}
break;
case CYLINDER:
hit=cylinder_intersect(&t0,&t1,x,y,z,vx,vy,vz,radius,yheight);
/*terminate neutrons which miss the component*/
if(!hit)ABSORB;
/*If we do hit - propagate to the start of the field unless the neutron is already there.*/
if(t0>FLT_EPSILON) {
PROP_DT(t0);
t1-=t0;
}else if (t0<-FLT_EPSILON && concentric){
PROP_DT(t1);
}
break;
case WINDOW:
PROP_Z0;
if (2*x>xwidth || 2*x<-xwidth || 2*y>yheight || 2*y<-yheight){
/*terminate neutrons which miss the component*/
ABSORB;
}
break;
default:
PROP_Z0;
}
mcmagnet_push(_particle, field_type,&(ROT_A_CURRENT_COMP),&(POS_A_CURRENT_COMP),0,Bprms);
#ifdef MCDEBUG
mcmagnet_print_stack();
#endif
/*does the field "close/stop" itself*/
if (!concentric){
switch (shape){
case BOX:
PROP_DT(t1);
break;
case CYLINDER:
PROP_DT(t1);
break;
case WINDOW:
PROP_Z0;
/*terminate neutrons which miss the exit window*/
if (2*x>xwidth || 2*x<-xwidth || 2*y>yheight || 2*y<-yheight){
ABSORB;
}
break;
default:
PROP_Z0;
}
mcmagnet_pop(_particle);
}
%}
MCDISPLAY
%{
enum shapes {NONE=0, BOX, WINDOW, CYLINDER, SPHERE, ANY};
switch (shape){
case BOX:
box(0,0,0,xwidth,yheight,zdepth);
break;
case CYLINDER:
circle("xz",0, yheight/2.0,0,radius);
circle("xz",0,-yheight/2.0,0,radius);
line(-radius,-yheight/2.0,0,-radius,yheight/2.0,0);
line( radius,-yheight/2.0,0, radius,yheight/2.0,0);
line(0,-yheight/2.0,-radius,0,yheight/2.0,-radius);
line(0,-yheight/2.0, radius,0,yheight/2.0, radius);
break;
case SPHERE:
circle("xz",0,0,0,radius);
circle("xy",0,0,0,radius);
circle("yz",0,0,0,radius);
break;
case WINDOW:
rectangle("xy",0,0,0,xwidth,yheight);
break;
}
%}
END
-------------- next part --------------
/* vim: set filetype=c:
*/
DEFINE INSTRUMENT rf_flipper
(
)
DECLARE
%{
#define GAMMA_N 1.83247185e+08 /* rad/(Ts), value from: boost/units/systems/si/codata/neutron_constants.hpp */
double lam;
double dlam;
double vel;
double B_rf;
double B0;
double freq;
double coillen;
double omega;
%}
INITIALIZE
%{
coillen = 0.02;
lam = 6;
dlam = 0.6;
freq = 450;
omega = freq * 2*M_PI;
B0 = omega /(GAMMA_N);
vel = K2V*(2*PI/lam);
B_rf =-PI * vel /(GAMMA_N*coillen);
printf("coil: B_0=%g\n B_rf=%g\n freq=%g\n vel=%g\n", B0, B_rf, freq,vel);
%}
TRACE
COMPONENT Origin = Progress_bar()
AT (0,0,0) ABSOLUTE
COMPONENT Source = Source_simple(
yheight = 0.01, xwidth = 0.01,
focus_xw = 0.01, focus_yh = 0.01, dist = 5.0,
lambda0 = lam, dlambda = dlam, gauss = 0)
AT (0, 0, 0) RELATIVE Origin
EXTEND
%{
sx = 0.; //
sy = 0.; // setting the polarization to z
sz = -1.; //
%}
COMPONENT Pol_monitor_before = PolLambda_monitor(
xwidth = 0.1, yheight = 0.1,
nL = 200, npol = 61, mx = 0, my=0,mz=1,
restore_neutron = 1, filename = "pre.pol",
Lmin = lam-dlam, Lmax = lam+dlam)
AT (0, 0, 0.5) RELATIVE Origin
COMPONENT stationary_bfield = Pol_Bfield(
field_type=1, concentric = 1,
Bx=0,By=0,Bz=B0,
xwidth=0.05,yheight=0.05, zdepth=coillen)
AT (0, 0, 1) RELATIVE Origin
COMPONENT rf_bfield = Pol_Bfield(
field_type=7, concentric = 0,
Bx=B_rf,By=omega,Bz=0,
xwidth=0.05,yheight=0.05, zdepth=coillen)
AT (0, 0, 1) RELATIVE Origin
COMPONENT stationary_bfield_end = Pol_Bfield_stop()
AT (0,0,1+coillen) RELATIVE Origin
COMPONENT Pol_monitor_after = PolLambda_monitor(
xwidth = 0.1, yheight = 0.1,
nL = 200, npol = 201, mx = 0, my=0,mz=1,
restore_neutron = 1, filename = "post.pol",
Lmin = lam-dlam, Lmax = lam+dlam)
AT (0, 0, 2) RELATIVE Origin
FINALLY
%{
%}
END
More information about the mcstas-users
mailing list