Skip to content

Commit

Permalink
Merge branch 'master' of https://github.com/McStasMcXtrace/McCode
Browse files Browse the repository at this point in the history
  • Loading branch information
@farhi2
farhi2 committed Dec 14, 2015
2 parents 0e1e88c + 618a80f commit 0920892
Show file tree
Hide file tree
Showing 2 changed files with 140 additions and 63 deletions.
198 changes: 136 additions & 62 deletions mcxtrace-comps/optics/Capillary.comp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,20 +11,29 @@
*
* %D
*
* A Capillary tube allowing for reflections along the tube. A material coating can be applied. No multilayer yet.
* A Capillary tube allowing for reflections along the tube. A material coating can be applied. Multilayer
* coatings may be handled by generating a reflectivity file (e.g. by IMD) and setting rtable=1.
* Waviness is implemented using the model described in
* Wang et.al., J. Appl. Phys., 1996
* where the grazing incidence angle $\theta$ is altered as
* $\theta' = \theta + \delta \theta \in [-min(theta,\Delta\theta,\Delta\theta]$
* This ensures that reflected rays will never be scattered into the capillary.
* \Delta\theta is the value specified by the parameter waviness.
*
* %P
* radius: [m] Radius of curvature.
* length: [m] Length of the unbent mirror.
* coating: [] Name of file containing the material data (i.e. f1 and f2) for the coating
* R0: [ ] Fixed constant reflectivity
*
* rtable: [ ] if nonzero, the coating file contains an E,theta paramterized matrix of raw reflectivities.
* waviness: [rad] The momentaneous waviness is uniformly distributed in the range [-waviness,waviness].
* longw: [ ] If non-zero, waviness is purely longitudinal in its nature.
* %E
***********************************************************************/

DEFINE COMPONENT Capillary
DEFINITION PARAMETERS (string coating="Be.txt")
SETTING PARAMETERS ( radius=1, length=0.2, R0=0)
DEFINITION PARAMETERS (string coating="Be.txt", longw=1)
SETTING PARAMETERS ( radius=1, length=0.2, R0=0, rtable=0, waviness=0)
OUTPUT PARAMETERS(prms,prmsp)
/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */

Expand All @@ -35,6 +44,8 @@ DECLARE
struct {
int Z;
double At, rho;
double E_max,E_min,E_step;
double theta_max,theta_min,theta_step;
t_Table T;
} prms,*prmsp;
%}
Expand All @@ -44,24 +55,41 @@ INITIALIZE
int status;

if (coating && !R0){
if ( coating && (status=Table_Read(&(prms.T),coating,0))==-1){
fprintf(stderr,"Error(%s): Could not parse file \"%s\"\n",NAME_CURRENT_COMP,coating);
exit(-1);
}
char **header_parsed;
header_parsed=Table_ParseHeader(prms.T.header,"Z","A[r]","rho","Z/A","sigma[a]",NULL);
if(header_parsed[2]){prms.rho=strtod(header_parsed[2],NULL);}
if(header_parsed[0]){prms.Z=strtod(header_parsed[0],NULL);}
if(header_parsed[1]){prms.At=strtod(header_parsed[1],NULL);}
prmsp=&prms;
char **header_parsed;
prmsp=&prms;
if ( coating && (status=Table_Read(&(prms.T),coating,0))==-1){
fprintf(stderr,"Error(%s): Could not parse file \"%s\"\n",NAME_CURRENT_COMP,coating);
exit(-1);
}
if (!rtable){
/*the given coating file is to be interpreted as a material data file*/
header_parsed=Table_ParseHeader(prms.T.header,"Z","A[r]","rho","Z/A","sigma[a]",NULL);
if(header_parsed[2]){prms.rho=strtod(header_parsed[2],NULL);}
if(header_parsed[0]){prms.Z=strtod(header_parsed[0],NULL);}
if(header_parsed[1]){prms.At=strtod(header_parsed[1],NULL);}
prmsp=&prms;
}else{
header_parsed = Table_ParseHeader(prms.T.header,
"E_min=", "E_max=", "E_step=",
"theta_min=", "theta_max=", "theta_step=",
NULL);
if (header_parsed[0] && header_parsed[1] && header_parsed[2] &&
header_parsed[3] && header_parsed[4] && header_parsed[5]){
sscanf(header_parsed[0], "%lf", &prmsp->E_min);
sscanf(header_parsed[1], "%lf", &prmsp->E_max);
sscanf(header_parsed[2], "%lf", &prmsp->E_step);
sscanf(header_parsed[3], "%lf", &prmsp->theta_min);
sscanf(header_parsed[4], "%lf", &prmsp->theta_max);
sscanf(header_parsed[5], "%lf", &prmsp->theta_step);
}
}
}else{
if (R0<0 || R0>1){
fprintf(stderr,"Error(%s) reflectivity (%g) is specified but is not in [0:1]\n",NAME_CURRENT_COMP,R0);
exit(-1);
}
prmsp=NULL;
if (R0<0 || R0>1){
fprintf(stderr,"Error(%s) reflectivity (%g) is specified but is not in [0:1]\n",NAME_CURRENT_COMP,R0);
exit(-1);
}
prmsp=NULL;
}

%}

TRACE
Expand All @@ -72,6 +100,8 @@ TRACE
/*Do we hit the cylindrical aperture.*/
PROP_Z0;
if(x*x+y*y <radius*radius){
k=sqrt(scalar_prod(kx,ky,kz,kx,ky,kz));

hit=cylinder_intersect(&l0,&l1,y,z-length*0.5,x,ky,kz,kx,radius,length);
while (hit){
/*if x-ray exits through the cylinder top break from loop.*/
Expand All @@ -83,56 +113,100 @@ TRACE
nx=x;ny=y;nz=0;
NORM(nx,ny,nz);
s=scalar_prod(kx,ky,kz,nx,ny,0);

/*if we have waviness alter the normal vector slightly*/
if(waviness!=0){
double theta=M_PI_2-acos(s/k); /*pi_2 since theta is supposed to be the grazing angle*/
/*assuming theta to be small we might disregard atan*/
if(longw){
double dtheta;
if(theta<waviness){
dtheta=rand01()*(theta+waviness)-theta;
}else{
dtheta=randpm1()*waviness;
}
nz=atan(dtheta);
NORM(nx,ny,nz);
}else{
/*waviness is also transversal but anisotropic*/
double radius;
if(theta<waviness){
radius=atan(waviness);
randvec_target_circle(&nx,&ny,&nz,NULL,nx,ny,0,radius);
}else{
radius=(atan(theta)+atan(waviness))/2.0;
randvec_target_circle(&nx,&ny,&nz,NULL,nx,ny,radius-atan(theta),radius);
}
NORM(nx,ny,nz);
}

/*recompute the scalar prod s*/
s=scalar_prod(kx,ky,kz,nx,ny,nz);
}

if(s!=0){
kx-=2*s*nx;
ky-=2*s*ny;
kz-=2*s*nz;
}

if (!prmsp){
/*apply constant reflectivity*/
p*=R0;//pow(R0,scatterc);
}else if (!rtable){
/*compute reflectivity from material data*/

/*adjust p according to reflectivity*/
double Qc,Q,f1,f2,delta,beta,na,e,k2,rho;
/*length of wavevector transfer may be compute from s and n_ above*/
Q=fabs(2*s*sqrt(nx*nx+ny*ny+nz*nz));

/*interpolate in material data*/
e=K2E*k;
f1=Table_Value(prms.T,e,1);

/*the conversion factor in the end is to transform the atomic density from cm^-3 to AA^-3
-> therefore we get Q in AA^-1*/
na=NA*prms.rho/prms.At*1e-24;
Qc=4*sqrt(M_PI*na*RE*f1);
double R=1;
if (Q>Qc){
complex double qp;
double q,b_mu;

q=Q/Qc;
/*delta=na*r0*2*M_PI/k2*f1;*/
f2=Table_Value(prms.T,e,2);
/*beta=na*r0*2*M_PI/k2*f2;*/
/*b_mu=beta*(2*k)^2 / Qc^2*/
b_mu=4*M_PI*na*RE*f2/(Qc*Qc);
if(q==1){
qp=sqrt(b_mu)*(1+I);
}else {
qp=csqrt(q*q-1+2*I*b_mu);
}
/*and from this compute the reflectivity*/
R=cabs((q-qp)/(q+qp));
p*=R;//pow(R,scatterc);
/*now also set the phase*/
phi+=carg((q-qp)/(q+qp));//*scatterc;
}
//printf("E= %g R= %g Q= %g Qc= %g\n",e,R,Q,Qc);
}else{
double e,R,theta;
e=k*K2E;
theta=RAD2DEG*(M_PI_2-acos(s/k));
/*find reflectivity by interpolation*/
R=Table_Value2d(prmsp->T, (e-prmsp->E_min)/prmsp->E_step, (theta-prmsp->theta_min)/prmsp->theta_step);
/*apply interpolated value*/
//printf("E= %g R= %g theta= %g norm_e= %g norm_th= %g\n",e,R,theta,(e-prmsp->E_min)/prmsp->E_step, (theta-prmsp->theta_min)/prmsp->theta_step );
p*=R;//pow(R,scatterc);
}

scatterc++;
SCATTER;
hit=cylinder_intersect(&l0,&l1,y,z-length*0.5,x,ky,kz,kx,radius,length);
}
if (!prmsp){
/*apply constant reflectivity*/
p*=pow(R0,scatterc);
}else{
/*compute reflectivity from material data*/

/*adjust p according to reflectivity*/
double Qc,Q,f1,f2,delta,beta,na,e,k2,rho;
/*length of wavevector transfer may be compute from s and n_ above*/
Q=fabs(2*s*sqrt(nx*nx+ny*ny));

/*interpolate in material data*/
k2=scalar_prod(kx,ky,kz,kx,ky,kz);
e=K2E*sqrt(k2);
f1=Table_Value(prms.T,e,1);

/*the conversion factor in the end is to transform the atomic density from cm^-3 to AA^-3
-> therefore we get Q in AA^-1*/
na=NA*prms.rho/prms.At*1e-24;
Qc=4*sqrt(M_PI*na*RE*f1);
if (Q>Qc){
complex double qp;
double q,b_mu,R;

q=Q/Qc;
/*delta=na*r0*2*M_PI/k2*f1;*/
f2=Table_Value(prms.T,e,2);
/*beta=na*r0*2*M_PI/k2*f2;*/
/*b_mu=beta*(2*k)^2 / Qc^2*/
b_mu=4*M_PI*na*RE*f2/(Qc*Qc);
if(q==1){
qp=sqrt(b_mu)*(1+I);
}else {
qp=csqrt(q*q-1+2*I*b_mu);
}
/*and from this compute the reflectivity*/
R=cabs((q-qp)/(q+qp));
p*=pow(R,scatterc);
/*now also set the phase*/
phi+=carg((q-qp)/(q+qp))*scatterc;
}
}
}else{
ABSORB;
//RESTORE_XRAY(INDEX_CURRENT_COMP,x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p);
Expand Down
5 changes: 4 additions & 1 deletion tools/Python/mcplot/gnuplot/mcgnuplotter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -209,7 +209,10 @@ def __init__(self, input_file, noqt=False, log_scale=False):

elif file_ext == '.dat':
data_struct = get_monitor(input_file)
gpo = McGnuplot1D(data_struct['file'], data_struct, Gnuplot.Gnuplot(persist=gp_persist))
if re.search('array_2d.*', data_struct['type']):
gpo = McGnuplotPSD(data_struct['file'], data_struct, Gnuplot.Gnuplot(persist=gp_persist))
else:
gpo = McGnuplot1D(data_struct['file'], data_struct, Gnuplot.Gnuplot(persist=gp_persist))
self.__gnuplot_objs[gpo.key] = gpo
else:
raise Exception('McGnuPlotter: input file must be .sim or .dat')
Expand Down

0 comments on commit 0920892

Please sign in to comment.