/****************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2003, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Pol_bender
*
* %I
* Written by: Peter Christiansen
* Date: August 2006
* Origin: RISOE
*
* Polarising bender.
*
* %D
* Based on Guide_curved written by Ross Stewart.
* Models a rectangular curved guide tube with entrance centered on the Z axis.
* The entrance lies in the X-Y plane.  Draws a true depiction
* of the guide with multiple slits (but without spacers), and trajectories.
* It relies on similar physics as the Monochromator_pol.
* The reflec function and parameters are passed to this component to
* give a bigger freedom.
* The up direction is hardcoded to be along the y-axis (0, 1, 0)
*
* The guide is asummed to have half a spacer on each side:
*    slit1  slit2  slit3
*  |+     ++     ++     +|
*  |+     ++     ++     +|
*  <---------------------> xwidth
*  <------> xwidth/nslit (nslit=3)
*         <> d
*
* The reflection functions and parameters defaults as follows:
* Bot defaults to Top, Left defaults to Top, Right defaults to left.
* Down defaults to down and up defaults to up for all functions and
* Top(Up and Down) defaults to StdReflecFunc and {0.99,0.0219,6.07,2.0,0.003}
* which stands for {R0, Qc, alpha, m, W}.
*
* Example:
* Pol_bender(xwidth = 0.08, yheight = 0.08, length = 1.0, radius= 10.0,
* 	     nslit=5, d=0.0, endFlat=0, drawOption=2,
* 	     rTopUpPar={0.99, 0.0219, 6.07, 3.0, 0.003},
* 	     rTopDownPar={0.99, 0.0219, 6.07, 1.0, 0.003})
*
* See also the example instruments Test_Pol_Bender and
* Test_Pol_Bender_Vs_Guide_Curved (under tests).
*
* %BUGS
* This component has been against tested Guide_curved and found to
* give the same intensities. Gravity option has not been tested.
*
* GRAVITY: YES (when gravity is along y-axis)
*
* %P
* INPUT PARAMETERS:
*
* xwidth: [m]        Width at the guide entry 
* yheight: [m]       Height at the guide entry 
* length: [m]        length of guide along center 
* radius: [m]        Radius of curvature of the guide (+:curve left/-:right) 
* G: [m/s^2]         Gravitational constant
* nslit: [1]         Number of slits 
* d: [m]             Width of spacers (subdividing absorbing walls) 
* endFlat: [1]       If endflat>0 then entrance and exit planes are parallel. 
* rTopUpPar: [1]     Top mirror Parameters for spin up standard reflectivity function
* rTopDownPar: [1]   Top mirror Parameters for spin down standard reflectivity function
* rBotUpPar: [1]     Bottom mirror Parameters for spin up standard reflectivity function
* rBotDownPar: [1]   Bottom mirror Parameters for spin down standard reflectivity function
* rLeftUpPar: [1]    Left mirror Parameters for spin up standard reflectivity function
* rLeftDownPar: [1]  Left mirror Parameters for spin down standard reflectivity function
* rRightUpPar: [1]   Right mirror Parameters for spin up standard reflectivity function
* rRightDownPar: [1] Right mirror Parameters for spin down standard reflectivity function
* rTopUpData: [1]     Reflectivity file for top mirror, spin up 
* rTopDownData: [1]   Reflectivity file for top mirror, spin down 
* rBotUpData: [1]     Reflectivity file for bottom mirror, spin up
* rBotDownData: [1]   Reflectivity file for bottom mirror, spin down 
* rLeftUpData: [1]    Reflectivity file for left mirror, spin up
* rLeftDownData: [1]  Reflectivity file for left mirror, spin down 
* rRightUpData: [1]   Reflectivity file for right mirror, spin up
* rRightDownData: [1] Reflectivity file for right mirror, spin down 
* drawOption: [1]    1: fine(all slits/90 points per arc), 2: normal (max 20/40), 3: rough (max 5/10) 
* debug: [1]         if debug > 0 print out some internal parameters 
*
* CALCULATED PARAMETERS:
*
* localG: [m/s/s]    Gravity vector in guide reference system 
* normalXXX: [1]     Several normal vector used for defining the geometry 
* pointXXX: [1]      Several points used for defining the geometry 
* rXXXParPtr: []     Pointers to reflection parameters used with ref. functions.
*
* %L
*
* %E
*******************************************************************************/

DEFINE COMPONENT Pol_bender

SETTING PARAMETERS (xwidth, yheight, length, radius, G=9.8, int nslit=1, d=0.0, int debug=0, int endFlat=0,
vector rTopUpPar={0.99,0.0219,6.07,2.0,0.003},
vector rTopDownPar={0.99,0.0219,6.07,2.0,0.003},
vector rBotUpPar={0.99,0.0219,6.07,2.0,0.003},
vector rBotDownPar={0.99,0.0219,6.07,2.0,0.003},
vector rLeftUpPar={0.99,0.0219,6.07,2.0,0.003},
vector rLeftDownPar={0.99,0.0219,6.07,2.0,0.003},
vector rRightUpPar={0.99,0.0219,6.07,2.0,0.003},
vector rRightDownPar={0.99,0.0219,6.07,2.0,0.003},
string rTopUpData="", string rTopDownData="",string rBotUpData="",string rBotDownData="",
string rLeftUpData="", string rLeftDownData="",string rRightUpData="",string rRightDownData="",
int drawOption=1)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */

SHARE
%{
  %include "pol-lib"
  %include "ref-lib"
%}

DECLARE
%{
  Coords localG;
  Coords normTopBot;
  Coords normIn;
  Coords normOut;
  Coords pointTop;
  Coords pointBot;
  Coords pointIn;
  Coords pointOut;

  t_Table rTopUpTable;
  t_Table rTopDownTable;
  t_Table rBotUpTable;
  t_Table rBotDownTable;
  t_Table rLeftUpTable;
  t_Table rLeftDownTable;
  t_Table rRightUpTable;
  t_Table rRightDownTable;
  int useTables;
%}

INITIALIZE
%{
  double angle;

  if(strlen(rTopUpData) && strcmp(rTopUpData,"NULL")){
    useTables=1;
    /*if rUpTopData is set assume we'll be usning tabled data for all reflectivities*/
    if (Table_Read(&rTopUpTable, rTopUpData, 1) <= 0) {
      fprintf(stderr,"Pol_bender: %s: can not read file %s\n",NAME_CURRENT_COMP, rTopUpData);
      exit(1);
    }
    if (Table_Read(&rTopDownTable, rTopDownData, 1) <= 0) {
      fprintf(stderr,"Pol_bender: %s: can not read file %s\n",NAME_CURRENT_COMP, rTopDownData);
      exit(1);
    }
    if (Table_Read(&rBotUpTable, rBotUpData, 1) <= 0) {
      fprintf(stderr,"Pol_bender: %s: can not read file %s\n",NAME_CURRENT_COMP, rBotUpData);
      exit(1);
    }
    if (Table_Read(&rBotDownTable, rBotDownData, 1) <= 0) {
      fprintf(stderr,"Pol_bender: %s: can not read file %s\n",NAME_CURRENT_COMP, rBotDownData);
      exit(1);
    }
    if (Table_Read(&rLeftUpTable, rLeftUpData, 1) <= 0) {
      fprintf(stderr,"Pol_bender: %s: can not read file %s\n",NAME_CURRENT_COMP, rLeftUpData);
      exit(1);
    }
    if (Table_Read(&rLeftDownTable, rLeftDownData, 1) <= 0) {
      fprintf(stderr,"Pol_bender: %s: can not read file %s\n",NAME_CURRENT_COMP, rLeftDownData);
      exit(1);
    }
    if (Table_Read(&rRightUpTable, rRightUpData, 1) <= 0) {
      fprintf(stderr,"Pol_bender: %s: can not read file %s\n",NAME_CURRENT_COMP, rRightUpData);
      exit(1);
    }
    if (Table_Read(&rRightDownTable, rRightDownData, 1) <= 0) {
      fprintf(stderr,"Pol_bender: %s: can not read file %s\n",NAME_CURRENT_COMP, rRightDownData);
      exit(1);
    }
  }
  if ((xwidth<=0) || (yheight <= 0) || (length<=0) || (radius==0)) {
    fprintf(stderr, "Pol_bender: %s: NULL or negative length scale!\n"
	    "ERROR      (xwidth,yheight,length, radius). Exiting\n",
	    NAME_CURRENT_COMP);
    exit(1);
  }

  if (drawOption<1 || drawOption>3) {
    fprintf(stderr, "Pol_bender: %s: drawOption %ld not supported. Exiting.\n",
	    NAME_CURRENT_COMP, drawOption);
    exit(1);
  }

  if (mcgravitation) {

    localG = rot_apply(ROT_A_CURRENT_COMP, coords_set(0,-GRAVITY,0));
    fprintf(stdout,"Pol_bender %s: Gravity is on!\n",
	    NAME_CURRENT_COMP);
    if (localG.x!=0 || localG.z!=0)
      fprintf(stderr,"WARNING: Pol_Bender: %s: "
	      "This component only gives correct resulta with gravitation,\n"
	      "when gravity is strictly along the y-axis!\n",
	      NAME_CURRENT_COMP);

  } else
    localG = coords_set(0, 0, 0);

  // To be able to handle the situation properly where a component of
  // the gravity is along the z-axis we also define entrance (in) and
  // exit (out) planes

  angle = length/radius;
  normIn    = coords_set(0, 0, 1);
  if (endFlat)
    normOut   = coords_set(0, 0, 1);
  else
    normOut   = coords_set(sin(angle), 0, cos(angle));
  pointIn   = coords_set(0, 0, 0);
  pointOut  = coords_set(radius-radius*cos(angle), 0, radius*sin(angle));

  // Top and bot plane (+y dir) can be spanned by (1, 0, 0) & (0, 0, 1)
  // and the top point (0, yheight/2, 0) and bot point (0, -yheight/2, 0)
  // A normal vector is (0, 1, 0)
  normTopBot  = coords_set(0, 1, 0);
  pointTop = coords_set(0, yheight/2, 0);
  pointBot = coords_set(0, -yheight/2, 0);

%}

TRACE
%{
  const double whalf  = 0.5*xwidth; /* half width of guide */
  const double hhalf  = 0.5*yheight; /* half height of guide */
  const double z_off  = radius*sin(length/radius); /* z-comp of guide length */
  const double dThreshold = 1e-10; /* distance threshold */
  const double tThreshold = dThreshold/sqrt(vx*vx + vy*vy + vz*vz);
  double angle_z_vout;       /* angle between z-axis and v_out */

  // Variables used in the case of multiple slits
  const double slitWidth = xwidth/nslit; // slitwidth
  const double spacerhalf  = 0.5*d; /* half width of spacers */
  int slitHit;                 // decide which slit is hit
  double posInSlit;            // position in slit

  double t11, t12, t21, t22, theta, alpha, endtime, phi;
  int i_bounce;
  int nerr=0;

  // Pol variables
  double FN, FM, Rup, Rdown, weight;

  double Rleft;  /* radius of curvature of left mirror */
  double Rright; /* radius of curvature of right mirror */
  double absR = fabs(radius);
  double sign = 1;
  if(radius<0)
    sign = -1;

  /* Propagate neutron to guide entrance. */
  PROP_Z0;
  if (!inside_rectangle(x, y, xwidth, yheight))
    ABSORB;

  if(nslit>1) {
    // check if neutron is absorbed on a spacer
    posInSlit = fmod(x+whalf, slitWidth);
    if(posInSlit <= spacerhalf ||
       posInSlit >= slitWidth-spacerhalf)
      ABSORB;

    // check which slat is hit
    slitHit = (int)((x+whalf)/slitWidth);

    // Modify R1 and R2 according to which slat was hit
    Rleft  = absR + sign*whalf - sign*(slitHit+1)*slitWidth + sign*spacerhalf;
    Rright = absR + sign*whalf - sign*slitHit*slitWidth - sign*spacerhalf;

    if(debug>0)
      printf("\nslitHit: %d/%f, Rleft: %f, Rright: %f\n",
	     slitHit, (x+whalf)/slitWidth, Rleft, Rright);
  } else { // only 1 slit

    Rleft  = absR - sign*whalf;
    Rright = absR + sign*whalf;
  }

  for(;;) {

    double tLeft, tRight, tTop, tBot, tIn, tOut, tMirror;
    double tUp, tSide, time, endtime;
    double R, Q;
    Coords vVec, xVec;
    int isPolarising;
    double vel_xz;

    isPolarising = 0;

    xVec = coords_set(x, y, z);
    vVec = coords_set(vx, vy, vz);

    solve_2nd_order(&tTop, NULL, 0.5*coords_sp(normTopBot,localG),
		    coords_sp(normTopBot, vVec),
		    coords_sp(normTopBot, coords_sub(xVec, pointTop)));

    solve_2nd_order(&tBot, NULL, 0.5*coords_sp(normTopBot,localG),
		    coords_sp(normTopBot, vVec),
		    coords_sp(normTopBot, coords_sub(xVec, pointBot)));

    solve_2nd_order(&tIn, NULL, 0.5*coords_sp(normIn,localG),
		    coords_sp(normIn, vVec),
		    coords_sp(normIn, coords_sub(xVec, pointIn)));

    solve_2nd_order(&tOut, NULL, 0.5*coords_sp(normOut,localG),
		    coords_sp(normOut, vVec),
		    coords_sp(normOut, coords_sub(xVec, pointOut)));

    /* Find itersection points with inside and outside guide walls */
    if (!cylinder_intersect(&t11, &t12 ,x - radius, y, z, vx, vy, vz, Rleft,  2*yheight)){
      /*neutron did not hit the cylinder*/
      t11=t12=0;
    }
    if (!cylinder_intersect(&t21, &t22 ,x - radius, y, z, vx, vy, vz, Rright, 2*yheight)){
      /*neutron did not hit the cylinder*/
      t21=t22=0;
    }

    /* Choose appropriate reflection time */
    tLeft  = (t11 < tThreshold) ? t12 : t11;
    tRight = (t21 < tThreshold) ? t22 : t21;

    /* Choose appropriate reflection time */
    if (tTop>tThreshold && (tTop<tBot || tBot<=tThreshold))
      tUp=tTop;
    else
      tUp=tBot;

    if (tLeft>tThreshold && (tLeft<tRight || tRight<=tThreshold))
      tSide=tLeft;
    else
      tSide=tRight;

    if (tUp>tThreshold && (tUp<tSide || tSide<=tThreshold))
      time=tUp;
    else
      time=tSide;

    if (time<=tThreshold) {
      nerr++;
      if (nerr < 10) {
        fprintf(stdout, "tTop: %e, tBot:%e, tRight: %e, tLeft: %e\n"
	      "tUp: %e, tSide: %e, time: %e\n",
	      tTop, tBot, tRight, tLeft, tUp, tSide, time);
      } else {
        fprintf(stdout, "Found 10 propagation error for this neutron, terminating!\n");
        break;
      }
    }

    /* Has neutron left the guide? */
    if (tOut>tThreshold && (tOut<tIn || tIn<=tThreshold))
      endtime=tOut;
    else
      endtime=tIn;

    if (time > endtime)
      break;

    PROP_DT(time);
    SCATTER;

    /* Find reflection surface */
    if(time==tSide) {      /* Left or right side */

      if(time==tLeft)
	R = sign*Rleft;
      else
	R = sign*Rright;

      phi   = atan(vx/vz);        /* angle of neutron trajectory */
      alpha = asin(z/R);      /* angle of guide wall */
      theta = fabs(phi - alpha);    /* angle of reflection */
      angle_z_vout = 2.0*alpha - phi;

      vel_xz = sqrt(vx*vx + vz*vz);      /* in plane velocity */
      vz    = vel_xz*cos(angle_z_vout);
      vx    = vel_xz*sin(angle_z_vout);

    } else { /* Top or Bottom wall */
      theta = fabs(atan(vy/vz));
      vy    = -vy;
    }

    /* Now compute reflectivity. */
    Q = 2.0*sin(theta)*sqrt(vx*vx + vy*vy + vz*vz)*V2K;

    // calculate reflection probability
    if(time==tTop) {
      if(useTables){
        Rup=Table_Value(rTopUpTable,Q,1);
        Rdown=Table_Value(rTopDownTable,Q,1);
      }else{
        StdReflecFunc(Q, rTopUpPar, &Rup);
        StdReflecFunc(Q, rTopDownPar, &Rdown);
      }
      if(debug>0)
          fprintf(stdout, "\tTop hit:\n");
    } else if(time==tBot) {
      if(useTables){
        Rup=Table_Value(rBotUpTable,Q,1);
        Rdown=Table_Value(rBotDownTable,Q,1);
      }else{
        StdReflecFunc(Q, rBotUpPar, &Rup);
        StdReflecFunc(Q, rBotDownPar, &Rdown);
      }
      if(debug>0)
          fprintf(stdout, "\tBot hit:\n");
    } else if(time==tRight) {
      if(useTables){
        Rup=Table_Value(rRightUpTable,Q,1);
        Rdown=Table_Value(rRightDownTable,Q,1);
      }else{
        StdReflecFunc(Q, rRightUpPar, &Rup);
        StdReflecFunc(Q, rRightDownPar, &Rdown);
      }
      if(debug>0)
          fprintf(stdout, "\tRight hit:\n");
    } else if(time==tLeft) {
      if(useTables){
        Rup=Table_Value(rLeftUpTable,Q,1);
        Rdown=Table_Value(rLeftDownTable,Q,1);
      }else{
        StdReflecFunc(Q, rLeftUpPar, &Rup);
        StdReflecFunc(Q, rLeftDownPar, &Rdown);
      }
      if(debug>0)
          fprintf(stdout, "\tLeft hit:\n");
    }
    if(Rup != Rdown) {
      isPolarising = 1;
      GetMonoPolFNFM(Rup, Rdown, &FN, &FM);
      GetMonoPolRefProb(FN, FM, sy, &weight);
    } else
      weight = Rup;

    if(debug>0)
      printf("\tlambda: %.2f AA, Q: %.4f, Rup: %.4f, Rdown: %.4f,"
	     " weight: %.4f\n",
	     2*PI/(sqrt(vx*vx + vy*vy + vz*vz)*V2K), Q,
	     Rup, Rdown, weight);

    // check that refWeight is meaningfull
    if (weight <= 0) ABSORB;
    if (weight >  1) weight =1 ;

    if(isPolarising) {
      SetMonoPolRefOut(FN, FM, weight, &sx, &sy, &sz);
      if(sx*sx+sy*sy+sz*sz>1.000001)
        fprintf(stderr,"Pol_bender: %s: Warning: polarisation |s| = %g > 1\n",
              NAME_CURRENT_COMP, sx*sx+sy*sy+sz*sz); // check that polarisation is meaningfull
    }

    p *= weight;

    if(p==0) {
      ABSORB;
      break;
    }
  }

%}

MCDISPLAY
%{
  double x1, x2, z1, z2;
  const int n = 90;
  double *xplot=malloc(n*sizeof(double));
  double *zplot=malloc(n*sizeof(double));
  int ns = 0;
  int j = 1;
  const double lengthOfGuide = sin(length/radius)*radius;
  const double slitWidth = xwidth/nslit;
  double R = 0; /* radius of arc */
  int nSlitsMax = nslit;
  int nMax      = n;

  if (lengthOfGuide<=0)
    exit(fprintf(stdout,"Pol_bender: %s: Negative guide length ! lengthOfGuide=%g\n",
	    NAME_CURRENT_COMP, lengthOfGuide));

  if (drawOption==2) {

    if(nSlitsMax>20)
      nSlitsMax = 20;
    nMax = 40;
  } else if (drawOption==3) {

    if(nSlitsMax>5)
      nSlitsMax = 5;
    nMax = 10;
  }

  

  // draw opening
  rectangle("xy", 0, 0, 0, xwidth, yheight);

  for(ns=0; ns < nSlitsMax+1; ns++) {

    // to make sure the sides are drawn properly
    if(ns==nSlitsMax && nSlitsMax<nslit)
      ns=nslit;

    // calculate x for this R
    R = radius - 0.5*xwidth + ns*slitWidth;

    for(j=0; j<nMax; j++) {

      if(endFlat) {

	if(ns==0)  // only calculate once
	  zplot[j] = j*lengthOfGuide/(double)(nMax-1);
      } else
	zplot[j] = R*sin(length/radius * (double)j/(double)(nMax-1));

      if(radius>0)
	xplot[j] = radius - sqrt(R*R - zplot[j]*zplot[j]);
      else
	xplot[j] = radius + sqrt(R*R - zplot[j]*zplot[j]);
    }

    // To be able to draw end we store some of the point values
    if(ns==0) { // first wall

      x1 = xplot[nMax-1];
      z1 = zplot[nMax-1];
    } else if(ns==nslit) { //last wall

      x2 = xplot[nMax-1];
      z2 = zplot[nMax-1];
    }

    for(j=0; j<nMax-1; j++) {
      line(xplot[j],  0.5*yheight, zplot[j], xplot[j+1],  0.5*yheight, zplot[j+1]);
      line(xplot[j], -0.5*yheight, zplot[j], xplot[j+1], -0.5*yheight, zplot[j+1]);
    }
  }

  // draw end gap
  line(x1, 0.5*yheight, z1, x2, 0.5*yheight, z2);
  line(x1, 0.5*yheight, z1, x1,-0.5*yheight, z1);
  line(x2,-0.5*yheight, z2, x2, 0.5*yheight, z2);
  line(x1,-0.5*yheight, z1, x2,-0.5*yheight, z2);
  free(xplot);free(zplot);
  %}

END