/*******************************************************************************
*
* McXtrace, X-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: Abs_objects
*
* %Identification
* Written by: Erik Knudsen
* Date: Jan 24, 2011
* Version: $Revision$
* Origin: DTU Physics
* Release: McXtrace 1.1
*
* Blocks of attenuating material in off format
*
* %Description
*
* This component is a model of 1 or more off-shaped blocks attenuating the x-ray beam.
* Which shapes are present and their relative positions are described in a file of the follwing format:
* #objects
* Material-filename  OFF-filename x y z xwidth yheight zdepth
* ...
*
* An example is;
* 2
* Be.txt  cube.off  0 0.01 0    0 0 0
* Rh.txt  chess.off 0 -0.01 0   0 0 0
*
* A xwidth etc of zero means use whatever dimensions are in the off/ply file. 
* If xwidth (or yheight or zdepth) is nonzero, the component scales the object 
* to fill that dimension. The xyz coordinates indicate the object shift.
*
* Example: Abs_objects(objects="input_abs_objects_template.dat")
*
* %Parameters
* INPUT PARAMETERS
* objects:    [str] Input file where the off-shapes are defined.
* refraction: [str] Flag to enable refraction at interfaces.
*
* %End
*******************************************************************************/
DEFINE COMPONENT Abs_objects

SETTING PARAMETERS (int refraction=1, string objects="input_abs_objects_template.dat")

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

SHARE
%{
#ifndef MX_ABS_OBJECTS_MAX_INTERSECTS
#define MX_ABS_OBJECTS_MAX_INTERSECTS 256
#endif

  %include "read_table-lib"
  %include "interoff-lib"

  struct _abs_object_data_struct {
    int mu_c;
    double Z,Ar, rho;
    double delta_prefactor,betafact;
    double xx,obj_xwidth,yy,obj_yheight,zz,obj_zdepth;
  };

#pragma acc routine
  int abs_object_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);
    }
    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 mu_c[128];
  double Z[128];
  double Ar[128];
  double rho[128];
  double delta_prefactor[128];
  double betafact[128];
  double xx[128];
  double obj_xwidth[128];
  double yy[128];
  double obj_yheight[128];
  double zz[128];
  double obj_zdepth[128];
  int objectsc;
  struct _abs_object_data_struct prms[128];
  t_Table mat_table[128];
  off_struct offdata[128];
%}


INITIALIZE
%{
  FILE *fp=Open_File(objects,"r",NULL);
  if (fp == NULL) {
    fprintf(stderr,"Error: %s: Could not open file \"%s\".\n", NAME_CURRENT_COMP, objects);
    exit(-1);
  }
  char line[512];
  fgets(line,512,fp);
  objectsc=atoi(line);

  int i;
  for (i=0;i<objectsc;i++){
    fgets(line,512,fp);
    char off_fn[512],mat_fn[512];
    float xx1,yy1,zz1,xwidth,yheight,zdepth;
    int status;
    //sscanf(line,"%s %s \n", mat_fn,off_fn);
    status=sscanf(line,"%s %s %g %g %g %g %g %g\n",
      mat_fn,off_fn, &xx1,&yy1,&zz1,&xwidth,&yheight,&zdepth);
    
    if (status < 8) {
      fprintf(stderr,"Error: %s: Invalid line %i in file \"%s\".\n",
        NAME_CURRENT_COMP, i+1, objects);
      exit(-1);
    }
    xx[i]=xx1;         yy[i]=yy1;           zz[i]=zz1;
    obj_xwidth[i]=xwidth; obj_yheight[i]=yheight; obj_zdepth[i]=zdepth;

    if ( (status=Table_Read(&(mat_table[i]),mat_fn,0))==-1){
      fprintf(stderr,"Error: %s: Could not parse file \"%s\".\n",
        NAME_CURRENT_COMP, mat_fn);
      exit(-1);
    }
    char **header_parsed;
    header_parsed=Table_ParseHeader(mat_table[i].header,"Z","A[r]","rho","Z/A","sigma[a]",NULL);
    if (header_parsed[0]){Z[i]=strtod(header_parsed[0],NULL);}
    if (header_parsed[1]){Ar[i]=strtod(header_parsed[1],NULL);}
    if (header_parsed[2]){rho[i]=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",mat_fn);rho[i]=1;}
    /*which columns holds the mus*/
    mu_c[i]=5;
    if (mat_table[i].columns==3) {
      /*format for compound material (XCOM) is assumed. For now set delta_prefactor=0 => delta=0 (no refraction).*/
      mu_c[i]=1;
      delta_prefactor[i]=0;
    }else{
      delta_prefactor[i]= NA*(rho[i]*1e-24)/Ar[i] * 2.0*M_PI*RE;
    }
    /*read in the off/ply data from file. The literal 1 is to not center the object in local (0,0,0)*/
    if (!off_init(off_fn, obj_xwidth[i], obj_yheight[i], obj_zdepth[i], 1, &(offdata[i]))) {
      fprintf(stderr,"Error: %s: Could not parse file \"%s\".\n",
        NAME_CURRENT_COMP, off_fn);
      exit(-1);
    }
    printf("%s: Volume[%d], material='%s' geometry='%s'\n", 
      NAME_CURRENT_COMP, i+1,mat_fn,off_fn);
  }
  fclose(fp);
%}

TRACE
%{

  double alpha,e,k,kx_vac,ky_vac,kz_vac,mu,mu0,delta,beta,f;
  double l0,l1,nx,ny,nz;
  Coords c0,c1;
  int i=1,j,m,status=0,exit=0;
  int intersect=0,entry_exit=0;
  double lmin;
  int idx;

  struct intersection_data_struct {
    double l0,dl;
    double mu0,beta,delta;
    double _kx,_ky,_kz;
    int idx;
    int entry_exit;
  } intersections[MX_ABS_OBJECTS_MAX_INTERSECTS],intersection_tmp;

  k=sqrt(kx*kx+ky*ky+kz*kz);
  e=k*K2E;
  kx_vac=kx;
  ky_vac=ky;
  kz_vac=kz;
  for (i=0;i<MX_ABS_OBJECTS_MAX_INTERSECTS;i++){
    intersections[i].dl=0;
  }

  int intersectionc=0;
  for (i=0;i<objectsc;i++){
    l0=0;
    l1=0;
    status = off_x_intersect(&l0, &l1, &c0,&c1, x-xx[i], y-yy[i], z-zz[i], kx, ky, kz, offdata[i] );
    if (status>1){/*ignore singlet intersections*/
      if (l0<0 && l1 >0){
        /*we're already inside some thing. Set l0 to 0;*/
        l0=0;
      }
      if (l0>=0) {
        int j=0;
        intersectionc++;

	if(intersectionc>=MX_ABS_OBJECTS_MAX_INTERSECTS){
	  fprintf(stderr, "Warning: %s: Intersection storage overrun. Rest of scene will be ignored. Please examine results carefully and/or recompile with \"-DMX_ABS_OBJECTS_MAX_INTERSECTS=N\" where N>256",NAME_CURRENT_COMP);
	  break;
	}
        /*a dl of 0 marks an unused slot*/
        while ( intersections[j].dl!=0 && (intersections[j].l0<l0 || (intersections[j].l0==l0 && intersections[j].dl<(l1-l0)) ) ) {
          j++;
        }

        int m=objectsc*2;//MX_ABS_OBJECTS_MAX_INTERSECTS-1;
        while (m>j){
          if(intersections[m-1].dl!=0){
            intersections[m]=intersections[m-1];
          }
	  m--;
        }
        intersections[j].l0=l0;
        intersections[j].dl=l1-l0;
	intersections[j].idx=i;

        /*precompute some material constants*/
        double f;
        f=Table_Value(mat_table[i],e,1);
        intersections[j].delta = f/(k*k) * delta_prefactor[i];
        intersections[j].mu0=Table_Value(mat_table[i],e,mu_c[i])*rho[i]*1e2;
        intersections[j].beta = mu0/(2*k);
        /*change k we're now inside some material*/
        intersections[j]._kx=kx_vac*(1.0-intersections[j].delta);
        intersections[j]._ky=ky_vac*(1.0-intersections[j].delta);
        intersections[j]._kz=kz_vac*(1.0-intersections[j].delta);
      }
    }
  }
  /*now we can loop through the intersection times and do propagation accordingly*/
  double lsum=0;
  for (j=0;j<intersectionc;j++){
    /*set k*/
    /*if l0-lsum>0 propagate to the first intersection in vacuum. I.e. prop. to object entry.*/
    double dl=intersections[j].l0-lsum;
    if (dl>0){
      PROP_DL(dl);
      lsum+=dl;
      SCATTER;
    }
    //kx=intersections[j].kx;
    //ky=intersections[j].ky;
    //kz=intersections[j].kz;
    dl=intersections[j].dl;
    p*=exp(-intersections[j].mu0*dl);
    PROP_DL(dl);
    lsum+=dl;
    SCATTER;/*we set scatterpoints at entry and exit*/
    //kx=kx_vac;
    //ky=ky_vac;
    //kz=kz_vac;

    /*we may reenter this object at a later stage so check for intersection again*/
    l0=-1;l1=-1;status=0;
    i=intersections[j].idx;
    status = off_x_intersect(&l0, &l1, &c0, &c1, x-xx[i], y-yy[i], z-zz[i], kx, ky, kz, offdata[i] );
    if (status>1 && l0>FLT_EPSILON && l1>FLT_EPSILON){/*ignore singlet intersections*/
      /*we do in fact hit the object again*/
      intersection_tmp=intersections[j];/*store the computed material values for later*/
      intersection_tmp.dl=l1-l0;
      if (l0<0 && l1 >0){
        /*already inside something. Set l0 to 0;*/
        l0=0;
      }
      if (l0>=0) {
        l0=l0+lsum;/*correct for the fact that we have travelled in the component already*/
	int m=j;/*start at the present object and go on to sort in this intersection*/
        /*a dl of 0 marks an unused slot*/
	while ( intersections[m+1].dl!=0 && (intersections[m+1].l0<l0 || (intersections[m+1].l0==l0 && intersections[m+1].dl<(l1-l0)) ) ) {
          intersections[m]=intersections[m+1];
	  m++;
        }
	intersections[m]=intersection_tmp;/*reinsert values but override new intersection lengths*/
        intersections[m].l0=l0;
      }
      /*we have added another intersection so don't progress j*/
      j--;
    }
  }
%}

FINALLY
%{
  int i;
  for (i=0;i<objectsc;i++){
    Table_Free(&(mat_table[i]));
  }
%}

MCDISPLAY
%{
  int i;
  for (i=0;i<objectsc;i++){
    off_display(offdata[i]);
  }
%}

END