/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright (C) 1997-2007, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Monochromator_curved
*
* %I
*
* Written by: Emmanuel Farhi, Kim, Lefmann, Peter Link
* Date: Aug. 24th 2001
* Origin: ILL
* Modified by: EF, Aug. 24th 2001: From Mosaic_anisotropic and Mon_2foc
* Modified by: EF, Feb 13th 2002: Read reflectivity table
* Modified by: EF, Oct 24th 2002: Read transmission table
*
* Double bent multiple crystal slabs with anisotropic gaussian mosaic.
*
* %Description
* Double bent infinitely thin mosaic crystal, useful as a monochromator or
* analyzer. which uses a small-mosaicity approximation and taking into account
* higher order scattering. The mosaic is anisotropic gaussian, with different
* FWHMs in the Y and Z directions. The scattering vector is perpendicular to the
* surface. For an unrotated monochromator component, the crystal plane lies in
* the y-z plane (ie. parallel to the beam). The component works in reflection, but
* also transmits the non-diffracted beam. Reflectivity and transmission files may
* be used. The slabs are positioned in the vertical plane (not on a
* cylinder/sphere), and are rotated according to the curvature radius.
* When curvatures are set to 0, the monochromator is flat.
* The curvatures approximation for parallel beam focusing to distance L, with
* monochromator rotation angle A1 are:
* RV = 2*L*sin(DEG2RAD*A1);
* RH = 2*L/sin(DEG2RAD*A1);
*
* When you rotate the component by A1 = asin(Q/2/Ki)*RAD2DEG, do not forget to
* rotate the following components by A2=2*A1 (for 1st order) !
*
* Example: Monochromator_curved(zwidth=0.01, yheight=0.01, gap=0.0005, NH=11, NV=11, mosaich=30.0, mosaicv=30.0, r0=0.7, Q=1.8734)
*
* Monochromator lattice parameter
* PG 002 DM=3.355 AA (Highly Oriented Pyrolythic Graphite)
* PG 004 DM=1.677 AA
* Heusler 111 DM=3.362 AA (Cu2MnAl)
* CoFe DM=1.771 AA (Co0.92Fe0.08)
* Ge 111 DM=3.266 AA
* Ge 311 DM=1.714 AA
* Ge 511 DM=1.089 AA
* Ge 533 DM=0.863 AA
* Si 111 DM=3.135 AA
* Cu 111 DM=2.087 AA
* Cu 002 DM=1.807 AA
* Cu 220 DM=1.278 AA
* Cu 111 DM=2.095 AA
*
* %Parameters
* INPUT PARAMETERS:
*
* zwidth: [m] horizontal width of an individual slab
* yheight: [m] vertical height of an individual slab
* gap: [m] typical gap between adjacent slabs
* NH: [int] number of slabs horizontal
* NV: [int] number of slabs vertical
* mosaich: [arc minutes] Horisontal mosaic FWHM
* mosaicv: [arc minutes] Vertical mosaic FWHM
* r0: [1] Maximum reflectivity. O unactivates component
* Q: [AA-1] Scattering vector
* RV: [m] radius of vertical focussing. flat for 0
* RH: [m] radius of horizontal focussing. flat for 0
*
* optional parameters
* DM: [AA] monochromator d-spacing instead of Q=2*pi/DM
* mosaic: [arc minutes] sets mosaich=mosaicv
* width: [m] total width of monochromator, along Z
* height: [m] total height of monochromator, along Y
* verbose: [0/1] verbosity level
* reflect: [str] reflectivity file name of text file as 2 columns [k, R]
* transmit: [str] transmission file name of text file as 2 columns [k, T]
* t0: [1] transmission efficiency
* order: [1] specify the diffraction order, 1 is usually prefered. Use 0 for all
*
* %Link
* Additional note from Peter Link.
* %L
* Obsolete Mosaic_anisotropic by Kristian Nielsen
* %L
* Contributed Monochromator_2foc by Peter Link
*
* %End
*******************************************************************************/
DEFINE COMPONENT Monochromator_curved
DEFINITION PARAMETERS ()
SETTING PARAMETERS (string reflect="NULL", string transmit="NULL",
zwidth=0.01, yheight=0.01,gap=0.0005,int NH=11, int NV=11, mosaich=30.0, mosaicv=30.0, r0=0.7, t0=1.0,Q=1.8734,
RV=0, RH=0, DM=0, mosaic=0, width=0, height=0, verbose=0, order=0)
OUTPUT PARAMETERS ()
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
#pragma acc routine
double GAUSS_monocurved(double x, double mean, double rms){
return (exp(-((x)-(mean))*((x)-(mean))/(2*(rms)*(rms)))/(sqrt(2*PI)*(rms)));
}
%include "read_table-lib"
%}
DECLARE
%{
double mos_rms_y; /* root-mean-square of mosaic, in radians */
double mos_rms_z;
double mos_rms_max;
double mono_Q;
double SlabWidth;
double SlabHeight;
t_Table rTable;
t_Table tTable;
int rTableFlag;
int tTableFlag;
double *tiltH;
double *tiltV;
%}
INITIALIZE
%{
int i;
if (mosaic != 0) {
mos_rms_y = MIN2RAD*mosaic/sqrt(8*log(2));
mos_rms_z = mos_rms_y; }
else {
mos_rms_y = MIN2RAD*mosaich/sqrt(8*log(2));
mos_rms_z = MIN2RAD*mosaicv/sqrt(8*log(2)); }
mos_rms_max = mos_rms_y > mos_rms_z ? mos_rms_y : mos_rms_z;
mono_Q = Q;
if (DM != 0) mono_Q = 2*PI/DM;
if (mono_Q <= 0) { fprintf(stderr,"Monochromator_curved: %s: Error scattering vector Q = 0\n", NAME_CURRENT_COMP); exit(-1); }
if (r0 < 0) { fprintf(stderr,"Monochromator_curved: %s: Error reflectivity r0 is negative\n", NAME_CURRENT_COMP); exit(-1); }
if (r0 == 0) { fprintf(stderr,"Monochromator_curved: %s: Reflectivity r0 is null. Ignoring component.\n", NAME_CURRENT_COMP); }
if (NH*NV == 0) { fprintf(stderr,"Monochromator_curved: %s: no slabs ??? (NH or NV=0)\n", NAME_CURRENT_COMP); exit(-1); }
if (verbose && r0)
{
printf("Monochromator_curved: component %s Q=%.3g Angs-1 (DM=%.4g Angs)\n", NAME_CURRENT_COMP, mono_Q, 2*PI/mono_Q);
if (NH*NV == 1) printf(" flat.\n");
else
{ if (NH > 1)
{ printf(" horizontal: %i blades", (int)NH);
if (RH != 0) printf(" focusing with RH=%.3g [m]", RH);
printf("\n");
}
if (NV > 1)
{ printf(" vertical: %i blades", (int)NV);
if (RV != 0) printf(" focusing with RV=%.3g [m]", RV);
printf("\n");
}
}
}
if (reflect != NULL && r0 && strlen(reflect) && strcmp(reflect,"NULL") && strcmp(reflect,"0"))
{
if (verbose) fprintf(stdout, "Monochromator_curved: %s: Reflectivity data (k, R) from %s\n", NAME_CURRENT_COMP, reflect);
Table_Read(&rTable, reflect, 1); /* read 1st block data from file into rTable */
Table_Rebin(&rTable); /* rebin as evenly, increasing array */
if (rTable.rows < 2) Table_Free(&rTable);
if (verbose) Table_Info(rTable);
rTableFlag = 1;
} else {
rTableFlag = 0;
}
if (transmit != NULL && strlen(transmit) && strcmp(transmit,"NULL") && strcmp(transmit,"0"))
{
if (verbose) fprintf(stdout, "Monochromator_curved: %s: Transmission data (k, T) from %s\n", NAME_CURRENT_COMP, transmit);
Table_Read(&tTable, transmit, 1); /* read 1st block data from file into rTable */
Table_Rebin(&tTable); /* rebin as evenly, increasing array */
if (tTable.rows < 2) Table_Free(&tTable);
if (verbose) Table_Info(tTable);
tTableFlag = 1;
} else {
tTableFlag = 0;
}
if (width == 0) SlabWidth = zwidth;
else SlabWidth = (width+gap)/NH - gap;
if (height == 0) SlabHeight = yheight;
else SlabHeight = (height+gap)/NV - gap;
tiltH=calloc((int)2*(NH+1),sizeof(double));
tiltV=calloc((int)2*(NV+1),sizeof(double));
if (!tiltH) printf("Monochromator_curved: %s: Warning: not enough memory to allocate tilts (NH=%g).\n", NAME_CURRENT_COMP, NH);
else if (RH) { /* pre-compute tilts */
for (i=0;i<=NH;i++){
tiltH[i]=asin((i-(NH+1)/2.0)*(SlabWidth+gap)/RH);
}
}
if (!tiltV) printf("Monochromator_curved: %s: Warning: not enough memory to allocate tilts (NV=%g).\n", NAME_CURRENT_COMP, NV);
else if (RV) {
for (i=0;i<=NV;i++){
tiltV[i]=-asin((i-(NV+1)/2.0)*(SlabHeight+gap)/RV);
}
}
%}
TRACE
%{
double dt;
double Gauss_X[] = {-0.987992518020485, -0.937273392400706, -0.848206583410427,
-0.724417731360170, -0.570972172608539, -0.394151347077563,
-0.201194093997435, 0, 0.201194093997435,
0.394151347077563, 0.570972172608539, 0.724417731360170,
0.848206583410427, 0.937273392400706, 0.987992518020485};
double Gauss_W[] = {0.030753241996117, 0.070366047488108, 0.107159220467172,
0.139570677926154, 0.166269205816994, 0.186161000115562,
0.198431485327111, 0.202578241925561, 0.198431485327111,
0.186161000115562, 0.166269205816994, 0.139570677926154,
0.107159220467172, 0.070366047488108, 0.030753241996117};
if(vx != 0.0 && (dt = -x/vx) >= 0.0 && r0)
{ /* Moving towards crystal? */
double zmin,zmax, ymin,ymax;
double yy,zz;
zmax = ((NH*(SlabWidth+gap))-gap)/2;
zmin = -zmax;
ymax = ((NV*(SlabHeight+gap))-gap)/2;
ymin = -ymax;
/* Test-propagate to crystal plane */
zz=z+vz*dt;
yy=y+vy*dt;
if (zz>zmin && zzymin && yy 2*k/mono_Q) Q_order--;
if((!order && Q_order > 0) || (Q_order == fabs(order) && order)) { /* Bragg scattering possible? */
double q0, q0x, theta, delta, p_reflect, my_r0;
q0 = Q_order*mono_Q;
q0x = ratio < 0 ? -q0 : q0;
theta = asin(q0/(2*k)); /* Actual bragg angle */
/* Make MC choice: reflect or transmit? */
delta = asin(fabs(kux)) - theta;
if (rTableFlag)
{
my_r0 = r0*Table_Value(rTable, k, 1); /* 2nd column */
}
else my_r0 = r0;
if (my_r0 > 1)
{
if (my_r0 > 1.01 && verbose)
fprintf(stdout, "Warning: Monochromator_curved : lowered reflectivity from %f to 1 (k=%f)\n", my_r0, k);
my_r0=0.999;
}
if (my_r0 < 0)
{
if (verbose)
fprintf(stdout, "Warning: Monochromator_curved : raised reflectivity from %f to 0 (k=%f)\n", my_r0, k);
my_r0=0;
}
p_reflect = fabs(my_r0)*exp(-kiz*kiz/(kiy*kiy + kiz*kiz)*(delta*delta)/
(2*mos_rms_y*mos_rms_y))*
exp(-kiy*kiy/(kiy*kiy + kiz*kiz)*(delta*delta)/
(2*mos_rms_z*mos_rms_z));
double rr=rand01();
if(rr <= p_reflect) { /* Reflect */
double bx,by,bz,ax,ay,az,phi;
double cos_2theta,k_sin_2theta,cos_phi,sin_phi,q_x,q_y,q_z;
double total,c1x,c1y,c1z,w,mos_sample;
int i=0;
cos_2theta = cos(2*theta);
k_sin_2theta = k*sin(2*theta);
/* Get unit normal to plane containing ki and most probable kf */
vec_prod(bx, by, bz, kix, kiy, kiz, q0x, 0, 0);
NORM(bx,by,bz);
bx = bx * k_sin_2theta;
by = by * k_sin_2theta;
bz = bz * k_sin_2theta;
/* Get unit vector normal to ki and b */
vec_prod(ax, ay, az, bx, by, bz, kux, kuy, kuz);
/* Compute the total scattering probability at this ki */
total = 0;
/* Choose width of Gaussian distribution to sample the angle
* phi on the Debye-Scherrer cone for the scattered neutron.
* The radius of the Debye-Scherrer cone is smaller by a
* factor 1/cos(theta) than the radius of the (partial) sphere
* describing the possible orientations of Q due to mosaicity, so we
* start with a width 1/cos(theta) greater than the largest of
* the two mosaics. */
mos_sample = mos_rms_max/cos(theta);
c1x = kix*(cos_2theta-1);
c1y = kiy*(cos_2theta-1);
c1z = kiz*(cos_2theta-1);
/* Loop, repeatedly reducing the sample width until it is small
* enough to avoid sampling scattering directions with
* ridiculously low scattering probability.
* Use a cut-off at 5 times the gauss width for considering
* scattering probability as well as for integration limits
* when integrating the sampled distribution below. */
for(i=0; i<100; i++) {
w = 5*mos_sample;
cos_phi = cos(w);
sin_phi = sin(w);
q_x = c1x + cos_phi*ax + sin_phi*bx;
q_y = (c1y + cos_phi*ay + sin_phi*by)/mos_rms_z;
q_z = (c1z + cos_phi*az + sin_phi*bz)/mos_rms_y;
/* Stop when we get near a factor of 25=5^2. */
if(q_z*q_z + q_y*q_y < (25/(2.0/3.0))*(q_x*q_x))
break;
mos_sample *= (2.0/3.0);
}
/* Now integrate the chosen sampling distribution, using a
* cut-off at five times sigma. */
for(i = 0; i < (sizeof(Gauss_X)/sizeof(double)); i++)
{
phi = w*Gauss_X[i];
cos_phi = cos(phi);
sin_phi = sin(phi);
q_x = c1x + cos_phi*ax + sin_phi*bx;
q_y = c1y + cos_phi*ay + sin_phi*by;
q_z = c1z + cos_phi*az + sin_phi*bz;
p_reflect = GAUSS_monocurved((q_z/q_x),0,mos_rms_y)*
GAUSS_monocurved((q_y/q_x),0,mos_rms_z);
total += Gauss_W[i]*p_reflect;
}
total *= w;
/* Choose point on Debye-Scherrer cone. Sample from a Gaussian of
* width 1/cos(theta) greater than the mosaic and correct for any
* error by adjusting the neutron weight later. */
phi = mos_sample*randnorm();
/* Compute final wave vector kf and scattering vector q = ki - kf */
cos_phi = cos(phi);
sin_phi = sin(phi);
q_x = c1x + cos_phi*ax + sin_phi*bx;
q_y = c1y + cos_phi*ay + sin_phi*by;
q_z = c1z + cos_phi*az + sin_phi*bz;
p_reflect = GAUSS_monocurved((q_z/q_x),0,mos_rms_y)*
GAUSS_monocurved((q_y/q_x),0,mos_rms_z);
vx = K2V*(kix+q_x);
vy = K2V*(kiy+q_y);
vz = K2V*(kiz+q_z);
p_reflect /= total*GAUSS_monocurved(phi,0,mos_sample);
if (p_reflect <= 0) ABSORB;
if (p_reflect > 1) p_reflect = 1;
p = p * p_reflect;
} /* End MC choice to reflect or transmit neutron (if tmp 1)
{
if (my_t0 > 1.01 && verbose)
fprintf(stdout, "Warning: Monochromator_curved : lowered transmission from %f to 1 (k=%f)\n", my_t0, k);
my_t0=0.999;
}
if (my_t0 > 0) p = p * my_t0;
else ABSORB;
}
} /* end if not in gap */
/* rotate back in component frame */
Rotation TT;
rot_transpose(T, TT);
/* now make the coordinate system change */
mccoordschange_polarisation(TT, &vx, &vy, &vz);
coords_get(rot_apply(TT,coords_set(x,y,z)),&x,&y,&z);
y= y+center_y;
z= z+center_z;
/* Visualise scattering point in proper, component frame
- but only if the neutron is reflected, that is none of:
* transmitted
* falling outside the slab material */
if(!do_transmit) SCATTER;
/* mccoordschange_polarisation(tt, &sx, &sy, &sz); */
} /* End intersect the crystal (if z) */
else {
/* restore neutron state when no interaction */
RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);
}
} /* End neutron moving towards crystal (if vx)*/
%}
FINALLY
%{
if(rTableFlag){
Table_Free(&rTable);
}
if(tTableFlag){
Table_Free(&tTable);
}
if (tiltH) free(tiltH);
if (tiltV) free(tiltV);
%}
MCDISPLAY
%{
int ih;
for(ih = 0; ih < NH; ih++)
{
int iv;
for(iv = 0; iv < NV; iv++)
{
double zmin,zmax,ymin,ymax;
double xt, yt;
zmin = (SlabWidth+gap)*(ih-NH/2.0)+gap/2;
zmax = zmin+SlabWidth;
ymin = (SlabHeight+gap)*(iv-NV/2.0)+gap/2;
ymax = ymin+SlabHeight;
if (RH) xt = -(zmax*zmax - zmin*zmin)/RH/2;
else xt = 0;
if (RV) yt = -(ymax*ymax - ymin*ymin)/RV/2;
else yt = 0;
multiline(5, xt+yt, (double)ymin, (double)zmin,
xt-yt, (double)ymax, (double)zmin,
-xt-yt, (double)ymax, (double)zmax,
-xt+yt, (double)ymin, (double)zmax,
xt+yt, (double)ymin, (double)zmin);
}
}
%}
END