/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2011, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component Spherical_Backscattering_Analyser
*
* %I
* Written by: Nikolaos Tsapatsaris (with help from Peter Willendrup, Ruep Lechner, Heloisa Bordallo)
* Date: 2015
* Origin: Niels Bohr Institute
* 
* %D
* Spherical backscattering analyser with variable reflectivity, and mosaic gaussian response.
*
* Based on earlier developments by Miguel A. Gonzalez 2000, Tilo Seydel 2007, Henrik Jacobsen 2011.
*
* Example:
* COMPONENT doppler = Spherical_Backscattering_Analyser(
*   xmin=0, xmax=0, ymin=0, ymax=0, mosaic=21.0, dspread=0.00035, Q=2.003886241, DM=0, radius=2.5, f_doppler=0, A_doppler=0, R0=1, debug=0)
*
* %P
* xmin:          [m] Minimum absolute value of x =~ 0.5 * Width of the monochromator
* xmax:          [m] Maximum absolute value of x =~ 0.5 * Width of the monochromator
* ymin:          [m] Minimum absolute value of x =~ 0.5 * Width of the monochromator
* ymax:          [m] Maximum absolute value of y =~ 0.5 * Height of the monochromator
* radius:        [m] Curvature radius 
* dspread:       [1] Relative d-sprea, FWHM
* mosaic:  [arc min] Mosaicity, FWHM
* Q:          [AA-1] Magnitude of scattering vector 
* R0:            [1] Scalar, maximum reflectivity of neutrons
* DM:           [AA] monochromator d-spacing instead of Q = 2*pi/DM
* f_doppler:    [Hz] Frequency of doppler drive
* A_doppler:     [m] Amplitude of doppler drive
* debug:         [1] if debug > 0 print out some info about the calculations
*
*******************************************************************************/

DEFINE COMPONENT Spherical_Backscattering_Analyser

SETTING PARAMETERS (xmin=0, xmax=0, ymin=0, ymax=0, mosaic=21.0, dspread=0.00035, Q=2.003886241, DM=0, radius=2.5, f_doppler=0, A_doppler=0, R0=1, debug=0)


SHARE
%{
  double z_sphere(double xin, double yin, double rin) {
    return -(rin - sqrt(rin*rin - xin*xin - yin*yin));
  }

%}

DECLARE
%{
  double d_rms;
  double mos_rms;
  double mono_Q;
%}

INITIALIZE
%{
  mos_rms = MIN2RAD*mosaic/sqrt(8*log(2));
  
  mono_Q = Q;
  if (DM != 0) 
    mono_Q = 2*PI/DM;
  
  DM = 2*PI/mono_Q;
  d_rms = dspread*DM/sqrt(8*log(2));
%}

TRACE
%{  
  double vel;
    double sinTheta, lambdaBragg, lambda, dLambda2, sigmaLambda2;

  double old_x, old_y, old_z, old_t, x0, y0, z0, xi, yi, zi;
  double a, b, c, dt1, dt2, dt;
  double nx, ny, nz, nmod, v;
  double kproj, ydarwin, dkz, p_reflect; 
  double q0mod, q0x, q0y, q0z, kix, kiy, kiz, kfx, kfy, kfz;
        
  double omega_doppler;
  double v_doppler_inst;
    
  omega_doppler=2*PI*f_doppler;
   old_x=x; old_y=y; old_z=z; old_t=t;
    
  // Time interval necessary for the neutron to reach the sphere: solve equation: neutron path (line) crosses monochromator (sphere). Center of sphere is at (0,0,-radius).

//point on neutron path satisfies: xpath= x+vx*t, ypath=y+vy*t, zpath=z+vz*t;
//point on monochromator sphere satisfies: radius^2=x^2+y^2+(z+r)^2
//settings these equal gives an equation for t:
  a = vx*vx + vy*vy + vz*vz;
  b = 2.0 * ( vx*x + vy*y + vz*(z+radius) );  
  c = x*x+y*y+z*z + 2*z*radius;
  if ((b*b-4*a*c) < 0)
     {
      printf("Imaginary solutions. Something has gone wrong. Unphysical.\n");
      ABSORB;
     } 
  dt1 = (-b + sqrt(b*b-4*a*c)) / (2.0*a);
  dt2 = (-b - sqrt(b*b-4*a*c)) / (2.0*a);
  if (dt1 > 0)
     dt = dt1;
  else
     dt = dt2;   
  // propagates the neutron the time dt, so it arrives to the sphere

  PROP_DT(dt);
  
// ABSORB if neutron hits outside sphere
  if (x<xmin || x>xmax || y<ymin || y>ymax)
    ABSORB;

//n is a vector pointing along the radius of the sphere
  nx =x;
  ny= y;
  nz=(z+radius);
//modulus of the vector  n
nmod = sqrt(nx*nx+ny*ny+nz*nz);

// change to moving coordinates (doppler motion parallel to Z) 
  v_doppler_inst=A_doppler*omega_doppler*cos(omega_doppler*t);
  kix = vx*V2K; kiy = vy*V2K ; kiz = vz*V2K + v_doppler_inst*V2K;

// projection of the incident vector on the normal 
   kproj = (kix*nx + kiy*ny + kiz*nz) / nmod;
  
  vel=sqrt(a);  
  
  sinTheta = fabs(K2V*kproj)/vel;  
    
    // calculate lambdaBragg
    lambdaBragg = 2.0*DM*sinTheta;
    
    // calculate lambda of neutron
    lambda = 2*PI/kproj;
        
    // calculate deltaLambda squared and sigmaLambda squared
    dLambda2 = (lambda-lambdaBragg)*(lambda-lambdaBragg);
    // The sigmaLambda is propagated by differentiating the bragg 
    // condition: Lambda = 2*d*sinTheta
  
    sigmaLambda2 = 2.0*2.0 * sinTheta*sinTheta * d_rms*d_rms+2.0*2.0 * DM*DM * (1.0-sinTheta*sinTheta) * mos_rms*mos_rms;
      
    p_reflect = R0*exp(-dLambda2/(2.0*sigmaLambda2));

      if (p_reflect < 1e-5)
        {
          ABSORB;
        }
        else
        {
    // reflection: kf = ki - Q0 (the projection): q points along the normal and to scatter elastically, the component of ki along the normal must change sign, i.e. q0mod=2kproj
    q0mod = 2.0 * kproj;
    q0x = (nx/nmod)*q0mod; q0y = (ny/nmod)*q0mod;  q0z = (nz/nmod)*q0mod;
    kfx = kix - q0x;
    kfy = kiy - q0y;
    kfz = kiz - q0z;
    
      /* change to static coordinates */
      kfz = kfz - v_doppler_inst*V2K;
    
      vx = K2V*kfx;       
    vy = K2V*kfy;
      vz = K2V*kfz;

    p *= p_reflect; 
    SCATTER;
  }
  
    if(debug > 0) {
      printf("\n Lambda: %f, Lambda_Bragg: %f\n", lambda, lambdaBragg);
      printf("sigmaLambda: %f, R0: %f, p_reflect: %f\n", 
       sqrt(sigmaLambda2), R0, p_reflect);}
%}

MCDISPLAY
%{

  /*  dashed_line(xmin,ymin,0,xmin,ymax,0,5);
  dashed_line(xmin,ymax,0,xmax,ymax,0,5);
  dashed_line(xmax,ymax,0,xmax,ymin,0,5);
  dashed_line(xmax,ymin,0,xmin,ymin,0,5);*/

/* Step across the sphere in 11x11 steps to show the analyzer geometry */
  int i,j,N;
  double yv1,yv2,zv1,zv2;
  double xh1,xh2,zh1,zh2;
  double xd1,xd2,yd1,yd2,xd3,xd4,yd3,yd4,zd1,zd2,zd3,zd4;
  double dx,dy,dd;
  N=11;
  dx = (xmax-xmin)/(N-1);
  dy = (ymax-ymin)/(N-1);
  /* Horizontal, vertical, diagonal lines... */
  xh1=xmin;
  yv1=ymin;
  xd1=xmin;
  yd1=ymin;
  xd3=xmin;
  yd3=ymax;
  for (i=0; i<N-1; i++) {
    /* Calculate z-coordinates of 1st line-point(s)*/
    zh1=z_sphere(xh1,ymin,radius);
    zv1=z_sphere(xmin,yv1,radius);
    zd1=z_sphere(xd1,yd1,radius);
    zd3=z_sphere(xd3,yd3,radius);
    /* 2nd point x-y values */
    xh2=xh1+dx;
    yv2=yv1+dy;

    xd2=xd1+dx;
    yd2=yd1+dy;

    xd4=xd3+dx;
    yd4=yd3-dy;
    /* Calculate z-coordinates of 2nd set of line-point(s)*/
    zh2=z_sphere(xh1,ymin,radius);
    zv2=z_sphere(xmin,yv2,radius);
    zd2=z_sphere(xd2,yd2,radius);
    zd4=z_sphere(xd4,yd4,radius);
    /* Horizontal: */
    line(xh1,ymin,zh1,xh2,ymin,zh2);
    line(xh1,ymax,zh1,xh2,ymax,zh2);
    /* Vertical: */
    line(xmin,yv1,zv1,xmin,yv2,zv2);
    line(xmax,yv1,zv1,xmax,yv2,zv2);
    /* Diagonal: */
    line(xd1,yd1,zd1,xd2,yd2,zd2);
    line(xd3,yd3,zd3,xd4,yd4,zd4);
    /* Shift to next set of points ... */
    xh1=xh2; yv1=yv2; xd1=xd2; yd1=yd2; xd3=xd4; yd3=yd4;
  }
  
%}

END