/*******************************************************************************
*
* McXtrace, x-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: Filter
*
* %Identification
* Written by: Erik Knudsen
* Date: Jan 24, 2011
* Origin: DTU Physics
* Release: McXtrace 1.1
*
* Block of an attenuating material
* 
* %Description
* A chunk of attenuating material. Attenuation is computed through
* the effective length travelled within the material. 
* No scattering is modelled at present.
*
* Filter shape may be a cylinder, a sphere, a box or any other shape.
*   box/plate:       xwidth x yheight x zdepth 
*   cylinder:        radius x yheight (along Y axis)
*   sphere:          radius
*   any shape:       geometry=OFF/PLY_file
*
* Example: Filter(material_datafile="Ge.txt",
*   geometry="wire.ply",xwidth=0.02,yheight=0,zdepth=0)
*
* %Parameters
* INPUT PARAMETERS
*
* xwidth:     [m] Width of block.
* yheight:    [m] Height of block.
* zdepth:     [m] Thickness of block.
* radius:     [m] Radius of cylinder or sphere.
* material_datafile: [str] File where the material parameters for the filter may be found. Format is similar to what may be found off the NIST website. [Be.txt]
* geometry:   [str] File containing the polygon definition of a general shape object. When xwidth is also given, the object is rescaled accordingly (OFF/PLY)
* fixed_delta:[0/1] Use a fixed delta to compute refraction - useful for debugging.
* refraction: [0/1] If nonzero, refraction is enabled.
*
* %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>
* %Linkink
* Powercrust https://www.cs.ucdavis.edu/~amenta/powercrust.html
* %Link
* material datafile obtained from http://physics.nist.gov/cgi-bin/ffast/ffast.pl
* %End
*******************************************************************************/
DEFINE COMPONENT Filter

SETTING PARAMETERS (refraction=1,fixed_delta=0,string material_datafile="Be.txt",
    string geometry=0,xwidth=0,yheight=0,zdepth=0,radius=0)

/* 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
  
  #pragma acc routine
  int filter_refract(double *kx, double *ky, double *kz, double nx, double ny, double nz, double delta0, double delta1){
    const double n2=scalar_prod(nx,ny,nz, nx,ny,nz);
    const double k=sqrt(scalar_prod(*kx,*ky,*kz, *kx,*ky,*kz));
    const double nr=(1.0-delta0)/(1.0-delta1);
    double kxi=*kx;
    double kyi=*ky;
    double kzi=*kz;
    double s;

    NORM(kxi,kyi,kzi);
    if(n2!=1){
      NORM(nx,ny,nz);
    }
    //kinv=1.0/sqrt(k2);
    s=scalar_prod(nx,ny,nz,kxi,kyi,kzi);
    if(s>0){
      /*n points in the direction of k - i.e. into material 1, so use -n instead*/
      double sinth2=nr*nr*(1.0-(s)*(s));
      *kx=nr* (kxi) - (nr*(-s)+sqrt(1.0-sinth2))*(-nx);
      *ky=nr* (kyi) - (nr*(-s)+sqrt(1.0-sinth2))*(-ny);
      *kz=nr* (kzi) - (nr*(-s)+sqrt(1.0-sinth2))*(-nz);
    } else {
      /*n points oppsite to k - i.e. out of material 1, into mat. 0*/ 
      double sinth2=nr*nr*(1.0-s*s);
      *kx=nr* (kxi) - (nr*s+sqrt(1.0-sinth2))*nx;
      *ky=nr* (kyi) - (nr*s+sqrt(1.0-sinth2))*ny;
      *kz=nr* (kzi) - (nr*s+sqrt(1.0-sinth2))*nz;
    }
    *kx *=k;
    *ky *=k;
    *kz *=k;
    return 1;
  }

%}

DECLARE
%{
  int Z;
  int mu_c;
  double Ar;
  double rho;
  double delta_prefactor;
  double betafact;
  t_Table table;
  int shape;
  off_struct offdata;
  int filter_row;
%}

INITIALIZE
%{
  int status=0;
 
  /*if geometry is specified - try to read the off/ply file and set shape*/
  if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) {
    if (off_init(geometry, xwidth, yheight, zdepth, 0, &(offdata)))
      shape=ANY; 
  } else if (radius && yheight){
    shape=CYLINDER;
  } else if (radius){
    shape=SPHERE;
  } else if (xwidth && yheight && zdepth){
    shape=CUBE;
  } else {
    fprintf(stderr,"Error (%s): Filter has zero effective volume\n",NAME_CURRENT_COMP);
    exit(-1);
  }
  
  if ( (status=Table_Read(&(table),material_datafile,0))==-1){
    fprintf(stderr,"Error (%s): Could not parse file \"%s\".\n",NAME_CURRENT_COMP,material_datafile);
    exit(-1);
  }
  char **header_parsed;
  header_parsed=Table_ParseHeader(table.header,"Z","A[r]","rho","Z/A","sigma[a]",NULL);
  if (header_parsed[0]){Z=strtol(header_parsed[0],NULL,10);}
  if (header_parsed[1]){Ar=strtod(header_parsed[1],NULL);}
  if (header_parsed[2]){rho=strtod(header_parsed[2],NULL);}
  else{fprintf(stderr,"Warning(%s): %s not found in header of %s, set to 1\n",NAME_CURRENT_COMP,"rho",material_datafile);rho=1;}
 
  /*If we are using a three column format assume only attenuation
   I.e. set the delta prefactor to 0*/
  if (table.columns==3) {
    mu_c=1;
    delta_prefactor=0;
  } else {
    mu_c=5;
    delta_prefactor= NA*(rho*1e-24)/Ar * 2.0*M_PI*RE;
  }
  if(10*(double)mcget_ncount()>(table.max_x-table.min_x)/(table.step_x)){
    /*it is probably worth Rebinning*/
    printf("%s: Rebinning material data (%s) for faster tracing...",NAME_CURRENT_COMP,material_datafile); 
    Table_Rebin(&(table));
    printf("done\n");
  }
  
%}

TRACE
%{
  double alpha,e,k,mu,mu0,delta,beta,f;
  double l0,l1,nx,ny,nz;
  int i,status=0;
  if (shape==CYLINDER) {
    status=cylinder_intersect(&l0,&l1,x,y,z,kx,ky,kz,radius,yheight);
  } else if (shape==SPHERE) {
    status=sphere_intersect(&l0,&l1,x,y,z,kx,ky,kz,radius);
  } else if (shape==CUBE) {
    status=box_intersect(&l0,&l1,x,y,z,kx,ky,kz,xwidth,yheight,zdepth);
  } else if (shape==ANY) {
    status = off_x_intersect(&l0, &l1, NULL, NULL, x, y, z, kx, ky, kz, offdata );
  }

  if(status){
    PROP_DL(l0);
    l1-=l0;

    k=sqrt(kx*kx+ky*ky+kz*kz);
    e=k*K2E;
    /*Material's Number Density of Electrons [e/A^3] incl f' scattering length correction*/
    f=Table_Value(table,e,1);
    delta = (fixed_delta ? fixed_delta : f/(k*k) * delta_prefactor);
    mu0=Table_Value(table,e,mu_c)*rho*1e2;
    beta = mu0/(2*k);

    /*change direction of ray due to refraction*/
    if(refraction){
      switch (shape){
        case CUBE:
          if( fabs(x-xwidth*0.5)<FLT_EPSILON || fabs(x+xwidth*0.5)<FLT_EPSILON ){
            nx=1.0;ny=0.0;nz=0.0;
          }else if( fabs(y-yheight*0.5)<FLT_EPSILON || fabs(y+yheight*0.5)<FLT_EPSILON ){ 
            nx=0.0;ny=1.0;nz=0.0;
          }else if( fabs(z-zdepth*0.5)<FLT_EPSILON || fabs(z+zdepth*0.5)<FLT_EPSILON ){ 
            nx=0.0;ny=0.0;nz=1.0;
          };
          break;
        default:
          1;
      }
      status=filter_refract(&kx,&ky,&kz,nx,ny,nz,0.0,delta);
      SCATTER;

      /*change of direction - we must now intersect again*/
      if (shape==CUBE) {
          status=box_intersect(&l0,&l1,x,y,z,kx,ky,kz,xwidth,yheight,zdepth);
      }else if (shape==ANY){
          status = off_x_intersect(&l0, &l1, NULL, NULL, x, y, z, kx, ky, kz, offdata );
      }
    }

    /*change k ' we're now inside some material*/
    kx*=1-delta;
    ky*=1-delta;
    kz*=1-delta;

    /*propagate to next intersection with filter (assumed to be exit)*/
    p*=exp(-mu0*l1);
    SCATTER;
    PROP_DL(l1);
    
    /*change direction of ray due to refraction*/
    if(refraction){
      switch (shape){
        case CUBE:
          if( fabs(x-xwidth*0.5)<FLT_EPSILON || fabs(x+xwidth*0.5)<FLT_EPSILON ){
            nx=1.0;ny=0.0;nz=0.0;
          }else if( fabs(y-yheight*0.5)<FLT_EPSILON || fabs(y+yheight*0.5)<FLT_EPSILON ){ 
            nx=0.0;ny=1.0;nz=0.0;
          }else if( fabs(z-zdepth*0.5)<FLT_EPSILON || fabs(z+zdepth*0.5)<FLT_EPSILON ){ 
            nx=0.0;ny=0.0;nz=1.0;
          };
          break;
        default:
          1;
      }
      SCATTER;
      status=filter_refract(&kx,&ky,&kz,nx,ny,nz,delta,0.0);
    }

    /*correct for optical path length through the material*/
    /*change k back*/
    kx/=1-delta;
    ky/=1-delta;
    kz/=1-delta;
  }
%}

FINALLY
%{
  Table_Free(&(table));
%}

MCDISPLAY
%{
  
  if (shape==CUBE){
    box(0,0,0,xwidth,yheight,zdepth,0, 0, 1, 0);
  } else if (shape==CYLINDER){
    circle("xy",0, yheight/2.0,0,radius);
    circle("xy",0,-yheight/2.0,0,radius);
    line( radius,yheight/2.0,0, radius,-yheight/2.0,0);
    line(-radius,yheight/2.0,0,-radius,-yheight/2.0,0);
    line(0,yheight/2.0, radius, 0,-yheight/2.0, radius);
    line(0,yheight/2.0,-radius, 0,-yheight/2.0,-radius);
  } else if (shape==SPHERE){
    circle("xy",0,0,0,radius);
    circle("xz",0,0,0,radius);
    circle("yz",0,0,0,radius);
  } else if (shape == ANY) {	/* OFF file */
    off_display(offdata);
  }
  
%}

END