/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2011, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: curved_mirror
*
* %I
* Written by: Henrik Jacobsen
* Date: April 2012
* Origin: RNBI
*
* Single mirror plate that is curved and fits into an elliptic guide. 
*
* %D
*
*
* %P
* INPUT PARAMETERS:
*
* focus_start_w: [m]           Horizontal position of first focal point of ellipse in ABSOLUTE COORDINATES
* focus_start_h: [m]           Vertical position of first focal point of ellipse in ABSOLUTE COORDINATES
* focus_end_w: [m]             Horizontal position of second focal point of ellipse in ABSOLUTE COORDINATES
* focus_end_h: [m]             VErtical position of second focal point of ellipse in ABSOLUTE COORDINATES
* mirror_start: [m]            start of mirror in ABSOLUTE COORDINATES
* guide_start: [m]             start of guide in ABSOLUTE COORDINATES
* smallaxis_w: [m]             Small-axis dimension of horizontal ellipse
* smallaxis_h: [m]             Small-axis dimension of vertical ellipse
* yheight: [m]                 the height of the mirror when not in the guide
* smallaxis: [m]               the smallaxis of the guide
* length: [m]                  the length of the mirror
* R0=0.99:		[]	low angle reflectivity
* Qc: [AA-1]                   critical scattering vector
* alpha: [AA]                  slope of reflectivity
* m:			[]	m-value of material
* W: [AA-1]                    width of supermirror cutoff
* transmit: [0/1]              0: non reflected neutrons are absorbed. 1: non reflected neutrons pass through
* substrate_thickness: [m]     thickness of substrate (for absorption)
* coating_thickness: [m]       thickness of coating (for absorption)
* theta_1: [deg]               angle of mirror wrt propagation direction at start of mirror
* theta_2: [deg]               angle of mirror wrt propagation direction at center of mirror
* theta_3: [deg]               angle of mirror wrt propagation direction at end of mirror
* reflect: [q(Angs-1) R(0-1)]  (str)  Name of relfectivity file. Format 
*
*
*
* %D
*
* %E
*******************************************************************************/

DEFINE COMPONENT Mirror_Elliptic_Bispectral

SETTING PARAMETERS (string reflect=0,focus_start_w,focus_end_w, focus_start_h, focus_end_h, mirror_start, m,  smallaxis_w, smallaxis_h, length, transmit=0, substrate_thickness=0.0005, coating_thickness=10e-6)


SHARE
%{
%include "read_table-lib"
%}

DECLARE
%{
  t_Table pTable;
%}

INITIALIZE
%{
  if (reflect && strlen(reflect) && strcmp(reflect,"NULL") && strcmp(reflect,"0")) {
    if (Table_Read(&pTable, reflect, 1) <= 0) /* read 1st block data from file into pTable */
      exit(fprintf(stderr,"Mirror: %s: can not read file %s\n", NAME_CURRENT_COMP, reflect));
  }
%}

TRACE
%{
double f; //half of distance between focal points
double asquared;
double a; //half of ellipse length
double b; //half of ellipse width

double xprime; //x in coordinates with center of ellipse at (xprime,zprime)=(0,0)
double ymirror; //height of the mirror


//Defining the mirror
double a1;
double b1;
double c1;

//solving the time the neutron hits the sample
double A, B, C, D, E, P, Q, R, U, V, I, J, K;

//finding rotation of mirror
double alpha1, beta1, gamma1;
double theta_m;
double xhi, zeta;


double b_w;
double f_w;
double asquared_w;
double A_w;
double B_w;
double C_w;
double determinant_w;


double b_h;
double f_h;
double asquared_h;
double xprime_h;


double xprime_w;
double zprime_w;
double asquare_z;
double bsquare_x;


double v_n; //speed of neutron perpendicular to surface

double Ref; //reflectivity

double dt;
double q;
  int intersect;

double discriminant;

double xprime_start_w;
double zprime_start_w;

double z_test;
double x_test;
double z_prime_test;

int hit_back_flag;

int prop_case;


double x_2;
double y_2;
double z_2;
double t_2;
int x_hit;
int x_hit_2;

double xprime_h_2;
double ymirror_2;
int intersect_2;

intersect=0;
x_2=0;
y_2=0;
z_2=0;
t_2=-1;
prop_case=0;
   double old_x = x, old_y = y, old_z = z, old_t=t, old_vx=vx, old_vz=vz, old_vy=vy;

// Check if neutron hits mirror. First find which z,x coordinates it hits.



//define the ellipse
b_w=smallaxis_w/2;
f_w=(focus_end_w-focus_start_w)*0.5;
asquared_w=f_w*f_w+b_w*b_w;

//in coordinate system of mirror: xprime_w(t)=xprime_w+vx*t. z-value is zprime(t)=b*sqrt(1-x'^2/a^2)+old_z+vz*t. This gives equation for t: A_w*t^2+B_w*t+C_w=0;

xprime_start_w=old_x-f_w-focus_start_w+mirror_start;
zprime_start_w=z;

A_w=b_w*b_w*vx*vx+asquared_w*vz*vz;
B_w=2*b_w*b_w*xprime_start_w*vx+2*asquared_w*old_z*vz;
C_w=b_w*b_w*xprime_start_w*xprime_start_w+asquared_w*old_z*old_z-asquared_w*b_w*b_w;

//this equation must now be solved for t;

if (A_w!=0){
determinant_w=B_w*B_w-4.0*A_w*C_w;
if (determinant_w>=0){
I=(-B_w-sqrt(determinant_w))/(2.0*A_w);
J=(-B_w+sqrt(determinant_w))/(2.0*A_w);

if (I<=0 && J<=0){dt = -1.0;} //both times are negative.
if (I<=0 && J>0 ){dt = J;} //set dt to only positive value.
if (I>0  && J<=0){dt = I;} //set dt to only positive value.

if (I>0  && J>0 ){prop_case=1; if (I>J) {dt=J;}else{dt=I;}} //set dt to smallest positive value  


} else {dt=-1.0;} //only complex solutions: set dt negative so no scattering

}else{ //end if (A!=0)  
if (B_w!=0){
dt=-C/B_w;
}else{ //end if (B!)=0
 printf("warning: A_w=B_w=C_w=0. Neutron is ignored\n"); }
} //end if (A!=0){}else{
//now intersection time has been found.

//printf("dt=%f\n",dt);

if (dt>0) { //if time is positive, propagate neutron to where it hits mirror. This is done without gravity.

    x += vx*dt;
    y += vy*dt;
    z += vz*dt;
    t += dt;


if (prop_case>0) //also check if neutron can hit mirror at second solution - it might not be in y-range for first solution
{
    x_2=x+vx*fabs(J-I); 
    y_2=y+vy*fabs(J-I);
    z_2=z+vz*fabs(J-I); 
    t_2=t+fabs(J-I); 

}else{
x_2=x;
y_2=y;
z_2=z;
t_2=t;
}

x_hit=(x>=0 &&x<=length);
x_hit_2=(x_2>=0 &&x_2<=length);

// printf("x=%f,y=%f,z=%f\n",x,y,z);
//if (x >=0 && x<=length){ //check if neutron is within x limits of the mirror. If so, check if it is within y limits.
if (x_hit || x_hit_2){

//define the ellipse
b_h=smallaxis_h/2;

f_h=(focus_end_h-focus_start_h)*0.5;

 asquared_h=f_h*f_h+b_h*b_h;

xprime_h=-f_h-focus_start_h+mirror_start+x; //xprime is the x-coordinate in a coordinate system centered at the center of the ellipse

ymirror=b_h*sqrt(1-xprime_h*xprime_h/asquared_h);


xprime_h_2=-f_h-focus_start_h+mirror_start+x_2; //xprime is the x-coordinate in a coordinate system centered at the center of the ellipse

ymirror_2=b_h*sqrt(1-xprime_h_2*xprime_h_2/asquared_h);


intersect = ( y>=-ymirror && y<=ymirror && x>=0 && x<=length && zprime_start_w+vz*dt>=0);

intersect_2 = ( y_2>=-ymirror_2 && y_2<=ymirror_2 && x_2>=0 && x_2<=length && zprime_start_w+vz*(dt+fabs(J-I))>=0);

if (!intersect && intersect_2){ //if neutron doesn't hit mirror with smallest t, but hits with largest t, propagte to largest t
intersect=intersect_2;
x=x_2;
y=y_2;
z=z_2;
t=t_2;
dt=t_2-old_t;
//printf("x=%f,y=%f,z=%f\n",x,y,z);
}

    if (intersect) { //if neutron is within ylimits of the mirror handle reflection/transmission


//now perform reflection. 

//First find out if neutron hits front or back of mirror: propagate backwards a bit and check if neutron is outside ellipse: if so, it hits back of mirror


b_w=smallaxis_w/2.0;

f_w=(focus_end_w-focus_start_w)*0.5;
 asquared_w=f_w*f_w+b_w*b_w;

z_test=zprime_start_w+vz*(dt-1e-6);
x_test=xprime_start_w+vx*(dt-1e-6);
z_prime_test=b_w*sqrt(1-x_test*x_test/asquared_w);

//find velocity in q direction.

xprime_w=xprime_start_w+vx*dt;
zprime_w=zprime_start_w+vz*dt;
asquare_z=asquared_w*zprime_w;
bsquare_x=b_w*b_w*xprime_w;


zeta=(asquare_z)/(sqrt(asquare_z*asquare_z+bsquare_x*bsquare_x));
xhi=-(bsquare_x)/(sqrt(asquare_z*asquare_z+bsquare_x*bsquare_x));

//printf("z_test=%f, z_prime_test=%f\n",z_test,z_prime_test);

if (z_test>z_prime_test) {
hit_back_flag=1;
}
//printf("vx=%f, vz=%f, vy=%f, xhi=%f, zeta=%f, prop_case=%d\n",vx,vz,vy,xhi,zeta,prop_case);
v_n=-xhi*vx+zeta*vz;

q=fabs(2.0*v_n*V2Q);


 //Reflectivity parameters calculated from SWISS neutronics data.
double R0=0.99;
double Qc=0.0217;
double m_value=m*0.9853+0.1978;
double W=-0.0002*m_value+0.0022;
double alpha=0.1204*m_value+5.0944;
double beta=-7.6251*m_value+68.1137;

if (m_value<=3)
{alpha=m_value;
beta=0;}


      /* Reflectivity (see component Guide). */
      if(m == 0)
        ABSORB;
      if (reflect && strlen(reflect) && strcmp(reflect,"NULL") && strcmp(reflect,"0"))
         Ref=Table_Value(pTable, q, 1);
      else {
          Ref = R0;
          if(q > Qc)
          {
            double arg = (q-m_value*Qc)/W;
            if(arg < 10)
              Ref *= .5*(1-tanh(arg))*(1-alpha*(q-Qc)+beta*(q-Qc)*(q-Qc)); //matches data from Swiss Neutronics
            else  Ref=0;
          }
      }
      if (Ref < 0) Ref=0;
      else if (Ref > 1) Ref=1;


//Now comes actual reflection/transmission
      if (!transmit) { //all neutrons are reflected
//printf("v_n=%f,q=%f, Ref=%f, lambda=%f, theta=%f, 1p_before=%f",v_n,q, Ref, (2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz),asin((2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz)*q/(4*PI))*RAD2DEG,p);
        if (!Ref) ABSORB;
        p *= Ref;
//printf("p_after=%f\n",p);
//handle reflection: change v_n -->-v_n

vx=old_vx*(zeta*zeta-xhi*xhi)+old_vz*(2*zeta*xhi);
vz=+old_vx*(2*zeta*xhi)+old_vz*(xhi*xhi-zeta*zeta);


        SCATTER; //after transmission or reflection

      } else { //if neutrons can be transmitted



//calculate absorption.
// substrate
double lambda=(2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz);
double sin_theta=lambda*q/(4*PI);
double substrate_path_length=substrate_thickness/sin_theta;
double mu_substrate=0.0318/lambda+0.0055*lambda-0.0050; //unit: cm^-1
mu_substrate=mu_substrate*100; //unit: m^-1;
               
//For nickel and titanium coating, the following formular is used:
// mu = rho/m(atom)*sigma_a,thermal*lambda/lambda_thermal

// lambda_thermal=1.798 Å

// rho_nickel=8.908g/cm^3
// m(atom)_nickel=58.6934*1.661*10^-27 kg
// sigma_a,thermal_nickel=4.49*10^-28 m^2

// rho_titanium=4.506g/cm^3
// m(atom)_titanium=47.867*1.661*10^-27 kg
// sigma_a,thermal_titanium=6.09*10^-28 m^2

double coating_path_length=coating_thickness/sin_theta;
double Ni_coefficient=22.8180; 
double Ti_coefficient=19.1961;

double mu_coating=(0.5*Ni_coefficient+0.5*Ti_coefficient)*lambda; //it is roughly 50% nickel and 50% titanium



        // transmit when rand > R
        if (Ref == 0 || rand01() >= Ref) { //transmit
if (substrate_thickness>0){ p=p*exp(-mu_substrate*substrate_path_length-mu_coating*coating_path_length);} //reduce weight of neutrons due to attenuation in the mirror

} else {//neutron is reflected
		if (hit_back_flag==1 && substrate_thickness>0) { //if neutron comes from behind the mirror
			p=p*exp(-2*mu_substrate*substrate_path_length-2*mu_coating*coating_path_length);} //reduce weight of neutrons due to attenuation in the mirror
//handle reflection: change v_n -->-v_n

vx=old_vx*(zeta*zeta-xhi*xhi)-old_vz*(2*zeta*xhi);
vz=-old_vx*(2*zeta*xhi)+old_vz*(xhi*xhi-zeta*zeta);
}

//printf("p_before=%f, q=%f",p,q);
        SCATTER; //after transmission or reflection
//printf("p_after=%f\n",p);
      } //end } else { after if (!transmit) {
    } 


 


} // end if (x >=-length/2 && x<=length/2)


   if (!intersect) {
      /* No intersection: restore neutron position. */
      x = old_x;
      y = old_y;

      z = old_z;
      t = old_t;


    }
  

} //end if (dt>0)


%}

MCDISPLAY
%{




double xstart;
double xend;
double xprime_start;
double ystart;

double xprime_end;
double yend;

double focus_2_h;
double focus_1_h;
double b_h;
double f_h;
double asquared_h;

double focus_2_w;
double focus_1_w;
double b_w;
double f_w;
double asquared_w;


int n_lines;
int j=0;
double xprimepos[51];
double ypos[51];
double x_plot[51];
double zpos[51];
double xstep;
double xprime_w[51];



focus_2_h=focus_end_h; //focus in local coordinates
focus_1_h=focus_start_h;

b_h=smallaxis_h/2.0;

f_h=(focus_2_h-focus_1_h)*0.5;
 asquared_h=f_h*f_h+b_h*b_h;



focus_2_w=focus_end_w; //focus in local coordinates
focus_1_w=focus_start_w;

b_w=smallaxis_w/2.0;

f_w=(focus_2_w-focus_1_w)*0.5;
 asquared_w=f_w*f_w+b_w*b_w;



xstart=0;
xend=length;

n_lines=50;

xstep=length/n_lines;













for (j=0; j<n_lines+1; j++)
{
xprimepos[j]=-f_h-focus_start_h+mirror_start+xstart+xstep*j;

ypos[j]=b_h*sqrt(1-xprimepos[j]*xprimepos[j]/asquared_h); //correct


x_plot[j]=xstart+xstep*j;



xprime_w[j]=-f_w-focus_start_w+mirror_start+xstart+xstep*j; //xprime is the x-coordinate in a coordinate system centered at the center of the ellipse

zpos[j]=b_w*sqrt(1-xprime_w[j]*xprime_w[j]/asquared_w);
//printf("xprime=%f, zpos=%f\n",xprime_w[j], zpos[j]);
}

for (j=0; j<n_lines; j++)
{
line(x_plot[j], -ypos[j], zpos[j], x_plot[j+1], -ypos[j+1],zpos[j+1]);
line(x_plot[j], ypos[j], zpos[j], x_plot[j+1], ypos[j+1],zpos[j+1]);
}


line(x_plot[0],-ypos[0],zpos[0],x_plot[0],ypos[0],zpos[0]);
line(x_plot[50],-ypos[50],zpos[50],x_plot[50],ypos[50],zpos[50]);




//printf("ypos0=%f xpos0=%f ypos50=%f xpos50=%f",ypos[0], x_plot[0], ypos[50], x_plot[50]);

/*  double xmax, xmin, ymax, ymin;
  

  if (center == 0) {
    xmax= x1; xmin=0;
    ymax= yheight; ymin=0;
  } else {
    xmax= x1/2; xmin=-xmax;
    ymax= yheight/2; ymin=-ymax;
  }
  multiline(5, (double)xmin, (double)ymin, 0.0,
               (double)xmax, (double)ymin, 0.0,
               (double)xmax, (double)ymax, 0.0,
               (double)xmin, (double)ymax, 0.0,
               (double)xmin, (double)ymin, 0.0);
*/
%}
END