/*******************************************************************************
*
* McXtrace, X-ray tracing package
*         Copyright, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*         University of Copenhagen, Copenhagen, Denmark
*
* Component: Multilayer_elliptic
*
* %I
*
* Written by: Jana Baltser, Peter Willendrup, Anette Vickery, Andrea Prodi, Erik Knudsen
* Date: February 2011
* Version: 1.0
* Origin: NBI
*
* Elliptic multilayer mirror (in XZ)
* 
* %Description
* Reads reflectivity values from a data input file (Ref.dat) for a Si/W multilayer.
* The multilayer code reflects ray in an ideal geometry, does not include surface imperfections
*
* The mirror is positioned such that the long axis of the mirror elliptical surface coincides with
* z-axis
* 
* The algorithm:
*  Incoming photon's coordinates and direction (k-vector) are transformed into an elliptical reference frame 
* (elliptical parameters are calculated according to the mirror's position and its focusing distances and the  
* incident angle), the intersection point is then defined. A new, reflected photon is then starting at the  
* point of intersection.
*
* Example: Multilayer_elliptic(
*   coating = "Ref_W_B4C.txt", theta = 1.2,
*   s1 = 1, s2 = 2, length = 0.1, width = 0.1, R0 = 1,
*   Emin=7, Emax=10, Estep=0.05)
*
* %Parameters
* Input parameters:
* theta: [deg] Design angle of incidence.
* s1:    [m]   Design distance from the source to the multilayer.
* s2:    [m]   Design focusing distance of the multilayer.
* zdepth:[m]   Length of the mirror along Z.
* xwidth:[m]   Width of the mirror along X-axis.
* Gamma: [ ]   High electron density fraction of bilayer (in kinematical appr.).
* Lambda:[m]   Thickness of bilayer (in kinematical appr.).
* rho_AB:[ ]   Number electron density constrast in bilayer (in kinematical appr.).
* N:     [1]   Number of bilayers (in kinematical appr.).
* coating: [str] Datafile containing reflectivity values as a function of q and E.
* Emin:  [keV] Lower limit of energy interval in datafile. Overrides what's written in the datafile header.
* Emax:  [keV] Upper limit of energy interval in datafile. Overrides what's written in the datafile header.
* Estep: [keV] Step between energy sample points in datafile. Overrides what's written in the datafile header.
* R0:    [1]   Maximal reflectivity
* length:  [m]   alternate name for zdepth (obsolete)
* width:   [m]   alternate name for xwidth (obsolete)
* %End
*******************************************************************************/

DEFINE COMPONENT Multilayer_elliptic

SETTING PARAMETERS (string coating="Ref_W_B4C.txt",
    theta=1.2,s1=0,s2=0,length=0.5,width=0.2,R0=1, 
    Emin=-1, Emax=-1, Estep=-1, Gamma=0, Lambda=0, rho_AB=0, int N=0,
    xwidth=0, zdepth=0)

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

SHARE 
%{
  #include <complex.h>
  %include "read_table-lib"
  %include "reflectivity-lib"
  /*something that would be relevant for ALL elliptical mirrors*/
  /* coordinate transformation McXtrace-Ellipse (ME) and Ellipse-McXtrace(EM) functions */
  #pragma acc routine
  void
  CoordTransME (double* x_el, double* y_el, double* z_el, double x0, double y0, double z0, double Zmir, double Ymir, double xi_mir) {
    *x_el = x0;
    *y_el = cos (xi_mir) * y0 + sin (xi_mir) * z0 + Ymir;
    *z_el = -sin (xi_mir) * y0 + cos (xi_mir) * z0 + Zmir;
  }

  #pragma acc routine
  void
  CoordTransEM (double* x_gen, double* y_gen, double* z_gen, double x0, double y0, double z0, double Zmir, double Ymir, double xi_mir) {
    *x_gen = x0;
    *y_gen = cos (xi_mir) * (y0 - Ymir) - sin (xi_mir) * (z0 - Zmir);
    *z_gen = sin (xi_mir) * (y0 - Ymir) + cos (xi_mir) * (z0 - Zmir);
  }
%}

DECLARE
%{
  double a;
  double b;
  double c;
  double M;
  double Z0;
  double Y0;
  double xi;
  double cost0;
  int kinematical;
  t_Reflec re;
%}

INITIALIZE
%{
  /* calculation of the elliptical parameters according to the input mirror parameters:
  ellipse major axis a/2, minor axis b/2, M-magnification factor, Z0&Y0 - position of the mirror centre in the elliptical coordinate system.*/
  double Theta = DEG2RAD * theta;

  if (xwidth)
    width = xwidth;
  if (zdepth)
    length = zdepth;

  M = s2 / s1;
  cost0 = (1 - M) / sqrt (1 - 2 * M + M * M + 4 * M * (cos (Theta) * cos (Theta)));
  a = (s1 * sqrt (1 - cost0 * cost0 + cos (Theta) * cos (Theta) * cost0 * cost0))
      / (cost0 * cos (Theta) + sqrt (1 - cost0 * cost0 + (cos (Theta) * cos (Theta)) * cost0 * cost0));
  c = a * cos (Theta) / sqrt (1 - cost0 * cost0 + (cos (Theta) * cos (Theta)) * cost0 * cost0);
  b = sqrt (a * a - c * c);
  Z0 = a * cost0;
  Y0 = b * sin (acos (cost0));
  xi = -atan ((Z0 * b * b) / (Y0 * a * a));

  int status = 0;
  if (!Gamma && !Lambda) {
    /*refrain from using kinematical approximation - instead use reflectivity datafile*/
    kinematical = 0;
    /* reflectivity datafile parsing COATING_UNDEFINED - means set the type according to what is found in the file*/
    status = reflec_Init (&re, COATING_UNDEFINED, coating, NULL);
  } else {
    kinematical = 1;
    re.type = KINEMATIC;
    status = reflec_Init_kinematic (&(re), N, Gamma, Lambda,
                                    rho_AB); /*number of layers, ratio of high e-density material in layer, thickness of layer, e-density contrast*/
  }
%}

TRACE
%{
  double K, vink;
  double x_el, y_el, z_el;    // beginning coordinates transformed into the ellipse system
  double kx_el, ky_el, kz_el; // kvector transformed into the ellipse system, hence

  double A, B, C, D, t0, t1;
  double x_int, y_int, z_int, dist; // intersection with the elliptical surface
  double nx, ny, nz;
  double kxn, kyn, kzn; // reflected ray's kvector

  /* get the photon's coordinates and kvector in the ellipse frame */
  K = sqrt (kx * kx + ky * ky + kz * kz);

  CoordTransME (&x_el, &y_el, &z_el, x, y, z, Z0, Y0, xi);
  CoordTransME (&kx_el, &ky_el, &kz_el, kx, ky, kz, 0, 0, xi);

  NORM (kx_el, ky_el, kz_el);

  /*intersection calculation*/
  A = b * b * kz_el * kz_el + a * a * ky_el * ky_el;
  B = 2.0 * (z_el * kz_el * b * b + y_el * ky_el * a * a);
  C = b * b * z_el * z_el + a * a * y_el * y_el - a * a * b * b;
  D = B * B - 4 * A * C;
  if (D >= 0) {
    t0 = (-B - sqrt (D)) / (2 * A);
    t1 = (-B + sqrt (D)) / (2 * A);
    if (t0 < 0 && t1 >= 0) {
      double ttmp = t0;
      t0 = t1;
      t1 = ttmp;
    }
    /* check whether our intersection lies within the boundaries of the mirror*/
    x_int = x_el + kx_el * t0;
    y_int = y_el + ky_el * t0;
    z_int = z_el + kz_el * t0;

    if (y_int >= 0 && fabs (x_int) <= width / 2) {
      dist = sqrt ((x_el - x_int) * (x_el - x_int) + (y_el - y_int) * (y_el - y_int) + (z_el - z_int) * (z_el - z_int));
      PROP_DL (dist);

      if (fabs (z) <= length / 2) { /*finally in business on the mirror! YAY! */
        nx = 0;
        if (fabs (z_int) == 0) {
          ny = 1;
          nz = 0;
        } else {
          ny = (a * a * y_int) / (b * b * z_int);
          nz = 1.0;
        }
        NORM (nx, ny, nz);
        vink = scalar_prod (nx, ny, nz, kx_el, ky_el, kz_el);
        kxn = kx_el - 2.0 * vink * nx;
        kyn = ky_el - 2.0 * vink * ny;
        kzn = kz_el - 2.0 * vink * nz;
        NORM (kxn, kyn, kzn);

        double kxo, kyo, kzo;
        kxo = kx;
        kyo = ky, kzo = kz;
        CoordTransEM (&kx, &ky, &kz, kxn, kyn, kzn, 0, 0, xi);

        kx = K * kx;
        ky = K * ky;
        kz = K * kz;

        double QQ, EE, Ref;
        QQ = sqrt ((kx - kxo) * (kx - kxo) + (ky - kyo) * (ky - kyo) + (kz - kzo) * (kz - kzo));
        EE = K * K2E;

        if (kinematical) {
          /*
           * \Lambda: thickness of bilayer - following notation in Als-Nielsen/McMorrow
           * \Gamma:  \Gamma*\Lambda thickness of high electron density material.
           * r1(zeta) = 2 i r_0 \rho_{AB} \left(\frac{\Lambda^2 \Gamma}{\zeta}\right) \frac{\sin\left(\pi\Gamma\zeta\right)}{\pi\Gamma\zeta);
           */
          Ref = reflecq (re, QQ, 0, 0, 0);
          if (Ref > 1) {
            /*Reflectivity can't be >1*/
            Ref = 1.0;
          }
        } else {
          /*interpolate in table*/
          Ref = reflecq (re, QQ, 0, 0, 0);
        }

        /* apply reflectivity */
        p *= Ref;
        SCATTER;

      } else {
        RESTORE_XRAY (INDEX_CURRENT_COMP, x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p);
      }
    }
  }
%}

MCDISPLAY
%{
  /*
  rectangle("xz",0,0,0,width,length); */
  int i, j, NN = 10;
  double x0, y0, z0;
  double x1, y1, z1, z_el, y_el;

  x0 = -width / 2.0;

  for (i = 0; i <= NN; i++) {
    z0 = -length / 2.0;
    z_el = cos (xi) * z0 + Z0; // transformation to EL reference frame
    y_el = b * sqrt (1.0 - ((z_el * z_el) / (a * a)));
    y0 = cos (xi) * (y_el - Y0) - sin (xi) * (z_el - Z0);
    line (-width / 2, y0, z0, width / 2, y0, z0);
    for (j = 0; j <= NN; j++) {
      z1 = z0 + length / NN;
      z_el = cos (xi) * z1 + Z0;
      y_el = b * sqrt (1.0 - ((z_el * z_el) / (a * a)));
      y1 = cos (xi) * (y_el - Y0) - sin (xi) * (z_el - Z0);
      line (x0, y0, z0, x0, y1, z1);
      y0 = y1;
      z0 = z1;
      line (-width / 2, y1, z1, width / 2, y1, z1);
    }
    x0 = x0 + width / NN;
  }
%}

END