/*******************************************************************************
*
* Component: FlatEllipse_finite_mirror
*
* %I
* Written by: Christoph Herb, TUM
* Version: 0.1
* Origin: TUM
* Date: 2022-2023
*
* NMO (nested mirror optic) modules 
*
* %D
* Simulates NMO (nested mirror optic) modules as concevied by Böni et al., see
* Christoph Herb et al., Nucl. Instrum. Meth. A 1040, 1671564 (1-18) 2022.
*
* The component relies on an updated version of conic.h from MIT.
*
* %P
* sourceDist: [m]   Distance used for calculating the spacing of the mirrors
* LStart:     [m]   Left focal point
* LEnd:       [m]   Right focal point
* lStart:     [m]   z-Value of the mirror start
* lEnd:       [m]   z-Value of the mirror end
* r_0: [m] distance to the mirror at lStart
* nummirror:  [1]   number of mirrors
* mf: [ ] mvalue of the inner side of the coating, m>10 results in perfect reflections
* mb: [ ] mvalue of the outer side of the coating, m>10 results in perfect reflections
* mirror_width: [m] width of the individual mirrors
* mirror_sidelength: [m] side length of the individual mirrors
* doubleReflections: [1] binary value determining whether the mirror backside is reflective
* rfront_inner_file: [str] file of distances to the optical axis of the individual mirrors
* R0: [1] Low-angle reflectivity
* Qc: [AA-1]                   critical scattering vector
* W: [AA-1]                    width of supermirror cutoff
* alpha: [AA] Slope of reflectivity
* %L
* Christoph Herb et al., Nucl. Instrum. Meth. A 1040, 1671564 (1-18) 2022.
* %E
*
*******************************************************************************/

DEFINE COMPONENT FlatEllipse_finite_mirror
SETTING PARAMETERS (
    sourceDist = 0,
    LStart=0.6,
    LEnd = 0.6,
    lStart = 0.0,
    lEnd = 0.0,
    r_0 = 0.02076,
    int nummirror= 9,
    mf = 4,
    mb = 0,
    R0 = 0.99,
    Qc = 0.021,
    W = 0.003,
    alpha = 6.07,
    mirror_width = 0.003,
    mirror_sidelength = 1,
    doubleReflections = 0,
    string rfront_inner_file = "NULL"
)


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

  /* Function originally defined in file "calciterativemirrors.h"

  /*! \brief Function to return an array of distances for a nested mirror assembly, see attached files. Also works reasonably well for parabolic mirrors
  see

  @param number number of entries in the array = number of mirrros
  @param z_0 z-coordinate of the initial point on the mirror
  @param r_0 r-coordinate of the initial point on the mirror
  @param z_extract z-coordinate at which the distances are extracted
  @param LStart z-coordinate of the left focal point
  @param LEnd z-coordinate of the right focal point
  @param lStart z-coordinate at which the mirrors begin
  @param lEnd z-coordinate at which the mirrros end
  @return pointer to array with number of distances
  */
  double*
  get_r_at_z0 (int number, double z_0, double r_0, double z_extract, double LStart, double LEnd, double lStart, double lEnd) {
    int n = number;
    double* r_zExtracts = malloc (n * sizeof (double_t)); /* n is an array of 10 integers */
    if (!r_zExtracts) {
      fprintf (stderr, "NMO comp function get_r_at_z0: malloc() failed. Exit! \n");
      exit (-1);
    }
    r_zExtracts[0] = r_0;
    // helper variables as in conic_finite_mirror.h and explained in swissneutronics_überlegungen
    double k1;
    double k2;
    double k3;
    double c;
    double u;
    double a;
    double r_lEnd;
    double r_lStart;
    // initial mirror is calculated from the initial point z0, r0
    c = (LEnd - LStart) / 2;
    u = (z_0 + c - LEnd);
    a = sqrt ((u * u + c * c + r_0 * r_0 + sqrt (pow (u * u + c * c + r_0 * r_0, 2) - 4 * c * c * u * u)) / 2);
    k3 = c * c / (a * a) - 1;
    k2 = 2 * k3 * (c - LEnd);
    k1 = k3 * (c - LEnd) * (c - LEnd) - c * c + a * a;
    printf ("k1 %f k2 %f k3 %f\n", k1, k2, k3);
    // next mirror will be calculated with the point on the surface being lStart, r_lStart
    for (int k = 0; k < number; ++k) {
      r_zExtracts[k] = sqrt (k1 + k2 * z_extract + k3 * z_extract * z_extract);
      r_lEnd = sqrt (k1 + k2 * lEnd + k3 * lEnd * lEnd);       // calculate the radius at the end
      r_lStart = r_lEnd * (lStart - LStart) / (lEnd - LStart); //

      c = (LEnd - LStart) / 2;
      u = (lStart + c - LEnd);
      a = sqrt ((u * u + c * c + r_lStart * r_lStart + sqrt (pow (u * u + c * c + r_lStart * r_lStart, 2) - 4 * c * c * u * u)) / 2);
      k3 = c * c / (a * a) - 1;
      k2 = 2 * k3 * (c - LEnd);
      k1 = k3 * (c - LEnd) * (c - LEnd) - c * c + a * a;
      printf ("k1 %f k2 %f k3 %f\n", k1, k2, k3);
      // r_lEnd = sqrt(k1+ k2*lEnd + k3*lEnd*lEnd);
      // r_lStart = r_lEnd*(lStart-LStart)/(lEnd-LStart);
    };
    return r_zExtracts;
  }
%}

DECLARE
%{
  // Scene where all geometry is added to
  Scene s;
  // point structure
  Point p1;
  // Function to handle Conic-Neutron collisions with reflectivity from McStas Tables
  double* rfront_inner; // all r-distances at lStart for all mirror surfaces
  int silicon; // +1: neutron in silicon, -1: neutron in air, 0: mirrorwidth is 0; neutron cannot be in silicon and also does not track mirror transitions
  t_Table rsTable;
%}

INITIALIZE
%{
  if (rfront_inner_file && strlen (rfront_inner_file) && strcmp (rfront_inner_file, "NULL") && strcmp (rfront_inner_file, "0")) {
    if (Table_Read (&rsTable, rfront_inner_file, 1) <= 0) { /* read 1st block data from file into pTable */
      exit (fprintf (stderr, "FlatEllipse_finite_mirror: %s: can not read file %s\n", NAME_CURRENT_COMP, rfront_inner_file));
    }
    // read the data from the file into an array and point rfron_inner to it
    nummirror = rsTable.rows;
    rfront_inner = malloc (sizeof (double) * nummirror);
    if (!rfront_inner) {
      fprintf (stderr, "Component %s: malloc() failed in INIT. Exit! \n", NAME_CURRENT_COMP);
      exit (-1);
    }
    for (int i = 0; i < nummirror; i++) {

      rfront_inner[i] = Table_Index (rsTable, i, 1); // reads the value of the second col where i sits in the first col
    }
  } else { // proceed as usual calculating the values from the outermost mirror and the number of mirrors
    printf ("automatic calulation\n");
    rfront_inner = get_r_at_z0 (nummirror, 0, r_0, lStart, sourceDist, LEnd, lStart, lEnd);
    // calculate the r-distances of all mirrors at the entry of the NMO, we will need this later
  }
  if (sourceDist == 0) { // obsolete?
    sourceDist = LStart;
  }
  silicon = (mirror_width == 0) ? 0 : -1; // neutron starts in air by default

  // Load Reflectivity Data File TODO
  // Make new scene
  s = makeScene ();

  // Set Scene to use custom trace function for conic
  // s.traceNeutronConic = traceNeutronConicWithTables;

  // Add Geometry Here
  for (int i = 0; i < nummirror; i++) {
    p1 = makePoint (rfront_inner[i], 0, lStart);
    addFlatEllipse (LStart, LEnd, p1, lStart, lEnd, -mirror_sidelength / 2, mirror_sidelength / 2, mf, R0, Qc, alpha, W, &s); // inner side of the mirror
    printf ("b[%d] = %f\n", i, rfront_inner[i]);
  }
  if (mirror_width > 0) {
    for (int i = 0; i < nummirror; i++) {
      p1 = makePoint (rfront_inner[i] + mirror_width, 0, lStart);
      addFlatEllipse (LStart, LEnd, p1, lStart, lEnd, -mirror_sidelength / 2, mirror_sidelength / 2, mb, R0, Qc, alpha, W,
                      &s); // backside of the above mirror shifted by mirror_width
    }
  }
  addEndDisk (lEnd, 0.0, 2000, &s); // neutrons will be propagated to the end of the assembly, important if they still have to move through silicon to be
                                    // refracted at the correct position addEllipsoid(-L, L,p1, -l,+l, 40,&s);
%}

TRACE
%{
  double dt;
  double x_check;

  dt = (-z + lStart) / vz;
  if (dt < 0) {
    printf ("negative time\n");
  }
  PROP_DT (dt); // propagate neutron to the entrance window of the NMO

  /* "_mctmp_a" defines a "silicon" state variable in underlying conic.h functions */
  _mctmp_a = silicon;

  if (mirror_width > 0) { // if the width of the mirrors is finite neutrons have to know whether they are in silicon or not
    x_check = fabs (x);   // lateral component of the neutron which determines whether the neutron arrives in silicon
    for (int i = 0; i < nummirror; i++) {
      dt = fabs (rfront_inner[i]);        // make sure the mirror distance to check against is positive, repeated use of same variable don't do this at home
      if (dt + mirror_width >= x_check) { // backside of the substrate further out than neutron
        if (dt <= x_check) {              // mirror itself closer to the optical axis than the neutrons, i.e., we arrive in silicon
          /* "_mctmp_a" defines a "silicon" state variable in underlying conic.h functions */
          _mctmp_a = 1;
          // First we have to refract at the entrance
          Vec nStart = makeVec (0, 0, 1); // surface normal is oriented in beam direction hopefully
          Vec init_vec = get_class_particleVel (*_particle);
          refractNeutronFlat (_particle, nStart, 0, 0.478); // m_{silicon} =  0.478 laut Peter
          break;
        }
      } else { // backside of the mirror is closer to optical axis than neutron; as all further mirrors are even closer we can break here
        break;
      }
    }
  }

  traceSingleNeutron (_particle, s);
  Vec nEnd = makeVec (0, 0, 1);
  if (_mctmp_a == 1) { // if the neutron arrives at the end of the mirror assembly while still in silicon, it will refract again at the end of the mirror
    refractNeutronFlat (_particle, nEnd, 0.478, 0); // TODO add functionality to put whatever critical angle
  }

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

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

MCDISPLAY//TODO this does not work as of now does not show the orientation of the flat conics
%{
  // PW 20260112 Suppressed. Cals to rConic make no sense since the flat mirrors are NOT conical...

  /* //Enlarge xy-plane when mcdisplay is ran with --zoom */
  /* magnify("xy"); */

  /* 	//Draw xy-axis contour for Conic Surfaces */
  /* 	int i; */
  /*   for (i = 0; i < s.num_f; i++) { */
  /*       double step = (s.f[i].ze-s.f[i].zs)/100; */
  /*       double cz; */
  /* 	    for (cz = s.f[i].zs+step; cz <= s.f[i].ze; cz+= step) { */
  /*           double rp = rConic(cz-step,s.f[i]); */
  /*           double rc = rConic(cz, s.f[i]); */

  /*           line(0,rp,cz-step,0,rc,cz); */
  /*           line(0,-rp,cz-step,0,-rc,cz); */

  /*           line(rp,0,cz-step,rc,0,cz); */
  /*           line(-rp,0,cz-step,-rc,0,cz); */
  /*       } */
  /*   } */

  /*   //Draw xy-axis cross hairs for Disks */
  /*   for (i = 0; i < s.num_di; i++) { */
  /*       line(s.di[i].r0, 0, s.di[i].z0, s.di[i].r1, 0, s.di[i].z0); */
  /*       line(-s.di[i].r0, 0, s.di[i].z0, -s.di[i].r1, 0, s.di[i].z0); */
  /*       line(0, s.di[i].r0, s.di[i].z0, 0, s.di[i].r1,s.di[i].z0); */
  /*       line(0, -s.di[i].r0, s.di[i].z0, 0, -s.di[i].r1,s.di[i].z0); */
  /*   } */
%}

END