/*******************************************************************************
*
* 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 B. 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