/****************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2003, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Pol_guide_vmirror
*
* %I
* Written by: Peter Christiansen and Erik B Knudsen
* Date: July 2006
* Origin: RISOE
*
* Polarising guide with a supermirror along its diagonal.
*
* %D
* Models a rectangular guide with entrance centered on the Z axis and
* with one supermirros sitting on the diagonal inside.
* The entrance lies in the X-Y plane.  Draws a true depiction
* of the guide with mirror and trajectories.
* The polarisation is handled similar to in Monochromator_pol.
* The reflec functions are handled similar to Pol_mirror.
* The up direction is hardcoded to be along the y-axis (0, 1, 0)
*
* Note that this component can also be used as a frame overlap-mirror
* if the up and down reflectivities are set equal. In this case the wall
* refletivity (rPar) should probably be set to absorb. 
*
* The parameters can either be
* double pointer initializations (e.g. {R0, Qc, alpha, m, W})
* or table names (e.g."supermirror_m2.rfl" AND useTables=1).
* NB! This might cause warnings by the compiler that can be ignored.
*
* GRAVITY: YES
*
* %BUGS
* No absorption by mirror.
*
* %P
* INPUT PARAMETERS:
*
* xwidth: [m]      Width at the guide entry
* yheight: [m]     Height at the guide entry
* length: [m]      length of guide
* rFunc: [1]       Guide Reflection function
* rPar: [1]        Guide Parameters for rFunc
* rUpFunc: [1]     Mirror Reflection function for spin up
* rUpPar: [1]      Mirror Parameters for rUpFunc
* rDownFunc: [1]   Mirror Reflection function for spin down
* rDownPar: [1]    Mirror Parameters for rDownFunc
* useTables: [1]   Parameters are 0: Values, 1: Table names
* debug: [1]       if debug > 0 print out some internal runtime parameters
*
* OUTPUT PARAMETERS:
*
* localG: [m/s/s]  Gravity vector in guide reference system 
* normalTop: [1]   One of several normal vectors used for defining the geometry 
* pointTop: [1]    One of several points used for defining the geometry 
* rParPtr:      One of several pointers to reflection parameters used with the ref. functions. []
* SCATTERED: []    is 1 for reflected, and 2 for transmitted neutrons
*
*
* %L
*
* %E
*******************************************************************************/

DEFINE COMPONENT Pol_guide_mirror
DEFINITION PARAMETERS (xwidth, yheight, length,
rFunc=StdReflecFunc, rUpFunc=StdReflecFunc,
rDownFunc=StdReflecFunc,
rPar    ={1.0, 0.0219, 4.07, 3.2, 0.003},
rUpPar  ={1.0, 0.0219, 4.07, 3.2, 0.003},
rDownPar={0.1, 0.0219, 4.07, 3.2, 0.003},
useTables=0)
SETTING PARAMETERS (int debug=0)
OUTPUT PARAMETERS (localG,
normalTop, normalBot, normalLeft, normalRight, normalInOut,
pointTop, pointBot, pointLeft, pointRight,
pointIn, pointOut,
rParPtr, rUpParPtr, rDownParPtr)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */

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

DECLARE
%{
Coords localG;
Coords normalTop, normalBot, normalLeft, normalRight, normalInOut;
Coords pointTop, pointBot, pointLeft, pointRight, pointIn, pointOut;

#if (useTables)
  t_Table *rParPtr     = 0;
  t_Table *rUpParPtr   = 0;
  t_Table *rDownParPtr = 0;
#else
  double rParPtr[]     = (double [])rPar;
  double rUpParPtr[]   = (double [])rUpPar;
  double rDownParPtr[] = (double [])rDownPar;
#endif
  %}

INITIALIZE
%{
#if (useTables)
  rParPtr     = (t_Table*) malloc(sizeof(t_Table));
  rUpParPtr   = (t_Table*) malloc(sizeof(t_Table));
  rDownParPtr = (t_Table*) malloc(sizeof(t_Table));
  if (Table_Read(rParPtr, rPar, 1) <= 0) {
    fprintf(stderr,"Pol_guide_vmirror: %s: can not read file %s\n",
	    NAME_CURRENT_COMP, rPar);
    exit(1);
  }
  if (Table_Read(rUpParPtr, rUpPar, 1) <= 0) {
    fprintf(stderr,"Pol_guide_vmirror: %s: can not read file %s\n",
	    NAME_CURRENT_COMP, rUpPar);
    exit(1);
  }
  if (Table_Read(rDownParPtr, rDownPar, 1) <= 0) {
    fprintf(stderr,"Pol_guide_vmirror: %s: can not read file %s\n",
	    NAME_CURRENT_COMP, rDownPar);
    exit(1);
  }
#endif

  if ((xwidth<=0) || (yheight<= 0) || (length<=0)) {
    fprintf(stderr, "Pol_guide_vmirror: %s: NULL or negative length scale!\n"
	                  "ERROR      (xwidth,yheight,length). Exiting\n",
	    NAME_CURRENT_COMP);
    exit(1);
  }

  if (mcgravitation) {

    localG = rot_apply(ROT_A_CURRENT_COMP, coords_set(0,-GRAVITY,0));
    fprintf(stdout,"Pol_guide_vmirror: %s: Gravity is on!\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
  // The entrance and exit plane are defined by the normal vector
  // (0, 0, 1)
  // and the two points pointIn=(0, 0, 0) and pointOut=(0, 0, length)

  normalInOut = coords_set(0, 0, 1);
  pointIn   = coords_set(0, 0, 0);
  pointOut  = coords_set(0, 0, length);

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

  // Bottom plane (-y dir) can be spanned by (1, 0, 0) & (0, 0, 1)
  // and the point (0, -yheight/2, 0)
  // A normal vector is (0, 1, 0)
  normalBot  = coords_set(0, 1, 0);
  pointBot = coords_set(0, -yheight/2, 0);

  // Left plane (+x dir) can be spanned by (0, 1, 0) & (0, 0, 1)
  // and the point (xwidth/2, 0, 0)
  // A normal vector is (1, 0, 0)
  normalLeft  = coords_set(1, 0, 0);
  pointLeft = coords_set(xwidth/2, 0, 0);

  // Right plane (-x dir) can be spanned by (0, 1, 0) & (0, 0, 1)
  // and the point (-xwidth/2, 0, 0)
  // A normal vector is (1, 0, 0)
  normalRight  = coords_set(1, 0, 0);
  pointRight = coords_set(-xwidth/2, 0, 0);

  %}

TRACE
%{
  /* time threshold */
  const double tThreshold = 1e-10/sqrt(vx*vx + vy*vy + vz*vz);
  const double xwhalf = xwidth/2;
  const double norm = 1.0/sqrt(xwidth*xwidth + length*length);
  double R;

  Coords normalMirror, pointMirror;
  Coords* normalPointer = 0;

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

  /* Propagate neutron to guide entrance. */
  PROP_Z0;

  if (!inside_rectangle(x, y, xwidth, yheight))
    ABSORB;

  normalMirror  = coords_set(norm*length, 0, -norm*xwidth);
  pointMirror = coords_set(-xwhalf, 0, 0);

  for(;;) {
    double tLeft, tRight, tTop, tBot, tIn, tOut, tMirror;
    double tUp, tSide, time, endtime;
    double Q; //, dummy1, dummy2, dummy3;
    Coords vVec, xVec;
    int mirrorReflect;

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

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

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

    solve_2nd_order(&tRight, NULL, 0.5*coords_sp(normalRight,localG),
		    coords_sp(normalRight, vVec),
		    coords_sp(normalRight, coords_sub(xVec, pointRight)));

    solve_2nd_order(&tLeft, NULL, 0.5*coords_sp(normalLeft,localG),
		    coords_sp(normalLeft, vVec),
		    coords_sp(normalLeft, coords_sub(xVec, pointLeft)));

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

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

    solve_2nd_order(&tMirror, NULL, 0.5*coords_sp(normalMirror,localG),
		    coords_sp(normalMirror, vVec),
		    coords_sp(normalMirror, coords_sub(xVec, pointMirror)));

    /* 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 (tMirror>tThreshold && tMirror<time) {

      time=tMirror;
      mirrorReflect = 1; // flag to show which reflection function to use
    }

    if (time<=tThreshold)
      fprintf(stdout, "Pol_guide_vmirror: %s: tTop: %f, tBot:%f, tRight: %f, tLeft: %f\n"
	      "tUp: %f, tSide: %f, time: %f\n", NAME_CURRENT_COMP,
	      tTop, tBot, tRight, tLeft, tUp, tSide, time);

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

    if (time > endtime)
      break;

    if(time <= tThreshold) {

      printf("Time below threshold!\n");
      fprintf(stdout, "Pol_guide_vmirror: %s: tTop: %f, tBot:%f, tRight: %f, tLeft: %f\n"
	      "tUp: %f, tSide: %f, time: %f\n", NAME_CURRENT_COMP,
	      tTop, tBot, tRight, tLeft, tUp, tSide, time);
      break;
    }

    if(debug>0 && time==tLeft) {

      fprintf(stdout, "\nPol_guide_vmirror: %s: Left side hit: x, v, normal, point, gravity\n", NAME_CURRENT_COMP);
      coords_print(xVec);
      coords_print(vVec);
      coords_print(normalLeft);
      coords_print(pointLeft);
      coords_print(localG);

      fprintf(stdout, "\nA: %f, B: %f, C: %f, tLeft: %f\n",
	      0.5*coords_sp(normalLeft,localG),coords_sp(normalLeft, vVec),
	      coords_sp(normalLeft, coords_sub(xVec, pointLeft)), tLeft);
    }

    if(debug>0)
      fprintf(stdout, "Pol_guide_vmirror: %s: tTop: %f, tBot:%f, tRight: %f, tLeft: %f\n"
	      "tUp: %f, tSide: %f, time: %f\n", NAME_CURRENT_COMP,
	      tTop, tBot, tRight, tLeft, tUp, tSide, time);

    if(debug>0)
      fprintf(stdout, "Pol_guide_vmirror: %s: Start v: (%f, %f, %f)\n", NAME_CURRENT_COMP, vx, vy, vz);

    PROP_DT(time);
    if (mcgravitation)
      vVec = coords_set(vx, vy, vz);
    SCATTER;

    if(time==tTop)
      normalPointer = &normalTop;
    else if(time==tBot)
      normalPointer = &normalBot;
    else if(time==tRight)
      normalPointer = &normalRight;
    else if(time==tLeft)
      normalPointer = &normalLeft;
    else if(time==tMirror)
      normalPointer = &normalMirror;
    else
      fprintf(stderr, "Pol_guide_vmirror: %s: This should never happen!!!!\n", NAME_CURRENT_COMP);

    Q = 2*coords_sp(vVec, *normalPointer)*V2K;

    if(!mirrorReflect) {
      // we have hit one of the sides. Always reflect.
      vVec = coords_add(vVec, coords_scale(*normalPointer, -Q*K2V));
      rFunc(fabs(Q), rParPtr, &refWeight);
      p *= refWeight;

    } else {
      // we have hit one of the mirrors
      rUpFunc(fabs(Q), rUpParPtr, &Rup);
      rDownFunc(fabs(Q), rDownParPtr, &Rdown);
      if (Rup <  0)   ABSORB;
      if (Rup >  1)   Rup =1 ;
      if (Rdown <  0) ABSORB;
      if (Rdown >  1) Rdown =1 ;

      GetMonoPolFNFM(Rup, Rdown, &FN, &FM);
      GetMonoPolRefProb(FN, FM, sy, &refWeight);
      // check that refWeight is meaningfull
      if (refWeight <  0) ABSORB;
      if (refWeight >  1) refWeight =1 ;

      if (rand01()<refWeight) {
	      //reflect: SCATTERED==1 for reflection

	      vVec = coords_add(vVec, coords_scale(*normalPointer, -Q*K2V));
	      SetMonoPolRefOut(FN, FM, refWeight, &sx, &sy, &sz);
      } else {
	      // transmit: SCATTERED==2 for transmission
	      SCATTER;
	      SetMonoPolTransOut(FN, FM, refWeight, &sx, &sy, &sz);
      }

      if (sx*sx+sy*sy+sz*sz>1.000001) { // check that polarisation is meaningfull
        fprintf(stderr, "Pol_guide_vmirror: %s: polarisation |s|=%g > 1 s=[%g,%g,%g]\n",
          NAME_CURRENT_COMP, sx*sx+sy*sy+sz*sz, sx, sy, sz);
      }
    }

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

    // set new velocity vector
    coords_get(vVec, &vx, &vy, &vz);

    if(debug>0)
      fprintf(stdout, "Pol_guide_vmirror: %s: End v: (%f, %f, %f)\n", NAME_CURRENT_COMP, vx, vy, vz);
  }

  %}

MCDISPLAY
%{
  int i, j;

  

  // draw box
  box(0, 0, length/2.0, xwidth, yheight, length);

  for(j = -1; j<=1; j+=2)
      line(-xwidth/2.0, j*yheight/2, 0, xwidth/2.0, j*yheight/2, length);

  %}

END