/*******************************************************************************
*
* McXtrace, X-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: Absorption_sample
*
* %Identification
*
* Written by:Erik B Knudsen
* Date: March 2011
* Version: 1.0
* Origin: Risoe
*
* Sample component with absorbing materials.
*
* %Description
* A sample component consisting of a volume of one material and volume of another 
* material inside. This is useful as a phantom for simulating tomography experiments.
* The inner material can be left unset (all 0), to only have one volume.
*
* Sample shape may be a cylinder, a sphere, a box or any other shape
*   box/plate:       xwidth x yheight x zdepth
*   cylinder:        radius x yheight
*   sphere:          radius (yheight=0)
*   any shape:       geometry=OFF/PLY file
*
* Example: Absorption_sample( material_datafile_o="Mn.txt", xwidth_o = 0.5, yheight_o = 0.5, zdepth_o = 0.0001, rho_o=7.15 )
*
* %P
* Input parameters: 
* radius_o:           [m] Radius of "outer" enclosing material cylinder (0)
* xwidth_o:           [m] Width of "outer" enclosing material box (yheight)
* yheight_o:          [m] Height of "outer enclosing material box (1) 
* zdepth_o:           [m] Thickness of outer enclosing material box (yheight)
* radius_i:           [m] Radius of "inner" enclosed material cylinder (0)
* xwidth_i:           [m] Width of "inner" enclosed material box (yheight)
* yheight_i:          [m] Height of "inner enclosed material (1) 
* zdepth_i:           [m] Thickness of inner enclosed material box (yheight)
* x_i:                [m] Center x-coordinate of "inner" object
* y_i:                [m] Center y-coordinate of "inner" object
* z_i:                [m] Center z-coordinate of "inner" object
* rho_i:              [g/cm^3] density of the enclosed material
* rho_o:              [g/cm^3] density of the enclosing material
* material_datafile_o:[str] Name of file containing material data for outer volume.
* material_datafile_i:[str] Name of file containing material data for inclusion.
* geometry_o:         [str] Name of the outer Object File Format (OFF) or PLY file for complex geometry. The OFF/PLY file may be generated from XYZ coordinates using qhull/powercrust
* geometry_i:         [str] Name of an inner Object File Format (OFF) or PLY file for complex geometry. The OFF/PLY file may be generated from XYZ coordinates using qhull/powercrust [str]
*
* %Link
* Meshlab https://www.meshlab.net/
* %Link
* Geomview and Object File Format (OFF) <http://www.geomview.org>
* %Link
* Java version of Geomview (display only) jroff.jar <http://www.holmes3d.net/graphics/roffview/> 
* %Link
* qhull <http://qhull.org>
* %Link
* Powercrust https://www.cs.ucdavis.edu/~amenta/powercrust.html
* %End
*******************************************************************************/

DEFINE COMPONENT Absorption_sample

SETTING PARAMETERS (string material_datafile_i="",string material_datafile_o="",
    radius_o=0,xwidth_o=0,yheight_o=1,zdepth_o=0,
    radius_i=0,xwidth_i=0,yheight_i=0.0,zdepth_i=0,
    x_i=0,y_i=0,z_i=0,rho_i=0,rho_o=0,
    string geometry_i="", string geometry_o="")
/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */ 

SHARE
%{
  %include "read_table-lib"
  %include "interoff-lib"

  #ifndef SHAPES_T
  #define SHAPES_T
  enum shapes_t { NONE = -1, SPHERE, CYLINDER, CUBE, ELLIPSOID, ANY };
  #endif
%}

DECLARE
%{
  double Z_A_o;
  double V_o;
  int E_c_o;
  int l_c_o;
  int mu_c_o;
  double Z_A_i;
  double V_i;
  int E_c_i;
  int l_c_i;
  int mu_c_i;
  t_Table t_o;
  t_Table t_i;
  off_struct offdata_o;
  off_struct offdata_i;
  int shape_i;
  int shape_o;
%}

INITIALIZE
%{
  int status;
  V_i = 1;

  /* Checking volumes */
  if (geometry_o && strlen (geometry_o) && strcmp (geometry_o, "NULL") && strcmp (geometry_o, "0")) {
    if (!off_init (geometry_o, xwidth_o, yheight_o, zdepth_o, 0, &offdata_o))
      exit (printf ("Absorption_sample: %s: FATAL: invalid OFF/PLY geometry specification for file '%s'.\n", NAME_CURRENT_COMP, geometry_o));
    shape_o = ANY;
  } else if (!yheight_o && !xwidth_o && !zdepth_o && radius_o)
    shape_o = SPHERE;
  else if (yheight_o && radius_o)
    shape_o = CYLINDER;
  else if (yheight_o && !radius_o) {
    if (!zdepth_o)
      zdepth_o = yheight_o;
    if (!xwidth_o)
      xwidth_o = yheight_o;
    shape_o = CUBE;
  } else {
    exit (fprintf (stderr, "ERROR: %s: Unmeaningful outer volume description.\n", NAME_CURRENT_COMP));
  }

  if (geometry_i && strlen (geometry_i) && strcmp (geometry_i, "NULL") && strcmp (geometry_i, "0")) {
    if (!off_init (geometry_i, xwidth_i, yheight_i, zdepth_i, 0, &offdata_i))
      exit (printf ("Absorption_sample: %s: FATAL: invalid OFF/PLY geometry specification for file '%s'.\n", NAME_CURRENT_COMP, geometry_i));
    shape_i = ANY;
  } else if (!yheight_i && !xwidth_i && !zdepth_i && radius_i)
    shape_i = SPHERE;
  else if (yheight_i && radius_i)
    shape_i = CYLINDER;
  else if (yheight_i && !radius_i) {
    if (!zdepth_i)
      zdepth_i = yheight_i;
    if (!xwidth_i)
      xwidth_i = yheight_i;
    shape_i = CUBE;
  } else {
    fprintf (stderr, "Warning: %s: Invalid inner volume specification. Ignoring.\n", NAME_CURRENT_COMP);
    V_i = 0;
    shape_i = NONE;
  }

  /* Loading datafiles */
  if (material_datafile_o && (status = Table_Read (&t_o, material_datafile_o, 0)) == -1) {
    fprintf (stderr, "Error: (%s) Could not parse file \"%s\".\n", NAME_CURRENT_COMP, material_datafile_o);
    exit (-1);
  }
  if (t_o.columns == 3) { /*which column is the energy in and which holds mu*/
    E_c_o = 0;
    mu_c_o = 1;
  } else {
    E_c_o = 0;
    mu_c_o = 5;
  }
  if (rho_o == 0) {
    char** header_parsed;
    header_parsed = Table_ParseHeader (t_o.header, "Z", "A[r]", "rho", "Z/A", NULL);
    if (header_parsed[3]) { /*assuming that a Z/A is given, i.e. Z and A[r] are redundant*/
      Z_A_o = strtod (header_parsed[3], NULL);
    } else if ((strlen (header_parsed[0])) && (strlen (header_parsed[1]))) {
      Z_A_o = strtod (header_parsed[0], NULL) / strtod (header_parsed[1], NULL);
    }
    if (strlen (header_parsed[2])) {
      rho_o = strtod (header_parsed[2], NULL);
      printf ("INFO: %s: Setting rho_o=%g from %s\n", NAME_CURRENT_COMP, rho_o, material_datafile_o);
    }
  }
  if (V_i) { /*if volume is zero - don't bother to read the file*/
    if (material_datafile_i && (status = Table_Read (&t_i, material_datafile_i, 0)) == -1) {
      fprintf (stderr, "Error: Could not parse file \"%s\" in COMP %s\n", material_datafile_i, NAME_CURRENT_COMP);
      exit (-1);
    }
    if (t_i.columns == 3) {
      E_c_i = 0;
      mu_c_i = 1;
    } else {
      E_c_i = 0;
      mu_c_i = 5;
    }
    if (rho_i == 0) { /*when density is not from input, try to read from tables */
      char** header_parsed;
      header_parsed = Table_ParseHeader (t_i.header, "Z", "A[r]", "rho", "Z/A", NULL);
      if (header_parsed[3]) { /*assuming that a Z/A is given, i.e. Z and A[r] are redundant*/
        Z_A_i = strtod (header_parsed[3], NULL);
      } else if ((strlen (header_parsed[0])) && (strlen (header_parsed[1]))) {
        Z_A_i = strtod (header_parsed[0], NULL) / strtod (header_parsed[1], NULL);
      }
      if (strlen (header_parsed[2])) {
        rho_i = strtod (header_parsed[2], NULL);
        printf ("INFO: %s: Setting rho_i=%g from %s\n", NAME_CURRENT_COMP, rho_i, material_datafile_i);
      }
    }
  }
%}

TRACE
%{
  double l0o = 0, l1o = 0, l0i = 0, l1i = 0, mu_o = 0, mu_i = 0;
  int status, status_i;

  if (shape_o == ANY) {
    status = off_x_intersect (&l0o, &l1o, NULL, NULL, x, y, z, kx, ky, kz, offdata_o);
  } else if (shape_o == CYLINDER) {
    status = cylinder_intersect (&l0o, &l1o, x, y, z, kx, ky, kz, radius_o, yheight_o);
  } else if (shape_o == CUBE) {
    status = box_intersect (&l0o, &l1o, x, y, z, kx, ky, kz, xwidth_o, yheight_o, zdepth_o);
  } else if (shape_o == SPHERE) {
    status = sphere_intersect (&l0o, &l1o, x, y, z, kx, ky, kz, radius_o);
  }

  if (status) { /*rays intersects the enclosing material*/
    PROP_DL (l0o);
    SCATTER;
    status_i = 0;
    double k = sqrt (scalar_prod (kx, ky, kz, kx, ky, kz));

    if (shape_i == ANY) {
      status_i = off_x_intersect (&l0i, &l1i, NULL, NULL, (x - x_i), (y - y_i), (z - z_i), kx, ky, kz, offdata_i);
    } else if (shape_i == CYLINDER) {
      status_i = cylinder_intersect (&l0i, &l1i, (x - x_i), (y - y_i), (z - z_i), kx, ky, kz, radius_i, yheight_i);
    } else if (shape_i == CUBE) {
      status_i = box_intersect (&l0i, &l1i, (x - x_i), (y - y_i), (z - z_i), kx, ky, kz, xwidth_i, yheight_i, zdepth_i);
    } else if (shape_i == SPHERE) {
      status_i = sphere_intersect (&l0i, &l1i, (x - x_i), (y - y_i), (z - z_i), kx, ky, kz, radius_i);
    }

    if (status_i) { /*rays intersect the inclusion*/
      PROP_DL (l0i);
      SCATTER;
      PROP_DL (l1i - l0i);
      SCATTER;
      /*now calculate the mu*/
      mu_i = Table_Value (t_i, k * K2E, mu_c_i) * rho_i;
      p *= exp (-(l1i - l0i) * mu_i * 1e2); /* 1e2 to have in m unit */
    }
    PROP_DL (l1o - l0o - l1i);
    SCATTER;
    mu_o = Table_Value (t_o, k * K2E, mu_c_o) * rho_o;
    p *= exp (-(l1o - l0o - (l1i - l0i)) * mu_o * 1e2);
  }
%}

MCDISPLAY
%{

  if (shape_o == CUBE) {
    box (0, 0, 0, xwidth_o, yheight_o, zdepth_o, 0, 0, 1, 0);
  } else if (shape_o == CYLINDER) {
    circle ("xy", 0, yheight_o / 2.0, 0, radius_o);
    circle ("xy", 0, -yheight_o / 2.0, 0, radius_o);
    line (radius_o, yheight_o / 2.0, 0, radius_o, -yheight_o / 2.0, 0);
    line (-radius_o, yheight_o / 2.0, 0, -radius_o, -yheight_o / 2.0, 0);
    line (0, yheight_o / 2.0, radius_o, 0, -yheight_o / 2.0, radius_o);
    line (0, yheight_o / 2.0, -radius_o, 0, -yheight_o / 2.0, -radius_o);
  } else if (shape_o == SPHERE) {
    circle ("xy", 0, 0, 0, radius_o);
    circle ("xz", 0, 0, 0, radius_o);
    circle ("yz", 0, 0, 0, radius_o);
  } else if (shape_o == ANY) { /* OFF file */
    off_display (offdata_o);
  }

  if (shape_i == CUBE) {
    box (0, 0, 0, xwidth_i, yheight_i, zdepth_i, 0, 0, 1, 0);
  } else if (shape_i == CYLINDER) {
    circle ("xy", 0, yheight_i / 2.0, 0, radius_i);
    circle ("xy", 0, -yheight_i / 2.0, 0, radius_i);
    line (radius_i, yheight_i / 2.0, 0, radius_i, -yheight_i / 2.0, 0);
    line (-radius_i, yheight_i / 2.0, 0, -radius_i, -yheight_i / 2.0, 0);
    line (0, yheight_i / 2.0, radius_i, 0, -yheight_i / 2.0, radius_i);
    line (0, yheight_i / 2.0, -radius_i, 0, -yheight_i / 2.0, -radius_i);
  } else if (shape_i == SPHERE) {
    circle ("xy", 0, 0, 0, radius_i);
    circle ("xz", 0, 0, 0, radius_i);
    circle ("yz", 0, 0, 0, radius_i);
  } else if (shape_i == ANY) { /* OFF file */
    off_display (offdata_i);
  }
%}

END