/*******************************************************************************
*
* Component: Conics_HE
*
* %Identification
* Written by: Peter Wilendrup and Erik Knudsen <br>(derived from Giacomo Resta skeleton-component)
* Origin: DTU
* Release: McStas 2.7
* Date: September 2021
*
* Hyperbolic-Elliptic shells for Wolter optic
*
* %Description
* Implements a set of nshells Wolter Ellipsoid/Hyperboloid pairs using conics.h from ConicTracer.
* 
* The component has two distinct modes of specifying the geometry:
* a) Via the radii vector, parametrized from largest to smallest radius with a length of 'nshells'
*
* b) By specifying radii rmax and rmin, between which a quadratic law distributes 'nshells' surfaces.
*
* The mirrors are assumed to be touching at the mid-optic plane, i.e. there is no gap between primary 
* and secondary mirror.
* By definition the ratio between primary and secondary mirror glancing angles is 1/3.
* At present a single m-value is used for all mirrors.  
*
* %Parameters
*
* Input parameters:
* nshells: [1]          Number of nested shells to expect
* rmin: [m]             Midoptic plane radius of innermost mirror pair.
* rmax: [m]             Midoptic plane radius of outermost mirror pair.
* radii: [m]            Optional vector of radii (length should match nshells)
* m: [1]                Critical angle of mirrors as multiples of Ni_c.
* focal_length_u: [m]   Focal length (upstream) of the mirror pairs.
* focal_length_d: [m]   Focal length (downstream) of the mirror pairs.
* le: [m]               Paraboloid mirror length.
* lh: [m]               Hyperboloid mirror length.
* disk: [ ]             Flag. If nonzero, insert a disk to block the central area within the innermost mirror.
* Qc: [AA^-1]            Critical scattering vector
* R0: [1]               Reflectivity at Low angles for reflectivity curve approximation
* alpha: [AA]           Slope of reflectivity for reflectivity curve approximation
* W: [AA^-1]             Width of supermirror cut-off
* transmit: [1]         Fraction of statistics to assign to transmitted beam - NOT YET IMPLEMENTED
* sigma_m: [rad]        RMS meridional slope (figure) error of the mirror surfaces; 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 of the mirror surfaces
* mirr_thick: [m]       Thickness of mirror shell surfaces - NOT YET IMPLEMENTED
* %End
*
*******************************************************************************/
     
DEFINE COMPONENT Conics_HE

SETTING PARAMETERS (
    rmin=0.0031416, rmax=0.05236,focal_length_u=10, focal_length_d=10, le=0.25, lh=0.25,
    int nshells=4, m=1, mirr_thick=0, int disk=1, vector radii=NULL,
    R0 = 0.99, Qc = 0.021, W = 0.003, alpha = 6.07, transmit = 0, sigma_m=1e-6, sigma_s=250e-6)

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

DECLARE
%{
  // Scene where all geometry is added to
  Scene s;
%}

INITIALIZE
%{
  ConicSurf* pm;
  double dr, rr, cH, theta_1, theta_2, theta_i;
  double constant, quadratic;
  int i;

  s = makeScene ();
  if (!radii) {
    /* Mode b, use quadratic law to distribute the shells. */
    dr = nshells > 1 ? (rmax - rmin) / (nshells - 1) : 0;

    quadratic = (rmax - rmin) / (rmax * rmax - rmin * rmin);
    constant = rmax - quadratic * rmax * rmax;
  }

  for (i = 0; i < nshells; i++) {
    /* Mode a uses the supplied radii; mode b uses the quadratic distribution. */
    if (radii) {
      rr = radii[i];
    } else {
      rr = rmax - dr * i;
      rr = constant + quadratic * rr * rr; // Quadratic distribution of radius covers angular space better
    }

    Point pi = makePoint (0, rr, 0);
    if (focal_length_u == focal_length_d) {
      /* The common focus is at infinity in this limit, so both surfaces are parabolic. */
      pm = addParaboloid (-focal_length_u, pi, -lh, 0, m, R0, Qc, alpha, W, &s);
      addParaboloid (focal_length_d, pi, 0, le, m, R0, Qc, alpha, W, &s);
    } else {
      theta_1 = atan (rr / focal_length_u);
      theta_2 = atan (rr / focal_length_d);
      theta_i = 0.25 * (theta_1 + theta_2);
      cH = fabs (0.5 * (rr / tan (theta_2 - 2.0 * theta_i) - focal_length_d));

      pm = addHyperboloid (-focal_length_u - 2.0 * cH, -focal_length_u, pi, -lh, 0, m, R0, Qc, alpha, W, &s);
      addEllipsoid (-focal_length_u - 2.0 * cH, focal_length_d, pi, 0, le, m, R0, Qc, alpha, W, &s);
    }
  }
  for (i = 0; i < s.num_c; i++) {
    s.c[i].figerr_m = sigma_m;
    s.c[i].figerr_s = sigma_s;
  }
  if (disk) {
    addDisk (pm->zs, 0.0, rConic (pm->ze, *pm), &s);
  }
%}                            

TRACE
%{
  /* "_mctmp_a" defines a "silicon" state variable in underlying conic.h functions */
  _mctmp_a = 0;
  traceSingleNeutron (_particle, &s);

  if (!_particle->_absorbed) {
    SCATTER;
  }
%}

FINALLY %{
    //Mainly Writes Inline Detector Data
    finishSimulation(&s);
%}

MCDISPLAY
%{
  double zz = 0;
  // Enlarge xy-plane when mcdisplay is ran with --zoom
  magnify ("xy");

  // Draw Conic Surfaces (full revolution about z, r(z) from conic.h) as polygon facets
  int i;
  for (i = 0; i < s.num_c; i++) {
    int nz = 10, nphi = 24, iz, j;
    double dz = (s.c[i].ze - s.c[i].zs) / nz;
    for (iz = 0; iz < nz; iz++) {
      double z0 = s.c[i].zs + iz * dz, z1 = z0 + dz;
      double r0 = rConic (z0, s.c[i]), r1 = rConic (z1, s.c[i]);
      for (j = 0; j < nphi; j++) {
        double p0 = 2.0 * PI * j / nphi, p1 = 2.0 * PI * (j + 1) / nphi;
        polygon (4, r0 * cos (p0), r0 * sin (p0), z0 - zz, r0 * cos (p1), r0 * sin (p1), z0 - zz, r1 * cos (p1), r1 * sin (p1), z1 - zz, r1 * cos (p0),
                 r1 * sin (p0), z1 - zz);
      }
    }
  }

  // Draw xy-axis cross hairs for Disks
  // Local variables to control maximal display-size of cross-hairs
  for (i = 0; i < s.num_di; i++) {
    double r0 = s.di[i].r0;
    double r1 = s.di[i].r1;
    double z0 = s.di[i].z0;
    if (r0 > 1.0)
      r0 = 1.0;
    if (r1 > 1.0)
      r1 = 1.0;
    line (r0, 0, z0 - zz, r1, 0, z0 - zz);
    line (-r0, 0, z0 - zz, -r1, 0, z0 - zz);
    line (0, r0, z0 - zz, 0, r1, z0 - zz);
    line (0, -r0, z0 - zz, 0, -r1, z0 - zz);
  }
%}
	
END