Capillary.comp

/************************************************
*
* McXtrace, x-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* %Identification
* Written by: Erik B Knudsen
* Date: July 2015
* Origin: DTU Physics
* Release: McXtrace 1.2
*
* A capillary tube
*
* %Description
*
* 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 <span class="latex">$\theta$</span> is altered as
* <br>
* <div class="latex">
* $\theta' = \theta + \delta \theta \in [-min(theta,\Delta\theta,\Delta\theta]$
* </div>
* This ensures that reflected rays will never be scattered into the capillary.
* \Delta\theta is the value specified by the parameter waviness.
*
* Example: Capillary(
*       radius=1e-4,length=0.1, R0=0, coating="Rh.txt")
* 
* %Parameters
* radius:   [m]   Radius of curvature.
* length:   [m]   Length of the unbent mirror.
* coating:  [str] Name of file containing the material data (i.e. f1 and f2) for the coating
* R0:       [0-1] Fixed constant reflectivity
* rtable:   [0/1] If nonzero, the coating file contains an E,theta parameterized matrix of raw reflectivities.
* waviness: [rad] The momentaneous waviness is uniformly distributed in the range [-waviness,waviness].
* longw:    [0/1] If non-zero, waviness is purely longitudinal in its nature.
* %End
***********************************************************************/

DEFINE COMPONENT Capillary

SETTING PARAMETERS ( string coating="Be.txt", longw=1, radius=1, length=0.2, R0=0, rtable=0, waviness=0)

/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */ 
SHARE
%{
#include <complex.h>
#include <math.h>
%include "read_table-lib"
%}
DECLARE
%{
  int Z;
  double At;
  double rho;
  double E_max;
  double E_min;
  double E_step;
  double theta_max;
  double theta_min;
  double theta_step;
  t_Table T;
  int constant_R;
%}

INITIALIZE
%{
  int status;

  if (coating && !R0){
      char **header_parsed;
      if ( coating && (status=Table_Read(&(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(T.header,"Z","A[r]","rho","Z/A","sigma[a]",NULL);
          if(header_parsed[2]){rho=strtod(header_parsed[2],NULL);}
          if(header_parsed[0]){Z=strtod(header_parsed[0],NULL);}
          if(header_parsed[1]){At=strtod(header_parsed[1],NULL);}
      }else{
          header_parsed = Table_ParseHeader(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", &E_min);
              sscanf(header_parsed[1], "%lf", &E_max);
              sscanf(header_parsed[2], "%lf", &E_step);
              sscanf(header_parsed[3], "%lf", &theta_min);
              sscanf(header_parsed[4], "%lf", &theta_max);
              sscanf(header_parsed[5], "%lf", &theta_step);
          }
      }
      constant_R=0;
  }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);
      }
      constant_R=1;
  }
%}

TRACE
%{
    double l0,l1,dl,alpha,n,nx,ny,nz,s,k,knx,knz;
    int hit,scatterc=0;

    /*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.*/
            if (hit & 010) {
                break;
            }
            PROP_DL(l1);
                
            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 scat_radius;
                    if(theta<waviness){
                        scat_radius=atan(waviness);
                        randvec_target_circle(&nx,&ny,&nz,NULL,nx,ny,0,scat_radius);
                    }else{
                        scat_radius=(atan(theta)+atan(waviness))/2.0;
                        randvec_target_circle(&nx,&ny,&nz,NULL,nx,ny,scat_radius-atan(theta),scat_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 (constant_R){
                /*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;
                /*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(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*rho/At*1e-24;
                Qc=4*sqrt(M_PI*na*RE*f1);
                double R=1;
                if (Q>Qc){
                    double complex qp,Rc;
                    double q,b_mu;

                    q=Q/Qc;
                    /*delta=na*r0*2*M_PI/k2*f1;*/
                    f2=Table_Value(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*/
                    Rc=(q-qp)/(q+qp);
		    R=cabs(Rc);
                    p*=R;//pow(R,scatterc);
                    /*now also set the phase - have to compute manually since carg is not GPU yet*/
                    double arg = atan2(cimag(Rc),creal(Rc));
		    phi= phi + arg;//*scatterc;
                }
            }else{
                double e,R,theta;
                e=k*K2E;
                theta=RAD2DEG*(M_PI_2-acos(s/k));
                /*find reflectivity by interpolation*/
                R=Table_Value2d(T, (e-E_min)/E_step, (theta-theta_min)/theta_step);
                /*apply interpolated value*/
                p*=R;
            }

            scatterc++;
            SCATTER;
            hit=cylinder_intersect(&l0,&l1,y,z-length*0.5,x,ky,kz,kx,radius,length);
        }
    }else{
        ABSORB;
        //RESTORE_XRAY(INDEX_CURRENT_COMP,x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p);
    }
%}

MCDISPLAY
%{
    circle("xy",0,0,0,radius);
    circle("xy",0,0,length,radius);
    line(radius,0,0,radius,0,length);
    line(0,radius,0,0,radius,length);
    line(-radius,0,0,-radius,0,length);
    line(0,-radius,0,0,-radius,length);
%}

END