mcxtrace-comps-1.2-src/sources/Source_lab.comp

/*******************************************************************************
*
* McXtrace, X-ray tracing package
*         Copyright, All rights reserved
*         Technical University of Denmark, Kgs. Lyngby, Denmark
*         Institut Laue Langevin, Grenoble, France
*         University of Copenhagen, Copenhagen, Denmark
*
* Component: Source_lab
*
* %I
* Written by: Erik Bergbaeck Knudsen
* Date: May 2012
* Version: 1.0
* Release: McXtrace 1.1
* Origin: Kgs. Lyngby
*
* Laboratory x-ray source.
*
* %D
* Model of a laboratory x-ray tube, generating x-rays by bombarding a target by electrons.
* Given a input energy E0 of the electron beam, x-rays are emitted from the accessible emission lines
* The geometry of the tube is assumed to be:
* # The electron beam hits a slab of surface material surface at a right angle illuminating an area of width by height,
* # where width is measured along the component X-axis.
* # The centre of the electron beam at the anode surface is the origin of the component.
* # The Z-axis of the component points at the centre of the exit window (focus_xw by focus yh)
* placed at a distance dist from the origin.
* # The angle between the Z-axis and the anode surface is the take_off angle.
* For a detailed sketch of the geometry see the componnent manual.
*
* The Bremsstrahlung emitted is modelled using the model of Kramer (1923) as restated in International
* Tables of Crystallography C 4.1
* Characteristic radiation is modelled by Gaussian energy profiles with line-energies from Bearden (1967), widths from
* Krause (1979) and intensity from Honkimäki (1990) and x-ray data booklet.
* Absoprtion of emitted x-rays while travelling through the target anode is included.
*
* Example: Source_lab(material_datafile="Cu.txt",Emin=1, E0=80)
*
* %P
* width:    [m] Width of electron beam impinging on the anode.
* height:    [m] Height of electron beam impinging on the anode.
* thickness: [m] Thickness of the anode material slab.
* take_off: [deg] Take off angle of beam centre.
* dist:  [m] Distance between centre of illuminated target and exit window.
* E0:    [kV] Acceleration voltage of xray tube.
* tube_current: [A] Electron beam current.
* Emax:  [keV] Maximum energy to sample. Default (Emax=0) is to set it to E0.
* Emin:  [keV] Minimum energy to sample.
* focus_xw: [m] Width of exit window.
* focus_yh:[m] Height of exit window.
* frac:  [] Fraction of statistic to use for Bremsstrahlung.
* material_datafile: [] Name of datafile which describes the target material.
*
* %E
*************************************************************************/

DEFINE COMPONENT Source_lab
DEFINITION PARAMETERS (string material_datafile="Cu.txt")
SETTING PARAMETERS (width=1e-3, height=1e-3, thickness=100e-6, E0=20, Emax=0, Emin=1, focus_xw=5e-3, focus_yh=5e-3,
      take_off=6, dist=1, tube_current=1e-3, frac=0.5)
OUTPUT PARAMETERS (R_xray_gen,R_xray_geni, O_xray_gen, prms,pmul_c,mu,p_continous)
/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */

SHARE
%{
  %include "read_table-lib"

  /*here are some material data- currently only for Cu, Mo and W*/
  const struct xray_em_data{
    int Z;/*atom number*/
    double Ek;/*ionazation energy*/
    double w_k;/*flourescence yield*/
    int linec;
    double e[6];/*line energy*/
    double w[6];/*natural width of line FWHM*/
    double i[6];/*relative intensity*/
  } xray_mat_data[]={
    {29,8.979 ,0.407,2,{8.02783,8.04778,0,0,0,0},{2.11e-3,2.17e-3,0,0,0,0},{0.51,1.0,0,0,0,0}},
    {31,10.367,0.0  ,2,{9.22482,9.25174,0,0,0,0},{2.59e-3,2.66e-3,0,0,0,0},{0.51,1.0,0,0,0,0}},
    {42,20.00 ,0.770,5,{17.3743,17.47934,19.5903, 19.6083, 19.965,0},{6.31e-3,6.49e-3,12e-3, 12e-3, 12e-3,0},{0.52,1.0,0.08,0.15,0.03,0}},
    {74,69.525,0.945,2,{57.9817,59.31824,0,0,0,0},{44.9e-3,45.2e-3,0,0,0,0},{0.58,1.0,0,0,0,0}},
    {0,0.0,0.0,0,{0,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,0}}
  };

%}

DECLARE
%{
  Rotation R_xray_gen,R_xray_geni;
  Coords O_xray_gen;
  const double BKRAMER=2e-6; /*photons /keV /electron*/
  struct {
    int Z;
    double At,rho;
    t_Table T;
    const struct xray_em_data *em_p;
    double Icont,Ichar;
  } prms;
  double p_continous;
  double mu=0.2e-6;
  double pmul_c;
%}

INITIALIZE
%{
  int status,ii;
  if (E0<=0){
    fprintf(stderr,"Error %s: Impinging electron energy (E0) must be >0, was %g\n",NAME_CURRENT_COMP, E0);
    exit(-1);
  }

  if (!Emax){/*if Emax is not set use the impinging electron energy*/
    Emax=E0;
  }
  if(Emin<=0){
    fprintf(stderr,"Error (%s): Emin must be > 0 (%g)\n",NAME_CURRENT_COMP,Emin);exit(-1);
  }
  if(Emax<Emin){
    fprintf(stderr,"Error (%s): Nonsensical emission energy interval [Emin,Emax]=[%g %g] at E0=%g\n",NAME_CURRENT_COMP,Emin,Emax,E0);
    exit(-1);
  }

  if ( (status=Table_Read(&(prms.T),material_datafile,0))==-1){
    fprintf(stderr,"Error %s: Could not parse file \"%s\"\n",NAME_CURRENT_COMP,material_datafile?material_datafile:"");
    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);}
  /*use the atom number to get at the right data structure*/
  prms.em_p=&(xray_mat_data[0]);
  while (prms.Z!=prms.em_p->Z){
    prms.em_p++;
    if ((prms.em_p->Z)==0){
      fprintf(stderr,"Error: %s (Z=%d) anode not implemented yet. Please contact the McXtrace team to fix this. Aborting.\n",material_datafile,prms.Z);
      exit(-1);
    }
  }

  /*Integrate the continuous spectrum and the characteristic so as to get the relative intenisities right*/
  prms.Icont=tube_current/CELE*BKRAMER*prms.Z*(E0*log(Emax)-E0*log(Emin) - Emax + Emin);
  /*check if E0 >Ek. If not, no characteristic emission can take place*/
  if (E0>prms.em_p->Ek){
    double Bk=1.2e-5*pow(prms.em_p->Ek,1.67)*exp(-0.077*prms.Z);
    prms.Ichar=tube_current/CELE*4*M_PI*(E0/prms.em_p->Ek-1)*Bk;
    double Ichar_tot=0;
    for (ii=0;ii<prms.em_p->linec;ii++){
      if(Emax>prms.em_p->e[ii]+5*prms.em_p->w[ii] && Emin<prms.em_p->e[ii]-5*prms.em_p->w[ii]){
        Ichar_tot+=prms.em_p->i[ii];
      }else{
        Ichar_tot+=prms.em_p->i[ii]*0.5*( erf(M_SQRT1_2*(Emax-prms.em_p->e[ii])/prms.em_p->w[ii]) - erf(M_SQRT1_2*(Emin-prms.em_p->e[ii])/prms.em_p->w[ii]) );
      }
    }
    p_continous=prms.Icont/(prms.Ichar*Ichar_tot+prms.Icont);
    printf("Bk,Ico,Ich,Icht,p_c=%g %g %g %g %g\n",Bk,prms.Icont,prms.Ichar,Ichar_tot,p_continous);
  }else{
    /*characteristic K-emission is not possible*/
    p_continous=1;
    frac=1;
  }
  O_xray_gen=coords_set(0,0,0);
  rot_set_rotation(R_xray_gen,-take_off*DEG2RAD,0,0);
  rot_set_rotation(R_xray_geni,take_off*DEG2RAD,0,0);

  pmul_c=1.0/(mcget_ncount());

%}

TRACE
%{
  double x1,y1,z1,x2,y2,z2,r,e,k,pdir,pmul;
  /* pick a point in the generating volume*/
  x1=rand01()*width-width/2.0;
  z1=rand01()*height-height/2.0;
  /* y is the absorption depth of the electron getting converted to an xray*/
  y1=log(rand01())*mu;

  double px,py,pz;
  Coords P;

  /* transform initial coords to ones relative to the exit window which is our reference point*/
  P=coords_set(x1,y1,z1);
  P=coords_add(rot_apply(R_xray_gen,P),O_xray_gen);
  coords_get(P,&x,&y,&z);

  //mccoordschange(O_xray_gen,R_xray_gen,&x1,&y1,&z1, 0,0,0, 0,0,0);

  /*set a scatter pt at the generation pt*/
  //x=x1;y=y1;z=z1;
  SCATTER;

  /*randvec_target_rect_real computes a target point and a solid angle correction factor, hence the k-vector has to be computed from
    generation point and target point. The (0,0,1) location of the target is due to a silent assumption in randvec() that
    the target cannot be situated in the origin.*/
  randvec_target_rect_real(&px,&py,&pz,&pdir,0,0,dist,focus_xw,focus_yh, R_xray_gen,x1,y1,z1,2);
  /*k is parallell to the line between generation and target points*/
  kx=px-x1;
  ky=py-y1;
  kz=pz-z1;

  /*Now for wavelength selection*/
  r=rand01();
  //printf("brems %g %g\n",r,p_continous);
  if(r<frac){
    //printf("brems %g \n",r);
    /*bremsstrahlung*/
    double e=rand01()*(Emax-Emin)+Emin;
    k=e*E2K;
    pmul=tube_current/CELE*BKRAMER*prms.Z*(E0/e-1);
    /*correct for not having the full E-window*/
    pmul*=(Emax-Emin)/E0;
    /*correct for monte-carlo statistics*/
    pmul*=p_continous/frac;
  }else{
    const struct xray_em_data *pt=prms.em_p;
    /*characteristic radiation*/
    /*first pick a possible line*/
    r=rand01()*pt->linec;
    int lineno=(int)floor(r);
    if (lineno==pt->linec) {
      lineno--;/*we might get overflow*/
    }
    pmul=pt->i[lineno]*prms.Ichar;
    k=E2K*(randnorm()*pt->w[lineno]+pt->e[lineno]);

    pmul*=(1-p_continous)/frac;
    /*correct for not having the full E-window if so*/

  }

  /*scale k accordingly*/
  NORM(kx,ky,kz);
  kx*=k;ky*=k;kz*=k;

  /*set the x-ray weight to whatever we computed just before and correct for only sampling the exit window, and correct for number of issued photons*/
  p=pmul*pmul_c;

  double mu_abs,l0;
  int ie;
  /*Correct for absorption*/
  if (take_off<0){
    double px,py,pz,nx,ny,nz;
    coords_get(rot_apply(R_xray_gen,coords_set(0,1,0)),&nx,&ny,&nz);
    coords_get(rot_apply(R_xray_gen,coords_set(0,-thickness,0)),&px,&py,&pz);
    /*this means we exit through the bottom of the anode slab*/
    ie=plane_intersect(&l0,x,y,z,kx,ky,kz,nx,ny,nz,px,py,pz);
  }else{
    double nx,ny,nz;
    coords_get(rot_apply(R_xray_gen,coords_set(0,1,0)),&nx,&ny,&nz);
    /*exit through top of anode slab*/
    ie=plane_intersect(&l0,x,y,z,kx,ky,kz,nx,ny,nz,0,0,0);
    if(l0<0){
    }
  }
  mu_abs=Table_Value(prms.T, k*K2E, 5)*prms.rho*1e2; /*mu_abs in m^-1*/
  p*=exp(-l0*mu_abs);
  //printf("%g %g %g %g\n",k*K2E,mu_abs,exp(-l0*mu_abs),Table_Value(prms.T, k*K2E, 5) );
  /*set a random phase - and propagate to the exit window*/
  phi=rand01()*2*M_PI;
%}

MCDISPLAY
%{
  magnify("xy");
  double x1,y1,z1,x2,y2,z2,width_2,height_2;
  double dx,dy,dz;
  /*these are just dummies*/
  double d1,d2,d3,d4,d5,d6;

  width_2=width/2.0;
  height_2=height/2.0;
  x1=-width_2;y1=0;z1=-height_2;
  x2=-width_2;y2=0;z2= height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&x1,&y1,&z1,&d1,&d2,&d3,&d4,&d5,&d6);
  mccoordschange(O_xray_gen,R_xray_gen,&x2,&y2,&z2,&d1,&d2,&d3,&d4,&d5,&d6);
  line(x1,y1,z1,x2,y2,z2);

  x1=width_2;y1=0;z1=height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&x1,&y1,&z1,&d1,&d2,&d3,&d4,&d5,&d6);
  line(x2,y2,z2,x1,y1,z1);

  x2=width_2;y2=0;z2=-height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&x2,&y2,&z2,&d1,&d2,&d3,&d4,&d5,&d6);
  line(x1,y1,z1,x2,y2,z2);

  x1=-width_2;y1=0;z1=-height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&x1,&y1,&z1,&d1,&d2,&d3,&d4,&d5,&d6);
  line(x2,y2,z2,x1,y1,z1);

  /*this is the mean penetration depth of electron that get converted to x-rays*/
  x1=-width_2;y1=-mu;z1=-height_2;
  x2=-width_2;y2=-mu;z2= height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&x1,&y1,&z1,&d1,&d2,&d3,&d4,&d5,&d6);
  mccoordschange(O_xray_gen,R_xray_gen,&x2,&y2,&z2,&d1,&d2,&d3,&d4,&d5,&d6);
  dashed_line(x1,y1,z1,x2,y2,z2,5);
  x1=width_2;y1=-mu;z1=height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&x1,&y1,&z1,&d1,&d2,&d3,&d4,&d5,&d6);
  dashed_line(x2,y2,z2,x1,y1,z1,5);
  x2= width_2;y2=-mu;z2=-height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&x2,&y2,&z2,&d1,&d2,&d3,&d4,&d5,&d6);
  dashed_line(x1,y1,z1,x2,y2,z2,5);
  x1=-width_2;y1=-mu;z1=-height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&x1,&y1,&z1,&d1,&d2,&d3,&d4,&d5,&d6);
  dashed_line(x2,y2,z2,x1,y1,z1,5);

  x1=-width_2;y1=-mu;z1=-height_2;
  x2=-width_2;y2=0;z2=-height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&x1,&y1,&z1,&d1,&d2,&d3,&d4,&d5,&d6);
  mccoordschange(O_xray_gen,R_xray_gen,&x2,&y2,&z2,&d1,&d2,&d3,&d4,&d5,&d6);
  line(x2,y2,z2,x1,y1,z1);
  x1=width_2;y1=-mu;z1=-height_2;
  x2=width_2;y2=0;z2=-height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&x1,&y1,&z1,&d1,&d2,&d3,&d4,&d5,&d6);
  mccoordschange(O_xray_gen,R_xray_gen,&x2,&y2,&z2,&d1,&d2,&d3,&d4,&d5,&d6);
  line(x2,y2,z2,x1,y1,z1);
  x1=-width_2;y1=-mu;z1=height_2;
  x2=-width_2;y2=0;z2=height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&x1,&y1,&z1,&d1,&d2,&d3,&d4,&d5,&d6);
  mccoordschange(O_xray_gen,R_xray_gen,&x2,&y2,&z2,&d1,&d2,&d3,&d4,&d5,&d6);
  line(x2,y2,z2,x1,y1,z1);
  x1=width_2;y1=-mu;z1=height_2;
  x2=width_2;y2=0;z2=height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&x1,&y1,&z1,&d1,&d2,&d3,&d4,&d5,&d6);
  mccoordschange(O_xray_gen,R_xray_gen,&x2,&y2,&z2,&d1,&d2,&d3,&d4,&d5,&d6);
  line(x2,y2,z2,x1,y1,z1);


  /*now draw "exit" window*/
  rectangle("xy",0,0,0,focus_xw,focus_yh);
%}

END