/*******************************************************************************
*
* McXtrace, X-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: Source_lab
*
* %Identification
* Written by: Erik Bergbaeck Knudsen 
* Date: May 2012
* Version: 1.0
* Origin: Kgs. Lyngby
*
* Laboratory x-ray source.
*
* %Description
* 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 Lorentzian (default) or 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)
*
* %Parameters
* width:      [m]    Width of electron beam impinging on the anode.
* height:     [m]    Height of electron beam impinging on the anode.
* xwidth:     [m]    Width of the anode material slab.
* yheight:    [m]    Height of the anode material slab.
* 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:       [0-1]  Fraction of statistic to use for Bremsstrahlung.
* material_datafile: [string] Name of datafile which describes the target material.
* lorentzian: [0/1]  If nonzero Lorentzian (more correct) line profiles are used.
* exit_window_refpt: [m] If set, the AT position and exit window will coincide (legacy behaviour).
*
* %End
*************************************************************************/

DEFINE COMPONENT Source_lab

SETTING PARAMETERS (string material_datafile="Cu.txt", 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.1, lorentzian=1, xwidth=0, yheight=0, exit_window_refpt=0 )

/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */ 

SHARE
%{
  %include "read_table-lib"
#ifndef MX_SOURCE_LAB
#define MX_SOURCE_LAB
  /*here are some material data- currently only for Cu, Mo, W, and Ag*/
  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*/
  };
  struct xray_em_data xray_mat_data[8]={
    {24, 5.989 ,0.265, 3, {5.41472,5.40551,5.94671,0,0,0}, {1.97e-3,2.39e-3,3.05e-3,0,0,0}, {100,50,15,0,0,0}},
    {27, 7.709 ,0.391, 3, {6.93032,6.9153,7.6494,0,0,0}, {2.26e-3,3.08e-3,4.36e-3,0,0,0}, {100,51,17,0,0,0}},
    {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}},
    {47, 25.52, 0.8,2,{21.99030,22.162917, 0,0,0,0},{9.32e-3,9.16e-3,0,0,0,0}, {0.53, 1, 0,0,0,0}}, /* Ag fluorscence yield just guessed */
    {74, 69.525,0.945,5,{57.9817,59.31824,66.9514,67.2443,69.067,0},{44.9e-3,45.2e-3, 51.1e-3, 50.8e-3 , 50.2,0},{0.58,1.0,0.11,0.22,0.08,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}}

  };
#pragma acc declare create(xray_mat_data)
#endif

%}

DECLARE
%{
  Rotation R_xray_gen;
  Rotation R_xray_geni;
  Coords O_xray_gen;
  int Z;
  double At;
  double rho;
  t_Table T;
  int em_idx;
  double Icont;
  double Ichar;
  int linemin;
  int linemax;
  double p_continous;
  double pmul_c;
  double mu_electron;

  double BKRAMER;
%}

INITIALIZE
%{
  BKRAMER=2e-6; /*photons /keV /electron*/
#pragma acc update device(xray_mat_data[0:6])

  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(&(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(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);}

  /*use the atom number to get at the right data structure*/  
  int idx=0;
  while (Z!=xray_mat_data[idx]._Z){
    idx++;
    if ((xray_mat_data[idx]._Z)==0){
      fprintf(stderr,"Error: %s (Z=%d) anode not implemented yet. Please contact the McXtrace team to fix this. Aborting.\n",material_datafile,Z);
      exit(-1);
    }
  }
  em_idx=idx;
  struct xray_em_data *em_p = &(xray_mat_data[em_idx]);
  
  /*Integrate the continuous spectrum and the characteristic so as to get the relative intenisities right*/
  Icont=tube_current/CELE*BKRAMER*Z*(E0*log(Emax)-E0*log(Emin) - Emax + Emin);
  /*check if E0 >Ek. If not, no characteristic emission can take place*/
  if (E0>em_p->Ek){
    double Bk=1.2e-5*pow(em_p->Ek,1.67)*exp(-0.077*Z);
    Ichar=tube_current/CELE*4*M_PI*(E0/em_p->Ek-1)*Bk;
    double Ichar_tot=0;
    int linec=0;
    linemin=0;
    linemax=em_p->linec;
    for (ii=0;ii<em_p->linec;ii++){
      /*if the interval [Emin,Emax] contains 5 sigma of the characteristic peak - use full peak.*/
      if (Emin>em_p->e[ii]+5*em_p->w[ii] ){
        /*way below of energy limit - do not use.*/
        linemin=ii;
      }else if ( Emax<em_p->e[ii]-5*em_p->w[ii]){
        /*way above of energy limit - do not use.*/
        linemax= linemax<ii?linemax:ii;
      }else{
        linec++;
      }
      if(Emax>em_p->e[ii]+5*em_p->w[ii] && Emin<em_p->e[ii]-5*em_p->w[ii]){
        Ichar_tot+=em_p->i[ii];
      }else{
        if(!lorentzian){
	  /*We only partially use the peak - so update the relative intensity to reflect that*/
          Ichar_tot+=em_p->i[ii]*0.5*( erf(Emax-em_p->e[ii]/em_p->w[ii]/M_SQRT2) - erf( Emin-em_p->e[ii]/em_p->w[ii]/M_SQRT2) );
          em_p->i[ii]= em_p->i[ii]*0.5*( erf((Emax-em_p->e[ii]/M_SQRT2)/em_p->w[ii]/M_SQRT2) - erf((Emin-em_p->e[ii])/em_p->w[ii]/M_SQRT2) );
        }else{
	  /* 1/pi atan(2x/w)) is integrated lorentzian of fwhm w, MHM, April 2015 */
          Ichar_tot+=em_p->i[ii]*(1.0/M_PI)*( atan(2*(Emax-em_p->e[ii])/em_p->w[ii]) - atan(2*(Emin-em_p->e[ii])/em_p->w[ii]) );
        }
      }
    }

    p_continous=Icont/(Ichar*Ichar_tot+Icont);
  }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());

  mu_electron=0;

%}

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_electron;

  double px,py,pz; 
  Coords P;

  /* transform initial coords to ones in a frame with origin at the center of e-beam with optical axis pointing towards exit window*/
  P=coords_set(x1,y1,z1);
  P=coords_add(rot_apply(R_xray_gen,P),O_xray_gen);
  coords_get(P,&x,&y,&z);

  /*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();
  if(r<frac){
    /*bremsstrahlung*/
    e=rand01()*(Emax-Emin)+Emin;
    k=e*E2K;
    pmul=tube_current/CELE*BKRAMER*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{
    struct xray_em_data *pt=&(xray_mat_data[em_idx]);
    /*characteristic radiation*/
    /*first pick a possible line*/
    r=rand01()*(linemax-linemin) + linemin;
    int lineno=(int)floor(r);
    if (lineno==pt->linec) {
      lineno--;/*we might get overflow*/
    }

    if(!lorentzian){
      pmul=pt->i[lineno]*Ichar;
      e=(randnorm()*pt->w[lineno]+pt->e[lineno]);
      /*this can be very inefficient*/
      while (e<Emin || e>Emax){
        e=(randnorm()*pt->w[lineno]+pt->e[lineno]);
      }
    }else{/* tan((rand-0.5)/pi) is Lorentzian with FWHM of 2, MHM April 2015 */
      /* compute upper and lower random range to map onto energy bounds such that a*tan(u)+b = energy. Note -0.5<u<0.5 */
      double umin=atan(2*(Emin-pt->e[lineno])/pt->w[lineno])/M_PI;
      double umax=atan(2*(Emax-pt->e[lineno])/pt->w[lineno])/M_PI;
      pmul=pt->i[lineno]*Ichar*(umax-umin); /* weight intensity for partial line strength */
      e=tan(((umax-umin)*rand01()+umin)*M_PI)*pt->w[lineno]/2+pt->e[lineno];
    }
    k=E2K*e;
    pmul*=(1-p_continous)/frac;
  }

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

  int ie;
  /*Correct for absorption*/
  double mu_abs,l0,l1,lx,ly,lz,klx,kly,klz,xw,yh;
  l0=0;l1=0;
  /*if dimensions are set the anode has a limited size - otherwise use a
    practically infinite slab*/
  if(!xwidth) xw=FLT_MAX;
  if(!yheight) yh=FLT_MAX;
  coords_get(rot_apply(R_xray_geni,coords_set(x,y,z)),&lx,&ly,&lz);
  coords_get(rot_apply(R_xray_geni,coords_set(kx,ky,kz)),&klx,&kly,&klz);
  /*find path length inside anode material*/
  ie=box_intersect(&l0,&l1,lx,ly+thickness/2.0,lz,klz,kly,klz,xw,thickness,yh);
  if(!ie || l0>0){
    /*photon is somehow outside the anode - this should not happen*/
    ABSORB;
  }

  mu_abs=Table_Value(T, k*K2E, 5)*rho*1e2; /*mu_abs in m^-1*/
  p*=exp(-l1*mu_abs);
  /*set a random phase*/
  phi=rand01()*2*M_PI;

  /*Finally, if exit_window_refpt is set revert to legacy behaviour where the centre of the exit window is the
    reference point of the component*/
  if(exit_window_refpt){
    z=z-dist;
  }
  /*set a scatter pt at the generation pt*/
  SCATTER;
%}

MCDISPLAY
%{
  _class_particle p1,p2;
  
  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;
  p1.x=-width_2;p1.y=0;p1.z=-height_2;
  p2.x=-width_2;p2.y=0;p2.z= height_2;

  mccoordschange(O_xray_gen,R_xray_gen,&p1);
  mccoordschange(O_xray_gen,R_xray_gen,&p2);
  line(p1.x,p1.y,p1.z,p2.x,p2.y,p2.z);

  p1.x=width_2;p1.y=0;p1.z=height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&p1);
  line(p2.x,p2.y,p2.z,p1.x,p1.y,p1.z);
  
  p2.x=width_2;p2.y=0;p2.z=-height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&p2);
  line(p1.x,p1.y,p1.z,p2.x,p2.y,p2.z);

  p1.x=-width_2;p1.y=0;p1.z=-height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&p1);
  line(p2.x,p2.y,p2.z,p1.x,p1.y,p1.z);

  /*this is the mean penetration depth of electron that get converted to x-rays*/  
  p1.x=-width_2;p1.y=-mu_electron;p1.z=-height_2;
  p2.x=-width_2;p2.y=-mu_electron;p2.z= height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&p1);
  mccoordschange(O_xray_gen,R_xray_gen,&p2);
  dashed_line(p1.x,p1.y,p1.z,p2.x,p2.y,p2.z,5);
  p1.x=width_2;p1.y=-mu_electron;p1.z=height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&p1);
  dashed_line(p2.x,p2.y,p2.z,p1.x,p1.y,p1.z,5);
  p2.x= width_2;p2.y=-mu_electron;p2.z=-height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&p2);
  dashed_line(p1.x,p1.y,p1.z,p2.x,p2.y,p2.z,5);
  p1.x=-width_2;p1.y=-mu_electron;p1.z=-height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&p1);
  dashed_line(p2.x,p2.y,p2.z,p1.x,p1.y,p1.z,5);
  
  p1.x=-width_2;p1.y=-mu_electron;p1.z=-height_2;
  p2.x=-width_2;p2.y=0;p2.z=-height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&p1);
  mccoordschange(O_xray_gen,R_xray_gen,&p2);
  line(p2.x,p2.y,p2.z,p1.x,p1.y,p1.z);
  p1.x=width_2;p1.y=-mu_electron;p1.z=-height_2;
  p2.x=width_2;p2.y=0;p2.z=-height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&p1);
  mccoordschange(O_xray_gen,R_xray_gen,&p2);
  line(p2.x,p2.y,p2.z,p1.x,p1.y,p1.z);
  p1.x=-width_2;p1.y=-mu_electron;p1.z=height_2;
  p2.x=-width_2;p2.y=0;p2.z=height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&p1);
  mccoordschange(O_xray_gen,R_xray_gen,&p2);
  line(p2.x,p2.y,p2.z,p1.x,p1.y,p1.z);
  p1.x=width_2;p1.y=-mu_electron;p1.z=height_2;
  p2.x=width_2;p2.y=0;p2.z=height_2;
  mccoordschange(O_xray_gen,R_xray_gen,&p1);
  mccoordschange(O_xray_gen,R_xray_gen,&p2);
  line(p2.x,p2.y,p2.z,p1.x,p1.y,p1.z);


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

END