/*******************************************************************************
*
* 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 allow 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 perfect elliptic guide is centered along the z-axis with the entrance
* and exit in the xy-plane. The horizontal and vertical ellipses defining
* the guide geometry is by default set by two focal points.
* These are placed a distance away from the guide openings along the z-axis;
* if distance given is positive, when the focal point is outside the guide.
*
* Multiple options for defining these ellipse exist including approximation of
* parabolas and half ellipses (mid point of the ellipse or one of the guide openings)
*
* The guide coating parameters can be set for each side of the guide.
* Furthermore the m-value can be specified for user defined segments
* of the guide.
*
* <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 defination of the guide geometry, will ignore w,h lin and lout parameters if this is nonzero. Length of the axis parallel to the z for the horizontal ellipse
* minorAxisxw: [m]          direct defination of the guide geometry, will ignore w,h lin and lout parameters if this is nonzero. Length of the axis Perpendicular to the z for the horizontal ellipse
* majorAxisyh: [m]          direct defination of the guide geometry, will ignore w,h lin and lout parameters if this is nonzero. Length of the axis parallel to the z for the vertical ellipse
* minorAxisyh: [m]          direct defination of the guide geometry, will ignore w,h lin and lout parameters if this is nonzero. Length of the axis Perpendicular to the z for the vertical ellipse
* majorAxisoffsetxw: [m]    direct defination of the guide geometry, distance between the center of the horizontal ellipse and the guide entrance
* majorAxisoffsetyh: [m]    direct defination 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,
	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 an 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;
	// 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 intern guide values from user input
	///////////////////////////////////////////////////////////////////////////

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

	if( majorAxisyh != 0 || minorAxisyh != 0
	 || majorAxisxw != 0 || minorAxisxw != 0)
	{
		directDefination = 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 && directDefination == 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 && directDefination == 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 && directDefination == 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: vector-specifcation of coating used, but nSegments=%i.\n  Please give provide nSegments = length of coating segment-arrays!\n", 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;
  SCATTER;

  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;


	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));

	// 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
		there the ellipse is a unit circle. And use the defination 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