/******************************************************************************* * * 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