1725 lines
63 KiB
Plaintext
1725 lines
63 KiB
Plaintext
/*******************************************************************************
|
|
*
|
|
* 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
|
|
* Version: 0.91
|
|
* Origin: NBI
|
|
* Release: McStas 2.1
|
|
*
|
|
* 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.
|
|
*
|
|
* %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
|
|
* he 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.
|
|
*
|
|
* OUTPUT PARAMETERS
|
|
* guideInfo: (1) structure containing internal variables
|
|
* Gx: x-component of (m/s^2) local gravity vector
|
|
* Gy: y-component of (m/s^2) local gravity vector
|
|
* Gz: z-component of (m/s^2) local gravity vector
|
|
* dynamicalSegLength: (1) automatic generate length of specific m-value segments
|
|
* if no seglength is given
|
|
*
|
|
* %D
|
|
*
|
|
* Example 1, Elliptical definition using focal points:
|
|
* 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
|
|
* )
|
|
*
|
|
* Example 2: Half elliptical definition:
|
|
* 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"
|
|
* )
|
|
*
|
|
* Example 3: Parabolic approximation:
|
|
* 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"
|
|
* )
|
|
*
|
|
* Example 4: Elliptical definition with varying m-values:
|
|
* 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];
|
|
*
|
|
* %E
|
|
*******************************************************************************/
|
|
|
|
DEFINE COMPONENT Elliptic_guide_gravity_custom
|
|
DEFINITION PARAMETERS(
|
|
mvaluesright=NULL, mvaluesleft=NULL, mvaluestop=NULL, mvaluesbottom=NULL,
|
|
seglength=NULL )
|
|
SETTING PARAMETERS (l,
|
|
xwidth = 0,yheight = 0,
|
|
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", // on or off
|
|
enableGravity = 1.0,
|
|
curvature=0 )
|
|
|
|
OUTPUT PARAMETERS (
|
|
guideInfo,latestParticleCollision,Gx,Gy,Gz,
|
|
Gx0, Gy0, Gz0,Circ,dynamicalSegLength
|
|
)
|
|
|
|
SHARE
|
|
%{
|
|
// m_local_refl is to be used in shielding applications and accessed by shielding logger
|
|
#ifndef MVALUE_LOCAL_IS_DEF
|
|
#define MVALUE_LOCAL_IS_DEF 1
|
|
double m_local_refl=-1;
|
|
#endif
|
|
|
|
%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]
|
|
*/
|
|
|
|
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;
|
|
}
|
|
|
|
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; }
|
|
}
|
|
|
|
|
|
|
|
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;
|
|
}
|
|
|
|
|
|
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;
|
|
}
|
|
|
|
/**
|
|
Check if the neutron is within the guide using the sign
|
|
of the crossproduct between the two points,
|
|
on each of the enclosing box surface and neutrons position.
|
|
|
|
@param x,y,z; position of the neutron. [m]
|
|
@param guideInfo; pointer to the guide infomation holding structure.
|
|
@return; return 1 if the neutron is inside the guide [boolean]
|
|
*/
|
|
|
|
/*
|
|
int guide_elliptical_InsideEnclosingBox(double x,double y,double z,struct SGI *guideInfo){
|
|
int guide_elliptical_IsPointInVolume(
|
|
double *x,double *y,double *z,
|
|
double px,double py,double pz){
|
|
int guide_elliptical_WhichSide( double p1x,double p1y,double p1z,
|
|
double p2x,double p2y,double p2z,
|
|
double p3x,double p3y,double p3z,
|
|
double px ,double py ,double pz ){
|
|
|
|
double v1x = p1x - p2x, v1y = p1y-p2y, v1z = p1z-p2z;
|
|
double v2x = p3x - p2x, v2y = p3y-p2y, v2z = p3z-p2z;
|
|
double v3x = v2y*v1z-v2z*v1y;
|
|
double v3y = v2z*v1x-v2x*v1z;
|
|
double v3z = v2x*v1y-v2y*v1x;
|
|
|
|
return 0 >= v3x*(px-p1x)+v3y*(py-p1y)+v3z*(pz-p1z);
|
|
}
|
|
|
|
if( //front
|
|
guide_elliptical_WhichSide(x[3],y[3],z[3],x[2],y[2],z[2],x[1],y[1],z[1],px,py,pz) &&
|
|
guide_elliptical_WhichSide(x[1],y[1],z[1],x[0],y[0],z[0],x[3],y[3],z[3],px,py,pz) &&
|
|
//back
|
|
guide_elliptical_WhichSide(x[5],y[5],z[5],x[6],y[6],z[6],x[7],y[7],z[7],px,py,pz) &&
|
|
guide_elliptical_WhichSide(x[7],y[7],z[7],x[4],y[4],z[4],x[5],y[5],z[5],px,py,pz) &&
|
|
//right
|
|
guide_elliptical_WhichSide(x[7],y[7],z[7],x[3],y[3],z[3],x[0],y[0],z[0],px,py,pz) &&
|
|
guide_elliptical_WhichSide(x[0],y[0],z[0],x[4],y[4],z[4],x[7],y[7],z[7],px,py,pz) &&
|
|
//left
|
|
guide_elliptical_WhichSide(x[1],y[1],z[1],x[2],y[2],z[2],x[6],y[6],z[6],px,py,pz) &&
|
|
guide_elliptical_WhichSide(x[6],y[6],z[6],x[5],y[5],z[5],x[1],y[1],z[1],px,py,pz) &&
|
|
//top
|
|
guide_elliptical_WhichSide(x[0],y[0],z[0],x[1],y[1],z[1],x[5],y[5],z[5],px,py,pz) &&
|
|
guide_elliptical_WhichSide(x[5],y[5],z[5],x[4],y[4],z[4],x[0],y[0],z[0],px,py,pz) &&
|
|
//bottom
|
|
guide_elliptical_WhichSide(x[6],y[6],z[6],x[2],y[2],z[2],x[3],y[3],z[3],px,py,pz) &&
|
|
guide_elliptical_WhichSide(x[3],y[3],z[3],x[7],y[7],z[7],x[6],y[6],z[6],px,py,pz) )
|
|
return 1;
|
|
else return 0;
|
|
}
|
|
return guide_elliptical_IsPointInVolume(
|
|
guideInfo->xArray,guideInfo->yArray,guideInfo->zArray,x,y,z);
|
|
}
|
|
*/
|
|
|
|
|
|
///////////////////////////////////////////////////////////////////////////
|
|
/////////////// 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
|
|
*/
|
|
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.
|
|
*/
|
|
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;
|
|
double combinedLength = 0;
|
|
int i;
|
|
for(i=0; i < sgi->numberOfSegments; i++){
|
|
combinedLength = combinedLength + sgi->segLength[i];
|
|
if(z0 < combinedLength) {
|
|
currentSegment = i; break;
|
|
}
|
|
}
|
|
|
|
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
|
|
struct Intersection latestParticleCollision; // Is changed duing trace
|
|
double Gx,Gy,Gz; // Local gravity vector, is set once in INITIALIZE
|
|
double Gx0, Gy0, 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.
|
|
|
|
|
|
// Sets the mirror type of the guides mirrors
|
|
// These variables are used in the collision functions
|
|
// to find the type of collision
|
|
|
|
// guideInfo.OuterSide[RightSide] =
|
|
// guide_elliptical_getMirrorTypeFromInput(outer_right_side_mirror,guideInfo.verboseSetting);
|
|
// guideInfo.OuterSide[TopSide] =
|
|
// guide_elliptical_getMirrorTypeFromInput(outer_top_side_mirror,guideInfo.verboseSetting);
|
|
// guideInfo.OuterSide[LeftSide] =
|
|
// guide_elliptical_getMirrorTypeFromInput(outer_left_side_mirror,guideInfo.verboseSetting);
|
|
// guideInfo.OuterSide[BottomSide] =
|
|
// guide_elliptical_getMirrorTypeFromInput(outer_bottom_side_mirror,guideInfo.verboseSetting);
|
|
|
|
// guideInfo.InnerSide[RightSide] =
|
|
// guide_elliptical_getMirrorTypeFromInput(inner_right_side_mirror,guideInfo.verboseSetting);
|
|
// guideInfo.InnerSide[TopSide] =
|
|
// guide_elliptical_getMirrorTypeFromInput(inner_top_side_mirror,guideInfo.verboseSetting);
|
|
// guideInfo.InnerSide[LeftSide] =
|
|
// guide_elliptical_getMirrorTypeFromInput(inner_left_side_mirror,guideInfo.verboseSetting);
|
|
// guideInfo.InnerSide[BottomSide] =
|
|
// guide_elliptical_getMirrorTypeFromInput(inner_bottom_side_mirror,guideInfo.verboseSetting);
|
|
|
|
// 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(directDefination == 0){ */
|
|
/* guideInfo.entranceHorizontalWidth = xwidth; */
|
|
/* guideInfo.entranceVerticalWidth = yheight; */
|
|
/* } */
|
|
|
|
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;
|
|
|
|
/*
|
|
double DefaultArray1[8] = { 1.0,-1.0,-1.0, 1.0, 1.0,-1.0,-1.0, 1.0};
|
|
double DefaultArray2[8] = { 1.0, 1.0,-1.0,-1.0, 1.0, 1.0,-1.0,-1.0};
|
|
double DefaultArray3[8] = {-1.0,-1.0,-1.0,-1.0, 1.0, 1.0, 1.0, 1.0};
|
|
|
|
guideInfo.EnclosingBoxOn = 0;
|
|
double *xinput;
|
|
if ( xInput != NULL ){ xinput = xInput; guideInfo.EnclosingBoxOn = 1; }
|
|
else { xinput = DefaultArray1; }
|
|
double *yinput;
|
|
if ( yInput != NULL ){ yinput = yInput; guideInfo.EnclosingBoxOn = 1;}
|
|
else { yinput = DefaultArray2; }
|
|
double *zinput;
|
|
if ( zInput != NULL ){ zinput = zInput; guideInfo.EnclosingBoxOn = 1;}
|
|
else { zinput = DefaultArray3; }
|
|
*/
|
|
|
|
/*
|
|
double xarray[8] ={ guideInfo.ellipseMinorAxis[0]*xinput[0],
|
|
guideInfo.ellipseMinorAxis[2]*xinput[1],
|
|
guideInfo.ellipseMinorAxis[2]*xinput[2],
|
|
guideInfo.ellipseMinorAxis[0]*xinput[3],
|
|
guideInfo.ellipseMinorAxis[0]*xinput[4],
|
|
guideInfo.ellipseMinorAxis[2]*xinput[5],
|
|
guideInfo.ellipseMinorAxis[2]*xinput[6],
|
|
guideInfo.ellipseMinorAxis[0]*xinput[7] };
|
|
double yarray[8] ={ guideInfo.ellipseMinorAxis[1]*yinput[0],
|
|
guideInfo.ellipseMinorAxis[1]*yinput[1],
|
|
guideInfo.ellipseMinorAxis[3]*yinput[2],
|
|
guideInfo.ellipseMinorAxis[3]*yinput[3],
|
|
guideInfo.ellipseMinorAxis[1]*yinput[4],
|
|
guideInfo.ellipseMinorAxis[1]*yinput[5],
|
|
guideInfo.ellipseMinorAxis[3]*yinput[6],
|
|
guideInfo.ellipseMinorAxis[3]*yinput[7] };
|
|
double zarray[8] ={ guideInfo.Length/2*zinput[0]+guideInfo.Length/2,
|
|
guideInfo.Length/2*zinput[1]+guideInfo.Length/2,
|
|
guideInfo.Length/2*zinput[2]+guideInfo.Length/2,
|
|
guideInfo.Length/2*zinput[3]+guideInfo.Length/2,
|
|
guideInfo.Length/2*zinput[4]+guideInfo.Length/2,
|
|
guideInfo.Length/2*zinput[5]+guideInfo.Length/2,
|
|
guideInfo.Length/2*zinput[6]+guideInfo.Length/2,
|
|
guideInfo.Length/2*zinput[7]+guideInfo.Length/2 };
|
|
int i = 0;
|
|
for(i = 0; i < 8; i++){
|
|
guideInfo.xArray[i] = xarray[i];
|
|
guideInfo.yArray[i] = yarray[i];
|
|
guideInfo.zArray[i] = zarray[i];
|
|
}
|
|
*/
|
|
|
|
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;
|
|
|
|
guideInfo.numberOfSegments = sizeof(mvaluesright)/sizeof(mvaluesright[0]);
|
|
|
|
//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;
|
|
|
|
if( guideInfo.numberOfSegments != sizeof(mvaluesright)/sizeof(mvaluesright[0])
|
|
|| guideInfo.numberOfSegments != sizeof(mvaluesleft)/sizeof(mvaluesleft[0])
|
|
|| guideInfo.numberOfSegments != sizeof(mvaluestop)/sizeof(mvaluestop[0])
|
|
|| guideInfo.numberOfSegments != sizeof(mvaluesbottom)/sizeof(mvaluesbottom[0])
|
|
|| (guideInfo.segLength == NULL
|
|
& guideInfo.numberOfSegments != sizeof(seglength)/sizeof(guideInfo.segLength[0])
|
|
) ) {
|
|
|
|
printf("Error in userinput inside %s, the length of the arrays"
|
|
" mvalues and seglength are not equal\n",NAME_CURRENT_COMP);
|
|
printf("The length of the arrays are: mValuesright is %lu,"
|
|
" mvaluesleft is %lu, mvaluestop is %lu, mvaluesbottom is"
|
|
" %lu and seglength is %lu and should be %d \n; Above assume that the arrays are using double \n",
|
|
sizeof(mvaluesright)/sizeof(double),
|
|
sizeof(mvaluesleft)/sizeof(double),
|
|
sizeof(mvaluestop)/sizeof(double),
|
|
sizeof(mvaluesbottom)/sizeof(double),
|
|
sizeof(guideInfo.segLength)/sizeof(double),
|
|
guideInfo.numberOfSegments
|
|
);
|
|
|
|
if ( guideInfo.verboseSetting ) {
|
|
int i;
|
|
|
|
printf("The Values of mvaluesright is: [");
|
|
for(i=0; i < sizeof(mvaluesright)/sizeof(mvaluesright[0]); i++)
|
|
printf("%e,",guideInfo.mValuesright[i] );
|
|
printf("]\n");
|
|
|
|
printf("The Values of mvaluesleft is: [");
|
|
for(i=0; i < sizeof(mvaluesleft)/sizeof(mvaluesleft[0]); i++)
|
|
printf("%e,",guideInfo.mValuesleft[i] );
|
|
printf("]\n");
|
|
|
|
printf("The Values of mvaluestop is: [");
|
|
for(i=0; i < sizeof(mvaluestop)/sizeof(mvaluestop[0]); i++)
|
|
printf("%e,",guideInfo.mValuestop[i] );
|
|
printf("]\n");
|
|
|
|
printf("The Values of mvaluesbottom is: [");
|
|
for(i=0; i < sizeof(mvaluesbottom)/sizeof(mvaluesbottom[0]); i++)
|
|
printf("%e,",guideInfo.mValuesbottom[i] );
|
|
printf("]\n");
|
|
|
|
printf("The Values of seglength is: [");
|
|
for(i=0; i < sizeof(guideInfo.segLength)/sizeof(guideInfo.segLength[0]); i++)
|
|
printf("%e,",guideInfo.segLength[i]);
|
|
printf("]\n");
|
|
}
|
|
exit( printf("Exit due to critical error in userinput for the"
|
|
" component %s, consider having a look at the input"
|
|
" for following: mvaluesright,mvaluesleft,mvaluestop,"
|
|
"mvaluesbottom and/or seglength.",NAME_CURRENT_COMP) );
|
|
}
|
|
//
|
|
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
|
|
%{
|
|
PROP_Z0;
|
|
m_local_refl=-1; // RK - added. Comment: neutron entered the guide and hasn't reached the coating yet.
|
|
SCATTER;
|
|
|
|
double Gloc;
|
|
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 ) { m_local_refl=-1; // RK - added.
|
|
// RK - Comment: neutron doesn't make it to the guide opening, doesn't reach coating, setting m_local_refl=-1
|
|
ABSORB;}
|
|
|
|
|
|
int bounces = 0;
|
|
for(bounces = 0; bounces <= 1000; bounces++){
|
|
//RK - Comment: now have entered the guide, will assign m_local_refl a value of the coating at the place of the hit.
|
|
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.
|
|
//RK -- Comment: this happens when neutron hits precisely the corner of the guide. Assume no coating there
|
|
if(boolean == 0) { m_local_refl=-1; ABSORB;} //if(boolean == 0) ABSORB; // RK -- modified line
|
|
//RK -- Comment: this happens when neutron is travelling parallel to the wall. Assume absorption in the upper coating layer which amounts to setting m=1
|
|
if(timeToCollision < 1e-15) { m_local_refl = 1; ABSORB;} //if(timeToCollision < 1e-15) ABSORB; // RK -- modified line
|
|
|
|
// 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);
|
|
m_local_refl=-1; // RK - added. Comment: exiting guide, no scattering on the coating
|
|
SCATTER;
|
|
break;
|
|
}
|
|
|
|
// Move the neutron and handle the reflection.
|
|
PROP_GRAV_DT(timeToCollision,Gx,Gy,Gz);
|
|
/*** RK now we are at the point of reflection, have to figure out what the m-value is here. ***/
|
|
if(!guideInfo.enableSegments){
|
|
m_local_refl=guideInfo.mArr[latestParticleCollision.side];
|
|
}
|
|
else{
|
|
int currentSegment;
|
|
double combinedLength = 0;
|
|
int i;
|
|
for(i=0; i < guideInfo.numberOfSegments; i++){
|
|
combinedLength = combinedLength + guideInfo.segLength[i];
|
|
if(z < combinedLength) {
|
|
currentSegment = i; break;
|
|
}
|
|
}
|
|
|
|
if(latestParticleCollision.side == RightSide)
|
|
m_local_refl = guideInfo.mValuesright[currentSegment];
|
|
if(latestParticleCollision.side == LeftSide)
|
|
m_local_refl=guideInfo.mValuesleft[currentSegment];
|
|
if(latestParticleCollision.side == TopSide)
|
|
m_local_refl = guideInfo.mValuestop[currentSegment];
|
|
if(latestParticleCollision.side == BottomSide)
|
|
m_local_refl = guideInfo.mValuesbottom[currentSegment];
|
|
}
|
|
/*****/
|
|
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;
|
|
}
|
|
}
|
|
// RK -- this happens when neutron didn't exit after 1000 bounces which is unrealistic situation. m_local_refl remains set to the m-value at last revlection.
|
|
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[ApproximationMirrors+1];
|
|
double xminus[ApproximationMirrors+1];
|
|
double yplus[ApproximationMirrors+1];
|
|
double yminus[ApproximationMirrors+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
|
|
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 );
|
|
}
|
|
%}
|
|
|
|
END
|