/*******************************************************************************
*
* McXtrace, x-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: Ring_c
*
* %Identification
*
* Written by: Erik B Knudsen and Desiree D. M. Ferreira
* Date: Feb. 2016
* Version: 1.0
* Release: McXtrace 1.2
* Origin: DTU Physics, DTU Space
* Modified by: Søren Jeppesen
* Date: Feb. 2017
* Version: 1.1
* Release: McXtrace 1.2
* Origin: DTU Physics, DTU Space
*
* Stack of conical shells as part of a Wolter optic.
*
* %Description
* A stack of conical shells is simulated. 
* To intersect the Wolter I plates we take advatage of the azimuthal symmetry and only consider the radial component
* of the photon's wavevector.
* 
* %Parameters
* Input parameters:
* radius_m: [m] Ring radius of the upper (reflecting) plate of the shell at the optic centre. 
* yheight:  [m] Height of the shell.
* chamferwidth:  [m] Width of side walls.
* gap:      [m] Gap between the plate and the intersection plane with the hyperbolic section. (currently ignored)
* Z0:       [m] Distance between optics centre plane and focal spot (essentially focal length).
* mirror_reflec: [ ] Data file containing reflectivities of the reflector surface (TOP).
* bottom_reflec: [ ]  Data file containing reflectivities of the bottom surface (BOTTOM).
* R_d:      [ ] Default reflectivity value to use if no reflectivity file is given. Useful f.i. is one surface is reflecting and the others absorbing.
* primary:  [ ] If non-zero, the shell is considered a primary reflector, and extends towards negative z. I.e. the entry plane is behind the z=0-plane. If zero, the shell is considered secondary
*               and extends from the z=0-plane and towards positive z.
* dalpha:   [deg] Offset to the alpha angle computed from the focal length. Useful for targeting the modified conical geometry (currently ignored).
* waviness: [rad] Waviness of the shell reflecting surface. The slope error is assumed to be uniformly distributed in the interval [-waviness:waviness].
* longw:    [ ] If non-zero, waviness is 1D and along the shell axis.
* %End
*******************************************************************************/

DEFINE COMPONENT Ring_c
SETTING PARAMETERS (int ring_nr, radius_m, Z0, xwidth, yheight, gap=0, chamferwidth=0, length=0, string mirror_reflec="", string bottom_reflec="", R_d=1, primary=1, dalpha=0, waviness=0, longw=0)

SHARE
%{
  #ifndef MCSPO_INTERSECT_CONE
  #define MCSPO_INTERSECT_CONE 1

  int
  intersect_cone (double* l0, double x, double y, double z, double kx, double ky, double kz, double alpha, double radius, double* nx, double* ny, double* nz) {
    double kxn = kx, kyn = ky, kzn = kz;
    NORM (kxn, kyn, kzn);
    double c = tan (alpha);
    double z0 = radius / c;
    double c2 = c * c;
    double A, B, C;
    A = kxn * kxn + kyn * kyn - c2 * kzn * kzn;
    B = 2 * (kxn * x + kyn * y - c2 * kzn * (z - z0));
    C = x * x + y * y - c2 * (z - z0) * (z - z0);

    int status;
    double l1;
    if ((status = solve_2nd_order (l0, &l1, A, B, C)) == 0) {
      /*note that if l1->NULL only the smallest positive solution is returned*/
      fprintf (stderr, "Error(%s): No solution to second order eq.\n", "Pore_c");
    }
    /*compute normal vector here*/
    x += kxn * (*l0);
    y += kyn * (*l0);
    z += kzn * (*l0);

    double vn = sqrt (x * x + y * y);
    *nx = x / vn;
    *ny = y / vn;

    *nz = 1;

    *nx *= cos (alpha);
    *ny *= cos (alpha);
    *nz *= sin (alpha);

    return status;
  }
  #endif

  #ifndef PROP_Z
   #define PROP_Z(zz)\
      mcPROP_Z(zz)
  #endif

  #ifndef mcPROP_Z
   #define mcPROP_Z(zz) \
    do { \
      MCNUM mc_dl,mc_k; \
      if(kz == 0) { ABSORB; }; \
      mc_k=sqrt(scalar_prod(kx,ky,kz,kx,ky,kz)); \
      mc_dl= ((zz)-z) * mc_k / kz; \
      if(mc_dl<0 && allow_backprop==0) { ABSORB; };\
      PROP_DL(mc_dl); \
    } while(0)
  #endif
%}

DECLARE
%{
  double nExit[3];
  double wExit[3];
  double nEntry[3];
  double wEntry[3];
  double nTop[3];
  double nBottom[3];
  double radius_1;
  double radius_2;
  double e_min[2];
  double e_step[2];
  double e_max[2];
  double theta_min[2];
  double theta_step[2];
  double theta_max[2];

  double zentry;
  double zexit;
  double* zentry_vec;
  double* zexit_vec;

  t_Table reflec_top_table;
  t_Table reflec_bottom_table;
  t_Table size_table;
%}

INITIALIZE
%{
  /*read data from files into tables using read_table-lib*/
  char* filenames[2] = { mirror_reflec, bottom_reflec };
  t_Table* ref_tables[2] = { &reflec_top_table, &reflec_bottom_table };
  int i;

  /*read data from files into tables using read_table-lib*/
  for (i = 0; i < 2; i++) {
    char* reflec = filenames[i];
    t_Table* tp = ref_tables[i];
    if (reflec && strlen (reflec)) {
      char** header_parsed;
      /* read 1st block data from file into tp */
      if (Table_Read (tp, reflec, 1) <= 0) {
        exit (fprintf (stderr, "Error(%s): can not read file %s\n", NAME_CURRENT_COMP, reflec));
      }
      header_parsed = Table_ParseHeader (tp->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]) {
        e_min[i] = strtod (header_parsed[0], NULL);
        e_max[i] = strtod (header_parsed[1], NULL);
        e_step[i] = strtod (header_parsed[2], NULL);
        theta_min[i] = strtod (header_parsed[3], NULL);
        theta_max[i] = strtod (header_parsed[4], NULL);
        theta_step[i] = strtod (header_parsed[5], NULL);
      } else {
        exit (fprintf (stderr, "Error (%s): wrong/missing header line(s) in file %s\n", NAME_CURRENT_COMP, reflec));
      }
      if (!((int)(e_max[i] - e_min[i]) == (int)((tp->rows - 1) * e_step[i]))) {
        exit (fprintf (stderr, "Error (%s): e_step does not match e_min and e_max in file %s\n", NAME_CURRENT_COMP, reflec));
      }
      if (!((int)(theta_max[i] - theta_min[i]) == (int)((tp->columns - 1) * theta_step[i]))) {
        exit (fprintf (stderr, "Error (%s): theta_step does not match theta_min and theta_max in file %s\n", NAME_CURRENT_COMP, reflec));
      }
    } else {
      /*mark the table as unread by setting "rows" to -1
        This will trigger the default reflectivity.*/
      tp->rows = -1;
    }
  }

  /* Read table with the individual shell size parameters */

  if (size_file) {
    if (Table_Read (&(size_table), size_file, ring_nr) <= 0)
      exit (fprintf (stderr, "Error (%s): Could not read %s. Aborting.\n", NAME_CURRENT_COMP, size_file));
  }

  zentry_vec = calloc (size_table.rows, sizeof (double));
  zexit_vec = calloc (size_table.rows, sizeof (double));

  double alpha, thetap, thetah, P, d, e, C0, Z;
  /*There are in general 68 reflecting planes*/
  for (i = 0; i < size_table.rows; i++) {
    if (primary) {
      Z0p = radius_m / tan (alpha);
      D = sqrt (radius_m * radius_m + Z0p * Z0p);
      zentry = Z0p * (1 - (D + length) / D);
      zexit = 0;
      radius_1 = (D + length) / D * radius_m;
      radius_2 = radius_m;
    } else {
      alpha *= 3.0;
      Z0p = radius_m / tan (alpha);
      D = sqrt (radius_m * radius_m + Z0p * Z0p);
      zentry = 0;
      zexit = Z0p * (1 - (D - length) / D);
      radius_1 = radius_m;
      radius_2 = (D - length) / D * radius_m;
    }
    zentry_vec[i] = zentry;
    zexit_vec[i] = zexit;
  }

  /*find minimum zentry - i.e. the "first" entry plane.*/
  zentry = zentry_vec[0];
  zexit = zexit_vec[0];
  for (i = 0; i < size_table.rows - 1; i++) {
    if (zentry_vec[i] < zentry_vec[i + 1]) {
      zentry = zentry_vec[i];
    }
    if (zexit_vec[i] < zexit_vec[i + 1]) {
      zexit = zexit_vec[i];
    }
  }

  nEntry[0] = 0;
  nEntry[1] = 0;
  nEntry[2] = -1;
  wEntry[0] = wEntry[1] = 0;
  wEntry[2] = zentry;

  nExit[0] = 0;
  nExit[1] = 0;
  nExit[2] = 1;
  wExit[0] = wExit[1] = wExit[2] = zexit;
%}

TRACE
%{
  double r_p_min, r_p_max, r2;
  int hit = 0;
  int hit_chamfer = 0;

  /*assuming the table to be sorted*/
  ALLOW_BACKPROP;
  PROP_Z (zentry_vec[0]);
  r_p_min = Table_Index (size_table, 0, 4);
  r_p_max = Table_Index (size_table, size_table.rows - 1, 4);
  r2 = (x * x + y * y);

  /*TODO this should also check for plate width*/
  hit = ((r2 > (r_p_min - yheight) * (r_p_min - yheight)) && (r2 < (r_p_max + pore_th) * (r_p_max + pore_th)));
  if (hit) {
    /*we have a chance of hitting a ring - we might still miss due to beam divergence though*/
    hit = 0;
    int ii = 0;
    int hit_chamfer;

    while (!hit && ii < size_table.rows) {
      ALLOW_BACKPROP;
      PROP_Z (zentry_vec[ii]);
      radius_p = Table_Index (size_table, ii, 4);
      radius_p_2 = Table_Index (size_table, ii - 1, 4);

      radius_m = Table_Index (size_table, ii, 3);
      radius_m_2 = Table_Index (size_table, ii - 1, 3);

      hit = ((x * x + y * y < radius_p * radius_p) && (x * x + y * y > (radius_p - yheight) * (radius_p - yheight)));
      ii++;
    }
    /*ii-1 is now the numbeer of the hit plate - unless it is outside all of them,
      radius_p and radius_m are the relevant radii for this plate.*/

    /*TODO figure out which pore we are at and set the side walls accordingly*/
    hit_chamfer = 0;

    enum { LEFT, RIGHT, TOP, BOTTOM, EXIT, NONE } wall;
    t_Table* reflec_table = NULL;
    double R;

    if (hit) {
      SCATTER;
      int exit = 0;
      int intersections[5];
      int i_small;
      double l[5];
      double l_small;

      double nx, ny, nz;

      while (!exit) {
        l_small = DBL_MAX;
        wall = NONE;
        double nx, ny, nz;
        double wx, wy, wz;
        int prm_idx;
        /*exit plane*/
        intersections[EXIT] = plane_intersect (l + EXIT, x, y, z, kx, ky, kz, nExit[0], nExit[1], nExit[2], wExit[0], wExit[1], wExit[2]);
        if (intersections[EXIT] && l[EXIT] > DBL_EPSILON && l[EXIT] < l_small) {
          l_small = l[EXIT];
          i_small = intersections[EXIT];
          wall = EXIT;
        }
        /*top surface - the real reflecting surface*/
        intersections[TOP] = intersect_cone ((l + TOP), x, y + radius_m, z, kx, ky, kz, alpha, radius_m, &(nTop[0]), &(nTop[1]), &(nTop[2]));
        if (intersections[TOP] && l[TOP] > DBL_EPSILON && l[TOP] < l_small) {
          l_small = l[TOP];
          i_small = intersections[TOP];
          wall = TOP;
        }
        /*bottom surface*/
        intersections[BOTTOM]
            = intersect_cone ((l + BOTTOM), x, y + radius_m, z, kx, ky, kz, alpha, radius_m - yheight, &(nBottom[0]), &(nBottom[1]), &(nBottom[2]));
        if (intersections[BOTTOM] && l[BOTTOM] > DBL_EPSILON && l[BOTTOM] < l_small) {
          l_small = l[BOTTOM];
          i_small = intersections[BOTTOM];
          wall = BOTTOM;
        }

        /*sort intersections to find the smallest positive one*/
        switch (wall) {
        case TOP:
          /*handle top wall reflection*/
          reflec_table = &reflec_top_table;
          nx = nTop[0];
          ny = nTop[1];
          nz = nTop[2];
          prm_idx = 0;
          break;
        case BOTTOM:
          /*handle bottom wall "reflection"*/
          reflec_table = &reflec_bottom_table;
          nx = nBottom[0];
          ny = nBottom[1];
          nz = nBottom[2];
          prm_idx = 1;
          break;
        case EXIT:
          /*photon will exit pore*/
          exit = 1;
          break;
        }
        if (exit) {
          continue;
        }

        PROP_DL (l_small);

        double kix = kx, kiy = ky, kiz = kz;
        double k = sqrt (kx * kx + ky * ky + kz * kz);
        double e = K2E * k;
        double s = scalar_prod (kx, ky, kz, nx, ny, nz);
        double theta = RAD2DEG * (M_PI_2 - acos (s / k)); /*pi_2 since theta is supposed to be the grazing angle*/

        /*if we have waviness alter the normal vector slightly*/
        if (waviness != 0) {
          /*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;
            }
            double tx, ty, tz;
            vec_prod (tx, ty, tz, 0, 0, 1, nx, ny, nz);
            rotate (nx, ny, nz, nx, ny, nz, dtheta, tx, ty, tz);
          } else {
            /*waviness is also transversal but isotropic*/
            double radius;
            if (theta < waviness) {
              radius = atan (waviness);
              randvec_target_circle (&nx, &ny, &nz, NULL, nx, ny, nx, radius);
            } else {
              radius = (atan (theta) + atan (waviness)) / 2.0;
              randvec_target_circle (&nx, &ny, &nz, NULL, nx, ny, nx + radius - atan (theta), radius);
            }
            NORM (nx, ny, nz);
          } /*reflect the photon through the surface normal*/
          /*recompute theta*/
          theta = RAD2DEG * 0.5 * acos (scalar_prod (kx, ky, kz, kix, kiy, kiz) / k / k);
        }
        /*reflect the photon through the surface normal*/
        if (s != 0) {
          kx -= 2 * s * nx;
          ky -= 2 * s * ny;
          kz -= 2 * s * nz;
        }
        SCATTER;

        if (reflec_table == NULL || reflec_table->rows == -1) {
          R = R_d;
        } else {
          R = Table_Value2d (*reflec_table, (e - e_min[prm_idx]) / e_step[prm_idx], (theta - theta_min[prm_idx]) / theta_step[prm_idx]);
        }
        p *= R;
      }
    } else if (hit_chamfer) {
      ABSORB;
    } else {
      /*no hit*/
      ABSORB;
    }
  }
%}

FINALLY
%{
  free (zentry_vec);
%}

MCDISPLAY
%{
  circle ("xy", 0, 0, zentry, radius_p);
  circle ("xy", 0, 0, zentry, radius_p - yheight);

  circle ("xy", 0, 0, 0, radius_m);
  circle ("xy", 0, 0, 0, radius_m - yheight);
%}

END