1/******************************************************************************* 2* 3* McStas, neutron ray-tracing package 4* Copyright (C) 1997-2011, All rights reserved 5* Risoe National Laboratory, Roskilde, Denmark 6* Institut Laue Langevin, Grenoble, France 7* 8* Component: Elliptic_guide_gravity 9* 10* %I 11* Written by: Henrik Bo Hoffmann Carlsen and Mads Bertelsen 12* Date: 27 Aug 2012 13* Origin: NBI 14* 15* Perfect elliptic guide which allow for simulations with gravity. 16* The guide mirrors can be divided into segments with individual m-values. 17* Parabolic guide components can also be simulated. 18* 19* %D 20* 21* The perfect elliptic guide is centered along the z-axis with the entrance 22* and exit in the xy-plane. The horizontal and vertical ellipses defining 23* the guide geometry is by default set by two focal points. 24* These are placed a distance away from the guide openings along the z-axis; 25* if distance given is positive, when the focal point is outside the guide. 26* 27* Multiple options for defining these ellipse exist including approximation of 28* parabolas and half ellipses (mid point of the ellipse or one of the guide openings) 29* 30* The guide coating parameters can be set for each side of the guide. 31* Furthermore the m-value can be specified for user defined segments 32* of the guide. 33* 34* %P 35* INPUT PARAMETERS: 36* mvaluesright: [pointer] Pointer to array of m-values, right mirror 37* mvaluesleft: [pointer] - same, left mirror 38* mvaluestop: [pointer] - same, top mirror 39* mvaluesbottom: [pointer] - same, bottom mirror 40* seglength: [pointer] Pointer to array of segment lengths for discrete mirror description 41* l: [m] length of the guide 42* 43* linxw: [m] distance from 1st focal point to guide entrance 44* - left and right horizontal mirrors 45* loutxw: [m] distance from 2nd focal point to guide exit 46* - left and right horizontal mirrors 47* linyh: [m] distance from 1st focal point to guide entrance 48* - top and bottom vertical mirrors 49* loutyh: [m] distance from 2nd focal point to guide exit 50* - top and bottom vertical mirrors 51* xwidth: [m] width at the guide entry, mid or exit (see dimensionsAt) 52* yheight: [m] height at the guide entry, mid or exit (see dimensionsAt) 53* R0: [1] Low-angle reflectivity 54* Qc: [AA-1] Critical scattering vector 55* alpha: [AA] Slope of reflectivity 56* m: [1] m-value of material for all mirrors, 57* zero means complete absorption. 58* W: [AA-1] Width of supermirror cut-off 59* 60* alphatop: [AA] Slope of reflectivity for top horizontal mirror, 61* overwrites alpha 62* mtop: [1] m-value of material for top horizontal mirror, 63* overwrites m 64* 65* alphabottom: [AA] Slope of reflectivity for bottom horizontal mirror 66* mbottom: [1] m-value of material for bottom horizontal mirror 67* 68* alpharight: [AA] Slope of reflectivity for right vertical mirror 69* mright: [1] m-value of material for right vertical mirror 70* 71* alphaleft: [AA] Slope of reflectivity for left vertical mirror 72* mleft: [1] m-value of material for left vertical mirror 73* 74* option: [string] options are 'ellipse' and 'halfEllipse'. Ellipse is defined by 75* both the focal points, while halfEllipse locked the center 76* of the ellipse either the entrance or exit of the guide, and 77* use the focal point of the other end to define the ellipse 78* 79* dimensionsAt: [string] define whether xwidth and yheight sets the size of 80* he opening, minor axis or the end of the guide. 81* 82* majorAxisxw: [m] direct defination of the guide geometry, will ignore 83* w,h lin and lout parameters if this is nonzero. Length of 84* the axis parallel to the z for the horizontal ellipse 85* minorAxisxw: [m] direct defination of the guide geometry, will ignore 86* w,h lin and lout parameters if this is nonzero. Length of 87* the axis Perpendicular to the z for the horizontal ellipse 88* majorAxisyh: [m] direct defination of the guide geometry, will ignore 89* w,h lin and lout parameters if this is nonzero. Length of 90* the axis parallel to the z for the vertical ellipse 91* minorAxisyh: [m] direct defination of the guide geometry, will ignore 92* w,h lin and lout parameters if this is nonzero. Length of 93* the axis Perpendicular to the z for the vertical ellipse 94* 95* majorAxisoffsetxw: [m] direct defination of the guide geometry, 96* distance between the center of the horizontal 97* ellipse and the guide entrance 98* majorAxisoffsetyh: [m] direct defination of the guide geometry, 99* distance between the center of the vertical 100* ellipse and the guide entrance 101* verbose: [bool] Give extra information about calculations 102* 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. 103* 104* OUTPUT PARAMETERS 105* guideInfo: [1] structure containing internal variables 106* Gx: [] x-component of (m/s^2) local gravity vector 107* Gy: [] y-component of (m/s^2) local gravity vector 108* Gz: [] z-component of (m/s^2) local gravity vector 109* dynamicalSegLength: [1] automatic generate length of specific m-value segments 110* if no seglength is given 111* 112* %D 113* 114* Example 1, Elliptical definition using focal points: 115* Elliptic_guide_gravity( 116* l=50, 117* linxw=5,linyh=5,loutxw=10,loutyh=10, 118* xwidth=0.05,yheight=0.05, 119* R0=0.99,Qc=0.0219,alpha=6.07,m=1.0,W=0.003 120* ) 121* 122* Example 2: Half elliptical definition: 123* Elliptic_guide_gravity( 124* l=50, 125* linxw=5,linyh=5,loutxw=10,loutyh=10, 126* xwidth=0.1,yheight=0.1, 127* R0=0.99,Qc=0.0219,alpha=6.07,m=1.0,W=0.003, 128* option = "halfEllipse", 129* dimensionsAt = "entrance" 130* ) 131* 132* Example 3: Parabolic approximation: 133* Elliptic_guide_gravity( 134* l=50, 135* linxw=5,linyh=5,loutxw=1e6,loutyh=1e6, // values larger than 1e8 may cause erroneous results 136* xwidth=0.1,yheight=0.1, 137* R0 = 0.99,Qc=0.0219,alpha=6.07,m=1.0,W=0.003, 138* dimensionsAt = "exit" 139* ) 140* 141* Example 4: Elliptical definition with varying m-values: 142* Elliptic_guide_gravity( 143* l=50, 144* linxw=5,linyh=5,loutxw=10,loutyh=10, 145* xwidth=0.1,yheight=0.1, 146* R0 = 0.99,Qc=0.0219,alpha=6.07,m=1.0,W=0.003, 147* mvaluesright=marray,mvaluesleft=marray,mvaluestop=marray,mvaluesbottom=marray 148* ) 149* where marray is initialized as 150* for(iter=0; iter < 50; iter++){ marray[iter] = 2; } 151* And Declared as 152* double mValues[50]; 153* 154* %E 155*******************************************************************************/ 156 157DEFINE COMPONENT Elliptic_guide_gravity 158DEFINITION PARAMETERS( 159mvaluesright=NULL, mvaluesleft=NULL, mvaluestop=NULL, mvaluesbottom=NULL, 160seglength=NULL ) 161SETTING PARAMETERS (l, 162xwidth = 0,yheight = 0, 163linxw = 0,loutxw = 0, 164linyh = 0,loutyh = 0, 165majorAxisxw = 0,minorAxisxw = 0, 166majorAxisyh = 0,minorAxisyh = 0, 167majorAxisoffsetxw = 0, 168majorAxisoffsetyh = 0, 169string dimensionsAt = "entrance", 170string option = "ellipse", 171R0=0.99, Qc=0.0218, alpha=6.07, m=2, W=0.003, 172alpharight=-1, mright=-1, 173alphaleft=-1, mleft=-1, 174alphatop=-1, mtop=-1, 175alphabottom=-1, mbottom=-1, 176string verbose = "on", // on or off 177enableGravity = 1.0, 178curvature=0 ) 179 180OUTPUT PARAMETERS ( 181guideInfo,latestParticleCollision,Gx,Gy,Gz, 182Gx0, Gy0, Gz0,Circ,dynamicalSegLength 183) 184 185SHARE 186%{ 187%include "ref-lib" 188 189/////////////////////////////////////////////////////////////////////////// 190/////////////// local structs and enums 191/////////////////////////////////////////////////////////////////////////// 192 193/** 194Sides of the guide 195*/ 196enum Side {RightSide,TopSide,LeftSide,BottomSide,None}; 197 198/** 199The type of the collision is set in the collision function 200and decide the functions called in trace() 201 Reflex (TODO change this name) calls the reflection function 202 Absorb calls the built in ABSORB funtion. 203 LeaveGuide calls break and end the calculations in this component 204 EnterGuide does nothing 205*/ 206enum CollisionType {Reflex,Absorb,LeaveGuide,EnterGuide}; 207 208/** 209 The Mirror type sets the CollisionType of particles colliding on the mirror 210*/ 211enum MirrorType {MirrorTypeReflection,MirrorTypeTransparent,MirrorTypeabsorption}; 212 213// enum IntersectionType {Reflex,Absorb,Transparent,Leave,Enter}; 214 215/** 216 Collision between guide and the particle 217 contain infomation on the time to the next collision, 218 which side of the guide it is on and whether this part of the guide 219 is a perfect or approximated ellipse. 220*/ 221struct Intersection 222{ 223 double delta_time_to_next_collision; 224 enum Side side; // A number from 0 to 4 (4 being an error warning) 225 int ApproxOn; 226 enum CollisionType collisionType; 227}; 228 229/** 230 Static Guide information (SGI) 231 contain information on the guide, the ellipses and the mirrors on all sides 232*/ 233struct SGI 234{ 235 // guide infomation 236 double Length; 237 double entranceHorizontalWidth, entranceVerticalWidth; 238 double exitHorizontalWidth, exitVerticalWidth; 239 240 // ellipses infomation 241 double ellipseMajorAxis[4], ellipseMinorAxis[4]; 242 double ellipseMajorOffset[4], ellipseMinorOffset[4]; 243 244 // mirror infomation 245 double R0Arr[4]; 246 double QcArr[4]; 247 double alphaArr[4]; 248 double mArr[4]; 249 double WArr[4]; 250 251 // mirror type 252 enum MirrorType InnerSide[4]; 253 enum MirrorType OuterSide[4]; 254 255 // selene 256 int EnclosingBoxOn; 257 double xArray[8]; 258 double yArray[8]; 259 double zArray[8]; 260 261 // segmentation 262 int numberOfSegments; 263 int enableSegments; 264 double *mValuesright; 265 double *mValuesleft; 266 double *mValuestop; 267 double *mValuesbottom; 268 double *segLength; 269 270 int verboseSetting; 271}; 272 273 274/////////////////////////////////////////////////////////////////////////// 275/////////////// Error Handling Functions 276/////////////////////////////////////////////////////////////////////////// 277 278/** 279 If a user input is less than zero and hence doesn't allow for a well 280 define geomtric of the guide or physical values for mirrors 281 @param var is the input varible there the error occurred [text] 282*/ 283int guide_elliptical_illegalInputLessThanZero(char var[],int verbose){ 284 if (verbose) 285 printf("The user defined variable %s in %s has an illegal value" 286 " less than zero\n",var,"Elliptic_guide_gravity"); 287 return 1; 288} 289 290/** 291 The first focal point is in and the second is out. 292 If -in-out > L then they would change position as the 293 first and second focal points. This is 294 @param in,out is the input varible there the error occurred [text] 295*/ 296int guide_elliptical_illegalInputFocalPointsHyperbola( 297 char in[],char out[], 298 double inValue,double outValue, int verbose){ 299 if (verbose){ 300 printf("The user defined length of the guide, length \ 301 and the focal points %s and %s does not result \ 302 in an well defined ellipse. swap the focal points \ 303 or increase L, %s or %s to fix this problem\n", 304 in,out,in,out); 305 printf("The mininum length of the should be around %e\n", 306 inValue+outValue+0.000001); 307 } 308 return 1; 309} 310 311/** 312 Gives a warning if a part of the code is called that 313 should not be accessible if the algoritmes are working correctly 314 Most likely errors are floating points and ill-defined cases 315*/ 316void guide_elliptical_callCriticalWarning(char func[],int verbose){ 317 if (verbose) 318 printf("A CRITICAL WARNING has been called inside %s by function %s." 319 "This is most likely due to a programming error \ 320 inside the component. \n", 321 "Elliptic_guide_gravity",func); 322 } 323 324/////////////////////////////////////////////////////////////////////////// 325/////////////// Collision handling functions 326/////////////////////////////////////////////////////////////////////////// 327 328int guide_elliptical_getMirrorTypeFromInput(char * input,int verbose){ 329 int type = -1; 330 char* r1 = "reflection"; char* r2 = "reflect"; char* r3 = "r"; 331 char* a1 = "absorption"; char* a2 = "absorb"; char* a3 = "a"; 332 char* t1 = "transparant";char* t2 = "trans"; char* t3 = "t"; 333 if (strcmp (input, r1) == 0 334 || strcmp (input, r2) == 0 335 || strcmp (input, r3) == 0) 336 type = MirrorTypeReflection; 337 if (strcmp (input, a1) == 0 338 || strcmp (input, a2) == 0 339 || strcmp (input, a3) == 0) 340 type = MirrorTypeabsorption; 341 if (strcmp (input, t1) == 0 342 || strcmp (input, t2) == 0 343 || strcmp (input, t3) == 0) 344 type = MirrorTypeTransparent; 345 if ( type == -1 && verbose) 346 printf( "Following string is not a valid type of a mirror: %s," 347 "use reflection,absorption or transparant. \n" ,input); 348 349 return type; 350 } 351 352/////////////////////////////////////////////////////////////////////////// 353/////////////// Collision functions 354/////////////////////////////////////////////////////////////////////////// 355 356/** 357 Find the intersection between the neutron and the ellipse using newton method. 358 As there is up to 4 solution to this problem, and only the 359 smallest positive root is the physical solution. Using the tuning points 360 it is possible to look the only the potential roots to speed up calculations. 361 362 @param coef; A pointer to the array holding the coeffecients 363 for the 4th order polynomial. 364 @param startPosition, The default starting point for newton method. [s] 365 @param limit; A point after all the roots of the polynial. [s] 366 @param solution A pointer which will hold the physical solution 367 if this function return true. 368 @return; return 1 if the physical solution is found. [boolean] 369*/ 370 371double guide_elliptical_foverdf(double *coefficients,double currentPoint){ 372 double numerator= coefficients[0]*currentPoint*currentPoint*currentPoint*currentPoint 373 + coefficients[1]*currentPoint*currentPoint*currentPoint 374 + coefficients[2]*currentPoint*currentPoint 375 + coefficients[3]*currentPoint 376 + coefficients[4]; 377 double denominator=4*coefficients[0]*currentPoint*currentPoint*currentPoint 378 + 3*coefficients[1]*currentPoint*currentPoint 379 + 2*coefficients[2]*currentPoint 380 + coefficients[3]; 381 return numerator/denominator; 382} 383 384int guide_elliptical_newtonRapsonsMethod4thOrder( 385 double *coefficients,double *solution,double startingPoint, 386 double tolerance,double max_iterations){ 387 388 double numerator; 389 double denominator; 390 double t_previous; 391 double t = startingPoint; 392 int iteration = 0; 393 394 do { 395 t_previous = t; 396 t = t_previous - guide_elliptical_foverdf(coefficients,t); 397 iteration++; 398 } while( fabs(t-t_previous) > tolerance && iteration < max_iterations ); 399 if( iteration == max_iterations ) { return 0; } 400 else { *solution = t; return 1; } 401} 402 403 404 405int guide_elliptical_findNeutronEllipseIntersection( 406 double *coef,double startPosition, 407 double limit,double *solution){ 408 409 // in the case of no gravity 410 if(coef[0] == 0 & coef[1] == 0){ 411 double t1=0; 412 double t2=0; 413 int boolean = solve_2nd_order(&t1,&t2,coef[2],coef[3],coef[4]); 414 415 if ( t1 > startPosition ){ *solution = t1; } 416 if ( t2 > startPosition ){ *solution = t2; } 417 return boolean; 418 } 419 420 double tol = 1e-15; 421 double max_iter = 1e3; 422 double turningP1,turningP2; 423 424 double sp = startPosition; 425 int inside; 426 if ( coef[0]*sp*sp*sp*sp 427 +coef[1]*sp*sp*sp 428 +coef[2]*sp*sp 429 +coef[3]*sp 430 +coef[4] < 0) 431 inside = 1; 432 else inside = 0; 433 434 int boolean = solve_2nd_order( 435 &turningP1,&turningP2, 436 12*coef[0],6*coef[1],2*coef[2]); 437 438 double t1=0,t2=0; 439 double ss=100; 440 441 if( inside ){ 442 if(boolean) guide_elliptical_newtonRapsonsMethod4thOrder(coef,&t1,turningP1,tol,max_iter); 443 guide_elliptical_newtonRapsonsMethod4thOrder(coef,&t2,limit,tol,max_iter); 444 } 445 else{ 446 if(boolean) guide_elliptical_newtonRapsonsMethod4thOrder(coef,&t1,turningP2,tol,max_iter); 447 guide_elliptical_newtonRapsonsMethod4thOrder(coef,&t2,startPosition,tol,max_iter); 448 } 449 450 if (ss > t1 && t1 > 1e-15) ss = t1; 451 if (ss > t2 && t2 > 1e-15) ss = t2; 452 *solution = ss; 453 454 return 1; 455} 456 457 458int guide_elliptical_handleGuideIntersection( 459 double x, double y, double z, 460 double vx,double vy,double vz, 461 double Gx,double Gy,double Gz, 462 struct SGI *guideInfo, 463 struct Intersection *currentCollision){ 464 // 465 double horExS = 1/( guideInfo->ellipseMinorAxis[RightSide] 466 *guideInfo->ellipseMinorAxis[RightSide]); 467 double horEzS = 1/( guideInfo->ellipseMajorAxis[RightSide] 468 *guideInfo->ellipseMajorAxis[RightSide]); 469 double hordiffx = x-guideInfo->ellipseMinorOffset[RightSide]; 470 double hordiffz = z-guideInfo->ellipseMajorOffset[RightSide]; 471 472 double horAlpha = ( Gx*Gx*horExS + Gz*Gz*horEzS )/4; 473 double horBeta = ( Gx*vx*horExS + Gz*vz*horEzS ); 474 double horGamma = horExS*vx*vx + horEzS*vz*vz 475 + horExS*Gx*hordiffx + horEzS*Gz*hordiffz; 476 double horDelta = 2*horExS*vx*hordiffx + 2*horEzS*vz*hordiffz; 477 double horEpsilon = horExS*hordiffx*hordiffx + horEzS*hordiffz*hordiffz - 1; 478 479 double horCoefficients[5] = {horAlpha,horBeta,horGamma,horDelta,horEpsilon}; 480 481 double verEyS = 1/( guideInfo->ellipseMinorAxis[TopSide] 482 *guideInfo->ellipseMinorAxis[TopSide]); 483 double verEzS = 1/( guideInfo->ellipseMajorAxis[TopSide] 484 *guideInfo->ellipseMajorAxis[TopSide]); 485 double verdiffy = y-guideInfo->ellipseMinorOffset[TopSide]; 486 double verdiffz = z-guideInfo->ellipseMajorOffset[TopSide]; 487 488 double verAlpha = ( Gy*Gy*verEyS + Gz*Gz*verEzS )/4; 489 double verBeta = ( Gy*vy*verEyS + Gz*vz*verEzS ); 490 double verGamma = verEyS*vy*vy + verEzS*vz*vz 491 + verEyS*Gy*verdiffy + verEzS*Gz*verdiffz; 492 double verDelta = 2*verEyS*vy*verdiffy + 2*verEzS*vz*verdiffz; 493 double verEpsilon = verEyS*verdiffy*verdiffy + verEzS*verdiffz*verdiffz - 1; 494 495 double verCoefficients[5] = {verAlpha,verBeta,verGamma,verDelta,verEpsilon}; 496 497 498 double upperlimit; 499 double startingPoint = 1e-15; 500 501 int boolean; 502 // Horizontal 503 double solutionH = 0; 504 solve_2nd_order( 505 &upperlimit,NULL, 506 -0.5*Gz,-vz,2*guideInfo->ellipseMajorAxis[RightSide]-z); 507 int booleanH = guide_elliptical_findNeutronEllipseIntersection( 508 horCoefficients,startingPoint,upperlimit,&solutionH); 509 // Vertical 510 double solutionV = 0; 511 solve_2nd_order( 512 &upperlimit,NULL, 513 -0.5*Gz,-vz,2*guideInfo->ellipseMajorAxis[TopSide]-z); 514 int booleanV = guide_elliptical_findNeutronEllipseIntersection( 515 verCoefficients,startingPoint,upperlimit,&solutionV); 516 517 if (solutionH <= 0) 518 currentCollision->delta_time_to_next_collision = solutionV; 519 else if (solutionV <= 0) 520 currentCollision->delta_time_to_next_collision = solutionH; 521 else if (fabs(solutionH - solutionV) < 1e-12) return 0; 522 else if (solutionH < solutionV){ 523 currentCollision->delta_time_to_next_collision = solutionH; 524 boolean = booleanH; 525 } 526 else{ 527 currentCollision->delta_time_to_next_collision = solutionV; 528 boolean = booleanV; 529 } 530 531 double tside = currentCollision->delta_time_to_next_collision; 532 double xside = x + vx*tside + 0.5*Gx*tside*tside; 533 double yside = y + vy*tside + 0.5*Gy*tside*tside; 534 double zside = z + vz*tside + 0.5*Gz*tside*tside; 535 536 double xfactor = 537 2*sqrt(1 - ( (zside-guideInfo->ellipseMajorOffset[RightSide]) 538 *(zside-guideInfo->ellipseMajorOffset[RightSide]) 539 )/(guideInfo->ellipseMajorAxis[RightSide] 540 *guideInfo->ellipseMajorAxis[RightSide] ) 541 )*guideInfo->ellipseMinorAxis[RightSide]; 542 543 double yfactor = 544 2*sqrt(1 - ( (zside-guideInfo->ellipseMajorOffset[BottomSide]) 545 *(zside-guideInfo->ellipseMajorOffset[BottomSide]) 546 )/(guideInfo->ellipseMajorAxis[BottomSide] 547 *guideInfo->ellipseMajorAxis[BottomSide] ) 548 )*guideInfo->ellipseMinorAxis[BottomSide]; 549 550 xside = xside/xfactor; 551 yside = yside/yfactor; 552 if( fabs(yside) >= fabs(xside) ){ 553 if(y > 0) currentCollision->side = TopSide; 554 else currentCollision->side = BottomSide; 555 } 556 else{ 557 if(x < 0) currentCollision->side = RightSide; 558 else currentCollision->side = LeftSide; 559 } 560 if (tside < 1e-15) printf("low time is: %e\n",tside); 561 562 return boolean; 563} 564 565/** 566 Check if the neutron is within the guide using the sign 567 of the crossproduct between the two points, 568 on each of the enclosing box surface and neutrons position. 569 570 @param x,y,z; position of the neutron. [m] 571 @param guideInfo; pointer to the guide infomation holding structure. 572 @return; return 1 if the neutron is inside the guide [boolean] 573*/ 574 575/* 576int guide_elliptical_InsideEnclosingBox(double x,double y,double z,struct SGI *guideInfo){ 577 int guide_elliptical_IsPointInVolume( 578 double *x,double *y,double *z, 579 double px,double py,double pz){ 580 int guide_elliptical_WhichSide( double p1x,double p1y,double p1z, 581 double p2x,double p2y,double p2z, 582 double p3x,double p3y,double p3z, 583 double px ,double py ,double pz ){ 584 585 double v1x = p1x - p2x, v1y = p1y-p2y, v1z = p1z-p2z; 586 double v2x = p3x - p2x, v2y = p3y-p2y, v2z = p3z-p2z; 587 double v3x = v2y*v1z-v2z*v1y; 588 double v3y = v2z*v1x-v2x*v1z; 589 double v3z = v2x*v1y-v2y*v1x; 590 591 return 0 >= v3x*(px-p1x)+v3y*(py-p1y)+v3z*(pz-p1z); 592 } 593 594 if( //front 595 guide_elliptical_WhichSide(x[3],y[3],z[3],x[2],y[2],z[2],x[1],y[1],z[1],px,py,pz) && 596 guide_elliptical_WhichSide(x[1],y[1],z[1],x[0],y[0],z[0],x[3],y[3],z[3],px,py,pz) && 597 //back 598 guide_elliptical_WhichSide(x[5],y[5],z[5],x[6],y[6],z[6],x[7],y[7],z[7],px,py,pz) && 599 guide_elliptical_WhichSide(x[7],y[7],z[7],x[4],y[4],z[4],x[5],y[5],z[5],px,py,pz) && 600 //right 601 guide_elliptical_WhichSide(x[7],y[7],z[7],x[3],y[3],z[3],x[0],y[0],z[0],px,py,pz) && 602 guide_elliptical_WhichSide(x[0],y[0],z[0],x[4],y[4],z[4],x[7],y[7],z[7],px,py,pz) && 603 //left 604 guide_elliptical_WhichSide(x[1],y[1],z[1],x[2],y[2],z[2],x[6],y[6],z[6],px,py,pz) && 605 guide_elliptical_WhichSide(x[6],y[6],z[6],x[5],y[5],z[5],x[1],y[1],z[1],px,py,pz) && 606 //top 607 guide_elliptical_WhichSide(x[0],y[0],z[0],x[1],y[1],z[1],x[5],y[5],z[5],px,py,pz) && 608 guide_elliptical_WhichSide(x[5],y[5],z[5],x[4],y[4],z[4],x[0],y[0],z[0],px,py,pz) && 609 //bottom 610 guide_elliptical_WhichSide(x[6],y[6],z[6],x[2],y[2],z[2],x[3],y[3],z[3],px,py,pz) && 611 guide_elliptical_WhichSide(x[3],y[3],z[3],x[7],y[7],z[7],x[6],y[6],z[6],px,py,pz) ) 612 return 1; 613 else return 0; 614 } 615 return guide_elliptical_IsPointInVolume( 616 guideInfo->xArray,guideInfo->yArray,guideInfo->zArray,x,y,z); 617} 618*/ 619 620 621/////////////////////////////////////////////////////////////////////////// 622/////////////// reflection functions 623/////////////////////////////////////////////////////////////////////////// 624 625 626/** 627 Calculate the new velocity vector for the particle colliding on 628 the inner side of the elliptic mirror and returns the loss-factor (TODO) 629 630 @param pos_V0,pos_W0 Is the 2d position vector of the particle, 631 assumed to be a point on the ellipse. [m] 632 @param pvel_V0,pvel_W0 Is the 2d velocity vector of the particle. [m/s] 633 @param ellipse_V_axis_squared,ellipse_W_axis_squared 634 are the axes of the ellipse. [m] 635 @param ellipse_V_offset,ellipse_W_offset Is the 2d vector difference 636 between the ellipse coordinate system (center of the ellipse) 637 and the guide coordinate system [m] 638 @param R0, Mvalue, Qc, W, Alpha #TODO 639 slaa beskrivelse af disse variabler i andre dokumenter 640 og hold dig til standarden. 641 @return the new wieght of the package 642*/ 643double guide_elliptical_ReflectionOnEllipticSurface( 644 double pos_V,double pos_W, 645 double *pvel_V,double *pvel_W, 646 double ellipse_V_axis,double ellipse_W_axis, 647 double ellipse_V_offset,double ellipse_W_offset, 648 double R0, double Qc, double alpha, double Mvalue, double W) 649{ 650 651 // Turns the velocity vector (vel_V0,vel_W0) into a local value 652 double vel_V = *pvel_V; 653 double vel_W = *pvel_W; 654 655 // Galilean transformation of the particles start position 656 // to the ellipse coordinate system 657 pos_V=pos_V-ellipse_V_offset; 658 pos_W=pos_W-ellipse_W_offset; 659 660 /* 661 * If we reflect the velocity vector in the normal 662 * to the ellipse in the point of intersection 663 * The resulting vector will be -f2, do to conservation of momentum. 664 * this result in the following equation 665 * f2 = -f1 + 2(f1 dot nhat)nhat 666 * which is equal to f2 = f1 - 2(f1 dot n)n/nlength^2 667 */ 668 669 // The normal vector to the point of intersection 670 double normVec_V = - pos_W*ellipse_V_axis/ellipse_W_axis; 671 double normVec_W = pos_V*ellipse_W_axis/ellipse_V_axis; 672 673 double normVec_length_squared = normVec_V*normVec_V + normVec_W*normVec_W; 674 675 // Dot product of (vel_V0,vel_W0) and the normal vector 676 double Vel_dot_NV = vel_V*normVec_V+vel_W*normVec_W; 677 678 // Calculate f2 679 double vel_V_2 = -vel_V + 2*Vel_dot_NV*normVec_V/normVec_length_squared; 680 double vel_W_2 = -vel_W + 2*Vel_dot_NV*normVec_W/normVec_length_squared; 681 682 // Apply the new velocity vector to the particle globally 683 *pvel_V=vel_V_2; 684 *pvel_W=vel_W_2; 685 686 // Calculate q and the weighting of the neutron package 687 // q=f1-f2 688 double delta_vel_V = vel_V-vel_V_2; 689 double delta_vel_W = vel_W-vel_W_2; 690 double q = V2Q*sqrt( delta_vel_V*delta_vel_V+delta_vel_W*delta_vel_W ); 691 692 // Calculate the loss of neutrons due to the reflection 693 double mirrorPar[] = {R0, Qc, alpha, Mvalue, W}; 694 double weight = 1.0; 695 StdReflecFunc(q, mirrorPar, &weight); 696 697 return weight; 698} 699 700/** 701 Use the found side of Intersection to call guide_elliptical_ReflectionOnEllipticSurface with 702 the parameters of that side. 703*/ 704double guide_elliptical_handleReflection(double x0, double y0, double z0, 705 double *vx_p,double *vy_p,double *vz_p, 706 struct SGI *sgi, 707 struct Intersection *cc) 708{ 709 710 if(!sgi->enableSegments){ 711 if(cc->side == RightSide || cc->side == LeftSide) 712 return guide_elliptical_ReflectionOnEllipticSurface(x0,z0,vx_p,vz_p, 713 sgi->ellipseMinorAxis[cc->side], 714 sgi->ellipseMajorAxis[cc->side], 715 sgi->ellipseMinorOffset[cc->side], 716 sgi->ellipseMajorOffset[cc->side], 717 sgi->R0Arr[cc->side], 718 sgi->QcArr[cc->side], 719 sgi->alphaArr[cc->side], 720 sgi->mArr[cc->side], 721 sgi->WArr[cc->side] 722 ); 723 if(cc->side == TopSide || cc->side == BottomSide) 724 return guide_elliptical_ReflectionOnEllipticSurface(y0,z0,vy_p,vz_p, 725 sgi->ellipseMinorAxis[cc->side], 726 sgi->ellipseMajorAxis[cc->side], 727 sgi->ellipseMinorOffset[cc->side], 728 sgi->ellipseMajorOffset[cc->side], 729 sgi->R0Arr[cc->side], 730 sgi->QcArr[cc->side], 731 sgi->alphaArr[cc->side], 732 sgi->mArr[cc->side], 733 sgi->WArr[cc->side] 734 ); 735 } 736 else{ 737 int currentSegment; 738 double combinedLength = 0; 739 int i; 740 for(i=0; i < sgi->numberOfSegments; i++){ 741 combinedLength = combinedLength + sgi->segLength[i]; 742 if(z0 < combinedLength) { 743 currentSegment = i; break; 744 } 745 } 746 747 if(cc->side == RightSide) 748 return guide_elliptical_ReflectionOnEllipticSurface(x0,z0,vx_p,vz_p, 749 sgi->ellipseMinorAxis[cc->side], 750 sgi->ellipseMajorAxis[cc->side], 751 sgi->ellipseMinorOffset[cc->side], 752 sgi->ellipseMajorOffset[cc->side], 753 sgi->R0Arr[cc->side], 754 sgi->QcArr[cc->side], 755 sgi->alphaArr[cc->side], 756 sgi->mValuesright[currentSegment], 757 sgi->WArr[cc->side] ); 758 if(cc->side == LeftSide) 759 return guide_elliptical_ReflectionOnEllipticSurface(x0,z0,vx_p,vz_p, 760 sgi->ellipseMinorAxis[cc->side], 761 sgi->ellipseMajorAxis[cc->side], 762 sgi->ellipseMinorOffset[cc->side], 763 sgi->ellipseMajorOffset[cc->side], 764 sgi->R0Arr[cc->side], 765 sgi->QcArr[cc->side], 766 sgi->alphaArr[cc->side], 767 sgi->mValuesleft[currentSegment], 768 sgi->WArr[cc->side] ); 769 if(cc->side == TopSide) 770 return guide_elliptical_ReflectionOnEllipticSurface(y0,z0,vy_p,vz_p, 771 sgi->ellipseMinorAxis[cc->side], 772 sgi->ellipseMajorAxis[cc->side], 773 sgi->ellipseMinorOffset[cc->side], 774 sgi->ellipseMajorOffset[cc->side], 775 sgi->R0Arr[cc->side], 776 sgi->QcArr[cc->side], 777 sgi->alphaArr[cc->side], 778 sgi->mValuestop[currentSegment], 779 sgi->WArr[cc->side] ); 780 if(cc->side == BottomSide) 781 return guide_elliptical_ReflectionOnEllipticSurface(y0,z0,vy_p,vz_p, 782 sgi->ellipseMinorAxis[cc->side], 783 sgi->ellipseMajorAxis[cc->side], 784 sgi->ellipseMinorOffset[cc->side], 785 sgi->ellipseMajorOffset[cc->side], 786 sgi->R0Arr[cc->side], 787 sgi->QcArr[cc->side], 788 sgi->alphaArr[cc->side], 789 sgi->mValuesbottom[currentSegment], 790 sgi->WArr[cc->side] ); 791 } 792 return 0; 793} 794 795/////////////////////////////////////////////////////////////////////////// 796/////////////// End of functions 797/////////////////////////////////////////////////////////////////////////// 798%} 799 800 801DECLARE 802%{ 803 /** 804 All variables below is locally declared 805 and hence accessible through OUTPUT PARAMETERS. 806 */ 807 struct SGI guideInfo; // Static Guide information, is set in INITIALIZE 808 struct Intersection latestParticleCollision; // Is changed duing trace 809 double Gx,Gy,Gz; // Local gravity vector, is set once in INITIALIZE 810 double Gx0, Gy0, Gz0; 811 double Circ; 812 double *dynamicalSegLength; 813 814%} 815 816INITIALIZE 817%{ 818 /////////////////////////////////////////////////////////////////////////// 819 /////////////// Test user input 820 /////////////////////////////////////////////////////////////////////////// 821 822 if(strcmp(verbose,"on") == 0) 823 guideInfo.verboseSetting = 1; 824 else guideInfo.verboseSetting = 0; 825 826 827 guideInfo.R0Arr[RightSide] = R0; 828 guideInfo.QcArr[RightSide] = Qc; 829 guideInfo.alphaArr[RightSide] = alpharight; 830 guideInfo.mArr[RightSide] = mright; 831 guideInfo.WArr[RightSide] = W; 832 833 guideInfo.R0Arr[TopSide] = R0; 834 guideInfo.QcArr[TopSide] = Qc; 835 guideInfo.alphaArr[TopSide] = alphatop; 836 guideInfo.mArr[TopSide] = mtop; 837 guideInfo.WArr[TopSide] = W; 838 839 guideInfo.R0Arr[LeftSide] = R0; 840 guideInfo.QcArr[LeftSide] = Qc; 841 guideInfo.alphaArr[LeftSide] = alphaleft; 842 guideInfo.mArr[LeftSide] = mleft; 843 guideInfo.WArr[LeftSide] = W; 844 845 guideInfo.R0Arr[BottomSide] = R0; 846 guideInfo.QcArr[BottomSide] = Qc; 847 guideInfo.alphaArr[BottomSide] = alphabottom; 848 guideInfo.mArr[BottomSide] = mbottom; 849 guideInfo.WArr[BottomSide] = W; 850 851 int sides; 852 for (sides = RightSide; sides <= BottomSide; sides++){ 853 if (guideInfo.alphaArr[sides] == -1) guideInfo.alphaArr[sides] = alpha; 854 if (guideInfo.mArr[sides] == -1) guideInfo.mArr[sides] = m; 855 } 856 857 // Test user input for illegal values 858 int inputErrors = 0; 859 // Lower or equal to zero 860 if(l <= 0) inputErrors += 861 guide_elliptical_illegalInputLessThanZero("length",guideInfo.verboseSetting); 862 if(guideInfo.alphaArr[TopSide] < 0) inputErrors += 863 guide_elliptical_illegalInputLessThanZero("alphatop",guideInfo.verboseSetting); 864 if(guideInfo.mArr[TopSide] < 0) inputErrors += 865 guide_elliptical_illegalInputLessThanZero("mtop",guideInfo.verboseSetting); 866 867 if(guideInfo.alphaArr[BottomSide] < 0) inputErrors += 868 guide_elliptical_illegalInputLessThanZero("alphabottom",guideInfo.verboseSetting); 869 if(guideInfo.mArr[BottomSide] < 0) inputErrors += 870 guide_elliptical_illegalInputLessThanZero("mbottom",guideInfo.verboseSetting); 871 872 if(guideInfo.alphaArr[RightSide] < 0) inputErrors += 873 guide_elliptical_illegalInputLessThanZero("alpharight",guideInfo.verboseSetting); 874 if(guideInfo.mArr[RightSide] < 0) inputErrors += 875 guide_elliptical_illegalInputLessThanZero("mright",guideInfo.verboseSetting); 876 877 if(guideInfo.alphaArr[LeftSide] < 0) inputErrors += 878 guide_elliptical_illegalInputLessThanZero("alphaleft",guideInfo.verboseSetting); 879 if(guideInfo.mArr[LeftSide] < 0) inputErrors += 880 guide_elliptical_illegalInputLessThanZero("mleft",guideInfo.verboseSetting); 881 882 // Focal points result in hyperbola instead of an ellipse 883 if(l <= -linxw-loutxw) inputErrors += guide_elliptical_illegalInputFocalPointsHyperbola( 884 "linw","loutw",linxw,loutxw,guideInfo.verboseSetting); 885 if(l <= -linyh-loutyh) inputErrors += guide_elliptical_illegalInputFocalPointsHyperbola( 886 "linh","louth",linyh,loutyh,guideInfo.verboseSetting); 887 888 if( strcmp(dimensionsAt,"entrance") != 0 889 && strcmp(dimensionsAt,"mid") != 0 890 && strcmp(dimensionsAt,"exit") != 0){ 891 inputErrors += 1; 892 printf("dimensionsAt were given an incorrect input." 893 "Input must be string containing \"entrance\",\"mid\" or \"exit\" \n"); 894 } 895 896 897 // Terminate program if any input errors occurred 898 if(inputErrors != 0 ){ 899 exit(printf("\nCRITICAL ERROR(S) IN COMPONENT %s" 900 " CONSIDER CHECKING USER INPUT AS %d INPUT ERRORS WAS FOUND.\n", 901 NAME_CURRENT_COMP,inputErrors) ); 902 } 903 904 905 /////////////////////////////////////////////////////////////////////////// 906 /////////////// Calculate intern guide values from user input 907 /////////////////////////////////////////////////////////////////////////// 908 909 /* Calculate the foci line for the ellipses. 910 These can be used to calculate the axes of the ellipses 911 using pyth and defination of the ellipse that says distance 912 between the foci and every point on the ellipse is constant. 913 */ 914 int directDefination = 0; 915 916 if( majorAxisyh != 0 || minorAxisyh != 0 917 || majorAxisxw != 0 || minorAxisxw != 0) 918 { 919 directDefination = 1; 920 guideInfo.Length = l; 921 922 guideInfo.ellipseMajorAxis[RightSide] = majorAxisxw; 923 guideInfo.ellipseMinorAxis[RightSide] = minorAxisxw; 924 guideInfo.ellipseMajorOffset[RightSide] = majorAxisoffsetxw; 925 guideInfo.ellipseMinorOffset[RightSide] = 0; 926 927 guideInfo.ellipseMajorAxis[TopSide] = majorAxisyh; 928 guideInfo.ellipseMinorAxis[TopSide] = minorAxisyh; 929 guideInfo.ellipseMajorOffset[TopSide] = majorAxisoffsetyh; 930 guideInfo.ellipseMinorOffset[TopSide] = 0; 931 932 guideInfo.ellipseMajorAxis[LeftSide] = majorAxisxw; 933 guideInfo.ellipseMinorAxis[LeftSide] = minorAxisxw; 934 guideInfo.ellipseMajorOffset[LeftSide] = majorAxisoffsetxw; 935 guideInfo.ellipseMinorOffset[LeftSide] = 0; 936 937 guideInfo.ellipseMajorAxis[BottomSide] = majorAxisyh; 938 guideInfo.ellipseMinorAxis[BottomSide] = minorAxisyh; 939 guideInfo.ellipseMajorOffset[BottomSide] = majorAxisoffsetyh; 940 guideInfo.ellipseMinorOffset[BottomSide] = 0; 941 942 guideInfo.entranceHorizontalWidth = 943 2*sqrt(1 - (majorAxisoffsetyh*majorAxisoffsetyh) 944 /(majorAxisyh*majorAxisyh) )*minorAxisyh; 945 guideInfo.entranceVerticalWidth = 946 2*sqrt(1 - (majorAxisoffsetxw*majorAxisoffsetxw) 947 /(majorAxisxw*majorAxisxw) )*minorAxisxw; 948 } 949 950 if ( strcmp(option,"ellipse") == 0 && directDefination == 0) 951 { 952 if ( strcmp(dimensionsAt,"entrance") == 0 ){ 953 double lofbs_horizontal = 954 sqrt( linxw*linxw + xwidth*xwidth*0.25) 955 + sqrt( (l + loutxw)*(l + loutxw) + xwidth*xwidth*0.25); 956 957 double lofbs_vertical = 958 sqrt( linyh*linyh + yheight*yheight*0.25) 959 + sqrt( (l + loutyh)*(l + loutyh) + yheight*yheight*0.25); 960 961 guideInfo.Length = l; 962 963 guideInfo.ellipseMajorAxis[RightSide] = lofbs_horizontal/2; 964 guideInfo.ellipseMinorAxis[RightSide] = 965 sqrt(0.25*lofbs_horizontal*lofbs_horizontal 966 -0.25*(l+linxw+loutxw)*(l+linxw+loutxw) ); 967 968 guideInfo.ellipseMajorOffset[RightSide] = (l+linxw+loutxw)/2-linxw; 969 guideInfo.ellipseMinorOffset[RightSide] = 0; 970 971 guideInfo.ellipseMajorAxis[LeftSide] = 972 guideInfo.ellipseMajorAxis[RightSide]; 973 guideInfo.ellipseMinorAxis[LeftSide] = 974 guideInfo.ellipseMinorAxis[RightSide]; 975 guideInfo.ellipseMajorOffset[LeftSide] = 976 guideInfo.ellipseMajorOffset[RightSide]; 977 guideInfo.ellipseMinorOffset[LeftSide] = 978 guideInfo.ellipseMinorOffset[RightSide]; 979 980 guideInfo.ellipseMajorAxis[TopSide] = lofbs_vertical/2; 981 982 guideInfo.ellipseMinorAxis[TopSide] = 983 sqrt(0.25*lofbs_vertical*lofbs_vertical 984 -0.25*(l+linyh+loutyh)*(l+linyh+loutyh) ); 985 986 guideInfo.ellipseMajorOffset[TopSide] = (l+linyh+loutyh)/2-linyh; 987 guideInfo.ellipseMinorOffset[TopSide] = 0; 988 989 guideInfo.ellipseMajorAxis[BottomSide] = 990 guideInfo.ellipseMajorAxis[TopSide]; 991 guideInfo.ellipseMinorAxis[BottomSide] = 992 guideInfo.ellipseMinorAxis[TopSide]; 993 guideInfo.ellipseMajorOffset[BottomSide] = 994 guideInfo.ellipseMajorOffset[TopSide]; 995 guideInfo.ellipseMinorOffset[BottomSide] = 996 guideInfo.ellipseMinorOffset[TopSide]; 997 } 998 if ( strcmp(dimensionsAt,"exit") == 0 ){ 999 double lofbs_horizontal = 1000 sqrt( loutxw*loutxw + xwidth*xwidth*0.25) 1001 + sqrt( (l + linxw)*(l + linxw) + xwidth*xwidth*0.25); 1002 1003 double lofbs_vertical = 1004 sqrt( loutyh*loutyh + yheight*yheight*0.25) 1005 + sqrt( (l + linyh)*(l + linyh) + yheight*yheight*0.25); 1006 1007 guideInfo.Length = l; 1008 1009 guideInfo.ellipseMajorAxis[RightSide] = lofbs_horizontal/2; 1010 guideInfo.ellipseMinorAxis[RightSide] = 1011 sqrt(0.25*lofbs_horizontal*lofbs_horizontal 1012 -0.25*(l+linxw+loutxw)*(l+linxw+loutxw) ); 1013 1014 guideInfo.ellipseMajorOffset[RightSide] =(l+linxw+loutxw)/2-linxw; 1015 guideInfo.ellipseMinorOffset[RightSide] = 0; 1016 1017 guideInfo.ellipseMajorAxis[LeftSide] = 1018 guideInfo.ellipseMajorAxis[RightSide]; 1019 guideInfo.ellipseMinorAxis[LeftSide] = 1020 guideInfo.ellipseMinorAxis[RightSide]; 1021 guideInfo.ellipseMajorOffset[LeftSide] = 1022 guideInfo.ellipseMajorOffset[RightSide]; 1023 guideInfo.ellipseMinorOffset[LeftSide] = 1024 guideInfo.ellipseMinorOffset[RightSide]; 1025 1026 guideInfo.ellipseMajorAxis[TopSide] = lofbs_vertical/2; 1027 1028 guideInfo.ellipseMinorAxis[TopSide] = 1029 sqrt(0.25*lofbs_vertical*lofbs_vertical 1030 -0.25*(l+linyh+loutyh)*(l+linyh+loutyh) ); 1031 1032 guideInfo.ellipseMajorOffset[TopSide] = (l+linyh+loutyh)/2-linyh; 1033 guideInfo.ellipseMinorOffset[TopSide] = 0; 1034 1035 guideInfo.ellipseMajorAxis[BottomSide] = 1036 guideInfo.ellipseMajorAxis[TopSide]; 1037 guideInfo.ellipseMinorAxis[BottomSide] = 1038 guideInfo.ellipseMinorAxis[TopSide]; 1039 guideInfo.ellipseMajorOffset[BottomSide] = 1040 guideInfo.ellipseMajorOffset[TopSide]; 1041 guideInfo.ellipseMinorOffset[BottomSide] = 1042 guideInfo.ellipseMinorOffset[TopSide]; 1043 } 1044 if ( strcmp(dimensionsAt,"mid") == 0 ){ 1045 1046 guideInfo.Length = l; 1047 1048 guideInfo.ellipseMajorAxis[RightSide] = 1049 sqrt( (linxw+l+loutxw)*(linxw+l+loutxw)/4+xwidth*xwidth/4); 1050 guideInfo.ellipseMinorAxis[RightSide] = xwidth/2; 1051 1052 guideInfo.ellipseMajorOffset[RightSide] = (l+linxw+loutxw)/2-linxw; 1053 guideInfo.ellipseMinorOffset[RightSide] = 0; 1054 1055 guideInfo.ellipseMajorAxis[LeftSide] = 1056 guideInfo.ellipseMajorAxis[RightSide]; 1057 guideInfo.ellipseMinorAxis[LeftSide] = 1058 guideInfo.ellipseMinorAxis[RightSide]; 1059 guideInfo.ellipseMajorOffset[LeftSide] = 1060 guideInfo.ellipseMajorOffset[RightSide]; 1061 guideInfo.ellipseMinorOffset[LeftSide] = 1062 guideInfo.ellipseMinorOffset[RightSide]; 1063 1064 guideInfo.ellipseMajorAxis[TopSide] = 1065 sqrt( (linyh+l+loutyh)*(linyh+l+loutyh)/4+yheight*yheight/4); 1066 guideInfo.ellipseMinorAxis[TopSide] = yheight/2; 1067 1068 guideInfo.ellipseMajorOffset[TopSide] = (l+linyh+loutyh)/2-linyh; 1069 guideInfo.ellipseMinorOffset[TopSide] = 0; 1070 1071 guideInfo.ellipseMajorAxis[BottomSide] = 1072 guideInfo.ellipseMajorAxis[TopSide]; 1073 guideInfo.ellipseMinorAxis[BottomSide] = 1074 guideInfo.ellipseMinorAxis[TopSide]; 1075 guideInfo.ellipseMajorOffset[BottomSide] = 1076 guideInfo.ellipseMajorOffset[TopSide]; 1077 guideInfo.ellipseMinorOffset[BottomSide] = 1078 guideInfo.ellipseMinorOffset[TopSide]; 1079 1080 } 1081 } 1082 1083 guideInfo.entranceHorizontalWidth = 2*sqrt( 1084 1 - guideInfo.ellipseMajorOffset[RightSide] 1085 *guideInfo.ellipseMajorOffset[RightSide] 1086 /(guideInfo.ellipseMajorAxis[RightSide] 1087 *guideInfo.ellipseMajorAxis[RightSide] ) ) 1088 *guideInfo.ellipseMinorAxis[RightSide]; 1089 guideInfo.entranceVerticalWidth = 2*sqrt( 1090 1 - guideInfo.ellipseMajorOffset[TopSide] 1091 *guideInfo.ellipseMajorOffset[TopSide] 1092 /(guideInfo.ellipseMajorAxis[TopSide] 1093 *guideInfo.ellipseMajorAxis[TopSide] ) ) 1094 *guideInfo.ellipseMinorAxis[TopSide]; 1095 1096 1097 if ( strcmp(option,"halfellipse") == 0 && directDefination == 0 ){ 1098 exit( printf("Critical error in %s; the option for option = halfellipse is currently disabled.",NAME_CURRENT_COMP) ); 1099 1100 double used_focal_vertical; 1101 double used_focal_horizontal; 1102 double major_offset_horizontal = 0; 1103 double major_offset_vertical = 0; 1104 1105 if ( strcmp(dimensionsAt,"entrance") == 0 ){ 1106 used_focal_vertical = sqrt( (yheight*yheight)/4 1107 + (l+linyh)*(l+linyh) ); 1108 used_focal_horizontal = sqrt( (xwidth*xwidth)/4 1109 + (l+linxw)*(l+linxw) ); 1110 major_offset_vertical = l; 1111 major_offset_horizontal = l; 1112 } 1113 else { 1114 used_focal_vertical = sqrt( (yheight*yheight)/4 1115 + (l+loutyh)*(l+loutyh) ); 1116 used_focal_horizontal = sqrt( (xwidth*xwidth)/4 1117 + (l+loutxw)*(l+loutxw) ); 1118 } 1119 1120 guideInfo.Length = l; 1121 1122 guideInfo.ellipseMajorAxis[RightSide] = used_focal_horizontal; 1123 guideInfo.ellipseMinorAxis[RightSide] = xwidth/2; 1124 1125 guideInfo.ellipseMajorOffset[RightSide] = major_offset_horizontal; 1126 guideInfo.ellipseMinorOffset[RightSide] = 0; 1127 1128 guideInfo.ellipseMajorAxis[LeftSide] = 1129 guideInfo.ellipseMajorAxis[RightSide]; 1130 guideInfo.ellipseMinorAxis[LeftSide] = 1131 guideInfo.ellipseMinorAxis[RightSide]; 1132 guideInfo.ellipseMajorOffset[LeftSide] = 1133 guideInfo.ellipseMajorOffset[RightSide]; 1134 guideInfo.ellipseMinorOffset[LeftSide] = 1135 guideInfo.ellipseMinorOffset[RightSide]; 1136 1137 guideInfo.ellipseMajorAxis[TopSide] = used_focal_vertical; 1138 guideInfo.ellipseMinorAxis[TopSide] = yheight/2; 1139 1140 guideInfo.ellipseMajorOffset[TopSide] = major_offset_vertical; 1141 guideInfo.ellipseMinorOffset[TopSide] = 0; 1142 1143 guideInfo.ellipseMajorAxis[BottomSide] = 1144 guideInfo.ellipseMajorAxis[TopSide]; 1145 guideInfo.ellipseMinorAxis[BottomSide] = 1146 guideInfo.ellipseMinorAxis[TopSide]; 1147 guideInfo.ellipseMajorOffset[BottomSide] = 1148 guideInfo.ellipseMajorOffset[TopSide]; 1149 guideInfo.ellipseMinorOffset[BottomSide] = 1150 guideInfo.ellipseMinorOffset[TopSide]; 1151 } 1152 1153 // Applies the properties of the mirrors in the guide given by the user. 1154 // These variables are used in the reflection functions. 1155 1156 1157 // Sets the mirror type of the guides mirrors 1158 // These variables are used in the collision functions 1159 // to find the type of collision 1160 1161 // guideInfo.OuterSide[RightSide] = 1162 // guide_elliptical_getMirrorTypeFromInput(outer_right_side_mirror,guideInfo.verboseSetting); 1163 // guideInfo.OuterSide[TopSide] = 1164 // guide_elliptical_getMirrorTypeFromInput(outer_top_side_mirror,guideInfo.verboseSetting); 1165 // guideInfo.OuterSide[LeftSide] = 1166 // guide_elliptical_getMirrorTypeFromInput(outer_left_side_mirror,guideInfo.verboseSetting); 1167 // guideInfo.OuterSide[BottomSide] = 1168 // guide_elliptical_getMirrorTypeFromInput(outer_bottom_side_mirror,guideInfo.verboseSetting); 1169 1170 // guideInfo.InnerSide[RightSide] = 1171 // guide_elliptical_getMirrorTypeFromInput(inner_right_side_mirror,guideInfo.verboseSetting); 1172 // guideInfo.InnerSide[TopSide] = 1173 // guide_elliptical_getMirrorTypeFromInput(inner_top_side_mirror,guideInfo.verboseSetting); 1174 // guideInfo.InnerSide[LeftSide] = 1175 // guide_elliptical_getMirrorTypeFromInput(inner_left_side_mirror,guideInfo.verboseSetting); 1176 // guideInfo.InnerSide[BottomSide] = 1177 // guide_elliptical_getMirrorTypeFromInput(inner_bottom_side_mirror,guideInfo.verboseSetting); 1178 1179 // Give a warning if all side of the guide is turned off, 1180 // as the guide is essentially turned off 1181 if( guideInfo.OuterSide[RightSide] == 1 1182 && guideInfo.OuterSide[TopSide] == 1 1183 && guideInfo.OuterSide[LeftSide] == 1 1184 && guideInfo.OuterSide[BottomSide] == 1 1185 && guideInfo.InnerSide[RightSide] == 1 1186 && guideInfo.InnerSide[TopSide] == 1 1187 && guideInfo.InnerSide[LeftSide] == 1 1188 && guideInfo.InnerSide[BottomSide] == 1) 1189 printf("Warning: In %s all the sides of the guide has been disabled," 1190 " so it not possible for any particle" 1191 " to collide with the guide, consider" 1192 " disabling this component",NAME_CURRENT_COMP); 1193 1194 if(guideInfo.mArr[RightSide] <= 0) guideInfo.InnerSide[RightSide] = 1195 MirrorTypeabsorption; 1196 if(guideInfo.mArr[TopSide] <= 0) guideInfo.InnerSide[TopSide] = 1197 MirrorTypeabsorption; 1198 if(guideInfo.mArr[LeftSide] <= 0) guideInfo.InnerSide[LeftSide] = 1199 MirrorTypeabsorption; 1200 if(guideInfo.mArr[BottomSide] <= 0) guideInfo.InnerSide[BottomSide] = 1201 MirrorTypeabsorption; 1202 /* if(directDefination == 0){ */ 1203 /* guideInfo.entranceHorizontalWidth = xwidth; */ 1204 /* guideInfo.entranceVerticalWidth = yheight; */ 1205 /* } */ 1206 1207 if( strcmp(option,"halfellipse") == 0 && directDefination == 0 ){ 1208 guideInfo.entranceHorizontalWidth = 1209 (guideInfo.ellipseMinorAxis[RightSide] 1210 * sqrt(1 - ( guideInfo.ellipseMajorOffset[RightSide] 1211 *guideInfo.ellipseMajorOffset[RightSide] ) 1212 /( guideInfo.ellipseMajorAxis[RightSide] 1213 *guideInfo.ellipseMajorAxis[RightSide] ) ) 1214 + guideInfo.ellipseMinorOffset[RightSide] )*2; 1215 guideInfo.entranceVerticalWidth = 1216 (guideInfo.ellipseMinorAxis[TopSide] 1217 * sqrt(1 - ( guideInfo.ellipseMajorOffset[TopSide] 1218 *guideInfo.ellipseMajorOffset[TopSide] ) 1219 /( guideInfo.ellipseMajorAxis[TopSide] 1220 *guideInfo.ellipseMajorAxis[TopSide] ) ) 1221 + guideInfo.ellipseMinorOffset[TopSide] )*2; 1222 } 1223 1224 1225 guideInfo.EnclosingBoxOn = 0; 1226 1227 /* 1228 double DefaultArray1[8] = { 1.0,-1.0,-1.0, 1.0, 1.0,-1.0,-1.0, 1.0}; 1229 double DefaultArray2[8] = { 1.0, 1.0,-1.0,-1.0, 1.0, 1.0,-1.0,-1.0}; 1230 double DefaultArray3[8] = {-1.0,-1.0,-1.0,-1.0, 1.0, 1.0, 1.0, 1.0}; 1231 1232 guideInfo.EnclosingBoxOn = 0; 1233 double *xinput; 1234 if ( xInput != NULL ){ xinput = xInput; guideInfo.EnclosingBoxOn = 1; } 1235 else { xinput = DefaultArray1; } 1236 double *yinput; 1237 if ( yInput != NULL ){ yinput = yInput; guideInfo.EnclosingBoxOn = 1;} 1238 else { yinput = DefaultArray2; } 1239 double *zinput; 1240 if ( zInput != NULL ){ zinput = zInput; guideInfo.EnclosingBoxOn = 1;} 1241 else { zinput = DefaultArray3; } 1242 */ 1243 1244 /* 1245 double xarray[8] ={ guideInfo.ellipseMinorAxis[0]*xinput[0], 1246 guideInfo.ellipseMinorAxis[2]*xinput[1], 1247 guideInfo.ellipseMinorAxis[2]*xinput[2], 1248 guideInfo.ellipseMinorAxis[0]*xinput[3], 1249 guideInfo.ellipseMinorAxis[0]*xinput[4], 1250 guideInfo.ellipseMinorAxis[2]*xinput[5], 1251 guideInfo.ellipseMinorAxis[2]*xinput[6], 1252 guideInfo.ellipseMinorAxis[0]*xinput[7] }; 1253 double yarray[8] ={ guideInfo.ellipseMinorAxis[1]*yinput[0], 1254 guideInfo.ellipseMinorAxis[1]*yinput[1], 1255 guideInfo.ellipseMinorAxis[3]*yinput[2], 1256 guideInfo.ellipseMinorAxis[3]*yinput[3], 1257 guideInfo.ellipseMinorAxis[1]*yinput[4], 1258 guideInfo.ellipseMinorAxis[1]*yinput[5], 1259 guideInfo.ellipseMinorAxis[3]*yinput[6], 1260 guideInfo.ellipseMinorAxis[3]*yinput[7] }; 1261 double zarray[8] ={ guideInfo.Length/2*zinput[0]+guideInfo.Length/2, 1262 guideInfo.Length/2*zinput[1]+guideInfo.Length/2, 1263 guideInfo.Length/2*zinput[2]+guideInfo.Length/2, 1264 guideInfo.Length/2*zinput[3]+guideInfo.Length/2, 1265 guideInfo.Length/2*zinput[4]+guideInfo.Length/2, 1266 guideInfo.Length/2*zinput[5]+guideInfo.Length/2, 1267 guideInfo.Length/2*zinput[6]+guideInfo.Length/2, 1268 guideInfo.Length/2*zinput[7]+guideInfo.Length/2 }; 1269 int i = 0; 1270 for(i = 0; i < 8; i++){ 1271 guideInfo.xArray[i] = xarray[i]; 1272 guideInfo.yArray[i] = yarray[i]; 1273 guideInfo.zArray[i] = zarray[i]; 1274 } 1275 */ 1276 1277 guideInfo.exitVerticalWidth = 1278 2*sqrt(1 - ( (guideInfo.Length-guideInfo.ellipseMajorOffset[BottomSide]) 1279 *(guideInfo.Length-guideInfo.ellipseMajorOffset[BottomSide]) 1280 )/(guideInfo.ellipseMajorAxis[BottomSide] 1281 *guideInfo.ellipseMajorAxis[BottomSide] ) 1282 )*guideInfo.ellipseMinorAxis[BottomSide]; 1283 1284 guideInfo.exitHorizontalWidth = 1285 2*sqrt(1 - ( (guideInfo.Length-guideInfo.ellipseMajorOffset[RightSide]) 1286 *(guideInfo.Length-guideInfo.ellipseMajorOffset[RightSide]) 1287 )/(guideInfo.ellipseMajorAxis[RightSide] 1288 *guideInfo.ellipseMajorAxis[RightSide] ) 1289 )*guideInfo.ellipseMinorAxis[RightSide]; 1290 1291 //////////////////segmentation of m values 1292 1293 // Are the arrays empty? 1294 if(mvaluesright != NULL || mvaluesleft != NULL 1295 || mvaluestop != NULL || mvaluesbottom != NULL) 1296 { 1297 guideInfo.enableSegments = 1; 1298 1299 guideInfo.numberOfSegments = sizeof(mvaluesright)/sizeof(mvaluesright[0]); 1300 1301 //printf("Length is %i\n",guideInfo.numberOfSegments); 1302 guideInfo.mValuesright = mvaluesright; 1303 guideInfo.mValuesleft = mvaluesleft; 1304 guideInfo.mValuestop = mvaluestop; 1305 guideInfo.mValuesbottom = mvaluesbottom; 1306 //printf("Seglength ... %f %f %f\n",seglength[0],seglength[1],seglength[2]); 1307 1308 // Are the arrays of equal length? 1309 if(seglength == NULL){ 1310 dynamicalSegLength = 1311 realloc(dynamicalSegLength, 1312 guideInfo.numberOfSegments*sizeof(double) ); 1313 int i; 1314 for (i = 0; i < guideInfo.numberOfSegments; ++i){ 1315 dynamicalSegLength[i] = 1316 guideInfo.Length/guideInfo.numberOfSegments; 1317 } 1318 guideInfo.segLength = dynamicalSegLength; 1319 } 1320 else guideInfo.segLength = seglength; 1321 1322 if( guideInfo.numberOfSegments != sizeof(mvaluesright)/sizeof(mvaluesright[0]) 1323 || guideInfo.numberOfSegments != sizeof(mvaluesleft)/sizeof(mvaluesleft[0]) 1324 || guideInfo.numberOfSegments != sizeof(mvaluestop)/sizeof(mvaluestop[0]) 1325 || guideInfo.numberOfSegments != sizeof(mvaluesbottom)/sizeof(mvaluesbottom[0]) 1326 || (guideInfo.segLength == NULL 1327 & guideInfo.numberOfSegments != sizeof(seglength)/sizeof(guideInfo.segLength[0]) 1328 ) ) { 1329 1330 printf("Error in userinput inside %s, the length of the arrays" 1331 " mvalues and seglength are not equal\n",NAME_CURRENT_COMP); 1332 printf("The length of the arrays are: mValuesright is %lu," 1333 " mvaluesleft is %lu, mvaluestop is %lu, mvaluesbottom is" 1334 " %lu and seglength is %lu and should be %d \n; Above assume that the arrays are using double \n", 1335 sizeof(mvaluesright)/sizeof(double), 1336 sizeof(mvaluesleft)/sizeof(double), 1337 sizeof(mvaluestop)/sizeof(double), 1338 sizeof(mvaluesbottom)/sizeof(double), 1339 sizeof(guideInfo.segLength)/sizeof(double), 1340 guideInfo.numberOfSegments 1341 ); 1342 1343 if ( guideInfo.verboseSetting ) { 1344 int i; 1345 1346 printf("The Values of mvaluesright is: ["); 1347 for(i=0; i < sizeof(mvaluesright)/sizeof(mvaluesright[0]); i++) 1348 printf("%e,",guideInfo.mValuesright[i] ); 1349 printf("]\n"); 1350 1351 printf("The Values of mvaluesleft is: ["); 1352 for(i=0; i < sizeof(mvaluesleft)/sizeof(mvaluesleft[0]); i++) 1353 printf("%e,",guideInfo.mValuesleft[i] ); 1354 printf("]\n"); 1355 1356 printf("The Values of mvaluestop is: ["); 1357 for(i=0; i < sizeof(mvaluestop)/sizeof(mvaluestop[0]); i++) 1358 printf("%e,",guideInfo.mValuestop[i] ); 1359 printf("]\n"); 1360 1361 printf("The Values of mvaluesbottom is: ["); 1362 for(i=0; i < sizeof(mvaluesbottom)/sizeof(mvaluesbottom[0]); i++) 1363 printf("%e,",guideInfo.mValuesbottom[i] ); 1364 printf("]\n"); 1365 1366 printf("The Values of seglength is: ["); 1367 for(i=0; i < sizeof(guideInfo.segLength)/sizeof(guideInfo.segLength[0]); i++) 1368 printf("%e,",guideInfo.segLength[i]); 1369 printf("]\n"); 1370 } 1371 exit( printf("Exit due to critical error in userinput for the" 1372 " component %s, consider having a look at the input" 1373 " for following: mvaluesright,mvaluesleft,mvaluestop," 1374 "mvaluesbottom and/or seglength.",NAME_CURRENT_COMP) ); 1375 } 1376 // 1377 double sumOfelements=0; 1378 int i; 1379 for(i=0;i< guideInfo.numberOfSegments; i++) { 1380 sumOfelements += guideInfo.segLength[i]; 1381 } 1382 if ( guideInfo.verboseSetting 1383 && fabs(sumOfelements-guideInfo.Length) > 1e-9 ) 1384 printf("Error in userinput inside %s, the difference between" 1385 " guidelength and elements of the seglength array is:" 1386 "%e consider changes the parameters l or seglength \n", 1387 NAME_CURRENT_COMP,sumOfelements-guideInfo.Length); 1388 } 1389 else guideInfo.enableSegments = 0; 1390 1391 1392/////////////////////////////////////////////////////////////////////////// 1393/////////////// Calculate gravity vector in the guides coordinatesystem 1394/////////////////////////////////////////////////////////////////////////// 1395 1396/* 1397 Sets the local gravity vector equal to the global gravity vector (0,-g,0) 1398 and when apply the same rotation matrix as applied to guide. 1399*/ 1400 if (enableGravity != 0){ 1401 Gx0=0, Gy0=-GRAVITY*enableGravity, Gz0=0; 1402 Coords mcLocG; 1403 mcLocG = rot_apply(ROT_A_CURRENT_COMP, coords_set(0,Gy0,0)); 1404 coords_get(mcLocG, &Gx0, &Gy0, &Gz0); 1405 } 1406 Circ=2*PI*curvature; 1407%} 1408 1409/** 1410 Finds the next collision between the particle and the guide. 1411 Decided by the type of collision, the particle will then be either 1412 reflected on the guide, leaving the guide, or absorbed. 1413 If reflected then the particle will have its velocity vector and p 1414 changed appropriable (TODO). 1415 In the case of a absorbed particle the ABSORB function is called. 1416 and if the particle is found to have left the guide the break command 1417 is called and the end of trace is reached 1418*/ 1419TRACE 1420%{ 1421 PROP_Z0; 1422 SCATTER; 1423 1424 double Gloc; 1425 if (curvature) { 1426 Gloc=(vx*vx+vy*vy+vz*vz)/curvature; 1427 } else { 1428 Gloc=0; 1429 } 1430 1431 1432 if( !guideInfo.EnclosingBoxOn ) 1433 if( fabs(x) > guideInfo.entranceHorizontalWidth/2.0 1434 || fabs(y) > guideInfo.entranceVerticalWidth/2.0 ) 1435 ABSORB; 1436 1437 1438 int bounces = 0; 1439 for(bounces = 0; bounces <= 1000; bounces++){ 1440 1441 Gx=Gx0; Gy=Gy0; Gz=Gz0; 1442 if (curvature) { 1443 // Add velocity-dependent, location-dependent approximation to centripetal force for curvature... 1444 Gx=Gx0+Gloc*cos(2*PI*z/Circ);Gz=Gz0+Gloc*sin(2*PI*z/Circ); 1445 } 1446 1447 // Find the next intersection between the guide and the neutron. 1448 int boolean = guide_elliptical_handleGuideIntersection( 1449 x,y,z,vx,vy,vz,Gx,Gy,Gz, 1450 &guideInfo,&latestParticleCollision); 1451 1452 double timeToCollision = 1453 latestParticleCollision.delta_time_to_next_collision; 1454 1455 // Handle special cases. 1456 if(boolean == 0) ABSORB; 1457 if(timeToCollision < 1e-15) ABSORB; 1458 1459 // If the neutron reach the end of the guide, when move 1460 // the neutron to the end of guide and leave this component. 1461 if( z+vz*timeToCollision+0.5*Gz*timeToCollision*timeToCollision 1462 >= guideInfo.Length ) 1463 { 1464 double timeToExit = 0; 1465 solve_2nd_order( 1466 &timeToExit,NULL, 1467 -0.5*Gz,-vz,guideInfo.Length-z-1e-9); 1468 PROP_GRAV_DT(timeToExit,Gx,Gy,Gz); 1469 SCATTER; 1470 break; 1471 } 1472 1473 // Move the neutron and handle the reflection. 1474 PROP_GRAV_DT(timeToCollision,Gx,Gy,Gz); 1475 if(latestParticleCollision.collisionType == Absorb){ ABSORB; } 1476 if(latestParticleCollision.collisionType == Reflex){ 1477 p *= guide_elliptical_handleReflection(x,y,z,&vx,&vy,&vz, 1478 &guideInfo,&latestParticleCollision); 1479 SCATTER; 1480 if(p == 0) ABSORB; 1481 } 1482 } 1483 1484 if( fabs(x) > guideInfo.exitHorizontalWidth/2 1485 || fabs(y) > guideInfo.exitVerticalWidth/2 ) 1486 ABSORB; 1487 1488%} 1489 1490 1491FINALLY 1492%{ 1493%} 1494 1495MCDISPLAY 1496%{ 1497 1498 1499 // Calculate the points need to draw approximation of the ellipses 1500 // defining the guide 1501 1502 // the number of lines used to draw one side of the guide 1503 int ApproximationMirrors = 500; 1504 1505 // The start of the guide 1506 double zvalue=0; 1507 1508 // The the different in z between point used to draw the ellipse 1509 double zdelta = guideInfo.Length/(1.0*ApproximationMirrors); 1510 1511 // The vector used to store the points defining the lines 1512 double xplus[ApproximationMirrors+1]; 1513 double xminus[ApproximationMirrors+1]; 1514 double yplus[ApproximationMirrors+1]; 1515 double yminus[ApproximationMirrors+1]; 1516 1517 // Temperary values for the loop 1518 double tempx; 1519 double tempy; 1520 1521 /* 1522 Calculate the second coordinates to the points on the ellipse with z_i 1523 as the first coordinate. We transform the point to the coordinate system 1524 there the ellipse is a unit circle. And use the defination of this circle 1525 to find second coordinate (x^2+z^2 = 1) 1526 */ 1527 1528 ///////////////////////////////////////////////////////// 1529 double Length; 1530 double entranceHorizontalWidth; 1531 double entranceVerticalWidth; 1532 1533 // ellipses infomation 1534 double ellipseMajorAxis[4], ellipseMinorAxis[4]; 1535 double ellipseMajorOffset[4], ellipseMinorOffset[4]; 1536 1537 enum Side {RightSide,TopSide,LeftSide,BottomSide,None}; 1538 ///////////////////////////////////////////////////////// 1539 1540 int i = 0; 1541 double tempz = 0; 1542 for(i=0;i<ApproximationMirrors+1;i++){ 1543 1544 tempx = sqrt( 1545 guideInfo.ellipseMinorAxis[RightSide] 1546 *guideInfo.ellipseMinorAxis[RightSide] 1547 -( guideInfo.ellipseMinorAxis[RightSide] 1548 *guideInfo.ellipseMinorAxis[RightSide] ) 1549 /( guideInfo.ellipseMajorAxis[RightSide] 1550 *guideInfo.ellipseMajorAxis[RightSide] ) 1551 *( zvalue+zdelta*i-guideInfo.ellipseMajorOffset[RightSide] ) 1552 *( zvalue+zdelta*i-guideInfo.ellipseMajorOffset[RightSide] ) 1553 ); 1554 1555 xplus[i] = tempx + guideInfo.ellipseMinorOffset[RightSide]; 1556 xminus[i]= -tempx + guideInfo.ellipseMinorOffset[RightSide]; 1557 1558 tempy = sqrt( 1559 guideInfo.ellipseMinorAxis[TopSide] 1560 *guideInfo.ellipseMinorAxis[TopSide] 1561 -( guideInfo.ellipseMinorAxis[TopSide] 1562 *guideInfo.ellipseMinorAxis[TopSide] ) 1563 /( guideInfo.ellipseMajorAxis[TopSide] 1564 *guideInfo.ellipseMajorAxis[TopSide] ) 1565 *( zvalue+zdelta*i-guideInfo.ellipseMajorOffset[TopSide] ) 1566 *( zvalue+zdelta*i-guideInfo.ellipseMajorOffset[TopSide] ) 1567 ); 1568 1569 yplus[i] = tempy + guideInfo.ellipseMinorOffset[TopSide]; 1570 yminus[i]= -tempy + guideInfo.ellipseMinorOffset[TopSide]; 1571 } 1572 1573 ///// Draw lines 1574 1575 // Drawing lines orthogonal with the z direction. 1576 // at both ends of the guide and at the boardest place at the guide 1577 1578 // These may not give correct result if one of the ends are closed 1579 1580 int j=0; 1581 1582 line( xplus[j], yplus[j], zvalue+j*zdelta,0 , yplus[j], zvalue+j*zdelta); 1583 line( 0 , yplus[j], zvalue+j*zdelta,xminus[j], yplus[j], zvalue+j*zdelta); 1584 line( xminus[j], yminus[j], zvalue+j*zdelta,0 , yminus[j], zvalue+j*zdelta); 1585 line( 0, yminus[j], zvalue+j*zdelta,xplus[j], yminus[j], zvalue+j*zdelta); 1586 line( xminus[j],yplus[j], zvalue+j*zdelta, xminus[j],0 , zvalue+j*zdelta); 1587 line( xminus[j],0 , zvalue+j*zdelta,xminus[j],yminus[j], zvalue+j*zdelta); 1588 line( xplus[j], 0, zvalue+j*zdelta, xplus[j], yplus[j], zvalue+j*zdelta); 1589 line( xplus[j], yminus[j], zvalue+j*zdelta,xplus[j],0 , zvalue+j*zdelta); 1590 1591 j=ApproximationMirrors; 1592 1593 line( xplus[j], yplus[j], zvalue+j*zdelta,0 , yplus[j], zvalue+j*zdelta); 1594 line( 0 , yplus[j],zvalue+j*zdelta,xminus[j] , yplus[j], zvalue+j*zdelta); 1595 line( xminus[j], yminus[j], zvalue+j*zdelta,0 , yminus[j], zvalue+j*zdelta); 1596 line( 0, yminus[j], zvalue+j*zdelta, xplus[j], yminus[j], zvalue+j*zdelta); 1597 line( xminus[j],yplus[j], zvalue+j*zdelta, xminus[j],0 , zvalue+j*zdelta); 1598 line( xminus[j],0 , zvalue+j*zdelta,xminus[j],yminus[j], zvalue+j*zdelta); 1599 line( xplus[j], 0, zvalue+j*zdelta, xplus[j], yplus[j], zvalue+j*zdelta); 1600 line( xplus[j], yminus[j], zvalue+j*zdelta,xplus[j], 0 , zvalue+j*zdelta); 1601 1602 1603 1604 // find boardest place on the guide and draw a band around the guide 1605 int m0; 1606 double boardestPlace = 0; 1607 int boardestPlaceNumber = 0; 1608 for(m0=0; m0<ApproximationMirrors; m0++){ 1609 if( boardestPlace <= fabs(yplus[m0]) ){ 1610 boardestPlace = fabs(yplus[m0]); 1611 boardestPlaceNumber = m0; 1612 } 1613 } 1614 j = boardestPlaceNumber; 1615 1616 line( xplus[j], yplus[j], zvalue+j*zdelta,0 , yplus[j], zvalue+j*zdelta); 1617 line( 0 , yplus[j], zvalue+j*zdelta,xminus[j], yplus[j], zvalue+j*zdelta); 1618 line( xminus[j], yminus[j], zvalue+j*zdelta,0 , yminus[j], zvalue+j*zdelta); 1619 line( 0, yminus[j], zvalue+j*zdelta, xplus[j], yminus[j], zvalue+j*zdelta); 1620 line( xminus[j],yplus[j], zvalue+j*zdelta, xminus[j],0 , zvalue+j*zdelta); 1621 line( xminus[j],0 , zvalue+j*zdelta, xminus[j],yminus[j], zvalue+j*zdelta); 1622 line( xplus[j], 0, zvalue+j*zdelta, xplus[j], yplus[j], zvalue+j*zdelta); 1623 line( xplus[j], yminus[j], zvalue+j*zdelta, xplus[j], 0 , zvalue+j*zdelta); 1624 1625 1626 1627 // Drawing lines parallel with the z direction 1628 1629 int k=0; 1630 for(k=0; k < ApproximationMirrors; k++){ 1631 1632 line( xplus[k], yplus[k], zvalue+k*zdelta,xplus[k+1], yplus[k+1], zvalue+(k+1)*zdelta); 1633 line( xminus[k],yplus[k], zvalue+k*zdelta,xminus[k+1],yplus[k+1], zvalue+(k+1)*zdelta); 1634 1635 line( xplus[k], yminus[k],zvalue+k*zdelta,xplus[k+1], yminus[k+1],zvalue+(k+1)*zdelta); 1636 1637 line( xminus[k],yminus[k],zvalue+k*zdelta,xminus[k+1],yminus[k+1],zvalue+(k+1)*zdelta); 1638 1639 line( xminus[k],0 , zvalue+k*zdelta, xminus[k+1],0 , zvalue+(k+1)*zdelta); 1640 line( xplus[k], 0 , zvalue+k*zdelta, xplus[k+1], 0 , zvalue+(k+1)*zdelta); 1641 1642 line( 0 , yminus[k], zvalue+k*zdelta, 0 ,yminus[k+1],zvalue+(k+1)*zdelta); 1643 line( 0 ,yplus[k] , zvalue+k*zdelta , 0 ,yplus[k] , zvalue+(k+1)*zdelta); 1644 1645 } 1646 1647 if(guideInfo.EnclosingBoxOn){ 1648 dashed_line( guideInfo.xArray[0],guideInfo.yArray[0],guideInfo.zArray[0], 1649 guideInfo.xArray[1],guideInfo.yArray[1],guideInfo.zArray[1],10 ); 1650 dashed_line( guideInfo.xArray[1],guideInfo.yArray[1],guideInfo.zArray[1], 1651 guideInfo.xArray[2],guideInfo.yArray[2],guideInfo.zArray[2],10 ); 1652 dashed_line( guideInfo.xArray[2],guideInfo.yArray[2],guideInfo.zArray[2], 1653 guideInfo.xArray[3],guideInfo.yArray[3],guideInfo.zArray[3],10 ); 1654 dashed_line( guideInfo.xArray[3],guideInfo.yArray[3],guideInfo.zArray[3], 1655 guideInfo.xArray[0],guideInfo.yArray[0],guideInfo.zArray[0],10 ); 1656 1657 dashed_line( guideInfo.xArray[4],guideInfo.yArray[4],guideInfo.zArray[4], 1658 guideInfo.xArray[5],guideInfo.yArray[5],guideInfo.zArray[5],10 ); 1659 dashed_line( guideInfo.xArray[5],guideInfo.yArray[5],guideInfo.zArray[5], 1660 guideInfo.xArray[6],guideInfo.yArray[6],guideInfo.zArray[6],10 ); 1661 dashed_line( guideInfo.xArray[6],guideInfo.yArray[6],guideInfo.zArray[6], 1662 guideInfo.xArray[7],guideInfo.yArray[7],guideInfo.zArray[7],10 ); 1663 dashed_line( guideInfo.xArray[7],guideInfo.yArray[7],guideInfo.zArray[7], 1664 guideInfo.xArray[4],guideInfo.yArray[4],guideInfo.zArray[4],10 ); 1665 1666 dashed_line( guideInfo.xArray[0],guideInfo.yArray[0],guideInfo.zArray[0], 1667 guideInfo.xArray[4],guideInfo.yArray[4],guideInfo.zArray[4],10 ); 1668 dashed_line( guideInfo.xArray[4],guideInfo.yArray[4],guideInfo.zArray[4], 1669 guideInfo.xArray[7],guideInfo.yArray[7],guideInfo.zArray[7],10 ); 1670 dashed_line( guideInfo.xArray[7],guideInfo.yArray[7],guideInfo.zArray[7], 1671 guideInfo.xArray[3],guideInfo.yArray[3],guideInfo.zArray[3],10 ); 1672 dashed_line( guideInfo.xArray[3],guideInfo.yArray[3],guideInfo.zArray[3], 1673 guideInfo.xArray[0],guideInfo.yArray[0],guideInfo.zArray[0],10 ); 1674 1675 dashed_line( guideInfo.xArray[1],guideInfo.yArray[1],guideInfo.zArray[1], 1676 guideInfo.xArray[5],guideInfo.yArray[5],guideInfo.zArray[5],10 ); 1677 dashed_line( guideInfo.xArray[5],guideInfo.yArray[5],guideInfo.zArray[5], 1678 guideInfo.xArray[6],guideInfo.yArray[6],guideInfo.zArray[6],10 ); 1679 dashed_line( guideInfo.xArray[6],guideInfo.yArray[6],guideInfo.zArray[6], 1680 guideInfo.xArray[2],guideInfo.yArray[2],guideInfo.zArray[2],10 ); 1681 dashed_line( guideInfo.xArray[2],guideInfo.yArray[2],guideInfo.zArray[2], 1682 guideInfo.xArray[1],guideInfo.yArray[1],guideInfo.zArray[1],10 ); 1683 } 1684%} 1685 1686END 1687