/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2008, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Mirror_Elliptic
*
* %I
* Written by: Sylvain Desert
* Date: 2007
* Origin: LLB
* Modified by: E. Farhi, uniformize parameter names (Jul 2008)
*
* Elliptical mirror.
*
* %D
* Models an elliptical mirror. The reflectivity profile is given by a 2-column reflectivity free
* text file with format [q(Angs-1) R(0-1)]. The component is centered on the ellipse.
*
* Example: Mirror_Elliptic(reflect="supermirror_m3.rfl", focus=6.6e-4,
* interfocus = 8.2, yheight = 0.0002, zmin=-3.24, zmax=-1.49)
*
*
* %P
* INPUT PARAMETERS:
* yheight: [m] height of the mirror
* focus: [m] focal length (m)
* interfocus: [m] Distance between the two elliptical focal points
* zmin: [m] z-axis coordinate of ellipse start
* zmax: [m] z-axis coordinate of ellipse end
* reflect: [q(Angs-1) R(0-1)] (str) Reflectivity file name. Format
* 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
*
* Example instrument file FocalisationMirrors.instr is available in the examples/ folder.
*
* %E
*******************************************************************************/
DEFINE COMPONENT Mirror_Elliptic
DEFINITION PARAMETERS ()
SETTING PARAMETERS (string reflect=0, focus=6.6e-4,interfocus=8.2, yheight=2e-4, zmin=0, zmax=0,
R0=0.99, Qc=0.0219, alpha=6.07, m=1.0, W=0.003)
OUTPUT PARAMETERS ()
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
#ifndef SIGN
#define SIGN(a) (a >= 0 ? (a == 0 ? 0 : 1) : -1)
#endif
%include "read_table-lib"
%include "ref-lib"
%}
DECLARE
%{
double beta1; /* ellipse parameters */
double alpha1;
double beta2; /* ellipse squared parameters */
double alpha2;
t_Table pTable;
int err;
%}
INITIALIZE
%{
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,"Mirror_Elliptic: %s: can not read file %s\n", NAME_CURRENT_COMP, reflect));
}
/* Calculation of ellipse parameters */
alpha1 = interfocus/2 +focus;
alpha2 = alpha1*alpha1;
beta2 = alpha2 - (interfocus*interfocus)/4;
beta1 = sqrt(beta2);
err = 0;
yheight/=2;
if(zmin==0&&zmax==0){
zmin = -alpha1;
zmax = alpha1;
}
else{
if(zmin>=zmax) exit(fprintf(stderr,"Mirror_Elliptic: %s: error definition zmin and zmax\n", NAME_CURRENT_COMP));
}
printf("Mirror_Elliptic: %s: alpha=%f alpha^2=%f beta=%f beta^2=%f\n",
NAME_CURRENT_COMP, alpha1,alpha2,beta1,beta2);
%}
TRACE
%{
double q, B;
double div,z1,x1,z2,x2;
double v;
double vx_2,vz_2;
int i=-1;
double oa,ob,ab,xa,za;
double angle;
double old_x;
double old_y;
double old_z;
double par[5] = {R0, Qc, alpha, m, W};
double a,b;
double delta;
/* First check if neutron has the right direction. */
if((vz != 0.0 && -z/vz >= 0) && x+beta1> 0)
{
i++;
old_z=z;
old_x=x;
old_y=y;
a=vx/vz;
b=x-a*z;
/* printf("\nx : %e / z : %f / y : %e\nvx : %e / vz : %e / vy : %e\na : %e / b : %f",x,z,y,vx,vz,a,b); */
/* Calculation of intersection with ellipse */
delta = sqrt(4*(a*a*b*b-(a*a+beta2/alpha2)*(b*b-beta2)));
/* printf("\nDELTA : %e",delta); */
z1 = (-2*a*b - delta)/(2*(a*a+beta2/alpha2));
z2 = (-2*a*b + delta)/(2*(a*a+beta2/alpha2));
x1 = a*z1+b;
x2 = a*z2+b;
/* printf("\nx1 : %f / z1 : %f\nx2 : %f / z2 : %f\n",x1,z1,x2,z2); */
/* Choose the right result */
if((z1>z2)&&(fabs(z1)0.001){
#pragma acc atomic
err = err + 1;
printf("Mirror_Elliptic: %s: x=%e z=%f X=%f (Absorb)",NAME_CURRENT_COMP,x,z,a*z+b);
ABSORB;
}
/* y calculation */
y+=vy*(z-old_z)/vz;
/*reflection*/
if(x<0 && fabs(y)<=yheight && z>=zmin && z<=zmax){
/*reflection angle in the plane xz*/
div = -atan(vx/vz);
angle = -atan((beta2*z)/(alpha2*x));
/*vx and vz calculation after reflection*/
v=sqrt(vx*vx+vz*vz);
vz = v*cos(2*angle+div);
vx = v*sin(2*angle+div);
/*
printf("reflection2D :\nv: %e / angle (tangeante) : %f / div : %f / incidence : %f\n",v,angle,div,2*angle+div);
printf("vx : %f /vz : %f\n",vx,vz);
*/
/*incidence angle in 3D*/
ob = sqrt((old_x-x)*(old_x-x)+(old_z-z)*(old_z-z));
xa = x-ob*cos(div+angle)*sin(angle);
za = z-ob*cos(div+angle)*cos(angle);
oa = sqrt((old_x-xa)*(old_x-xa)+(old_z-za)*(old_z-za));
ob = sqrt((old_x-x)*(old_x-x)+(old_y-y)*(old_y-y)+(old_z-z)*(old_z-z));
ab = sqrt((xa-x)*(xa-x)+(old_y-y)*(old_y-y)+(za-z)*(za-z));
angle = acos((-ab*ab-ob*ob+oa*oa)/(2*ab*ob));
/* printf("3D :\nxa : %f / za : %f\noa : %f / ob : %f / ab : %f\nangle : %f / v : %e\n",xa,za,oa,ob,ab,angle,v); */
v=sqrt(vx*vx+vy*vy+vz*vz);
q = fabs(2*sin(angle)*v*V2Q);
/* Reflectivity (see component Guide). */
if (reflect && strlen(reflect) && strcmp(reflect,"NULL") && strcmp(reflect,"0"))
TableReflecFunc(q, &pTable, &B);
else {
StdReflecFunc(q, par, &B);
}
if (B <= 0) { ABSORB; }
else p *= B;
}
else ABSORB;
SCATTER;
}
else{
ABSORB;
}
%}
FINALLY
%{
if(err!=0){
fprintf(stderr,"Mirror_Elliptic: %s: WARNING : %d neutrons absorbed for inadapted divergence !\n", NAME_CURRENT_COMP, err);
}
%}
MCDISPLAY
%{
double xi,zi,xf,zf,delta_z;
delta_z = (zmax-zmin)/99;
xi=-beta1*sqrt(1-zmin*zmin/alpha2);
line(xi,-yheight,zmin,xi,yheight,zmin);
zi=zmin;
/* printf("delta_z : %f / xi : %f / zi : %f\n",delta_z,xi,zi); */
do{
zf = zi + delta_z;
xf=-beta1*sqrt(1-zf*zf/alpha2);
line(xi,yheight,zi,xf,yheight,zf);
line(xf,yheight,zf,xf,-yheight,zf);
line(xf,-yheight,zf,xi,-yheight,zi);
xi=xf;
zi=zf;
} while(zf<=zmax);
%}
END