/****************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2003, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Pol_mirror
*
* %I
* Written by: Peter Christiansen
* Date: July 2006
* Origin: RISOE
*
* Polarising mirror.
*
* %D
* This component models a rectangular infinitely thin mirror.
* For an unrotated component, the mirror surface lies in the Y-Z
* plane (ie. parallel to the beam).
* It relies on similar physics as the Monochromator_pol.
* The reflectivity function (see e.g. share/ref-lib for examples) 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)
* IMPORTANT: At present the component only works correctly for polarization along the up/down
* direction and for completely unpolarized beams, i.e. sx=sy=sz=0 for all rays.
*
* For now we assume:
* P(Transmit|Q) = 1 - P(Reflect|Q)
* i.e. NO ABSORPTION!
*
*
* The component can both reflect and transmit neutrons with a respective proportion
* depending on the p_reflect parameter:
*   p_reflect=-1 Reflect and transmit (proportions given from reflectivity) [default]
*   p_reflect=1 Only handle reflected events
*   p_reflect=0 Only handle transmitted events (reduce weight)
*   p_reflect=0-1 Both transmit and reflect with fixed statistics proportions
*
* 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.
*
* Examples:
* Reflection function parametrization
* Pol_mirror(zwidth = 0.40, yheight = 0.40, rUpFunc=StdReflecFunc, rUpPar={1.0, 0.0219, 6.07, 2.0, 0.003})
*
* Table function
* Pol_mirror(zwidth = 0.40, yheight = 0.40, rUpFunc=TableReflecFunc, rUpPar="supermirror_m2.rfl", rDownFunc=TableReflecFunc, rDownPar="supermirror_m3.rfl", useTables=1)
*
* See also the example instrument Test_Pol_Mirror (under tests).
*
* GRAVITY: YES
*
* %BUGS
* NO ABSORPTION
*
* %P
* INPUT PARAMETERS:
*
* zwidth: [m]     Width of the mirror
* yheight: [m]    Height of the mirror
* rUpFunc: [1]    Reflection function for spin up (q, *par, *r)
* rUpPar: [1]     Parameters for rUpFunc
* rDownFunc: [1]  Reflection function for spin down (q, *par, *r)
* rDownPar: [1]   Parameters for rDownFunc
* rUpData: [str]   Mirror Reflectivity data file for spin up
* rDownData: [str] Mirror Reflectivity data file for spin down
* useTables: [1]  Parameters are 0: Values, 1: Table names
* p_reflect: [1]  Proportion of reflected events. Use 0 to only get the transmitted beam, and 1 to get only the reflected beam. Use -1 to use the mirror reflectivity. This value is purely computational and is not related to the actual reflectivity
*
* CALCULATED PARAMETERS:
* SCATTERED: []   is 1 for reflected, and 2 for transmitted neutrons
*
* %L
*
* %E
*******************************************************************************/

DEFINE COMPONENT Pol_mirror


SETTING PARAMETERS (
  vector rUpPar={0.99,0.0219,6.07,2.0,0.003},
  vector rDownPar={0.99,0.0219,6.07,2.0,0.003},
  string rUpData="", string rDownData="",
  p_reflect=-1,
  zwidth,
  yheight)

/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */

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

DECLARE
%{
  t_Table rUpTable;
  int rUpTableFlag;
  t_Table rDownTable;
  int rDownTableFlag;
%}

INITIALIZE
%{

  if (strlen (rUpData) && strcmp (rUpData, "NULL")) {
    if (Table_Read (&rUpTable, rUpData, 1) <= 0) {
      fprintf (stderr, "Pol_mirror: %s: can not read file %s\n", NAME_CURRENT_COMP, rUpData);
      exit (1);
    }
    rUpTableFlag = 1;
  } else {
    rUpTableFlag = 0;
  }
  if (strlen (rUpData) && strcmp (rUpData, "NULL")) {
    if (Table_Read (&rDownTable, rDownData, 1) <= 0) {
      fprintf (stderr, "Pol_mirror: %s: can not read file %s\n", NAME_CURRENT_COMP, rDownData);
      exit (1);
    }
    rDownTableFlag = 1;
  } else {
    rDownTableFlag = 0;
  }

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

TRACE
%{
  double Q, Rup, Rdown, FN, FM, refWeight;
  int reflect = -1;
  int isPolarising = 0;

  // propagate to mirror plane
  PROP_X0;

  if (inside_rectangle (z, y, zwidth, yheight)) { /* Intersect the crystal? */

    // calculate scattering vector magnitude
    Q = fabs (2 * vx * V2K);

    // calculate reflection probability

    // downgraded from defpars version
    // rUpFunc(Q, rUpParPtr, &Rup);
    if (rUpTableFlag) {
      Rup = Table_Value (rUpTable, Q, 1);
    } else {
      StdReflecFunc (Q, rUpPar, &Rup);
    }
    // downgraded from defpars version
    // rDownFunc(Q, rDownParPtr, &Rdown);
    if (rDownTableFlag) {
      Rdown = Table_Value (rDownTable, Q, 1);
    } else {
      StdReflecFunc (Q, rDownPar, &Rdown);
    }

    if (Rup != Rdown) {
      isPolarising = 1;
      GetMonoPolFNFM (Rup, Rdown, &FN, &FM);
      GetMonoPolRefProb (FN, FM, sy, &refWeight);
      /* Output of PW discussions with Hal Lee 2024/03/08
       We have now done our QM "measurement", thus
       forcing the spin to assume up/down: */
      sx = 0;
      sz = 0;
    } else
      refWeight = Rup;

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

    // find out if neutrons is reflected or transmitted
    if (p_reflect < 0 || p_reflect > 1) { // reflect OR transmit using mirror reflectivity

      if (rand01 () < refWeight) // reflect
        reflect = 1;
      else
        reflect = 0;

    } else { // reflect OR transmit using a specified weighting

      if (rand01 () < p_reflect) {
        reflect = 1; // reflect
        p *= refWeight / p_reflect;
      } else {
        reflect = 0; // transmit
        p *= (1 - refWeight) / (1 - p_reflect);
      }
    }

    // set outgoing velocity and polarisation
    if (reflect == 1) { // reflect: SCATTERED==1 for reflection

      vx = -vx;
      if (isPolarising)
        SetMonoPolRefOut (FN, FM, refWeight, &sx, &sy, &sz);

    } else { // transmit: SCATTERED==2 for transmission
      if (isPolarising)
        SetMonoPolTransOut (FN, FM, refWeight, &sx, &sy, &sz);
      SCATTER;
    }

    if (isPolarising)
      if (sx * sx + sy * sy + sz * sz > 1.000001)
        fprintf (stderr, "Pol_mirror: %s: Warning: polarisation |s| = %g > 1\n", NAME_CURRENT_COMP,
                 sx * sx + sy * sy + sz * sz); // check that polarisation is meaningfull

    SCATTER;

  } /* End intersect the mirror */
  else {
    /* neutron will miss the mirror, so don't propagate i.e restore it */
    RESTORE_NEUTRON (INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);
  }
%}

MCDISPLAY
%{

  rectangle ("yz", 0, 0, 0, zwidth, yheight);
%}

END