/*******************************************************************************
*
* 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: Multilayer_elliptic
*
* %I
*
* Written by: Jana Baltser, Peter Willendrup, Anette Vickery, Andrea Prodi, Erik Knudsen
* Date: February 2011
* Version: 1.0
* Origin: NBI
*
* Elliptic multilayer mirror (in XZ)
* 
* %Description
* Reads reflectivity values from a data input file (Ref.dat) for a Si/W multilayer.
* The multilayer code reflects ray in an ideal geometry, does not include surface imperfections
*
* The mirror is positioned such that the long axis of the mirror elliptical surface coincides with
* 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: Multilayer_elliptic(
*   coating = "Ref_W_B4C.txt", theta = 1.2,
*   s1 = 1, s2 = 2, length = 0.1, width = 0.1, R0 = 1,
*   Emin=7, Emax=10, Estep=0.05)
*
* %Parameters
* Input parameters:
* theta: [deg] Design angle of incidence.
* s1:    [m]   Design distance from the source to the multilayer.
* s2:    [m]   Design focusing distance of the multilayer.
* zdepth:[m]   Length of the mirror along Z.
* xwidth:[m]   Width of the mirror along X-axis.
* Gamma: [ ]   High electron density fraction of bilayer (in kinematical appr.).
* Lambda:[m]   Thickness of bilayer (in kinematical appr.).
* rho_AB:[ ]   Number electron density constrast in bilayer (in kinematical appr.).
* N:     [1]   Number of bilayers (in kinematical appr.).
* coating: [str] Datafile containing reflectivity values as a function of q and E.
* Emin:  [keV] Lower limit of energy interval in datafile. Overrides what's written in the datafile header.
* Emax:  [keV] Upper limit of energy interval in datafile. Overrides what's written in the datafile header.
* Estep: [keV] Step between energy sample points in datafile. Overrides what's written in the datafile header.
* R0:    [1]   Maximal reflectivity
* length:  [m]   alternate name for zdepth (obsolete)
* width:   [m]   alternate name for xwidth (obsolete)
* %End
*******************************************************************************/

DEFINE COMPONENT Multilayer_elliptic

SETTING PARAMETERS (string coating="Ref_W_B4C.txt",
    theta=1.2,s1=0,s2=0,length=0.5,width=0.2,R0=1, 
    Emin=-1, Emax=-1, Estep=-1, Gamma=0, Lambda=0, rho_AB=0, int N=0,
    xwidth=0, zdepth=0)

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

SHARE 
%{
#include <complex.h>
  %include "read_table-lib"
  %include "reflectivity-lib"
  /*something that would be relevant for ALL elliptical mirrors*/
  /* coordinate transformation McXtrace-Ellipse (ME) and Ellipse-McXtrace(EM) functions */
#pragma acc routine
  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;
  }

#pragma acc routine
  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;
  int kinematical;
  t_Reflec re;
%}

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;
  
  if (xwidth) width=xwidth;
  if (zdepth) length=zdepth;
  
  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)); 

  int status=0;
  if(!Gamma && !Lambda){
    /*refrain from using kinematical approximation - instead use reflectivity datafile*/
    kinematical = 0;
    /* reflectivity datafile parsing COATING_UNDEFINED - means set the type according to what is found in the file*/
    status=reflec_Init(&re,COATING_UNDEFINED,coating, NULL);
  }else{
    kinematical = 1;
    re.type=KINEMATIC;
    status=reflec_Init_kinematic(&(re), N, Gamma, Lambda, rho_AB);/*number of layers, ratio of high e-density material in layer, thickness of layer, e-density contrast*/
  }
%}

TRACE
%{
  double K,vink; 
  double x_el,y_el,z_el;	// beginning coordinates transformed into the ellipse system
  double kx_el,ky_el,kz_el;	// kvector transformed into the ellipse system, hence 
  
  double A,B,C,D,t0,t1;
  double x_int,y_int,z_int,dist;	// intersection with the elliptical surface
  double nx,ny,nz;
  double kxn,kyn,kzn;		// reflected ray's kvector
  
  /* get the photon's coordinates and kvector in the ellipse frame */
  K=sqrt(kx*kx+ky*ky+kz*kz);
  
  CoordTransME(&x_el,&y_el,&z_el,x,y,z,Z0,Y0,xi);
  CoordTransME(&kx_el,&ky_el,&kz_el,kx,ky,kz,0,0,xi);
    
  NORM(kx_el,ky_el,kz_el);
  
  /*intersection calculation*/
  A=b*b*kz_el*kz_el+a*a*ky_el*ky_el;
  B=2.0*(z_el*kz_el*b*b+y_el*ky_el*a*a);
  C=b*b*z_el*z_el+a*a*y_el*y_el-a*a*b*b;
  D=B*B-4*A*C;
  if (D>=0){
    t0=(-B-sqrt(D))/(2*A);
    t1=(-B+sqrt(D))/(2*A);
    if (t0<0 && t1>=0) {
	double ttmp=t0;
        t0=t1; t1=ttmp;
    }
    /* check whether our intersection lies within the boundaries of the mirror*/
    x_int=x_el+kx_el*t0;
    y_int=y_el+ky_el*t0;
    z_int=z_el+kz_el*t0;
    
    if (y_int>=0 && fabs(x_int)<=width/2){
	dist=sqrt((x_el-x_int)*(x_el-x_int)+(y_el-y_int)*(y_el-y_int)+(z_el-z_int)*(z_el-z_int));
	PROP_DL(dist); 
	
	if (fabs(z)<=length/2) { /*finally in business on the mirror! YAY! */
	  nx=0;
	  if (fabs(z_int)==0){
	      ny=1;
	      nz=0;
	  } else {
            ny=(a*a*y_int)/(b*b*z_int);
            nz=1.0;
	  }
	  NORM(nx,ny,nz);
	  vink=scalar_prod(nx,ny,nz,kx_el,ky_el,kz_el); 
	  kxn=kx_el-2.0*vink*nx;
	  kyn=ky_el-2.0*vink*ny;
	  kzn=kz_el-2.0*vink*nz;
	  NORM(kxn,kyn,kzn); 
	  
	  double kxo,kyo,kzo;
	  kxo=kx;kyo=ky,kzo=kz;
	  CoordTransEM(&kx,&ky,&kz,kxn,kyn,kzn,0,0,xi);
	 	  
	  kx=K*kx;
	  ky=K*ky;
	  kz=K*kz;
	  	  
	  double QQ,EE,Ref;
	  QQ=sqrt((kx-kxo)*(kx-kxo)+(ky-kyo)*(ky-kyo)+(kz-kzo)*(kz-kzo)); 
	  EE=K*K2E; 

          if(kinematical){
            /*
             * \Lambda: thickness of bilayer - following notation in Als-Nielsen/McMorrow
             * \Gamma:  \Gamma*\Lambda thickness of high electron density material.
             * r1(zeta) = 2 i r_0 \rho_{AB} \left(\frac{\Lambda^2 \Gamma}{\zeta}\right) \frac{\sin\left(\pi\Gamma\zeta\right)}{\pi\Gamma\zeta);
             */
            Ref=reflecq(re,QQ,0,0,0);
            if (Ref>1){
              /*Reflectivity can't be >1*/
              Ref=1.0;
            }
          }else{
            /*interpolate in table*/
            Ref=reflecq(re,QQ,0,0,0);
          }
	  
	  /* apply reflectivity */
	  p*=Ref; 
	  SCATTER;
	  
	} else {
	  RESTORE_XRAY(INDEX_CURRENT_COMP, x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p);
	} 
    }
  }
  
%}

MCDISPLAY
%{
  /*
  rectangle("xz",0,0,0,width,length); */
  int i,j,NN=10;
  double x0,y0,z0;
  double x1,y1,z1,z_el,y_el;
  
  x0=-width/2.0;
    
  for (i=0;i<=NN;i++){
    z0=-length/2.0; 
    z_el=cos(xi)*z0+Z0; //transformation to EL reference frame
    y_el=b*sqrt(1.0-((z_el*z_el)/(a*a)));
    y0=cos(xi)*(y_el-Y0)-sin(xi)*(z_el-Z0);
    line(-width/2,y0,z0,width/2,y0,z0);
    for (j=0;j<=NN;j++){
      z1=z0+length/NN;
      z_el=cos(xi)*z1+Z0;
      y_el=b*sqrt(1.0-((z_el*z_el)/(a*a))); 
      y1=cos(xi)*(y_el-Y0)-sin(xi)*(z_el-Z0); 
      line(x0,y0,z0,x0,y1,z1);
      y0=y1;
      z0=z1;
      line(-width/2,y1,z1,width/2,y1,z1);
    }
    x0=x0+width/NN;
  }   
 
%}

END