1/******************************************************************************* 2* 3* McXtrace, X-ray tracing package 4* Copyright, All rights reserved 5* Risoe National Laboratory, Roskilde, Denmark 6* Institut Laue Langevin, Grenoble, France 7* University of Copenhagen, Copenhagen, Denmark 8* 9* Component: Lens_parab 10* 11* % 12* 13* Written by: Jana Baltser and Erik Knudsen 14* 15* Date: August 2010, modified July 2011 16* Version: 1.0 17* Release: McXtrace 0.1 18* Origin: NBI 19* 20* %D 21* A simple X-ray compound refractive lens (CRL) with a profile of the parabola in rotation simulates the photons' movement on passing through it. The CRL focuses in 2D 22* 23* %P 24* Input parameters: 25* r - raduis of curvature (circular approximation at the tip of the profile) [mm]; 26* yheight - the CRL's dimensions along Y, aka aperture[mm]; 27* xwidth - the CRL's dimensions along X [mm]; 28* d - distance between two surfaces of the lens along the propagation axis; 29* N - amount of single lenses in a stack. 30* T - transmission of the lens 31* rough_xy [rms] - waviness along x and y 32* rough_z [rms] -waviness along z 33* 34* attenuation coefficient mu is taken from the NIST database and 35* Be.txt 36* 37* 38*******************************************************************************/ 39 40 41DEFINE COMPONENT Lens_parab_rough 42DEFINITION PARAMETERS (string material_datafile="Be.txt") 43SETTING PARAMETERS (r=0.5e-3,yheight=1.4e-3,xwidth=1.4e-3,d=.1e-3,T=.99,N=1,rough_z=0,rough_xy=0) 44OUTPUT PARAMETERS (prms,parab) 45/*STATE PARAMETERS (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p)*/ 46 47SHARE 48%{ 49 %include "read_table-lib" 50 struct datastruct{ 51 int Z; 52 double Ar, rho; 53 double *E, *mu, *f; 54 }; 55 typedef struct { 56 double coord[3]; 57 double k[3]; 58 } incom; 59typedef struct { 60 double constants[8]; 61} lens; 62 63 incom intersection_lens_parab_rough(incom a,lens b){ 64 incom result={a.coord[0],a.coord[1],a.coord[2],a.k[0],a.k[1],a.k[2]}; 65 int i; 66 double A,B,C,D,rr; 67 double t[2],p[3],knorm[3],k[3],pos1_tmp[3],pos_tmp[3]; 68 69 double nxn,nyn,nzn,Nx,Ny,Nz,NORM,Knorm; 70 double cos_theta,cos_theta1,Arg,Arg1,s,q; 71 double k_new[3],k_new1[3],M,Sign,dd; 72 double yh,xw; 73 74 double tx,ty,tz,tnorm,txn,tyn,tzn; 75 double v,w; 76 double rough_xy,rough_z; 77 78 for(i=0;i<=2;i++){ 79 k[i]=a.k[i]; 80 p[i]=a.coord[i]; 81 } 82 83 Knorm=sqrt(k[0]*k[0]+k[1]*k[1]+k[2]*k[2]); 84 knorm[0]=k[0]/Knorm; 85 knorm[1]=k[1]/Knorm; 86 knorm[2]=k[2]/Knorm; 87 88 rr=b.constants[0]; 89 yh=b.constants[1]; 90 xw=b.constants[2]; 91 dd=b.constants[3]; 92 M=b.constants[4]; 93 Sign=b.constants[5]; 94 95 rough_xy=b.constants[6]; 96 rough_z=b.constants[7]; 97 98 A=knorm[0]*knorm[0]+knorm[1]*knorm[1]; 99 B= 2.0*(p[0]*knorm[0]+p[1]*knorm[1] - Sign*rr*knorm[2]); 100 C=p[0]*p[0]+p[1]*p[1] - Sign*2*rr*p[2] + Sign*rr*2*dd; 101 D=B*B-4.0*A*C; 102 103 if (D<0) { /*ray does not intersect the parabola*/ 104 fprintf(stderr,"line 96: D<0 %s\n",NAME_CURRENT_COMP); 105 return result; 106 } 107 if (A==0){ /*incident k-vector is parallel (exactly) to the z-axis. Thus, the eq. becomes linear*/ 108 if(B==0){ 109 fprintf(stderr, "Division by zero in %s\n",NAME_CURRENT_COMP); 110 return result; 111 } 112 t[0]=-C/B; 113 for(i=0; i<=2; i++){ 114 result.coord[i]=p[i]+t[0]*knorm[i]; 115 } 116 } else { 117 double qq; 118 if (B<0){ 119 qq=-0.5*(B-sqrt(D)); 120 }else{ 121 qq=-0.5*(B+sqrt(D)); 122 } 123 t[0]=qq/A; 124 t[1]=C/qq; 125 126 for(i=0;i<=2;i++){ 127 pos_tmp[i]=p[i]+t[0]*knorm[i]; 128 pos1_tmp[i]=p[i]+t[1]*knorm[i]; 129 } 130 if ( fabs(pos_tmp[1])<=fabs(yh/2) && fabs(pos_tmp[0])<=fabs(xw/2) ){ 131 for(i=0;i<=2;i++){ 132 result.coord[i]=pos_tmp[i]; 133 } 134 } else if ( fabs(pos1_tmp[1])<=fabs(yh/2) && fabs(pos1_tmp[0])<=fabs(xw/2) ){ 135 for(i=0;i<=2;i++){ 136 result.coord[i]=pos1_tmp[i]; 137 } 138 } 139 else return result; 140 } 141 142 // introducing waviness into the code 143 double d_xy, d_z; 144 145 d_xy=rough_xy*randnorm(); 146 d_z=rough_z*randnorm(); 147 148 /* Calculating tangential vector */ 149 150 if (result.coord[0]==0 && result.coord[1]==0){ // incoming ray is along the axis, so it does not refract 151 k_new[0]=k[0]; 152 k_new[1]=k[1]; 153 k_new[2]=k[2]; 154 for(i=0;i<3;i++) { 155 result.k[i]=k_new[i]; 156 } 157 return result; 158 } 159 else if (result.coord[0]!=0 && result.coord[1]!=0){ 160 Nx=-Sign*(result.coord[0]/rr); // surface normal 161 Ny=-Sign*(result.coord[1]/rr); 162 Nz=1; 163 164 if (rough_xy) { 165 Nx+=d_xy; 166 Ny+=d_xy; 167 } 168 if (rough_z) { 169 Nz+=d_z; 170 } 171 172 173 NORM=sqrt(Nx*Nx+Ny*Ny+Nz*Nz); 174 nxn=Nx/NORM; 175 nyn=Ny/NORM; 176 nzn=Nz/NORM; 177 178 double cos_chi; 179 cos_chi=knorm[0]*nxn+knorm[1]*nyn+knorm[2]*nzn; 180 w=1/(sqrt(1-cos_chi*cos_chi)); // tangential vector 181 v=-w*cos_chi; 182 183 tx=v*nxn+w*knorm[0]; 184 ty=v*nyn+w*knorm[1]; 185 tz=v*nzn+w*knorm[2]; 186 187 } 188 else if (result.coord[0]==0){ 189 tx=0; 190 ty=Sign*(rr/result.coord[1]); 191 tz=1; 192 } 193 else if (result.coord[1]==0){ 194 tx=Sign*(rr/result.coord[0]); 195 ty=0; 196 tz=1; 197 } 198 199 tnorm=sqrt(tx*tx+ty*ty+tz*tz); 200 txn=tx/tnorm; 201 tyn=ty/tnorm; 202 tzn=tz/tnorm; 203 204 cos_theta=txn*knorm[0]+tyn*knorm[1]+tzn*knorm[2]; 205 cos_theta1=M*cos_theta; // Snell's law 206 207 /* new k vector */ 208 if ((1.0-cos_theta*cos_theta)==0) { 209 fprintf(stderr,"line 134: Division by zero\n"); 210 return result; 211 } 212 213 Arg=(1.0-cos_theta1*cos_theta1)/(1.0-cos_theta*cos_theta); 214 s=(1/M)*sqrt(Arg); 215 q=(Knorm/tnorm)*((1/M)*cos_theta1-s*cos_theta); 216 217 k_new[0]=q*tx+s*k[0]; 218 k_new[1]=q*ty+s*k[1]; 219 k_new[2]=q*tz+s*k[2]; 220 221 for(i=0;i<3;i++) { 222 result.k[i]=k_new[i]; 223 } 224 return result; 225 } 226 227 const double Re=2.8179402894e-5; /* Thomson Scattering length [Angstrom] */ 228 const double Na=6.02214179e23; /* Avogadro's number [atoms per gram-mole]*/ 229%} 230 231DECLARE 232%{ 233 struct datastruct *prms; 234%} 235 236INITIALIZE 237%{ 238 int status=0; 239 t_Table T; 240 if ( (status=Table_Read(&T,material_datafile,0))==-1){ 241 fprintf(stderr,"Error: Could not parse file \"%s\" in COMP %s\n",material_datafile,NAME_CURRENT_COMP); 242 exit(-1); 243 } 244 char **header_parsed; 245 header_parsed=Table_ParseHeader(T.header,"Z","A[r]","rho",NULL); 246 prms=calloc(1,sizeof(struct datastruct)); 247 if (!prms->Z) prms->Z=strtol(header_parsed[0],NULL,10); 248 if (!prms->Ar) prms->Ar=strtod(header_parsed[1],NULL); 249 if (!prms->rho) prms->rho=strtod(header_parsed[2],NULL); 250 prms->E=malloc(sizeof(double)*(T.rows+1)); 251 prms->f=malloc(sizeof(double)*(T.rows+1)); 252 prms->mu=malloc(sizeof(double)*(T.rows+1)); 253 int i; 254 for(i=0;i<T.rows;i++){ 255 prms->E[i]=T.data[i*T.columns]; 256 prms->mu[i]=T.data[5+i*T.columns]*prms->rho*1e2; /*mu is now in SI, [m^-1]*/ 257 prms->f[i]=T.data[1+i*T.columns]; 258 } 259 Table_Free(&T); 260%} 261 262 263 264 265TRACE 266%{ 267 incom incid,refr,outg; 268 lens parab; 269 270 double E,mu,f,rhoel,dl,e,k,delta,beta,Refractive_Index_Re,Refractive_Index_Im,w; 271 272 int i=0,nr; 273 274 parab.constants[0]=r; 275 parab.constants[1]=yheight; 276 parab.constants[2]=xwidth; 277 278 parab.constants[6]=rough_xy; 279 parab.constants[7]=rough_z; 280 281 w=(yheight*yheight)/(8.0*r); 282 k=sqrt(kx*kx+ky*ky+kz*kz); 283 e=K2E*k; /*Energy in KeV, same unit as datafile */ 284 285/*Interpolation of Table Values*/ 286 287 while (e>prms->E[i]){ 288 i++; 289 if (prms->E[i]==-1){ 290 fprintf(stderr,"Photon energy (%g keV) is outside of the lens' material data\n",e); ABSORB; 291 } 292 } 293 E=(e-prms->E[i-1])/(prms->E[i]-prms->E[i-1]); 294 mu=(1-E)*prms->mu[i-1]+E*prms->mu[i]; 295 //mu= 1e-10*mu; /*factor conversion from m^-1 to A^-1*/ 296 f=(1-E)*prms->f[i-1]+E*prms->f[i]; 297 298 /*Calculation of Refractive Index */ 299 300 rhoel= f*Na*(prms->rho*1e-24)/prms->Ar; /*Material's Number Density of Electrons [e/A^3] incl f' scattering length correction*/ 301 302 delta= 2.0*M_PI*Re*rhoel/(k*k); 303 //beta=mu/(2.0*k); /*mu and k in A^-1*/ 304 //printf("Delta=%g\n",delta); 305 Refractive_Index_Re = 1.0-delta; 306 //Refractive_Index_Im = beta; 307 308 /*Ray Tracing*/ 309 310 incid.k[0]=kx; 311 incid.k[1]=ky; 312 incid.k[2]=kz; 313 314 incid.coord[0]=x; 315 incid.coord[1]=y; 316 incid.coord[2]=z; 317 318 for (nr=0;nr<=(N-1);nr++){ 319 parab.constants[3]=nr*d+nr*2*w; // d constant 320 parab.constants[4]=1.0/Refractive_Index_Re; // M constant 321 parab.constants[5]=-1.0; // Sign constant 322 323 refr=intersection_lens_parab_rough(incid,parab); 324 325 if(refr.k[0]==0 && refr.k[1]==0 && refr.k[2]==0) continue; 326 327 dl=sqrt( (refr.coord[0]-x)*(refr.coord[0]-x) + (refr.coord[1]-y)*(refr.coord[1]-y) + (refr.coord[2]-z)*(refr.coord[2]-z) ); 328 PROP_DL(dl); 329 SCATTER; 330 331 kx=refr.k[0]; 332 ky=refr.k[1]; 333 kz=refr.k[2]; 334 //alter parabolic input to match second parabola 335 336 parab.constants[3]=(nr+1)*d+nr*2*w; 337 parab.constants[4]=Refractive_Index_Re; 338 parab.constants[5]=1.0; 339 340 outg=intersection_lens_parab_rough(refr,parab); 341 342 dl=sqrt( (outg.coord[0]-x)*(outg.coord[0]-x) + (outg.coord[1]-y)*(outg.coord[1]-y) + (outg.coord[2]-z)*(outg.coord[2]-z) ); 343 PROP_DL(dl); 344 SCATTER; 345 346 kx=outg.k[0]; ky=outg.k[1]; kz=outg.k[2]; 347 incid=outg; 348 } 349 // transmission calculation 350 double mu_rho, ap; 351 352 mu_rho=mu/(prms->rho*1e2); // mass absorption coefficient [cm2] 353 354 ap=mu_rho*((parab.constants[0]*N*delta*prms->Ar*k*k)/(M_PI*Na*Re*1e-10*(prms->Z+f)))*1e16; //1e16 -dimension coefficient 355 // ap - effective aperture 356 357 if (T==0) 358 ABSORB; 359 else 360 T=exp(-mu*N*d)*(1/(2*ap))*(1-exp(-2*ap)); 361 p*=T; 362 363%} 364 365MCDISPLAY 366%{ 367 magnify("xy"); 368 double z_c,zdepth,w; 369 w=(yheight*yheight)/(8*r); 370 zdepth=N*(2*w+d); 371 z_c=zdepth/2.0-w; 372 box(0,0,z_c,yheight/2,yheight/2,zdepth); 373%} 374 375END 376