/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2011, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Elliptic_guide_gravity
*
* %I
* Written by: Henrik Bo Hoffmann Carlsen and Mads Bertelsen
* Date: 27 Aug 2012
* Origin: NBI
*
* Perfect elliptic guide which allows for simulations with gravity.
* The guide mirrors can be divided into segments with individual m-values.
* Parabolic guide components can also be simulated.
*
* %D
*
* The Elliptic_guide_gravity component models a true elliptic guide, with separate
* ellipses in the vertical (top/bottom) and horizontal (left/right) directions,
* allowing gravity to be simulated.
*
* The guide is centered along the z-axis, with entry and exit planes parallel to the
* xy-plane. The guide entry is at z=0 (in line with other guide components).
* By default, the ellipses are defined by specifying their focal points as
* distances from the guide entry (linxw, loutxw, linyh, loutyh). In addition,
* it is possible to directly define the ellipse through the parameter set:
* majorAxisxw, minorAxisxw, majorAxisyh, minorAxisyh, together with
* majorAxisoffsetxw and majorAxisoffsetyh, where the offset parameters place the
* ellipse centers respectively.
* Further, it is possible to approximate the ellipses through parabolas, or as
* half ellipses.
*
* m-values for the guides may be specified independently for each of the four walls,
* as a scalar value for the entire wall, or as a vector distributing the m-values
* evenly (seglength left unset) or according to the vector seglength.
*
* <b>Example 1, Elliptical definition using focal points:</b>
*
* Elliptic_guide_gravity(
*                 l=50,
*                 linxw=5,linyh=5,loutxw=10,loutyh=10,
*                 xwidth=0.05,yheight=0.05,
*                 R0=0.99,Qc=0.0219,alpha=6.07,m=1.0,W=0.003
* )
*
* <b>Example 2: Half elliptical definition:</b>
*
* Elliptic_guide_gravity(
*                 l=50,
*                 linxw=5,linyh=5,loutxw=10,loutyh=10,
*                 xwidth=0.1,yheight=0.1,
*                 R0=0.99,Qc=0.0219,alpha=6.07,m=1.0,W=0.003,
*                 option = "halfEllipse",
*                 dimensionsAt = "entrance"
* )
*
* <b>Example 3: Parabolic approximation:</b>
*
* Elliptic_guide_gravity(
*                 l=50,
*                 linxw=5,linyh=5,loutxw=1e6,loutyh=1e6, // values larger than 1e8 may cause erroneous results
*                 xwidth=0.1,yheight=0.1,
*                 R0 = 0.99,Qc=0.0219,alpha=6.07,m=1.0,W=0.003,
*                 dimensionsAt = "exit"
* )
*
* <b>Example 4: Elliptical definition with varying m-values:</b>
*
* Elliptic_guide_gravity(
*                 l=50,
*                 linxw=5,linyh=5,loutxw=10,loutyh=10,
*                 xwidth=0.1,yheight=0.1,
*                 R0 = 0.99,Qc=0.0219,alpha=6.07,m=1.0,W=0.003,
*                 mvaluesright=marray,mvaluesleft=marray,mvaluestop=marray,mvaluesbottom=marray
* )
*
* where marray is initialized as
* for(iter=0; iter < 50; iter++){ marray[iter] = 2; }
* And Declared as
* double mValues[50];
* If you are using the array-based coating-specification, you **must** give nSegments a compatible value.
*
* %P
* INPUT PARAMETERS:
* mvaluesright: [pointer]   Pointer to array of m-values, right mirror
* mvaluesleft: [pointer]    - same, left mirror
* mvaluestop: [pointer]     - same, top mirror
* mvaluesbottom: [pointer]  - same, bottom mirror
* seglength: [pointer]      Pointer to array of segment lengths for discrete mirror description
* l: [m]                    length of the guide
* linxw: [m]                distance from 1st focal point to guide entrance - left and right horizontal mirrors
* loutxw: [m]               distance from 2nd focal point to guide exit - left and right horizontal mirrors
* linyh: [m]                distance from 1st focal point to guide entrance - top and bottom vertical mirrors
* loutyh: [m]               distance from 2nd focal point to guide exit - top and bottom vertical mirrors
* xwidth: [m]               width at the guide entry, mid or exit (see dimensionsAt)
* yheight: [m]              height at the guide entry, mid or exit (see dimensionsAt)
* R0: [1]                   Low-angle reflectivity
* Qc: [AA^-1]                Critical scattering vector
* alpha: [AA]               Slope of reflectivity
* m: [1]                    m-value of material for all mirrors, zero means complete absorption.
* W: [AA^-1]                 Width of supermirror cut-off
* alphatop: [AA]            Slope of reflectivity for top horizontal mirror, overwrites alpha
* mtop: [1]                 m-value of material for top horizontal mirror, overwrites m
* alphabottom: [AA]         Slope of reflectivity for bottom horizontal mirror
* mbottom: [1]              m-value of material for bottom horizontal mirror
* alpharight: [AA]          Slope of reflectivity for right vertical mirror
* mright: [1]               m-value of material for right vertical mirror
* alphaleft: [AA]           Slope of reflectivity for left vertical mirror
* mleft: [1]                m-value of material for left vertical mirror
* option: [string]          Options are 'ellipse' and 'halfEllipse'. Ellipse is defined by both the focal points, while halfEllipse locked the center of the ellipse either the entrance or exit of the guide, and use the focal point of the other end to define the ellipse
* dimensionsAt: [string]    Define whether xwidth and yheight sets the size of the opening, minor axis or the end of the guide.
** majorAxisxw: [m]         Direct definition of the guide geometry, will ignore linxw, loutxw, linyh, and loutyh if this is nonzero. Length of the axis parallel to z for the horizontal ellipse.
* minorAxisxw: [m]          Direct definition of the guide geometry, will ignore linxw, loutxw, linyh, and loutyh if this is nonzero. Length of the axis perpendicular to z for the horizontal ellipse.
* majorAxisyh: [m]          Direct definition of the guide geometry, will ignore linxw, loutxw, linyh, and loutyh if this is nonzero. Length of the axis parallel to z for the vertical ellipse.
* minorAxisyh: [m]          Direct definition of the guide geometry, will ignore linxw, loutxw, linyh, and loutyh if this is nonzero. Length of the axis perpendicular to z for the vertical ellipse.
* majorAxisoffsetxw: [m]    Direct definition of the guide geometry, distance between the center of the horizontal ellipse and the guide entrance.
* majorAxisoffsetyh: [m]    Direct definition of the guide geometry, distance between the center of the vertical ellipse and the guide entrance.
* verbose: [bool]           Give extra information about calculations
* curvature: [m]            Simulate horizontal radius of curvature by centripetal force added to the gravity. Note: Does not curve the guide in mcdisplay but "curves the neutron". Has opposite sign definition of Guide_curved.
* nSegments: [m]            Must be used to specify number of guide segments, i.e. when giving inputs mvaluesright ... etc.
* enableGravity: [m]        Flag to select propagation with gravity.
*
* OUTPUT PARAMETERS
*
* %E
*******************************************************************************/

DEFINE COMPONENT Elliptic_guide_gravity
SETTING PARAMETERS (xwidth = 0,yheight = 0,l,
	linxw = 0,loutxw = 0,
	linyh = 0,loutyh = 0,
	majorAxisxw = 0,minorAxisxw = 0,
	majorAxisyh = 0,minorAxisyh = 0,
	majorAxisoffsetxw = 0,
	majorAxisoffsetyh = 0,
	string dimensionsAt = "entrance",
	string option = "ellipse",
	R0=0.99, Qc=0.0218, alpha=6.07, m=2, W=0.003,
	alpharight=-1,  mright=-1,
	alphaleft=-1,   mleft=-1,
	alphatop=-1,    mtop=-1,
	alphabottom=-1, mbottom=-1,
	string verbose = "on",
	enableGravity = 1.0,
	curvature=0,
	int nSegments=-1,
	vector mvaluesright=NULL, vector mvaluesleft=NULL, vector mvaluestop=NULL, vector mvaluesbottom=NULL,
		    vector seglength=NULL)



SHARE
%{
  %include "ref-lib"

  ///////////////////////////////////////////////////////////////////////////
  /////////////// local structs and enums
  ///////////////////////////////////////////////////////////////////////////

  /**
  Sides of the guide
  */
  enum Side { RightSide, TopSide, LeftSide, BottomSide, None };

  /**
  The type of the collision is set in the collision function
  and decide the functions called in trace()
          Reflex (TODO change this name) calls the reflection function
          Absorb calls the built in ABSORB funtion.
          LeaveGuide calls break and end the calculations in this component
          EnterGuide does nothing
  */
  enum CollisionType { Reflex, Absorb, LeaveGuide, EnterGuide };

  /**
          The Mirror type sets the CollisionType of particles colliding on the mirror
  */
  enum MirrorType { MirrorTypeReflection, MirrorTypeTransparent, MirrorTypeabsorption };

  // enum IntersectionType {Reflex,Absorb,Transparent,Leave,Enter};

  /**
          Collision between guide and the particle
          contain infomation on the time to the next collision,
          which side of the guide it is on and whether this part of the guide
          is a perfect or approximated ellipse.
  */
  struct Intersection {
    double delta_time_to_next_collision;
    enum Side side; // A number from 0 to 4 (4 being an error warning)
    int ApproxOn;
    enum CollisionType collisionType;
  };

  /**
          Static Guide information (SGI)
          contain information on the guide, the ellipses and the mirrors on all sides
  */
  struct SGI {
    // guide infomation
    double Length;
    double entranceHorizontalWidth, entranceVerticalWidth;
    double exitHorizontalWidth, exitVerticalWidth;

    // ellipses infomation
    double ellipseMajorAxis[4], ellipseMinorAxis[4];
    double ellipseMajorOffset[4], ellipseMinorOffset[4];

    // mirror infomation
    double R0Arr[4];
    double QcArr[4];
    double alphaArr[4];
    double mArr[4];
    double WArr[4];

    // mirror type
    enum MirrorType InnerSide[4];
    enum MirrorType OuterSide[4];

    // selene
    int EnclosingBoxOn;
    double xArray[8];
    double yArray[8];
    double zArray[8];

    // segmentation
    int numberOfSegments;
    int enableSegments;
    double* mValuesright;
    double* mValuesleft;
    double* mValuestop;
    double* mValuesbottom;
    double* segLength;

    int verboseSetting;
  };

  ///////////////////////////////////////////////////////////////////////////
  /////////////// Error Handling Functions
  ///////////////////////////////////////////////////////////////////////////

  /**
          If a user input is less than zero and hence doesn't allow for a well
          define geomtric of the guide or physical values for mirrors
          @param var is the input varible there the error occurred [text]
  */
  int
  guide_elliptical_illegalInputLessThanZero (char* var, int verbose) {
    if (verbose)
      printf ("The user defined variable %s in %s has an illegal value"
              " less than zero\n",
              var, "Elliptic_guide_gravity");
    return 1;
  }

  /**
          The first focal point is in and the second is out.
          If -in-out > L then they would change position as the
          first and second focal points. This is
          @param in,out is the input varible there the error occurred [text]
  */
  int
  guide_elliptical_illegalInputFocalPointsHyperbola (char* in, char* out, double inValue, double outValue, int verbose) {
    if (verbose) {
      printf ("The user defined length of the guide, length \
            and the focal points %s and %s does not result \
            in a well defined ellipse. Swap the focal points \
            or increase L, %s, or %s to fix this problem\n",
              in, out, in, out);
      printf ("The mininum length of the should be around %e\n", inValue + outValue + 0.000001);
    }
    return 1;
  }

  /**
          Gives a warning if a part of the code is called that
          should not be accessible if the algoritmes are working correctly
          Most likely errors are floating points and ill-defined cases
  */
  void
  guide_elliptical_callCriticalWarning (char* func, int verbose) {
    if (verbose)
      printf ("A CRITICAL WARNING has been called inside %s by function %s."
              "This is most likely due to a programming error \
            inside the component. \n",
              "Elliptic_guide_gravity", func);
  }

  ///////////////////////////////////////////////////////////////////////////
  /////////////// Collision handling functions
  ///////////////////////////////////////////////////////////////////////////

  int
  guide_elliptical_getMirrorTypeFromInput (char* input, int verbose) {
    int type = -1;
    char* r1 = "reflection";
    char* r2 = "reflect";
    char* r3 = "r";
    char* a1 = "absorption";
    char* a2 = "absorb";
    char* a3 = "a";
    char* t1 = "transparant";
    char* t2 = "trans";
    char* t3 = "t";
    if (strcmp (input, r1) == 0 || strcmp (input, r2) == 0 || strcmp (input, r3) == 0)
      type = MirrorTypeReflection;
    if (strcmp (input, a1) == 0 || strcmp (input, a2) == 0 || strcmp (input, a3) == 0)
      type = MirrorTypeabsorption;
    if (strcmp (input, t1) == 0 || strcmp (input, t2) == 0 || strcmp (input, t3) == 0)
      type = MirrorTypeTransparent;
    if (type == -1 && verbose)
      printf ("Following string is not a valid type of a mirror: %s,"
              "use reflection,absorption or transparant. \n",
              input);

    return type;
  }

  ///////////////////////////////////////////////////////////////////////////
  /////////////// Collision functions
  ///////////////////////////////////////////////////////////////////////////

  /**
          Find the intersection between the neutron and the ellipse using newton method.
          As there is up to 4 solution to this problem, and only the
          smallest positive root is the physical solution. Using the tuning points
          it is possible to look the only the potential roots to speed up calculations.

          @param coef; A pointer to the array holding the coeffecients
                          for the 4th order polynomial.
          @param startPosition, The default starting point for newton method. [s]
          @param limit; A point after all the roots of the polynial. [s]
          @param solution A pointer which will hold the physical solution
                          if this function return true.
          @return; return 1 if the physical solution is found. [boolean]
  */
  #pragma acc routine seq
  double
  guide_elliptical_foverdf (double* coefficients, double currentPoint) {
    double numerator = coefficients[0] * currentPoint * currentPoint * currentPoint * currentPoint + coefficients[1] * currentPoint * currentPoint * currentPoint
                       + coefficients[2] * currentPoint * currentPoint + coefficients[3] * currentPoint + coefficients[4];
    double denominator = 4 * coefficients[0] * currentPoint * currentPoint * currentPoint + 3 * coefficients[1] * currentPoint * currentPoint
                         + 2 * coefficients[2] * currentPoint + coefficients[3];
    return numerator / denominator;
  }
  #pragma acc routine seq
  int
  guide_elliptical_newtonRapsonsMethod4thOrder (double* coefficients, double* solution, double startingPoint, double tolerance, double max_iterations) {

    double numerator;
    double denominator;
    double t_previous;
    double t = startingPoint;
    int iteration = 0;

    do {
      t_previous = t;
      t = t_previous - guide_elliptical_foverdf (coefficients, t);
      iteration++;
    } while (fabs (t - t_previous) > tolerance && iteration < max_iterations);
    if (iteration == max_iterations) {
      return 0;
    } else {
      *solution = t;
      return 1;
    }
  }

  #pragma acc routine seq
  int
  guide_elliptical_findNeutronEllipseIntersection (double* coef, double startPosition, double limit, double* solution) {

    // in the case of no gravity
    if (coef[0] == 0 & coef[1] == 0) {
      double t1 = 0;
      double t2 = 0;
      int boolean = solve_2nd_order (&t1, &t2, coef[2], coef[3], coef[4]);

      if (t1 > startPosition) {
        *solution = t1;
      }
      if (t2 > startPosition) {
        *solution = t2;
      }
      return boolean;
    }

    double tol = 1e-15;
    double max_iter = 1e3;
    double turningP1, turningP2;

    double sp = startPosition;
    int inside;
    if (coef[0] * sp * sp * sp * sp + coef[1] * sp * sp * sp + coef[2] * sp * sp + coef[3] * sp + coef[4] < 0)
      inside = 1;
    else
      inside = 0;

    int boolean = solve_2nd_order (&turningP1, &turningP2, 12 * coef[0], 6 * coef[1], 2 * coef[2]);

    double t1 = 0, t2 = 0;
    double ss = 100;

    if (inside) {
      if (boolean)
        guide_elliptical_newtonRapsonsMethod4thOrder (coef, &t1, turningP1, tol, max_iter);
      guide_elliptical_newtonRapsonsMethod4thOrder (coef, &t2, limit, tol, max_iter);
    } else {
      if (boolean)
        guide_elliptical_newtonRapsonsMethod4thOrder (coef, &t1, turningP2, tol, max_iter);
      guide_elliptical_newtonRapsonsMethod4thOrder (coef, &t2, startPosition, tol, max_iter);
    }

    if (ss > t1 && t1 > 1e-15)
      ss = t1;
    if (ss > t2 && t2 > 1e-15)
      ss = t2;
    *solution = ss;

    return 1;
  }

  #pragma acc routine seq
  int
  guide_elliptical_handleGuideIntersection (double x, double y, double z, double vx, double vy, double vz, double Gx, double Gy, double Gz, struct SGI* guideInfo,
                                            struct Intersection* currentCollision) {
    //
    double horExS = 1 / (guideInfo->ellipseMinorAxis[RightSide] * guideInfo->ellipseMinorAxis[RightSide]);
    double horEzS = 1 / (guideInfo->ellipseMajorAxis[RightSide] * guideInfo->ellipseMajorAxis[RightSide]);
    double hordiffx = x - guideInfo->ellipseMinorOffset[RightSide];
    double hordiffz = z - guideInfo->ellipseMajorOffset[RightSide];

    double horAlpha = (Gx * Gx * horExS + Gz * Gz * horEzS) / 4;
    double horBeta = (Gx * vx * horExS + Gz * vz * horEzS);
    double horGamma = horExS * vx * vx + horEzS * vz * vz + horExS * Gx * hordiffx + horEzS * Gz * hordiffz;
    double horDelta = 2 * horExS * vx * hordiffx + 2 * horEzS * vz * hordiffz;
    double horEpsilon = horExS * hordiffx * hordiffx + horEzS * hordiffz * hordiffz - 1;

    double horCoefficients[5] = { horAlpha, horBeta, horGamma, horDelta, horEpsilon };

    double verEyS = 1 / (guideInfo->ellipseMinorAxis[TopSide] * guideInfo->ellipseMinorAxis[TopSide]);
    double verEzS = 1 / (guideInfo->ellipseMajorAxis[TopSide] * guideInfo->ellipseMajorAxis[TopSide]);
    double verdiffy = y - guideInfo->ellipseMinorOffset[TopSide];
    double verdiffz = z - guideInfo->ellipseMajorOffset[TopSide];

    double verAlpha = (Gy * Gy * verEyS + Gz * Gz * verEzS) / 4;
    double verBeta = (Gy * vy * verEyS + Gz * vz * verEzS);
    double verGamma = verEyS * vy * vy + verEzS * vz * vz + verEyS * Gy * verdiffy + verEzS * Gz * verdiffz;
    double verDelta = 2 * verEyS * vy * verdiffy + 2 * verEzS * vz * verdiffz;
    double verEpsilon = verEyS * verdiffy * verdiffy + verEzS * verdiffz * verdiffz - 1;

    double verCoefficients[5] = { verAlpha, verBeta, verGamma, verDelta, verEpsilon };

    double upperlimit;
    double startingPoint = 1e-15;

    int boolean = -1;
    // Horizontal
    double solutionH = 0;
    solve_2nd_order (&upperlimit, NULL, -0.5 * Gz, -vz, 2 * guideInfo->ellipseMajorAxis[RightSide] - z);
    int booleanH = guide_elliptical_findNeutronEllipseIntersection (horCoefficients, startingPoint, upperlimit, &solutionH);
    // Vertical
    double solutionV = 0;
    solve_2nd_order (&upperlimit, NULL, -0.5 * Gz, -vz, 2 * guideInfo->ellipseMajorAxis[TopSide] - z);
    int booleanV = guide_elliptical_findNeutronEllipseIntersection (verCoefficients, startingPoint, upperlimit, &solutionV);

    if (solutionH <= 0)
      currentCollision->delta_time_to_next_collision = solutionV;
    else if (solutionV <= 0)
      currentCollision->delta_time_to_next_collision = solutionH;
    else if (fabs (solutionH - solutionV) < 1e-12)
      return 0;
    else if (solutionH < solutionV) {
      currentCollision->delta_time_to_next_collision = solutionH;
      boolean = booleanH;
    } else {
      currentCollision->delta_time_to_next_collision = solutionV;
      boolean = booleanV;
    }

    double tside = currentCollision->delta_time_to_next_collision;
    double xside = x + vx * tside + 0.5 * Gx * tside * tside;
    double yside = y + vy * tside + 0.5 * Gy * tside * tside;
    double zside = z + vz * tside + 0.5 * Gz * tside * tside;

    double xfactor = 2
                     * sqrt (1
                             - ((zside - guideInfo->ellipseMajorOffset[RightSide]) * (zside - guideInfo->ellipseMajorOffset[RightSide]))
                                   / (guideInfo->ellipseMajorAxis[RightSide] * guideInfo->ellipseMajorAxis[RightSide]))
                     * guideInfo->ellipseMinorAxis[RightSide];

    double yfactor = 2
                     * sqrt (1
                             - ((zside - guideInfo->ellipseMajorOffset[BottomSide]) * (zside - guideInfo->ellipseMajorOffset[BottomSide]))
                                   / (guideInfo->ellipseMajorAxis[BottomSide] * guideInfo->ellipseMajorAxis[BottomSide]))
                     * guideInfo->ellipseMinorAxis[BottomSide];

    xside = xside / xfactor;
    yside = yside / yfactor;
    if (fabs (yside) >= fabs (xside)) {
      if (y > 0)
        currentCollision->side = TopSide;
      else
        currentCollision->side = BottomSide;
    } else {
      if (x < 0)
        currentCollision->side = RightSide;
      else
        currentCollision->side = LeftSide;
    }
    if (tside < 1e-15)
      printf ("low time is: %e\n", tside);

    return boolean;
  }

  ///////////////////////////////////////////////////////////////////////////
  /////////////// reflection functions
  ///////////////////////////////////////////////////////////////////////////

  /**
          Calculate the new velocity vector for the particle colliding on
          the inner side of the elliptic mirror and returns the loss-factor (TODO)

          @param pos_V0,pos_W0 Is the 2d position vector of the particle,
                          assumed to be a point on the ellipse. [m]
          @param pvel_V0,pvel_W0 Is the 2d velocity vector of the particle. [m/s]
          @param ellipse_V_axis_squared,ellipse_W_axis_squared
                          are the axes of the ellipse. [m]
          @param ellipse_V_offset,ellipse_W_offset Is the 2d vector difference
                          between the ellipse coordinate system (center of the ellipse)
                          and the guide coordinate system [m]
          @param R0, Mvalue, Qc, W, Alpha #TODO
                  slaa beskrivelse af disse variabler i andre dokumenter
                  og hold dig til standarden.
          @return the new wieght of the package
  */
  #pragma acc routine seq
  double
  guide_elliptical_ReflectionOnEllipticSurface (double pos_V, double pos_W, double* pvel_V, double* pvel_W, double ellipse_V_axis, double ellipse_W_axis,
                                                double ellipse_V_offset, double ellipse_W_offset, double R0, double Qc, double alpha, double Mvalue, double W) {

    // Turns the velocity vector (vel_V0,vel_W0) into a local value
    double vel_V = *pvel_V;
    double vel_W = *pvel_W;

    // Galilean transformation of the particles start position
    // to the ellipse coordinate system
    pos_V = pos_V - ellipse_V_offset;
    pos_W = pos_W - ellipse_W_offset;

    /*
     *   If we reflect the velocity vector in the normal
     *	to the ellipse in the point of intersection
     *   The resulting vector will be -f2, do to conservation of momentum.
     *   this result in the following equation
     *   f2 = -f1 + 2(f1 dot nhat)nhat
     *   which is equal to f2 = f1 - 2(f1 dot n)n/nlength^2
     */

    // The normal vector to the point of intersection
    double normVec_V = -pos_W * ellipse_V_axis / ellipse_W_axis;
    double normVec_W = pos_V * ellipse_W_axis / ellipse_V_axis;

    double normVec_length_squared = normVec_V * normVec_V + normVec_W * normVec_W;

    // Dot product of (vel_V0,vel_W0) and the normal vector
    double Vel_dot_NV = vel_V * normVec_V + vel_W * normVec_W;

    // Calculate f2
    double vel_V_2 = -vel_V + 2 * Vel_dot_NV * normVec_V / normVec_length_squared;
    double vel_W_2 = -vel_W + 2 * Vel_dot_NV * normVec_W / normVec_length_squared;

    // Apply the new velocity vector to the particle globally
    *pvel_V = vel_V_2;
    *pvel_W = vel_W_2;

    // Calculate q and the weighting of the neutron package
    // q=f1-f2
    double delta_vel_V = vel_V - vel_V_2;
    double delta_vel_W = vel_W - vel_W_2;
    double q = V2Q * sqrt (delta_vel_V * delta_vel_V + delta_vel_W * delta_vel_W);

    // Calculate the loss of neutrons due to the reflection
    double mirrorPar[] = { R0, Qc, alpha, Mvalue, W };
    double weight = 1.0;
    StdReflecFunc (q, mirrorPar, &weight);

    return weight;
  }

  /**
          Use the found side of Intersection to call guide_elliptical_ReflectionOnEllipticSurface with
          the parameters of that side.
  */
  #pragma acc routine seq
  double
  guide_elliptical_handleReflection (double x0, double y0, double z0, double* vx_p, double* vy_p, double* vz_p, struct SGI* sgi, struct Intersection* cc) {

    if (!sgi->enableSegments) {
      if (cc->side == RightSide || cc->side == LeftSide)
        return guide_elliptical_ReflectionOnEllipticSurface (x0, z0, vx_p, vz_p, sgi->ellipseMinorAxis[cc->side], sgi->ellipseMajorAxis[cc->side],
                                                             sgi->ellipseMinorOffset[cc->side], sgi->ellipseMajorOffset[cc->side], sgi->R0Arr[cc->side],
                                                             sgi->QcArr[cc->side], sgi->alphaArr[cc->side], sgi->mArr[cc->side], sgi->WArr[cc->side]);
      if (cc->side == TopSide || cc->side == BottomSide)
        return guide_elliptical_ReflectionOnEllipticSurface (y0, z0, vy_p, vz_p, sgi->ellipseMinorAxis[cc->side], sgi->ellipseMajorAxis[cc->side],
                                                             sgi->ellipseMinorOffset[cc->side], sgi->ellipseMajorOffset[cc->side], sgi->R0Arr[cc->side],
                                                             sgi->QcArr[cc->side], sgi->alphaArr[cc->side], sgi->mArr[cc->side], sgi->WArr[cc->side]);
    } else {
      int currentSegment = -1;
      double combinedLength = 0;
      int i;
      for (i = 0; i < sgi->numberOfSegments; i++) {
        combinedLength = combinedLength + sgi->segLength[i];
        if (z0 < combinedLength) {
          currentSegment = i;
          break;
        }
      }
      if (currentSegment < 0) {
        printf ("Elliptic_guide_gravity: Error indexing guide segment\n");
        return 0;
      }

      if (cc->side == RightSide)
        return guide_elliptical_ReflectionOnEllipticSurface (x0, z0, vx_p, vz_p, sgi->ellipseMinorAxis[cc->side], sgi->ellipseMajorAxis[cc->side],
                                                             sgi->ellipseMinorOffset[cc->side], sgi->ellipseMajorOffset[cc->side], sgi->R0Arr[cc->side],
                                                             sgi->QcArr[cc->side], sgi->alphaArr[cc->side], sgi->mValuesright[currentSegment],
                                                             sgi->WArr[cc->side]);
      if (cc->side == LeftSide)
        return guide_elliptical_ReflectionOnEllipticSurface (x0, z0, vx_p, vz_p, sgi->ellipseMinorAxis[cc->side], sgi->ellipseMajorAxis[cc->side],
                                                             sgi->ellipseMinorOffset[cc->side], sgi->ellipseMajorOffset[cc->side], sgi->R0Arr[cc->side],
                                                             sgi->QcArr[cc->side], sgi->alphaArr[cc->side], sgi->mValuesleft[currentSegment],
                                                             sgi->WArr[cc->side]);
      if (cc->side == TopSide)
        return guide_elliptical_ReflectionOnEllipticSurface (y0, z0, vy_p, vz_p, sgi->ellipseMinorAxis[cc->side], sgi->ellipseMajorAxis[cc->side],
                                                             sgi->ellipseMinorOffset[cc->side], sgi->ellipseMajorOffset[cc->side], sgi->R0Arr[cc->side],
                                                             sgi->QcArr[cc->side], sgi->alphaArr[cc->side], sgi->mValuestop[currentSegment], sgi->WArr[cc->side]);
      if (cc->side == BottomSide)
        return guide_elliptical_ReflectionOnEllipticSurface (y0, z0, vy_p, vz_p, sgi->ellipseMinorAxis[cc->side], sgi->ellipseMajorAxis[cc->side],
                                                             sgi->ellipseMinorOffset[cc->side], sgi->ellipseMajorOffset[cc->side], sgi->R0Arr[cc->side],
                                                             sgi->QcArr[cc->side], sgi->alphaArr[cc->side], sgi->mValuesbottom[currentSegment],
                                                             sgi->WArr[cc->side]);
    }
    return 0;
  }

  ///////////////////////////////////////////////////////////////////////////
  /////////////// End of functions
  ///////////////////////////////////////////////////////////////////////////
%}


DECLARE
%{
  /**
          All variables below is locally declared
          and hence accessible through OUTPUT PARAMETERS.
  */
  struct SGI guideInfo; // Static Guide information is set in INITIALIZE

  double Gx0; // Local gravity vector, is set once in INITIALIZE
  double Gy0;
  double Gz0;
  double Circ;
  double* dynamicalSegLength;
%}

INITIALIZE
%{
  ///////////////////////////////////////////////////////////////////////////
  /////////////// Test user input
  ///////////////////////////////////////////////////////////////////////////

  if (strcmp (verbose, "on") == 0)
    guideInfo.verboseSetting = 1;
  else
    guideInfo.verboseSetting = 0;

  guideInfo.R0Arr[RightSide] = R0;
  guideInfo.QcArr[RightSide] = Qc;
  guideInfo.alphaArr[RightSide] = alpharight;
  guideInfo.mArr[RightSide] = mright;
  guideInfo.WArr[RightSide] = W;

  guideInfo.R0Arr[TopSide] = R0;
  guideInfo.QcArr[TopSide] = Qc;
  guideInfo.alphaArr[TopSide] = alphatop;
  guideInfo.mArr[TopSide] = mtop;
  guideInfo.WArr[TopSide] = W;

  guideInfo.R0Arr[LeftSide] = R0;
  guideInfo.QcArr[LeftSide] = Qc;
  guideInfo.alphaArr[LeftSide] = alphaleft;
  guideInfo.mArr[LeftSide] = mleft;
  guideInfo.WArr[LeftSide] = W;

  guideInfo.R0Arr[BottomSide] = R0;
  guideInfo.QcArr[BottomSide] = Qc;
  guideInfo.alphaArr[BottomSide] = alphabottom;
  guideInfo.mArr[BottomSide] = mbottom;
  guideInfo.WArr[BottomSide] = W;

  int sides;
  for (sides = RightSide; sides <= BottomSide; sides++) {
    if (guideInfo.alphaArr[sides] == -1)
      guideInfo.alphaArr[sides] = alpha;
    if (guideInfo.mArr[sides] == -1)
      guideInfo.mArr[sides] = m;
  }

  // Test user input for illegal values
  int inputErrors = 0;
  // Lower or equal to zero
  if (l <= 0)
    inputErrors += guide_elliptical_illegalInputLessThanZero ("length", guideInfo.verboseSetting);
  if (guideInfo.alphaArr[TopSide] < 0)
    inputErrors += guide_elliptical_illegalInputLessThanZero ("alphatop", guideInfo.verboseSetting);
  if (guideInfo.mArr[TopSide] < 0)
    inputErrors += guide_elliptical_illegalInputLessThanZero ("mtop", guideInfo.verboseSetting);

  if (guideInfo.alphaArr[BottomSide] < 0)
    inputErrors += guide_elliptical_illegalInputLessThanZero ("alphabottom", guideInfo.verboseSetting);
  if (guideInfo.mArr[BottomSide] < 0)
    inputErrors += guide_elliptical_illegalInputLessThanZero ("mbottom", guideInfo.verboseSetting);

  if (guideInfo.alphaArr[RightSide] < 0)
    inputErrors += guide_elliptical_illegalInputLessThanZero ("alpharight", guideInfo.verboseSetting);
  if (guideInfo.mArr[RightSide] < 0)
    inputErrors += guide_elliptical_illegalInputLessThanZero ("mright", guideInfo.verboseSetting);

  if (guideInfo.alphaArr[LeftSide] < 0)
    inputErrors += guide_elliptical_illegalInputLessThanZero ("alphaleft", guideInfo.verboseSetting);
  if (guideInfo.mArr[LeftSide] < 0)
    inputErrors += guide_elliptical_illegalInputLessThanZero ("mleft", guideInfo.verboseSetting);

  // Focal points result in hyperbola instead of an ellipse
  if (l <= -linxw - loutxw)
    inputErrors += guide_elliptical_illegalInputFocalPointsHyperbola ("linw", "loutw", linxw, loutxw, guideInfo.verboseSetting);
  if (l <= -linyh - loutyh)
    inputErrors += guide_elliptical_illegalInputFocalPointsHyperbola ("linh", "louth", linyh, loutyh, guideInfo.verboseSetting);

  if (strcmp (dimensionsAt, "entrance") != 0 && strcmp (dimensionsAt, "mid") != 0 && strcmp (dimensionsAt, "exit") != 0) {
    inputErrors += 1;
    printf ("dimensionsAt were given an incorrect input."
            "Input must be string containing \"entrance\",\"mid\" or \"exit\" \n");
  }

  // Terminate program if any input errors occurred
  if (inputErrors != 0) {
    exit (printf ("\nCRITICAL ERROR(S) IN COMPONENT %s"
                  " CONSIDER CHECKING USER INPUT AS %d INPUT ERRORS WAS FOUND.\n",
                  NAME_CURRENT_COMP, inputErrors));
  }

  ///////////////////////////////////////////////////////////////////////////
  /////////////// Calculate internal guide values from user input
  ///////////////////////////////////////////////////////////////////////////

  /*  Calculate the focal lines for the ellipses.
          These can be used to calculate the axes of the ellipses
          using pyth and definition of the ellipse that says distance
          between the foci and every point on the ellipse is constant.
  */
  int directDefinition = 0;

  if (majorAxisyh != 0 || minorAxisyh != 0 || majorAxisxw != 0 || minorAxisxw != 0) {
    directDefinition = 1;
    guideInfo.Length = l;

    guideInfo.ellipseMajorAxis[RightSide] = majorAxisxw;
    guideInfo.ellipseMinorAxis[RightSide] = minorAxisxw;
    guideInfo.ellipseMajorOffset[RightSide] = majorAxisoffsetxw;
    guideInfo.ellipseMinorOffset[RightSide] = 0;

    guideInfo.ellipseMajorAxis[TopSide] = majorAxisyh;
    guideInfo.ellipseMinorAxis[TopSide] = minorAxisyh;
    guideInfo.ellipseMajorOffset[TopSide] = majorAxisoffsetyh;
    guideInfo.ellipseMinorOffset[TopSide] = 0;

    guideInfo.ellipseMajorAxis[LeftSide] = majorAxisxw;
    guideInfo.ellipseMinorAxis[LeftSide] = minorAxisxw;
    guideInfo.ellipseMajorOffset[LeftSide] = majorAxisoffsetxw;
    guideInfo.ellipseMinorOffset[LeftSide] = 0;

    guideInfo.ellipseMajorAxis[BottomSide] = majorAxisyh;
    guideInfo.ellipseMinorAxis[BottomSide] = minorAxisyh;
    guideInfo.ellipseMajorOffset[BottomSide] = majorAxisoffsetyh;
    guideInfo.ellipseMinorOffset[BottomSide] = 0;

    guideInfo.entranceHorizontalWidth = 2 * sqrt (1 - (majorAxisoffsetyh * majorAxisoffsetyh) / (majorAxisyh * majorAxisyh)) * minorAxisyh;
    guideInfo.entranceVerticalWidth = 2 * sqrt (1 - (majorAxisoffsetxw * majorAxisoffsetxw) / (majorAxisxw * majorAxisxw)) * minorAxisxw;
  }

  if (strcmp (option, "ellipse") == 0 && directDefinition == 0) {
    if (strcmp (dimensionsAt, "entrance") == 0) {
      double lofbs_horizontal = sqrt (linxw * linxw + xwidth * xwidth * 0.25) + sqrt ((l + loutxw) * (l + loutxw) + xwidth * xwidth * 0.25);

      double lofbs_vertical = sqrt (linyh * linyh + yheight * yheight * 0.25) + sqrt ((l + loutyh) * (l + loutyh) + yheight * yheight * 0.25);

      guideInfo.Length = l;

      guideInfo.ellipseMajorAxis[RightSide] = lofbs_horizontal / 2;
      guideInfo.ellipseMinorAxis[RightSide] = sqrt (0.25 * lofbs_horizontal * lofbs_horizontal - 0.25 * (l + linxw + loutxw) * (l + linxw + loutxw));

      guideInfo.ellipseMajorOffset[RightSide] = (l + linxw + loutxw) / 2 - linxw;
      guideInfo.ellipseMinorOffset[RightSide] = 0;

      guideInfo.ellipseMajorAxis[LeftSide] = guideInfo.ellipseMajorAxis[RightSide];
      guideInfo.ellipseMinorAxis[LeftSide] = guideInfo.ellipseMinorAxis[RightSide];
      guideInfo.ellipseMajorOffset[LeftSide] = guideInfo.ellipseMajorOffset[RightSide];
      guideInfo.ellipseMinorOffset[LeftSide] = guideInfo.ellipseMinorOffset[RightSide];

      guideInfo.ellipseMajorAxis[TopSide] = lofbs_vertical / 2;

      guideInfo.ellipseMinorAxis[TopSide] = sqrt (0.25 * lofbs_vertical * lofbs_vertical - 0.25 * (l + linyh + loutyh) * (l + linyh + loutyh));

      guideInfo.ellipseMajorOffset[TopSide] = (l + linyh + loutyh) / 2 - linyh;
      guideInfo.ellipseMinorOffset[TopSide] = 0;

      guideInfo.ellipseMajorAxis[BottomSide] = guideInfo.ellipseMajorAxis[TopSide];
      guideInfo.ellipseMinorAxis[BottomSide] = guideInfo.ellipseMinorAxis[TopSide];
      guideInfo.ellipseMajorOffset[BottomSide] = guideInfo.ellipseMajorOffset[TopSide];
      guideInfo.ellipseMinorOffset[BottomSide] = guideInfo.ellipseMinorOffset[TopSide];
    }
    if (strcmp (dimensionsAt, "exit") == 0) {
      double lofbs_horizontal = sqrt (loutxw * loutxw + xwidth * xwidth * 0.25) + sqrt ((l + linxw) * (l + linxw) + xwidth * xwidth * 0.25);

      double lofbs_vertical = sqrt (loutyh * loutyh + yheight * yheight * 0.25) + sqrt ((l + linyh) * (l + linyh) + yheight * yheight * 0.25);

      guideInfo.Length = l;

      guideInfo.ellipseMajorAxis[RightSide] = lofbs_horizontal / 2;
      guideInfo.ellipseMinorAxis[RightSide] = sqrt (0.25 * lofbs_horizontal * lofbs_horizontal - 0.25 * (l + linxw + loutxw) * (l + linxw + loutxw));

      guideInfo.ellipseMajorOffset[RightSide] = (l + linxw + loutxw) / 2 - linxw;
      guideInfo.ellipseMinorOffset[RightSide] = 0;

      guideInfo.ellipseMajorAxis[LeftSide] = guideInfo.ellipseMajorAxis[RightSide];
      guideInfo.ellipseMinorAxis[LeftSide] = guideInfo.ellipseMinorAxis[RightSide];
      guideInfo.ellipseMajorOffset[LeftSide] = guideInfo.ellipseMajorOffset[RightSide];
      guideInfo.ellipseMinorOffset[LeftSide] = guideInfo.ellipseMinorOffset[RightSide];

      guideInfo.ellipseMajorAxis[TopSide] = lofbs_vertical / 2;

      guideInfo.ellipseMinorAxis[TopSide] = sqrt (0.25 * lofbs_vertical * lofbs_vertical - 0.25 * (l + linyh + loutyh) * (l + linyh + loutyh));

      guideInfo.ellipseMajorOffset[TopSide] = (l + linyh + loutyh) / 2 - linyh;
      guideInfo.ellipseMinorOffset[TopSide] = 0;

      guideInfo.ellipseMajorAxis[BottomSide] = guideInfo.ellipseMajorAxis[TopSide];
      guideInfo.ellipseMinorAxis[BottomSide] = guideInfo.ellipseMinorAxis[TopSide];
      guideInfo.ellipseMajorOffset[BottomSide] = guideInfo.ellipseMajorOffset[TopSide];
      guideInfo.ellipseMinorOffset[BottomSide] = guideInfo.ellipseMinorOffset[TopSide];
    }
    if (strcmp (dimensionsAt, "mid") == 0) {

      guideInfo.Length = l;

      guideInfo.ellipseMajorAxis[RightSide] = sqrt ((linxw + l + loutxw) * (linxw + l + loutxw) / 4 + xwidth * xwidth / 4);
      guideInfo.ellipseMinorAxis[RightSide] = xwidth / 2;

      guideInfo.ellipseMajorOffset[RightSide] = (l + linxw + loutxw) / 2 - linxw;
      guideInfo.ellipseMinorOffset[RightSide] = 0;

      guideInfo.ellipseMajorAxis[LeftSide] = guideInfo.ellipseMajorAxis[RightSide];
      guideInfo.ellipseMinorAxis[LeftSide] = guideInfo.ellipseMinorAxis[RightSide];
      guideInfo.ellipseMajorOffset[LeftSide] = guideInfo.ellipseMajorOffset[RightSide];
      guideInfo.ellipseMinorOffset[LeftSide] = guideInfo.ellipseMinorOffset[RightSide];

      guideInfo.ellipseMajorAxis[TopSide] = sqrt ((linyh + l + loutyh) * (linyh + l + loutyh) / 4 + yheight * yheight / 4);
      guideInfo.ellipseMinorAxis[TopSide] = yheight / 2;

      guideInfo.ellipseMajorOffset[TopSide] = (l + linyh + loutyh) / 2 - linyh;
      guideInfo.ellipseMinorOffset[TopSide] = 0;

      guideInfo.ellipseMajorAxis[BottomSide] = guideInfo.ellipseMajorAxis[TopSide];
      guideInfo.ellipseMinorAxis[BottomSide] = guideInfo.ellipseMinorAxis[TopSide];
      guideInfo.ellipseMajorOffset[BottomSide] = guideInfo.ellipseMajorOffset[TopSide];
      guideInfo.ellipseMinorOffset[BottomSide] = guideInfo.ellipseMinorOffset[TopSide];
    }
  }

  guideInfo.entranceHorizontalWidth = 2
                                      * sqrt (1
                                              - guideInfo.ellipseMajorOffset[RightSide] * guideInfo.ellipseMajorOffset[RightSide]
                                                    / (guideInfo.ellipseMajorAxis[RightSide] * guideInfo.ellipseMajorAxis[RightSide]))
                                      * guideInfo.ellipseMinorAxis[RightSide];
  guideInfo.entranceVerticalWidth = 2
                                    * sqrt (1
                                            - guideInfo.ellipseMajorOffset[TopSide] * guideInfo.ellipseMajorOffset[TopSide]
                                                  / (guideInfo.ellipseMajorAxis[TopSide] * guideInfo.ellipseMajorAxis[TopSide]))
                                    * guideInfo.ellipseMinorAxis[TopSide];

  if (strcmp (option, "halfellipse") == 0 && directDefinition == 0) {
    exit (printf ("Critical error in %s; the option for option = halfellipse is currently disabled.", NAME_CURRENT_COMP));

    double used_focal_vertical;
    double used_focal_horizontal;
    double major_offset_horizontal = 0;
    double major_offset_vertical = 0;

    if (strcmp (dimensionsAt, "entrance") == 0) {
      used_focal_vertical = sqrt ((yheight * yheight) / 4 + (l + linyh) * (l + linyh));
      used_focal_horizontal = sqrt ((xwidth * xwidth) / 4 + (l + linxw) * (l + linxw));
      major_offset_vertical = l;
      major_offset_horizontal = l;
    } else {
      used_focal_vertical = sqrt ((yheight * yheight) / 4 + (l + loutyh) * (l + loutyh));
      used_focal_horizontal = sqrt ((xwidth * xwidth) / 4 + (l + loutxw) * (l + loutxw));
    }

    guideInfo.Length = l;

    guideInfo.ellipseMajorAxis[RightSide] = used_focal_horizontal;
    guideInfo.ellipseMinorAxis[RightSide] = xwidth / 2;

    guideInfo.ellipseMajorOffset[RightSide] = major_offset_horizontal;
    guideInfo.ellipseMinorOffset[RightSide] = 0;

    guideInfo.ellipseMajorAxis[LeftSide] = guideInfo.ellipseMajorAxis[RightSide];
    guideInfo.ellipseMinorAxis[LeftSide] = guideInfo.ellipseMinorAxis[RightSide];
    guideInfo.ellipseMajorOffset[LeftSide] = guideInfo.ellipseMajorOffset[RightSide];
    guideInfo.ellipseMinorOffset[LeftSide] = guideInfo.ellipseMinorOffset[RightSide];

    guideInfo.ellipseMajorAxis[TopSide] = used_focal_vertical;
    guideInfo.ellipseMinorAxis[TopSide] = yheight / 2;

    guideInfo.ellipseMajorOffset[TopSide] = major_offset_vertical;
    guideInfo.ellipseMinorOffset[TopSide] = 0;

    guideInfo.ellipseMajorAxis[BottomSide] = guideInfo.ellipseMajorAxis[TopSide];
    guideInfo.ellipseMinorAxis[BottomSide] = guideInfo.ellipseMinorAxis[TopSide];
    guideInfo.ellipseMajorOffset[BottomSide] = guideInfo.ellipseMajorOffset[TopSide];
    guideInfo.ellipseMinorOffset[BottomSide] = guideInfo.ellipseMinorOffset[TopSide];
  }

  // Applies the properties of the mirrors in the guide given by the user.
  // These variables are used in the reflection functions.

  // Give a warning if all side of the guide is turned off,
  // as the guide is essentially turned off
  if (guideInfo.OuterSide[RightSide] == 1 && guideInfo.OuterSide[TopSide] == 1 && guideInfo.OuterSide[LeftSide] == 1 && guideInfo.OuterSide[BottomSide] == 1
      && guideInfo.InnerSide[RightSide] == 1 && guideInfo.InnerSide[TopSide] == 1 && guideInfo.InnerSide[LeftSide] == 1 && guideInfo.InnerSide[BottomSide] == 1)
    printf ("Warning: In %s all the sides of the guide has been disabled,"
            " so it not possible for any particle"
            " to collide with the guide, consider"
            " disabling this component",
            NAME_CURRENT_COMP);

  if (guideInfo.mArr[RightSide] <= 0)
    guideInfo.InnerSide[RightSide] = MirrorTypeabsorption;
  if (guideInfo.mArr[TopSide] <= 0)
    guideInfo.InnerSide[TopSide] = MirrorTypeabsorption;
  if (guideInfo.mArr[LeftSide] <= 0)
    guideInfo.InnerSide[LeftSide] = MirrorTypeabsorption;
  if (guideInfo.mArr[BottomSide] <= 0)
    guideInfo.InnerSide[BottomSide] = MirrorTypeabsorption;

  if (strcmp (option, "halfellipse") == 0 && directDefinition == 0) {
    guideInfo.entranceHorizontalWidth = (guideInfo.ellipseMinorAxis[RightSide]
                                             * sqrt (1
                                                     - (guideInfo.ellipseMajorOffset[RightSide] * guideInfo.ellipseMajorOffset[RightSide])
                                                           / (guideInfo.ellipseMajorAxis[RightSide] * guideInfo.ellipseMajorAxis[RightSide]))
                                         + guideInfo.ellipseMinorOffset[RightSide])
                                        * 2;
    guideInfo.entranceVerticalWidth = (guideInfo.ellipseMinorAxis[TopSide]
                                           * sqrt (1
                                                   - (guideInfo.ellipseMajorOffset[TopSide] * guideInfo.ellipseMajorOffset[TopSide])
                                                         / (guideInfo.ellipseMajorAxis[TopSide] * guideInfo.ellipseMajorAxis[TopSide]))
                                       + guideInfo.ellipseMinorOffset[TopSide])
                                      * 2;
  }

  guideInfo.EnclosingBoxOn = 0;

  guideInfo.exitVerticalWidth
      = 2
        * sqrt (1
                - ((guideInfo.Length - guideInfo.ellipseMajorOffset[BottomSide]) * (guideInfo.Length - guideInfo.ellipseMajorOffset[BottomSide]))
                      / (guideInfo.ellipseMajorAxis[BottomSide] * guideInfo.ellipseMajorAxis[BottomSide]))
        * guideInfo.ellipseMinorAxis[BottomSide];

  guideInfo.exitHorizontalWidth
      = 2
        * sqrt (1
                - ((guideInfo.Length - guideInfo.ellipseMajorOffset[RightSide]) * (guideInfo.Length - guideInfo.ellipseMajorOffset[RightSide]))
                      / (guideInfo.ellipseMajorAxis[RightSide] * guideInfo.ellipseMajorAxis[RightSide]))
        * guideInfo.ellipseMinorAxis[RightSide];

  //////////////////segmentation of m values

  // Are the arrays empty?
  if (mvaluesright != NULL || mvaluesleft != NULL || mvaluestop != NULL || mvaluesbottom != NULL) {
    guideInfo.enableSegments = 1;

    if (nSegments == -1) {
      printf ("\nError (%s): vector-specifcation of coating used, but nSegments=%i.\n  Please give provide nSegments = length of coating segment-arrays!\n",
              _comp->_name, nSegments);
      exit (-1);
    } else {
      guideInfo.numberOfSegments = (int)nSegments;
    }

    // printf("Length is %i\n",guideInfo.numberOfSegments);
    guideInfo.mValuesright = mvaluesright;
    guideInfo.mValuesleft = mvaluesleft;
    guideInfo.mValuestop = mvaluestop;
    guideInfo.mValuesbottom = mvaluesbottom;
    // printf("Seglength ... %f %f %f\n",seglength[0],seglength[1],seglength[2]);

    // Are the arrays of equal length?
    if (seglength == NULL) {
      dynamicalSegLength = realloc (dynamicalSegLength, guideInfo.numberOfSegments * sizeof (double));
      int i;
      for (i = 0; i < guideInfo.numberOfSegments; ++i) {
        dynamicalSegLength[i] = guideInfo.Length / guideInfo.numberOfSegments;
      }
      guideInfo.segLength = dynamicalSegLength;
    } else
      guideInfo.segLength = seglength;

    double sumOfelements = 0;
    int i;
    for (i = 0; i < guideInfo.numberOfSegments; i++) {
      sumOfelements += guideInfo.segLength[i];
    }
    if (guideInfo.verboseSetting && fabs (sumOfelements - guideInfo.Length) > 1e-9)
      printf ("Error in userinput inside %s, the difference between"
              " guidelength and elements of the seglength array is:"
              "%e consider changes the parameters l or seglength \n",
              NAME_CURRENT_COMP, sumOfelements - guideInfo.Length);
  } else
    guideInfo.enableSegments = 0;

  ///////////////////////////////////////////////////////////////////////////
  /////////////// Calculate gravity vector in the guides coordinatesystem
  ///////////////////////////////////////////////////////////////////////////

  /*
          Sets the local gravity vector equal to the global gravity vector (0,-g,0)
          and when apply the same rotation matrix as applied to guide.
  */
  if (enableGravity != 0) {
    Gx0 = 0, Gy0 = -GRAVITY * enableGravity, Gz0 = 0;
    Coords mcLocG;
    mcLocG = rot_apply (ROT_A_CURRENT_COMP, coords_set (0, Gy0, 0));
    coords_get (mcLocG, &Gx0, &Gy0, &Gz0);
  }
  Circ = 2 * PI * curvature;
%}

/**
	Finds the next collision between the particle and the guide.
	Decided by the type of collision, the particle will then be either
	reflected on the guide, leaving the guide, or absorbed.
	If reflected then the particle will have its velocity vector and p
	changed appropriable (TODO).
	In the case of a absorbed particle the ABSORB function is called.
	and if the particle is found to have left the guide the break command
	is called and the end of trace is reached
*/
TRACE
%{
  struct Intersection latestParticleCollision;
  latestParticleCollision.delta_time_to_next_collision = 0;
  latestParticleCollision.side = 0;
  latestParticleCollision.ApproxOn = 0;
  latestParticleCollision.collisionType = 0;

  PROP_Z0;

  double Gloc;
  double Gx, Gy, Gz;
  if (curvature) {
    Gloc = (vx * vx + vy * vy + vz * vz) / curvature;
  } else {
    Gloc = 0;
  }

  if (!guideInfo.EnclosingBoxOn)
    if (fabs (x) > guideInfo.entranceHorizontalWidth / 2.0 || fabs (y) > guideInfo.entranceVerticalWidth / 2.0)
      ABSORB;

  SCATTER;

  int bounces = 0;
  for (bounces = 0; bounces <= 1000; bounces++) {

    Gx = Gx0;
    Gy = Gy0;
    Gz = Gz0;
    if (curvature) {
      // Add velocity-dependent, location-dependent approximation to centripetal force for curvature...
      Gx = Gx0 + Gloc * cos (2 * PI * z / Circ);
      Gz = Gz0 + Gloc * sin (2 * PI * z / Circ);
    }

    // Find the next intersection between the guide and the neutron.
    int boolean = guide_elliptical_handleGuideIntersection (x, y, z, vx, vy, vz, Gx, Gy, Gz, &guideInfo, &latestParticleCollision);

    double timeToCollision = latestParticleCollision.delta_time_to_next_collision;

    // Handle special cases.
    if (boolean == 0)
      ABSORB;
    if (timeToCollision < 1e-15)
      ABSORB;

    // If the neutron reach the end of the guide, when move
    // the neutron to the end of guide and leave this component.
    if (z + vz * timeToCollision + 0.5 * Gz * timeToCollision * timeToCollision >= guideInfo.Length) {
      double timeToExit = 0;
      solve_2nd_order (&timeToExit, NULL, -0.5 * Gz, -vz, guideInfo.Length - z - 1e-9);
      PROP_GRAV_DT (timeToExit, Gx, Gy, Gz);
      SCATTER;
      break;
    }

    // Move the neutron and handle the reflection.
    PROP_GRAV_DT (timeToCollision, Gx, Gy, Gz);
    if (latestParticleCollision.collisionType == Absorb) {
      ABSORB;
    }
    if (latestParticleCollision.collisionType == Reflex) {
      p *= guide_elliptical_handleReflection (x, y, z, &vx, &vy, &vz, &guideInfo, &latestParticleCollision);
      SCATTER;
      if (p == 0)
        ABSORB;
    }
  }

  if (fabs (x) > guideInfo.exitHorizontalWidth / 2 || fabs (y) > guideInfo.exitVerticalWidth / 2)
    ABSORB;
%}


FINALLY
%{
%}

MCDISPLAY
%{


  // Calculate the points need to draw approximation of the ellipses
  // defining the guide

  // the number of lines used to draw one side of the guide
  int ApproximationMirrors = 500;

  // The start of the guide
  double zvalue = 0;

  // The the different in z between point used to draw the ellipse
  double zdelta = guideInfo.Length / (1.0 * ApproximationMirrors);

  // The vector used to store the points defining the lines
  double* xplus = malloc ((ApproximationMirrors + 1) * sizeof (double));
  double* xminus = malloc ((ApproximationMirrors + 1) * sizeof (double));
  double* yplus = malloc ((ApproximationMirrors + 1) * sizeof (double));
  double* yminus = malloc ((ApproximationMirrors + 1) * sizeof (double));

  if (!xplus || !xminus || !yplus || !yminus) {
    fprintf (stderr, "Elliptic_guide_gravity %s: ERROR allocating visualisation arrays. Exit!\n", NAME_CURRENT_COMP);
    exit (-1);
  }

  // Temperary values for the loop
  double tempx;
  double tempy;

  /*
          Calculate the second coordinates to the points on the ellipse with z_i
          as the first coordinate. We transform the point to the coordinate system
          where the ellipse is a unit circle. And use the definition of this circle
          to find second coordinate (x^2+z^2 = 1)
  */

  /////////////////////////////////////////////////////////
  double Length;
  double entranceHorizontalWidth;
  double entranceVerticalWidth;

  // ellipses infomation
  double ellipseMajorAxis[4], ellipseMinorAxis[4];
  double ellipseMajorOffset[4], ellipseMinorOffset[4];

  enum Side { RightSide, TopSide, LeftSide, BottomSide, None };
  /////////////////////////////////////////////////////////

  int i = 0;
  double tempz = 0;
  for (i = 0; i < ApproximationMirrors + 1; i++) {

    tempx = sqrt (guideInfo.ellipseMinorAxis[RightSide] * guideInfo.ellipseMinorAxis[RightSide]
                  - (guideInfo.ellipseMinorAxis[RightSide] * guideInfo.ellipseMinorAxis[RightSide])
                        / (guideInfo.ellipseMajorAxis[RightSide] * guideInfo.ellipseMajorAxis[RightSide])
                        * (zvalue + zdelta * i - guideInfo.ellipseMajorOffset[RightSide]) * (zvalue + zdelta * i - guideInfo.ellipseMajorOffset[RightSide]));

    xplus[i] = tempx + guideInfo.ellipseMinorOffset[RightSide];
    xminus[i] = -tempx + guideInfo.ellipseMinorOffset[RightSide];

    tempy = sqrt (guideInfo.ellipseMinorAxis[TopSide] * guideInfo.ellipseMinorAxis[TopSide]
                  - (guideInfo.ellipseMinorAxis[TopSide] * guideInfo.ellipseMinorAxis[TopSide])
                        / (guideInfo.ellipseMajorAxis[TopSide] * guideInfo.ellipseMajorAxis[TopSide])
                        * (zvalue + zdelta * i - guideInfo.ellipseMajorOffset[TopSide]) * (zvalue + zdelta * i - guideInfo.ellipseMajorOffset[TopSide]));

    yplus[i] = tempy + guideInfo.ellipseMinorOffset[TopSide];
    yminus[i] = -tempy + guideInfo.ellipseMinorOffset[TopSide];
  }

  ///// Draw lines

  // Drawing lines orthogonal with the z direction.
  // at both ends of the guide and at the boardest place at the guide

  // These may not give correct result if one of the ends are closed

  int j = 0;

  line (xplus[j], yplus[j], zvalue + j * zdelta, 0, yplus[j], zvalue + j * zdelta);
  line (0, yplus[j], zvalue + j * zdelta, xminus[j], yplus[j], zvalue + j * zdelta);
  line (xminus[j], yminus[j], zvalue + j * zdelta, 0, yminus[j], zvalue + j * zdelta);
  line (0, yminus[j], zvalue + j * zdelta, xplus[j], yminus[j], zvalue + j * zdelta);
  line (xminus[j], yplus[j], zvalue + j * zdelta, xminus[j], 0, zvalue + j * zdelta);
  line (xminus[j], 0, zvalue + j * zdelta, xminus[j], yminus[j], zvalue + j * zdelta);
  line (xplus[j], 0, zvalue + j * zdelta, xplus[j], yplus[j], zvalue + j * zdelta);
  line (xplus[j], yminus[j], zvalue + j * zdelta, xplus[j], 0, zvalue + j * zdelta);

  j = ApproximationMirrors;

  line (xplus[j], yplus[j], zvalue + j * zdelta, 0, yplus[j], zvalue + j * zdelta);
  line (0, yplus[j], zvalue + j * zdelta, xminus[j], yplus[j], zvalue + j * zdelta);
  line (xminus[j], yminus[j], zvalue + j * zdelta, 0, yminus[j], zvalue + j * zdelta);
  line (0, yminus[j], zvalue + j * zdelta, xplus[j], yminus[j], zvalue + j * zdelta);
  line (xminus[j], yplus[j], zvalue + j * zdelta, xminus[j], 0, zvalue + j * zdelta);
  line (xminus[j], 0, zvalue + j * zdelta, xminus[j], yminus[j], zvalue + j * zdelta);
  line (xplus[j], 0, zvalue + j * zdelta, xplus[j], yplus[j], zvalue + j * zdelta);
  line (xplus[j], yminus[j], zvalue + j * zdelta, xplus[j], 0, zvalue + j * zdelta);

  // find boardest place on the guide and draw a band around the guide
  int m0;
  double boardestPlace = 0;
  int boardestPlaceNumber = 0;
  for (m0 = 0; m0 < ApproximationMirrors; m0++) {
    if (boardestPlace <= fabs (yplus[m0])) {
      boardestPlace = fabs (yplus[m0]);
      boardestPlaceNumber = m0;
    }
  }
  j = boardestPlaceNumber;

  line (xplus[j], yplus[j], zvalue + j * zdelta, 0, yplus[j], zvalue + j * zdelta);
  line (0, yplus[j], zvalue + j * zdelta, xminus[j], yplus[j], zvalue + j * zdelta);
  line (xminus[j], yminus[j], zvalue + j * zdelta, 0, yminus[j], zvalue + j * zdelta);
  line (0, yminus[j], zvalue + j * zdelta, xplus[j], yminus[j], zvalue + j * zdelta);
  line (xminus[j], yplus[j], zvalue + j * zdelta, xminus[j], 0, zvalue + j * zdelta);
  line (xminus[j], 0, zvalue + j * zdelta, xminus[j], yminus[j], zvalue + j * zdelta);
  line (xplus[j], 0, zvalue + j * zdelta, xplus[j], yplus[j], zvalue + j * zdelta);
  line (xplus[j], yminus[j], zvalue + j * zdelta, xplus[j], 0, zvalue + j * zdelta);

  // Drawing lines parallel with the z direction

  int k = 0;
  for (k = 0; k < ApproximationMirrors; k++) {

    line (xplus[k], yplus[k], zvalue + k * zdelta, xplus[k + 1], yplus[k + 1], zvalue + (k + 1) * zdelta);
    line (xminus[k], yplus[k], zvalue + k * zdelta, xminus[k + 1], yplus[k + 1], zvalue + (k + 1) * zdelta);

    line (xplus[k], yminus[k], zvalue + k * zdelta, xplus[k + 1], yminus[k + 1], zvalue + (k + 1) * zdelta);

    line (xminus[k], yminus[k], zvalue + k * zdelta, xminus[k + 1], yminus[k + 1], zvalue + (k + 1) * zdelta);

    line (xminus[k], 0, zvalue + k * zdelta, xminus[k + 1], 0, zvalue + (k + 1) * zdelta);
    line (xplus[k], 0, zvalue + k * zdelta, xplus[k + 1], 0, zvalue + (k + 1) * zdelta);

    line (0, yminus[k], zvalue + k * zdelta, 0, yminus[k + 1], zvalue + (k + 1) * zdelta);
    line (0, yplus[k], zvalue + k * zdelta, 0, yplus[k], zvalue + (k + 1) * zdelta);
  }

  if (guideInfo.EnclosingBoxOn) {
    dashed_line (guideInfo.xArray[0], guideInfo.yArray[0], guideInfo.zArray[0], guideInfo.xArray[1], guideInfo.yArray[1], guideInfo.zArray[1], 10);
    dashed_line (guideInfo.xArray[1], guideInfo.yArray[1], guideInfo.zArray[1], guideInfo.xArray[2], guideInfo.yArray[2], guideInfo.zArray[2], 10);
    dashed_line (guideInfo.xArray[2], guideInfo.yArray[2], guideInfo.zArray[2], guideInfo.xArray[3], guideInfo.yArray[3], guideInfo.zArray[3], 10);
    dashed_line (guideInfo.xArray[3], guideInfo.yArray[3], guideInfo.zArray[3], guideInfo.xArray[0], guideInfo.yArray[0], guideInfo.zArray[0], 10);

    dashed_line (guideInfo.xArray[4], guideInfo.yArray[4], guideInfo.zArray[4], guideInfo.xArray[5], guideInfo.yArray[5], guideInfo.zArray[5], 10);
    dashed_line (guideInfo.xArray[5], guideInfo.yArray[5], guideInfo.zArray[5], guideInfo.xArray[6], guideInfo.yArray[6], guideInfo.zArray[6], 10);
    dashed_line (guideInfo.xArray[6], guideInfo.yArray[6], guideInfo.zArray[6], guideInfo.xArray[7], guideInfo.yArray[7], guideInfo.zArray[7], 10);
    dashed_line (guideInfo.xArray[7], guideInfo.yArray[7], guideInfo.zArray[7], guideInfo.xArray[4], guideInfo.yArray[4], guideInfo.zArray[4], 10);

    dashed_line (guideInfo.xArray[0], guideInfo.yArray[0], guideInfo.zArray[0], guideInfo.xArray[4], guideInfo.yArray[4], guideInfo.zArray[4], 10);
    dashed_line (guideInfo.xArray[4], guideInfo.yArray[4], guideInfo.zArray[4], guideInfo.xArray[7], guideInfo.yArray[7], guideInfo.zArray[7], 10);
    dashed_line (guideInfo.xArray[7], guideInfo.yArray[7], guideInfo.zArray[7], guideInfo.xArray[3], guideInfo.yArray[3], guideInfo.zArray[3], 10);
    dashed_line (guideInfo.xArray[3], guideInfo.yArray[3], guideInfo.zArray[3], guideInfo.xArray[0], guideInfo.yArray[0], guideInfo.zArray[0], 10);

    dashed_line (guideInfo.xArray[1], guideInfo.yArray[1], guideInfo.zArray[1], guideInfo.xArray[5], guideInfo.yArray[5], guideInfo.zArray[5], 10);
    dashed_line (guideInfo.xArray[5], guideInfo.yArray[5], guideInfo.zArray[5], guideInfo.xArray[6], guideInfo.yArray[6], guideInfo.zArray[6], 10);
    dashed_line (guideInfo.xArray[6], guideInfo.yArray[6], guideInfo.zArray[6], guideInfo.xArray[2], guideInfo.yArray[2], guideInfo.zArray[2], 10);
    dashed_line (guideInfo.xArray[2], guideInfo.yArray[2], guideInfo.zArray[2], guideInfo.xArray[1], guideInfo.yArray[1], guideInfo.zArray[1], 10);
  }
  free (xminus);
  free (yminus);
  free (xplus);
  free (yplus);
%}

END