/*******************************************************************************
*
* McXtrace, xray-tracing package
*         Copyright, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Mirror_elliptic
* 
* %I
* Written by: Erik Knudsen
* Date: Feb 11, 2010
* Version: $Revision$
* Origin: Risoe
* Modified by: Bac Stephane, Padovani Antoine, 09th May 2022 : melded Bragg_crystal_bent.comp into here.
*
* Idealized elliptic mirror.
* 
* %D
* Takes a reflectivity (default=1) as input and reflects rays in a ideal geometry
* elliptic mirror.
* The mirror is positioned such that the a-axis of the mirror ellipsoid is on the
* x-axis, the b-axis is along the y-axis and the c is along the z-axis.
* The reference point of the mirror is the ellipsoid centre, offset by one half-axis along the y-axis.
* (See the component manual for a drawing). 
* This means that to position the mirror correctly.
* the user positions the ellipsoid governing the mirror shape, not the mirror itself.
* 
* %P
* INPUT PARAMETERS
* x_a:     [m]    1st short half axis (along x). Commonly set to zero, which really implies infinite value, so crystal is an elliptic cylinder.
* y_b:     [m]    2nd short half axis (along y), which is also the presumed near-normal direction, reflection near the y-z plane.
* z_c:     [m]    long  half axis (along z). Commonly a=0. b=c, which creates a circular cylindrical surface.
* width:   [m]    x width of the crystal.
* length:  [m]    z depth (length) of the crystal.
* R0: [1]         Reflectivity of mirror.
* coating: [str] Datafile containing either mirror material constants or reflectivity  numbers.
*
* %E
*******************************************************************************/

DEFINE COMPONENT Mirror_elliptic
DEFINITION PARAMETERS ()
SETTING PARAMETERS (x_a=0, y_b=1.0, z_c=1.0, 
        length=0.05, width=0.02,   
        R0=1,string coating="")
OUTPUT PARAMETERS (a2inv,b2inv,c2inv)
DEPENDENCY "-std=c99"
/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */

SHARE
%{
  %include "perfect_crystals-lib"
  #include <complex.h>
  %include "read_table-lib" 
  %include "reflectivity-lib"  
%}

DECLARE
%{
  double a2inv,b2inv,c2inv; /* 1/r^2 for physical ellipse */
  t_Reflec re;
%}

INITIALIZE
%{
  int status;

  if (coating && strlen(coating)) {
    status=reflec_Init(&re,COATING_UNDEFINED,coating,NULL);
  }else{
    /*assume a constant reflectivity*/
    status=reflec_Init(&re,CONSTANT,NULL, &(R0));
  }

    a2inv=(x_a)?1/(x_a*x_a):0; /* 0 really means infinity for x direction */
    b2inv=(y_b)?1/(y_b*y_b):0;
    c2inv=(z_c)?1/(z_c*z_c):0;

    printf("Component %s: initialized\n",NAME_CURRENT_COMP);

    fflush(NULL); // put diagnostics in order!
%}

TRACE
%{
    double E;				// (keV) x-ray energy
    double K; 				// length of k-vector
    double kxu,kyu,kzu;			// unit vector in the direction of k-vector.
    double x_int,y_int,z_int;
    double R;
    
    /* get the photon's kvector and energy */
    K=sqrt(kx*kx+ky*ky+kz*kz);
    E = K2E*K; /* use built-in constants for consistency */
    /* make unit vector in the direction of k :*/
    kxu = kx; kyu = ky; kzu = kz;
    NORM(kxu,kyu,kzu);

    double A,B,C, xt, yt, zt;
    double t0,t1;
    /*an offset to the mirror parameters perhaps*/

    xt=x;
    zt=z;
    yt=y-y_b;

    C=xt*xt*a2inv + yt*yt*b2inv + zt*zt*c2inv -1;
    B=2*(kxu*xt*a2inv + kyu*yt*b2inv + kzu*zt*c2inv);
    A=kxu*kxu*a2inv + kyu*kyu*b2inv + kzu*kzu*c2inv;

    if(solve_2nd_order(&t0,&t1,A,B,C)){
        double xx0, xx1, yy0, yy1, zz0, zz1; /* we will have to tentatively propagate twice to see which surface we hit */
        xx0=x+kxu*t0; yy0=y+kyu*t0; zz0=z+kzu*t0;
        xx1=x+kxu*t1; yy1=y+kyu*t1; zz1=z+kzu*t1;
        
        /*Check if we hit the mirror and whether the hit it is in front of the ray.
         * This does not account for mirror curvature */
                              
        int hit0=(fabs(xx0)<width/2.0) && (fabs(zz0)<length/2.0) && t0>0;
        int hit1=(fabs(xx1)<width/2.0) && (fabs(zz1)<length/2.0) && t1>0;
        int doit=hit0 || hit1;
                
        if(hit0 && !hit1) PROP_DL(t0); /* only one intersection actually on mirror */
        else if (hit1 && !hit0) PROP_DL(t1); /* other intersection */
        else if (hit0 && hit1) { /* both, take first strike (which may be back of mirror) */
            PROP_DL(t0<t1?t0:t1);	
        } else {
            RESTORE_XRAY(INDEX_CURRENT_COMP, x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p);
        }

        if ( doit ){
            SCATTER;
            xt=x; yt=y-y_b; zt=z; /* update shifted coordinates to intersection point */
            double nx, ny, nz;
            /* grad is the inner normal to the surface at this point */
            nx=-xt*a2inv; ny=-yt*b2inv; nz=-zt*c2inv; NORM(nx, ny, nz);

            double kdotn;
            kdotn=-scalar_prod(kx,ky,kz,nx,ny,nz);
            double kfx,kfy,kfz;
            kfx=kx+2*kdotn*nx;
            kfy=ky+2*kdotn*ny;
            kfz=kz+2*kdotn*nz;

            kx=kfx;ky=kfy;kz=kfz;

            /*adjust weight of ray*/

            double q=2.0*kdotn;
            //double R=reflecq(re,q,0,K*K2E,fabs(90-acos(kdotn/K)*RAD2DEG));
            double R=reflecq(re,q,0,K,fabs(90-acos(kdotn/K)*RAD2DEG));
            p*=R;
            /*update phase - as an approximation turn by 180 deg.*/;
            phi+=M_PI; 

            /*catch dead rays*/

            //if (p==0) ABSORB;
        } else {
        RESTORE_XRAY(INDEX_CURRENT_COMP, x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p);
        }
    }   
%}

MCDISPLAY
%{
  int i,j,N=12;
  double t,xx,yy,zz,xx0,yy0,zz0;
  double xwidth;
  double zdepth;
  zdepth=length;
  xwidth=width;  
  double a2,b2,c2; 
  a2=1.0/a2inv;
  b2=1.0/b2inv;
  c2=1.0/c2inv; 
  if (xwidth && zdepth){
    for (i=0;i<N;i++){
      xx0=-xwidth/2.0;
      zz=i*zdepth/(N-1)-zdepth/2.0;
      yy0=-(1)*sqrt(b2*(1-xx0*xx0/a2-zz*zz/c2))+sqrt(b2);
      for (j=1;j<N;j++){
        xx=j*xwidth/(N-1)-xwidth/2.0;
        yy=-(1)*sqrt(b2*(1-xx*xx/a2-zz*zz/c2))+sqrt(b2);
        line(xx0,yy0,zz,xx,yy,zz);
        xx0=xx;
        yy0=yy; 
      }
    }
    for (i=0;i<N;i++){
      zz0=-zdepth/2.0;
      xx=i*xwidth/(N-1)-xwidth/2.0;
      yy0=-(1)*sqrt(b2*(1-xx*xx/a2-zz0*zz0/c2))+sqrt(b2);
      for (j=1;j<N;j++){
        zz=j*zdepth/(N-1)-zdepth/2.0;
        yy=-(1)*sqrt(b2*(1-xx*xx/a2-zz*zz/c2))+sqrt(b2);
        line(xx,yy0,zz0,xx,yy,zz);
        zz0=zz;
        yy0=yy; 
      }
    }
  }    
%}

END