/*******************************************************************************
*
* 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: TwinKB_ML
*
* %Identification
*
* Written by: Jana Baltser, Peter Willendrup, Anette Vickery, Andrea Prodi, Erik Knudsen, Jesper Buch Jensen
* Date: May 2012
* Version: 1.0
* Origin: NBI
*
* Montel optic model (aka side-by-side Kirkpatrick Baez)
* 
* %Description
* Models a Montel optic, or Twin Kirkpatrick Baez mirror optic (hence the component name).
* The mirror are fully abutting, i.e. there's is no gap between them, and perfectly elliptic.
*
* Reads reflectivity values from a data input file for a W/B4C multilayer.
* The multilayer code reflects ray in an ideal geometry, the reflectivity datafile accounts for surface roughness, sigma.
*
* The mirror is positioned such that the long axis of the mirror elliptical surface coincides with the 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: TwinKB_ML( theta=1.2, s1=.045 , s2=.9 , length=0.06 , width=0.2 , R0=0 , reflectivity_datafile="Ref_W_B4C.txt")
*
* %Parameters
* Input parameters:
* theta:  [deg] Incident angle
* s1:     [m]   Distance from the source to the multilayer
* s2:     [m]   Focusing distance of the multilayer
* length: [m]   Length of the mirrors
* width:  [m]   Width of the mirror along x-axis
* R0:     [0-1] Constant reflectivity, R0=1 for an ideal situation. If R0=0, the code reads the reflectivity from the datafile
* reflectivity_datafile: [str] File which contains reflectivities as a function of q.
*
* %End
*******************************************************************************/

DEFINE COMPONENT TwinKB_ML

SETTING PARAMETERS (string reflectivity_datafile="Ref.txt", theta=1.2,s1,s2,length=0.6,width=0.2,R0=0)

//STATE PARAMETERS (x,y,z,kx,ky,kz,phi,Ex,Ey,Ez,p)

SHARE 
%{
   %include "read_table-lib"
  
  struct prms_multilayer_elliptic{
     double *Q;
     double **R;
   };
   
  /*something that would be relevant for ALL elliptical mirrors*/
  /* coordinate transformation McXtrace-Ellipse (ME) and Ellipse-McXtrace(EM) functions */
  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;
  }
  
  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;
  Rotation Q1;
  Rotation Q2;
  t_Table ref_table;
%}

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; 

  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)); 
  
  // reflectivity datafile parsing
  int status=0;
   
  if ((status=Table_Read(&ref_table, reflectivity_datafile,0))==-1){
    fprintf(stderr,"Error: Could not parse the file \'%s\' in COMP %s\n",reflectivity_datafile,NAME_CURRENT_COMP);
    exit(-1);
  }
  
%}

TRACE
%{
  double K,vink; 
  double x_e1,y_e1,z_e1,kx_e1,ky_e1,kz_e1;	// beginning coordinates transformed into the ellipse system
  double x_e2,y_e2,z_e2,kx_e2,ky_e2,kz_e2;	// kvector transformed into the ellipse system, hence 
  
  double A,B,C,D,l01,l11,l02,l12;
  double x_test1,y_test1,z_test1,x_test2,y_test2,z_test2,dist;	// intersection with the elliptical surface
  double nx,ny,nz;
  double kxn,kyn,kzn;		// reflected ray's kvector
 
  int status1,status2,bounce; 
  
  /* get the photon's coordinates and kvector in the ellipse frame */
  K=sqrt(kx*kx+ky*ky+kz*kz); 
  
  bounce=CHAR_MAX;
  while (bounce){
    bounce=0;
    /*switch to the ellipsoid frames. Note that the order of x and y has been swapped in the second set of calls*/
    CoordTransME(&x_e1,&y_e1,&z_e1,x,y,z,Z0,Y0,xi);
    CoordTransME(&kx_e1,&ky_e1,&kz_e1,kx,ky,kz,0,0,xi);
    NORM(kx_e1,ky_e1,kz_e1);
    CoordTransME(&y_e2,&x_e2,&z_e2,y,x,z,Z0,Y0,xi);
    CoordTransME(&ky_e2,&kx_e2,&kz_e2,ky,kx,kz,0,0,xi);
    NORM(kx_e2,ky_e2,kz_e2);
#ifdef MCDEBUG
    printf("coord transform1: r=(%g %g %g) k=(%g %g %g) => r=(%g %g %g) k=(%g %g %g) Z0,Y0=(%g,%g)\n",x,y,z,kx,ky,kz,x_e1,y_e1,z_e1,kx_e1,ky_e1,kz_e1,Z0,Y0);
    printf("coord transform2: r=(%g %g %g) k=(%g %g %g) => r=(%g %g %g) k=(%g %g %g)\n",x,y,z,kx,ky,kz,x_e2,y_e2,z_e2,kx_e2,ky_e2,kz_e2);
#endif
    double QQ[3][3]={{1,0,0},{0,1,0},{0,0,1}};
    /*compute intersections with the ellipsoid surfaces that contain the mirror surfaces
      using 1e6 as a half-axis to emulate something flat in that dimension*/
    status1=ellipsoid_intersect(&l01,&l11,x_e1,y_e1,z_e1,kx_e1,ky_e1,kz_e1,1e6,b,a,QQ);
    status2=ellipsoid_intersect(&l02,&l12,x_e2,y_e2,z_e2,kx_e2,ky_e2,kz_e2,b,1e6,a,QQ);
#define SWAP(a,b) \
    do { \
      double tmp=(a); \
      (a)=(b);(b)=tmp;\
    } while(0)\



    if (status1) {
      if (l01>0){
        double dl=l01;
        double xx,yy,zz;
        xx=x_e1+kx_e1*dl; yy=y_e1+ky_e1*dl; zz=z_e1+kz_e1*dl;
        if ((yy)<=0 && xx>0 && xx<width && fabs(zz-Z0)<length/2.0){    
#ifdef MCDEBUG
          printf("hit for intersection 1,0 %g %g (%g %g %g)\n",l01,l11,xx,yy,zz);
#endif
          ABSORB;
        }
      }
      if(l11>0){
        double dl=l11;
        double xx,yy,zz;
        xx=x_e1+kx_e1*dl; yy=y_e1+ky_e1*dl; zz=z_e1+kz_e1*dl;
        if ((yy)<=0 && xx>0 && xx<width && fabs(zz-Z0)<length/2.0){    
          if(bounce!=1){
            bounce|=1;
          }
          x_test1=xx;y_test1=yy;z_test1=zz;
#ifdef MCDEBUG
          printf("hit for intersection 1,1 %g %g (%g %g %g)\n",l01,l11,xx,yy,zz);
#endif
        }
      }
    }
    
    if (status2) {
      if(l02>0){
        double dl=l02;
        double xx,yy,zz;
        xx=x_e2+kx_e2*dl; yy=y_e2+ky_e2*dl; zz=z_e2+kz_e2*dl;
        if ((xx)<=0 && yy>0 && yy<width && fabs(zz-Z0)<length/2.0){
#ifdef MCDEBUG
          printf("hit for intersection 2,0 %g %g (%g %g %g)\n",l02,l12,xx,yy,zz);
          ABSORB;
#endif
        } 
      }
      if (l12>0){
        double dl=l12;
        double xx,yy,zz;
        xx=x_e2+kx_e2*dl; yy=y_e2+ky_e2*dl; zz=z_e2+kz_e2*dl;
        if ((xx)<=0 && yy>0 && yy<width && fabs(zz-Z0)<length/2.0){    
          if (bounce!=2) bounce|=2;
#ifdef MCDEBUG
          printf("hit for intersection 2,1 %g %g (%g %g %g)\n",l02,l12,xx,yy,zz);
#endif
          x_test2=xx;y_test2=yy;z_test2=zz;
        }
      }
    }
    /*if we're about to hit both mirrors - pick the first one*/
    if (bounce==3){
      if (l11<l12){
        bounce=1;
      }else{
        bounce=2;
      }
    } else if (!bounce)
      continue;
  
    /*propagate to the selected mirror and reflect
      first store the old wavevector though*/
    double kxo,kyo,kzo;
    kxo=kx;kyo=ky,kzo=kz;
    if (bounce==1){
      PROP_DL(l11);
      SCATTER;
      nx=0;
      if (fabs(z_test1)<FLT_EPSILON){
        ny=-1;
        nz=0;
      } else {
        ny=(a*a*y_test1)/(b*b*z_test1);
        nz=1.0;
      }
      NORM(nx,ny,nz);
      vink=scalar_prod(nx,ny,nz,kx_e1,ky_e1,kz_e1); 
      kxn=kx_e1-2.0*vink*nx;
      kyn=ky_e1-2.0*vink*ny;
      kzn=kz_e1-2.0*vink*nz;
      NORM(kxn,kyn,kzn); 
#ifdef MCDEBUG
      printf("1: vink=%g k=(%g %g %g), k_e=(%g %g %g), n=(%g %g %g), kn=(%g %g %g) ",vink,kx,ky,kz,kx_e1,ky_e1,kz_e1,nx,ny,nz,kxn,kyn,kzn); 
#endif
      CoordTransEM(&kx,&ky,&kz,kxn,kyn,kzn,0,0,xi);
    }else if (bounce==2){
      PROP_DL(l12);
      SCATTER;
      ny=0;
      if (fabs(z_test2)<FLT_EPSILON){
        nx=1;
        nz=0;
      } else {
        nx=(a*a*x_test2)/(b*b*z_test2);
        nz=1.0;
      }
      NORM(nx,ny,nz);
      vink=scalar_prod(nx,ny,nz,kx_e2,ky_e2,kz_e2); 
      kxn=kx_e2-2.0*vink*nx;
      kyn=ky_e2-2.0*vink*ny;
      kzn=kz_e2-2.0*vink*nz;
      NORM(kxn,kyn,kzn); 
#ifdef MCDEBUG
      printf("2: vink=%g k=(%g %g %g), k_e=(%g %g %g), n=(%g %g %g), kn=(%g %g %g) ",vink,kx,ky,kz,kx_e2,ky_e2,kz_e2,nx,ny,nz,kxn,kyn,kzn); 
#endif
      CoordTransEM(&ky,&kx,&kz,kyn,kxn,kzn,0,0,xi);
    }
    
    kx=K*kx;
    ky=K*ky;
    kz=K*kz;
#ifdef MCDEBUG
      printf("ko=(%g %g %g)\n",kx,ky,kz);
#endif
    PROP_DL(FLT_EPSILON); 
    if (R0){
        p*=R0;
    }else{
        // Implementing the reflectivity.
        double dSingle,lambda;	// dSingle - thickness of one layer of the multilayer, AA
        double Q,Refl;
        double ex1=z+Z0, ex2=-(-b*b*a*a+b*b*ex1*ex1)/(a*a);	// ex1 & ex2 - expressions to simplify the equation dSingle(z).

        lambda=2*M_PI/K;   
        dSingle= (0.5*lambda)/cos( acos(sqrt(ex2)/sqrt(ex2 + c*c + 2*c*ex1 + ex1*ex1) ) - atan((ex1*b*b)/(sqrt(ex2)*a*a)) );
        Q=sqrt((kx-kxo)*(kx-kxo)+(ky-kyo)*(ky-kyo)+(kz-kzo)*(kz-kzo));
#ifdef MCDEBUG 
        printf("ML d-spacing dSingle=%g,z=%g,K=%g,lambda=%g, wave-vector transfer Q=%g\n",dSingle,z,K,lambda,Q);
#endif   
        Refl=Table_Value(ref_table, Q, 1);
        p*=Refl;
    } 
  }
%}

FINALLY
%{
    Table_Free(&ref_table);
%}

MCDISPLAY
%{
  /*
  rectangle("xz",0,0,0,width,length); */
  int i,j,N=10;
  double w_2=width/2.0;
  double l_2=length/2.0;
  const double x_e[]={w_2,0,0,2*w_2,2*w_2};
  double y_e[]={Y0,0,0,0,0};
  double z_e[]={Z0,0,0,0,0};
  if (s1<s2){
    y_e[1]=Y0-l_2*sin(xi);
    y_e[2]=y_e[3]=Y0+l_2*sin(xi);
    y_e[4]=y_e[1];
    z_e[1]=Z0-l_2*cos(xi);
    z_e[2]=z_e[3]=Z0+l_2*cos(xi);
    z_e[4]=z_e[1];
  }else if (s2<=s1){
    y_e[1]=Y0+l_2*sin(xi);
    y_e[2]=y_e[3]=Y0-l_2*sin(xi);
    y_e[4]=y_e[1];
    z_e[1]=Z0-l_2*cos(xi);
    z_e[2]=z_e[3]=Z0+l_2*cos(xi);
    z_e[4]=z_e[1];
  }
  double xx[5],yy[5],zz[5];
  for (i=0;i<5;i++){
    CoordTransEM(xx+i,yy+i,zz+i,x_e[i],y_e[i],z_e[i],Z0,Y0,xi);
    //printf("%d; %g %g %g => %g %g %g\n",i,x_e[i],y_e[i],z_e[i],xx[i],yy[i],zz[i]); 
  }

  multiline(5,xx[1],yy[1],zz[1],xx[2],yy[2],zz[2],xx[3],yy[3],zz[3],xx[4],yy[4],zz[4],xx[1],yy[1],zz[1]);
  multiline(3,xx[1],yy[1],zz[1],xx[0],yy[0],zz[0],xx[3],yy[3],zz[3]);
  multiline(3,xx[2],yy[2],zz[2],xx[0],yy[0],zz[0],xx[4],yy[4],zz[4]);

  for (i=0;i<5;i++){
    CoordTransEM(yy+i,xx+i,zz+i,x_e[i],y_e[i],z_e[i],Z0,Y0,xi);
    //printf("%d; %g %g %g => %g %g %g\n",i,y_e[i],x_e[i],z_e[i],xx[i],yy[i],zz[i]); 
  }

  multiline(5,xx[1],yy[1],zz[1],xx[2],yy[2],zz[2],xx[3],yy[3],zz[3],xx[4],yy[4],zz[4],xx[1],yy[1],zz[1]);
  multiline(3,xx[1],yy[1],zz[1],xx[0],yy[0],zz[0],xx[3],yy[3],zz[3]);
  multiline(3,xx[2],yy[2],zz[2],xx[0],yy[0],zz[0],xx[4],yy[4],zz[4]);
 
%}

END