/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2026, All rights reserved
*         DTU Physics, Kongens Lyngby, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Mirror_conic
*
* %I
* Written by: Peter Willendrup
* Date: October 2026
* Origin: DTU
*
* Single flat mirror plate using the conic.h reflection routines (test component)
*
* %D
* Test component for the reflection routines in conic.h (ConicTracer), which are
* used by e.g. Conics_*, FlatEllipse_finite_mirror and NMO. It models a single
* flat mirror so that the results can be compared directly with Mirror.comp.
*
* The mirror lies in the local y-z plane at x=0, extending from z=0 to z=length
* and from y=-yheight/2 to y=yheight/2. Only neutrons arriving from the x<0 side
* (vx>0) interact; everything else passes unaffected. Each neutron interacts at
* most once.
*
* conic=0: the mirror is a conic.h FlatSurf and is handled by reflectNeutronFlat(),
*          which reflects with probability R and otherwise transmits the neutron
*          unaffected. Compare with Mirror(transmit=1).
* conic=1: the mirror is a conic.h ConicSurf cylinder of (very large) radius 1 km,
*          handled by reflectNeutronConic(), which always reflects and weights by R.
*          Compare with Mirror(transmit=0). A Gaussian figure (slope) error can be
*          applied via sigma_m and sigma_s; the default is 0 here, whereas the
*          Conics_* components default to the conic.h values of 1 and 250 urad.
*
* Example: Mirror_conic(length=1, yheight=0.1, m=2)
*
* %P
* length: [m]      Length of the mirror along z
* yheight: [m]     Height of the mirror along y
* R0: [1]          Low-angle reflectivity
* Qc: [AA^-1]      Critical scattering vector
* alpha: [AA]      Slope of reflectivity
* m: [1]           m-value of material
* W: [AA^-1]       Width of supermirror cut-off
* conic: [1]       0: use FlatSurf/reflectNeutronFlat(), 1: use ConicSurf/reflectNeutronConic()
* sigma_m: [rad]   RMS meridional slope (figure) error, conic=1 only; for grazing angles below
*                  10*sigma_m the waviness model of Hansen et al., arXiv:1612.03763 is used
* sigma_s: [rad]   RMS sagittal slope (figure) error, conic=1 only
*
* %E
*******************************************************************************/

DEFINE COMPONENT Mirror_conic

SETTING PARAMETERS (length=1, yheight=0.1, R0=0.99, Qc=0.021, alpha=6.07, m=2, W=0.003, int conic=0, sigma_m=0, sigma_s=0)

SHARE
%{
  %include "ref-lib"
  %include "conic.h"

  /* Set up the flat and the cylindrical mirror surface |x| = r, z in [0, len] */
  void mirror_conic_init (FlatSurf* f, ConicSurf* c, double r, double len, double yh, double* par) {
    memset (f, 0, sizeof (*f));
    f->k1 = r * r;
    f->ze = len;
    f->ll = -yh / 2;
    f->rl = yh / 2;
    f->R0 = par[0];
    f->Qc = par[1];
    f->alpha = par[2];
    f->m = par[3];
    f->W = par[4];

    memset (c, 0, sizeof (*c));
    c->k1 = r * r;
    c->ze = len;
    c->R0 = par[0];
    c->Qc = par[1];
    c->alpha = par[2];
    c->m = par[3];
    c->W = par[4];
    c->figerr_m = par[5];
    c->figerr_s = par[6];
    c->type = PARA;
  }
%}

DECLARE
%{
  FlatSurf fs;
  ConicSurf cs;
  double xshift;
%}

INITIALIZE
%{
  /* Both surfaces are |x| = xshift (FlatSurf) or r = xshift (ConicSurf) in conic.h
     coordinates; the local mirror plane x=0 is mapped to x=xshift. */
  xshift = 1e3;

  double par[7] = { R0, Qc, alpha, m, W, sigma_m, sigma_s };
  mirror_conic_init (&fs, &cs, xshift, length, yheight, par);
%}

TRACE
%{
  double dt;

  if (vx > 0) {
    /* "_mctmp_a" is the "silicon" state in conic.h; 0 means no substrate tracking */
    _mctmp_a = 0;
    x += xshift;
    dt = conic ? getTimeOfFirstCollisionConic (*_particle, cs) : getTimeOfFirstCollisionFlat (*_particle, fs);
    x -= xshift;
    if (dt > 0) {
      PROP_DT (dt);
      if (fabs (y) <= yheight / 2) {
        SCATTER;
        x += xshift;
        if (conic)
          reflectNeutronConic (_particle, cs);
        else
          reflectNeutronFlat (_particle, fs);
        x -= xshift;
        if (_particle->_absorbed)
          ABSORB;
      }
    }
  }
%}

MCDISPLAY
%{
  polygon (4, 0.0, -yheight / 2, 0.0, 0.0, -yheight / 2, length, 0.0, yheight / 2, length, 0.0, yheight / 2, 0.0);
%}

END