source: trunk/source/geometry/solids/BREPS/src/G4CylindricalSurface.cc@ 1202

Last change on this file since 1202 was 1058, checked in by garnier, 17 years ago

file release beta
File size: 30.9 KB
RevLine 
[831]1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4CylindricalSurface.cc,v 1.8 2006/06/29 18:42:08 gunter Exp $
[1058]28// GEANT4 tag $Name: geant4-09-02-ref-02 $
[831]29//
30// ----------------------------------------------------------------------
31// GEANT 4 class source file
32//
33// G4CylindricalSurface.cc
34//
35// ----------------------------------------------------------------------
36
37#include "G4CylindricalSurface.hh"
38#include "G4Sort.hh"
39#include "G4Globals.hh"
40
41G4CylindricalSurface::G4CylindricalSurface() : G4Surface()
42{
43 // default constructor
44 // default axis is ( 1.0, 0.0, 0.0 ), default radius is 1.0
45 axis = G4Vector3D( 1.0, 0.0, 0.0 );
46 radius = 1.0;
47}
48
49
50G4CylindricalSurface::G4CylindricalSurface( const G4Vector3D& o,
51 const G4Vector3D& a,
52 G4double r ) //: G4Surface( o )
53{
54 // Normal constructor
55 // require axis to be a unit vector
56 G4double amag = a.mag();
57 if ( amag != 0.0 )
58 axis = a * (1/ amag); // this makes the axis a unit vector
59 else
60 {
61 G4cerr << "Error in G4CylindricalSurface::G4CylindricalSurface--axis "
62 <<"has zero length\n"
63 << "\tDefault axis ( 1.0, 0.0, 0.0 ) is used.\n";
64
65 axis = G4Vector3D( 1.0, 0.0, 0.0 );
66 }
67
68 // Require radius to be non-negative
69 if ( r >= 0.0 )
70 radius = r;
71 else
72 {
73 G4cerr << "Error in G4CylindricalSurface::G4CylindricalSurface"
74 << "--asked for negative radius\n"
75 << "\tDefault radius of 1.0 is used.\n";
76
77 radius = 1.0;
78 }
79
80 origin =o;
81}
82
83
84G4CylindricalSurface::~G4CylindricalSurface()
85{
86}
87
88/*
89G4CylindricalSurface::G4CylindricalSurface( const G4CylindricalSurface& c )
90 : G4Surface( c.origin )
91{
92 axis = c.axis; radius = c.radius;
93}
94*/
95
96const char* G4CylindricalSurface::NameOf() const
97{
98 return "G4CylindricalSurface";
99}
100
101void G4CylindricalSurface::PrintOn( std::ostream& os ) const
102{
103 // printing function using C++ std::ostream class
104 os << "G4CylindricalSurface surface with origin: " << origin << "\t"
105 << "radius: " << radius << "\tand axis " << axis << "\n";
106}
107
108//G4int G4Surface::Intersect(const G4Ray& ry)
109G4int G4CylindricalSurface::Intersect(const G4Ray& ry)
110{
111
112 // L. Broglia : copy of G4FCylindricalSurface::Intersect
113
114 // Distance along a Ray (straight line with G4ThreeVec) to leave or enter
115 // a G4CylindricalSurface. The input variable which_way should be set
116 // to +1 to indicate leaving a G4CylindricalSurface, -1 to indicate
117 // entering a G4CylindricalSurface.
118 // p is the point of intersection of the Ray with the G4CylindricalSurface.
119 // If the G4Vector3D of the Ray is opposite to that of the Normal to
120 // the G4CylindricalSurface at the intersection point, it will not leave
121 // the G4CylindricalSurface.
122 // Similarly, if the G4Vector3D of the Ray is along that of the Normal
123 // to the G4CylindricalSurface at the intersection point, it will not enter
124 // the G4CylindricalSurface.
125 // This method is called by all finite shapes sub-classed to
126 // G4CylindricalSurface.
127 // Use the virtual function table to check if the intersection point
128 // is within the boundary of the finite shape.
129 // A negative result means no intersection.
130 // If no valid intersection point is found, set the distance
131 // and intersection point to large numbers.
132
133 // G4int which_way = -1;
134 //Originally a parameter.Read explanation above.
135
136 G4int which_way=1;
137
138 if(!Inside(ry.GetStart()))
139 which_way = -1;
140
141 distance = FLT_MAXX;
142 G4Vector3D lv ( FLT_MAXX, FLT_MAXX, FLT_MAXX );
143
144 closest_hit = lv;
145
146 // Origin and G4Vector3D unit vector of Ray.
147 G4Vector3D x = ry.GetStart();
148 G4Vector3D dhat = ry.GetDir();
149
150 // Axis unit vector of the G4CylindricalSurface.
151 G4Vector3D ahat = GetAxis();
152 G4int isoln = 0,
153 maxsoln = 2;
154
155 // array of solutions in distance along the Ray
156 G4double s[2];
157 s[0] = -1.0;
158 s[1] = -1.0 ;
159
160 // calculate the two solutions (quadratic equation)
161 G4Vector3D d = x - GetOrigin();
162 G4double radiu = GetRadius();
163
164 //quit with no intersection if the radius of the G4CylindricalSurface is zero
165 // if ( radiu <= 0.0 )
166 // return 0;
167
168 G4double dsq = d * d;
169 G4double da = d * ahat;
170 G4double dasq = da * da;
171 G4double rsq = radiu * radiu;
172 G4double qsq = dsq - dasq;
173 G4double dira = dhat * ahat;
174 G4double a = 1.0 - dira * dira;
175
176 if ( a <= 0.0 )
177 return 0;
178
179 G4double b = 2. * ( d * dhat - da * dira );
180 G4double c = rsq - qsq;
181 G4double radical = b * b + 4. * a * c;
182
183 if ( radical < 0.0 )
184 return 0;
185
186 G4double root = std::sqrt( radical );
187 s[0] = ( - b + root ) / ( 2. * a );
188 s[1] = ( - b - root ) / ( 2. * a );
189
190 // order the possible solutions by increasing distance along the Ray
191 // (G4Sorting routines are in support/G4Sort.h)
192 sort_double( s, isoln, maxsoln-1 );
193
194 // now loop over each positive solution, keeping the first one (smallest
195 // distance along the Ray) which is within the boundary of the sub-shape
196 // and which also has the correct G4Vector3D with respect to the Normal to
197 // the G4CylindricalSurface at the intersection point
198 for ( isoln = 0; isoln < maxsoln; isoln++ )
199 {
200 if ( s[isoln] >= kCarTolerance*0.5 )
201 {
202 if ( s[isoln] >= FLT_MAXX ) // quit if too large
203 return 0;
204
205 distance = s[isoln];
206 closest_hit = ry.GetPoint( distance );
207 G4double tmp = dhat * (Normal( closest_hit ));
208
209 if ((tmp * which_way) >= 0.0 )
210 if ( WithinBoundary( closest_hit ) == 1 )
211 distance = distance*distance;
212
213 return 1;
214 }
215 }
216
217 // get here only if there was no solution within the boundary, Reset
218 // distance and intersection point to large numbers
219 distance = FLT_MAXX;
220 closest_hit = lv;
221
222
223 return 0;
224}
225
226
227
228
229G4double G4CylindricalSurface::HowNear( const G4Vector3D& x ) const
230{
231 // Distance from the point x to the infinite G4CylindricalSurface.
232 // The distance will be positive if the point is Inside the
233 // G4CylindricalSurface, negative if the point is outside.
234 // Note that this may not be correct for a bounded cylindrical object
235 // subclassed to G4CylindricalSurface.
236 G4Vector3D d = x - origin;
237 G4double dA = d * axis;
238 G4double rad = std::sqrt( d.mag2() - dA*dA );
239 G4double hownear = std::fabs( radius - rad );
240
241 return hownear;
242}
243
244
245/*
246G4double G4CylindricalSurface::distanceAlongRay( G4int which_way, const G4Ray* ry,
247 G4Vector3D& p ) const
248{ // Distance along a Ray (straight line with G4Vector3D) to leave or enter
249 // a G4CylindricalSurface. The input variable which_way should be set to +1
250 // to indicate leaving a G4CylindricalSurface, -1 to indicate entering a
251 // G4CylindricalSurface.
252 // p is the point of intersection of the Ray with the G4CylindricalSurface.
253 // If the G4Vector3D of the Ray is opposite to that of the Normal to
254 // the G4CylindricalSurface at the intersection point, it will not leave the
255 // G4CylindricalSurface.
256 // Similarly, if the G4Vector3D of the Ray is along that of the Normal
257 // to the G4CylindricalSurface at the intersection point, it will not enter
258 // the G4CylindricalSurface.
259 // This method is called by all finite shapes sub-classed to
260 // G4CylindricalSurface.
261 // Use the virtual function table to check if the intersection point
262 // is within the boundary of the finite shape.
263 // A negative result means no intersection.
264 // If no valid intersection point is found, set the distance
265 // and intersection point to large numbers.
266 G4double Dist = FLT_MAXX;
267 G4Vector3D lv ( FLT_MAXX, FLT_MAXX, FLT_MAXX );
268 p = lv;
269// Origin and G4Vector3D unit vector of Ray.
270 G4Vector3D x = ry->Position();
271 G4Vector3D dhat = ry->Direction( 0.0 );
272// Axis unit vector of the G4CylindricalSurface.
273 G4Vector3D ahat = GetAxis();
274 G4int isoln = 0, maxsoln = 2;
275// array of solutions in distance along the Ray
276// G4double s[2] = { -1.0, -1.0 };
277 G4double s[2];s[0] = -1.0; s[1]= -1.0 ;
278// calculate the two solutions (quadratic equation)
279 G4Vector3D d = x - GetOrigin();
280 G4double radius = GetRadius();
281// quit with no intersection if the radius of the G4CylindricalSurface is zero
282 if ( radius <= 0.0 )
283 return Dist;
284 G4double dsq = d * d;
285 G4double da = d * ahat;
286 G4double dasq = da * da;
287 G4double rsq = radius * radius;
288 G4double qsq = dsq - dasq;
289 G4double dira = dhat * ahat;
290 G4double a = 1.0 - dira * dira;
291 if ( a <= 0.0 )
292 return Dist;
293 G4double b = 2. * ( d * dhat - da * dira );
294 G4double c = rsq - qsq;
295 G4double radical = b * b + 4. * a * c;
296 if ( radical < 0.0 )
297 return Dist;
298 G4double root = std::sqrt( radical );
299 s[0] = ( - b + root ) / ( 2. * a );
300 s[1] = ( - b - root ) / ( 2. * a );
301// order the possible solutions by increasing distance along the Ray
302// (G4Sorting routines are in support/G4Sort.h)
303 G4Sort_double( s, isoln, maxsoln-1 );
304// now loop over each positive solution, keeping the first one (smallest
305// distance along the Ray) which is within the boundary of the sub-shape
306// and which also has the correct G4Vector3D with respect to the Normal to
307// the G4CylindricalSurface at the intersection point
308 for ( isoln = 0; isoln < maxsoln; isoln++ ) {
309 if ( s[isoln] >= 0.0 ) {
310 if ( s[isoln] >= FLT_MAXX ) // quit if too large
311 return Dist;
312 Dist = s[isoln];
313 p = ry->Position( Dist );
314 if ( ( ( dhat * Normal( p ) * which_way ) >= 0.0 )
315 && ( WithinBoundary( p ) == 1 ) )
316 return Dist;
317 }
318 }
319// get here only if there was no solution within the boundary, Reset
320// distance and intersection point to large numbers
321 p = lv;
322 return FLT_MAXX;
323}
324
325*/
326
327
328/*
329
330
331G4double G4CylindricalSurface::distanceAlongHelix( G4int which_way,
332 const Helix* hx,
333 G4Vector3D& p ) const
334{ // Distance along a Helix to leave or enter a G4CylindricalSurface.
335 // The input variable which_way should be set to +1 to
336 // indicate leaving a G4CylindricalSurface, -1 to indicate entering a
337 // G4CylindricalSurface.
338 // p is the point of intersection of the Helix with the G4CylindricalSurface.
339 // If the G4Vector3D of the Helix is opposite to that of the Normal to
340 // the G4CylindricalSurface at the intersection point, it will not leave the
341 // G4CylindricalSurface.
342 // Similarly, if the G4Vector3D of the Helix is along that of the Normal
343 // to the G4CylindricalSurface at the intersection point, it will not enter
344 // the G4CylindricalSurface.
345 // This method is called by all finite shapes sub-classed to
346 // G4CylindricalSurface.
347 // Use the virtual function table to check if the intersection point
348 // is within the boundary of the finite shape.
349 // If no valid intersection point is found, set the distance
350 // and intersection point to large numbers.
351 // Possible negative distance solutions are discarded.
352 G4double Dist = FLT_MAXX;
353 G4Vector3D lv ( FLT_MAXX, FLT_MAXX, FLT_MAXX );
354 G4Vector3D zerovec; // zero G4Vector3D
355 p = lv;
356 G4int isoln = 0, maxsoln = 4;
357// Array of solutions in turning angle
358// G4double s[4] = { -1.0, -1.0, -1.0, -1.0 };
359 G4double s[4];s[0]=-1.0;s[1]= -1.0;s[2]= -1.0;s[3]= -1.0;
360
361
362// Flag set to 1 if exact solution is found
363 G4int exact = 0;
364// Helix parameters
365 G4double rh = hx->GetRadius(); // radius of Helix
366 G4Vector3D ah = hx->GetAxis(); // axis of Helix
367 G4Vector3D oh = hx->position(); // origin of Helix
368 G4Vector3D dh = hx->direction( 0.0 ); // initial G4Vector3D of Helix
369 G4Vector3D prp = hx->getPerp(); // perpendicular vector
370 G4double prpmag = prp.Magnitude();
371 G4double rhp = rh / prpmag;
372// G4CylindricalSurface parameters
373 G4double rc = GetRadius(); // radius of G4CylindricalSurface
374 if ( rc == 0.0 ) // quit if zero radius
375 return Dist;
376 G4Vector3D oc = GetOrigin(); // origin of G4CylindricalSurface
377 G4Vector3D ac = GetAxis(); // axis of G4CylindricalSurface
378//
379// Calculate quantities of use later on.
380 G4Vector3D alpha = rhp * prp;
381 G4Vector3D beta = rhp * dh;
382 G4Vector3D gamma = oh - oc;
383// Declare variables used later on in several places.
384 G4double rcd2 = 0.0, alpha2 = 0.0;
385 G4double A = 0.0, B = 0.0, C = 0.0, F = 0.0, G = 0.0, H = 0.0;
386 G4double CoverB = 0.0, radical = 0.0, root = 0.0, s1 = 0.0, s2 = 0.0;
387 G4Vector3D ghat;
388//
389// Set flag for special cases
390 G4int special_case = 0; // 0 means general case
391//
392// Test to see if axes of Helix and G4CylindricalSurface are parallel, in which
393// case there are exact solutions.
394 if ( ( std::fabs( ah.AngleBetween(ac) ) < FLT_EPSILO )
395 || ( std::fabs( ah.AngleBetween(ac) - pi ) < FLT_EPSILO ) ) {
396 special_case = 1;
397// If, in addition, gamma is a zero vector or is parallel to the
398// G4CylindricalSurface axis, this simplifies the previous case.
399 if ( gamma == zerovec ) {
400 special_case = 3;
401 ghat = gamma;
402 }
403 else {
404 ghat = gamma / gamma.Magnitude();
405 if ( ( std::fabs( ghat.AngleBetween(ac) ) < FLT_EPSILO )
406 || ( std::fabs( ghat.AngleBetween(ac) - pi ) <
407 FLT_EPSILO ) )
408 special_case = 3;
409 }
410// Test to see if, in addition to the axes of the Helix and G4CylindricalSurface
411// being parallel, the axis of the G4CylindricalSurface is perpendicular to the
412// initial G4Vector3D of the Helix.
413 if ( std::fabs( ( ac * dh ) ) < FLT_EPSILO ) {
414// And, if, in addition to all this, the difference in origins of the Helix
415// and G4CylindricalSurface is perpendicular to the initial G4Vector3D of the
416// Helix, there is a separate special case.
417 if ( std::fabs( ( ghat * dh ) ) < FLT_EPSILO )
418 special_case = 4;
419 }
420 } // end of section with axes of Helix and G4CylindricalSurface parallel
421//
422// Another peculiar case occurs if the axis of the G4CylindricalSurface and the
423// initial G4Vector3D of the Helix line up and their origins are the same.
424// This will require a higher order approximation than the general case.
425 if ( ( ( std::fabs( dh.AngleBetween(ac) ) < FLT_EPSILO )
426 || ( std::fabs( dh.AngleBetween(ac) - pi ) < FLT_EPSILO ) )
427 && ( gamma == zerovec ) )
428 special_case = 2;
429//
430// Now all the special cases have been tagged, so solutions are found
431// for each case. Exact solutions are indicated by setting exact = 1.
432// [For some reason switch doesn't work here, so do series of if's.]
433 if ( special_case == 0 ) { // approximate quadratic solutions
434 A = beta * beta - ( beta * ac ) * ( beta * ac )
435 + gamma * alpha - ( gamma * ac ) * ( alpha * ac );
436 B = 2.0 * gamma * beta
437 - 2.0 * ( gamma * ac ) * ( beta * ac );
438 C = gamma * gamma
439 - ( gamma * ac ) * ( gamma * ac ) - rc * rc;
440 if ( std::fabs( A ) < FLT_EPSILO ) { // no quadratic term
441 if ( B == 0.0 ) // no intersection, quit
442 return Dist;
443 else // B != 0
444 s[0] = -C / B;
445 }
446 else { // A != 0, general quadratic solution
447 radical = B * B - 4.0 * A * C;
448 if ( radical < 0.0 ) // no solution, quit
449 return Dist;
450 root = std::sqrt( radical );
451 s[0] = ( -B + root ) / ( 2.0 * A );
452 s[1] = ( -B - root ) / ( 2.0 * A );
453 if ( rh < 0.0 ) {
454 s[0] = -s[0];
455 s[1] = -s[1];
456 }
457 s[2] = s[0] + 2.0 * pi;
458 s[3] = s[1] + 2.0 * pi;
459 }
460 }
461//
462 else if ( special_case == 1 ) { // exact solutions
463 exact = 1;
464 H = 2.0 * ( alpha * alpha + gamma * alpha );
465 F = gamma * gamma
466 - ( ( gamma * ac ) * ( gamma * ac ) )
467 - rc * rc + H;
468 G = 2.0 * rhp *
469 ( gamma * dh - ( gamma * ac ) * ( ac * dh ) );
470 A = G * G + H * H;
471 B = -2.0 * F * H;
472 C = F * F - G * G;
473 if ( std::fabs( A ) < FLT_EPSILO ) { // no quadratic term
474 if ( B == 0.0 ) // no intersection, quit
475 return Dist;
476 else { // B != 0
477 CoverB = -C / B;
478 if ( std::fabs( CoverB ) > 1.0 )
479 return Dist;
480 s[0] = std::acos( CoverB );
481 }
482 }
483 else { // A != 0, general quadratic solution
484// Try a different method of calculation using F, G, and H to avoid
485// precision problems.
486// radical = B * B - 4.0 * A * C;
487// if ( radical < 0.0 ) {
488 if ( std::fabs( H ) > FLT_EPSILO ) {
489 G4double r1 = G / H;
490 G4double r2 = F / H;
491 G4double radsq = 1.0 + r1*r1 - r2*r2;
492 if ( radsq < 0.0 )
493 return Dist;
494 root = G * std::sqrt( radsq );
495 G4double denominator = H * ( 1.0 + r1*r1 );
496 s1 = ( F + root ) / denominator;
497 s2 = ( F - root ) / denominator;
498 }
499 else
500 return Dist;
501// } // end radical < 0 condition
502// else {
503// root = std::sqrt( radical );
504// s1 = ( -B + root ) / ( 2.0 * A );
505// s2 = ( -B - root ) / ( 2.0 * A );
506// }
507 if ( std::fabs( s1 ) <= 1.0 ) {
508 s[0] = std::acos( s1 );
509 s[2] = 2.0 * pi - s[0];
510 }
511 if ( std::fabs( s2 ) <= 1.0 ) {
512 s[1] = std::acos( s2 );
513 s[3] = 2.0 * pi - s[1];
514 }
515// Must take only solutions which satisfy original unsquared equation:
516// Gsin(s) - Hcos(s) + F = 0. Take best solution of pair and set false
517// solutions to -1. Only do this if the result is significantly different
518// from zero.
519 G4double temp1 = 0.0, temp2 = 0.0;
520 G4double rsign = 1.0;
521 if ( rh < 0.0 ) rsign = -1.0;
522 if ( s[0] > 0.0 ) {
523 temp1 = G * rsign * std::sin( s[0] )
524 - H * std::cos( s[0] ) + F;
525 temp2 = G * rsign * std::sin( s[2] )
526 - H * std::cos( s[2] ) + F;
527 if ( std::fabs( temp1 ) > std::fabs( temp2 ) )
528 if ( std::fabs( temp1 ) > FLT_EPSILO )
529 s[0] = -1.0;
530 else
531 if ( std::fabs( temp2 ) > FLT_EPSILO )
532 s[2] = -1.0;
533 }
534 if ( s[1] > 0.0 ) {
535 temp1 = G * rsign * std::sin( s[1] )
536 - H * std::cos( s[1] ) + F;
537 temp2 = G * rsign * std::sin( s[3] )
538 - H * std::cos( s[3] ) + F;
539 if ( std::fabs( temp1 ) > std::fabs( temp2 ) )
540 if ( std::fabs( temp1 ) > FLT_EPSILO )
541 s[1] = -1.0;
542 else
543 if ( std::fabs( temp2 ) > FLT_EPSILO )
544 s[3] = -1.0;
545 }
546 }
547 }
548//
549 else if ( special_case == 2 ) { // approximate solution
550 G4Vector3D e = ah.cross( ac );
551 G4double re = std::fabs( rhp ) * e.Magnitude();
552 s[0] = std::sqrt( 2.0 * rc / re );
553 }
554//
555 else if ( special_case == 3 ) { // exact solutions
556 exact = 1;
557 alpha2 = alpha * alpha;
558 rcd2 = rhp * rhp * ( 1.0 - ( (ac*dh) * (ac*dh) ) );
559 A = alpha2 - rcd2;
560 B = - 2.0 * alpha2;
561 C = alpha2 + rcd2 - rc*rc;
562 if ( std::fabs( A ) < FLT_EPSILO ) { // no quadratic term
563 if ( B == 0.0 ) // no intersection, quit
564 return Dist;
565 else { // B != 0
566 CoverB = -C / B;
567 if ( std::fabs( CoverB ) > 1.0 )
568 return Dist;
569 s[0] = std::acos( CoverB );
570 }
571 }
572 else { // A != 0, general quadratic solution
573 radical = B * B - 4.0 * A * C;
574 if ( radical < 0.0 )
575 return Dist;
576 root = std::sqrt( radical );
577 s1 = ( -B + root ) / ( 2.0 * A );
578 s2 = ( -B - root ) / ( 2.0 * A );
579 if ( std::fabs( s1 ) <= 1.0 )
580 s[0] = std::acos( s1 );
581 if ( std::fabs( s2 ) <= 1.0 )
582 s[1] = std::acos( s2 );
583 }
584 }
585//
586 else if ( special_case == 4 ) { // exact solution
587 exact = 1;
588 F = gamma * gamma
589 - ( ( gamma * ac ) * ( gamma * ac ) )
590 - rc * rc;
591 G = 2.0 * ( rhp * rhp + gamma * alpha );
592 if ( G == 0.0 ) // no intersection, quit
593 return Dist;
594 G4double cs = 1.0 + ( F / G );
595 if ( std::fabs( cs ) > 1.0 ) // no intersection, quit
596 return Dist;
597 s[0] = std::acos( cs );
598 }
599//
600 else // shouldn't get here
601 return Dist;
602//
603// **************************************************************************
604//
605// Order the possible solutions by increasing turning angle
606// (G4Sorting routines are in support/G4Sort.h).
607 G4Sort_double( s, isoln, maxsoln-1 );
608//
609// Now loop over each positive solution, keeping the first one (smallest
610// distance along the Helix) which is within the boundary of the sub-shape.
611 for ( isoln = 0; isoln < maxsoln; isoln++ ) {
612 if ( s[isoln] >= 0.0 ) {
613// Calculate distance along Helix and position and G4Vector3D vectors.
614 Dist = s[isoln] * std::fabs( rhp );
615 p = hx->position( Dist );
616 G4Vector3D d = hx->direction( Dist );
617 if ( exact == 0 ) { // only for approximate solns
618// Now do approximation to get remaining distance to correct this solution
619// iterate it until the accuracy is below the user-set surface precision.
620 G4double delta = 0.0;
621 G4double delta0 = FLT_MAXX;
622 G4int dummy = 1;
623 G4int iter = 0;
624 G4int in0 = Inside( hx->position ( 0.0 ) );
625 G4int in1 = Inside( p );
626 G4double sc = Scale();
627 while ( dummy ) {
628 iter++;
629// Terminate loop after 50 iterations and Reset distance to large number,
630// indicating no intersection with G4CylindricalSurface.
631// This generally occurs if the Helix curls too tightly to Intersect it.
632 if ( iter > 50 ) {
633 Dist = FLT_MAXX;
634 p = lv;
635 break;
636 }
637// Find distance from the current point along the above-calculated
638// G4Vector3D using a Ray.
639// The G4Vector3D of the Ray and the Sign of the distance are determined
640// by whether the starting point of the Helix is Inside or outside of
641// the G4CylindricalSurface.
642 in1 = Inside( p );
643 if ( in1 ) { // current point Inside
644 if ( in0 ) { // starting point Inside
645 Ray* r = new Ray( p, d );
646 delta =
647 distanceAlongRay( 1, r, p );
648 delete r;
649 }
650 else { // starting point outside
651 Ray* r = new Ray( p, -d );
652 delta =
653 -distanceAlongRay( 1, r, p );
654 delete r;
655 }
656 }
657 else { // current point outside
658 if ( in0 ) { // starting point Inside
659 Ray* r = new Ray( p, -d );
660 delta =
661 -distanceAlongRay( -1, r, p );
662 delete r;
663 }
664 else { // starting point outside
665 Ray* r = new Ray( p, d );
666 delta =
667 distanceAlongRay( -1, r, p );
668 delete r;
669 }
670 }
671// Test if distance is less than the surface precision, if so Terminate loop.
672 if ( std::fabs( delta / sc ) <= SURFACE_PRECISION )
673 break;
674// If delta has not changed sufficiently from the previous iteration,
675// skip out of this loop.
676 if ( std::fabs( ( delta - delta0 ) / sc ) <=
677 SURFACE_PRECISION )
678 break;
679// If delta has increased in absolute value from the previous iteration
680// either the Helix doesn't Intersect the G4CylindricalSurface or the approximate
681// solution is too far from the real solution. Try groping for a solution.
682// If not found, Reset distance to large number, indicating no intersection with
683// the G4CylindricalSurface.
684 if ( std::fabs( delta ) > std::fabs( delta0 ) ) {
685 Dist = std::fabs( rhp ) *
686 gropeAlongHelix( hx );
687 if ( Dist < 0.0 ) {
688 Dist = FLT_MAXX;
689 p = lv;
690 }
691 else
692 p = hx->position( Dist );
693 break;
694 }
695// Set old delta to new one.
696 delta0 = delta;
697// Add distance to G4CylindricalSurface to distance along Helix.
698 Dist += delta;
699// Negative distance along Helix means Helix doesn't Intersect
700// G4CylindricalSurface.
701// Reset distance to large number, indicating no intersection with
702// G4CylindricalSurface.
703 if ( Dist < 0.0 ) {
704 Dist = FLT_MAXX;
705 p = lv;
706 break;
707 }
708// Recalculate point along Helix and the G4Vector3D.
709 p = hx->position( Dist );
710 d = hx->direction( Dist );
711 } // end of while loop
712 } // end of exact == 0 condition
713// Now have best value of distance along Helix and position for this
714// solution, so test if it is within the boundary of the sub-shape
715// and require that it point in the correct G4Vector3D with respect to
716// the Normal to the G4CylindricalSurface.
717 if ( ( Dist < FLT_MAXX ) &&
718 ( ( hx->direction( Dist ) * Normal( p ) *
719 which_way ) >= 0.0 ) &&
720 ( WithinBoundary( p ) == 1 ) )
721 return Dist;
722 } // end of if s[isoln] >= 0.0 condition
723 } // end of for loop over solutions
724// if one gets here, there is no solution, so set distance along Helix
725// and position to large numbers
726 Dist = FLT_MAXX;
727 p = lv;
728 return Dist;
729}
730*/
731
732
733G4Vector3D G4CylindricalSurface::Normal( const G4Vector3D& p ) const
734{
735 // return the Normal unit vector to the G4CylindricalSurface
736 // at a point p on (or nearly on) the G4CylindricalSurface
737
738 G4Vector3D n = ( p - origin ) - ( ( p - origin ) * axis ) * axis;
739 G4double nmag = n.mag();
740
741 if ( nmag != 0.0 )
742 n = n * (1/nmag);
743
744 return n;
745}
746
747
748G4Vector3D G4CylindricalSurface::SurfaceNormal( const G4Point3D& p ) const
749{
750 // return the Normal unit vector to the G4CylindricalSurface at a point
751 // p on (or nearly on) the G4CylindricalSurface
752
753 G4Vector3D n = ( p - origin ) - ( ( p - origin ) * axis ) * axis;
754 G4double nmag = n.mag();
755
756 if ( nmag != 0.0 )
757 n = n * (1/nmag);
758
759 return n;
760}
761
762
763G4int G4CylindricalSurface::Inside ( const G4Vector3D& x ) const
764{
765 // Return 0 if point x is outside G4CylindricalSurface, 1 if Inside.
766 // Outside means that the distance to the G4CylindricalSurface would
767 // be negative.
768 // Use the HowNear function to calculate this distance.
769 if ( HowNear( x ) >= -0.5*kCarTolerance )
770 return 1;
771 else
772 return 0;
773}
774
775
776G4int G4CylindricalSurface::WithinBoundary( const G4Vector3D& x ) const
777{
778 // return 1 if point x is on the G4CylindricalSurface, otherwise return zero
779 // base this on the surface precision factor set in support/globals.h
780 if ( std::fabs( HowNear( x ) / Scale() ) <= SURFACE_PRECISION )
781 return 1;
782 else
783 return 0;
784}
785
786
787G4double G4CylindricalSurface::Scale() const
788{
789 // Returns the radius of a G4CylindricalSurface unless it is zero, in which
790 // case returns the arbitrary number 1.0.
791 // This is ok since derived finite-sized classes will overwrite this.
792 // Used for Scale-invariant tests of surface thickness.
793 if ( radius == 0.0 )
794 return 1.0;
795 else
796 return radius;
797}
798
799
800//void G4CylindricalSurface::rotate( G4double alpha, G4double beta,
801// G4double gamma, G4ThreeMat& m, G4int inverse )
802// // rotate G4CylindricalSurface first about global x-axis by angle alpha,
803// second about global y-axis by angle beta,
804// and third about global z-axis by angle gamma
805// by creating and using G4ThreeMat objects in Surface::rotate
806// angles are assumed to be given in radians
807// if inverse is non-zero, the order of rotations is reversed
808// the axis is rotated here, the origin is rotated by calling
809// Surface::rotate
810// G4Surface::rotate( alpha, beta, gamma, m, inverse );
811// axis = m * axis;
812//}
813
814//void G4CylindricalSurface::rotate( G4double alpha, G4double beta,
815// G4double gamma, G4int inverse )
816//{ // rotate G4CylindricalSurface first about global x-axis by angle alpha,
817// second about global y-axis by angle beta,
818// and third about global z-axis by angle gamma
819// by creating and using G4ThreeMat objects in Surface::rotate
820// angles are assumed to be given in radians
821// if inverse is non-zero, the order of rotations is reversed
822// the axis is rotated here, the origin is rotated by calling
823// Surface::rotate
824// G4ThreeMat m;
825// G4Surface::rotate( alpha, beta, gamma, m, inverse );
826// axis = m * axis;
827//}
828
829
830void G4CylindricalSurface::SetRadius( G4double r )
831{
832 // Reset the radius of the G4CylindricalSurface
833 // Require radius to be non-negative
834 if ( r >= 0.0 )
835 radius = r;
836 // use old value (do not change radius) if out of the range,
837 // but Print message
838 else
839 {
840 G4cerr << "Error in G4CylindricalSurface::SetRadius"
841 << "--asked for negative radius\n"
842 << "\tDefault radius of " << radius << " is used.\n";
843 }
844}
845
846
847/*
848G4double G4CylindricalSurface::gropeAlongHelix( const Helix* hx ) const
849{ // Grope for a solution of a Helix intersecting a G4CylindricalSurface.
850 // This function returns the turning angle (in radians) where the
851 // intersection occurs with only positive values allowed, or -1.0 if
852 // no intersection is found.
853// The idea is to start at the beginning of the Helix, then take steps
854// of some fraction of a turn. If at the end of a Step, the current position
855// along the Helix and the previous position are on opposite sides of the
856// G4CylindricalSurface, then the solution must lie somewhere in between.
857 G4int one_over_f = 8; // one over fraction of a turn to go in each Step
858 G4double turn_angle = 0.0;
859 G4double dist_along = 0.0;
860 G4double d_new;
861 G4double fk = 1.0 / G4double( one_over_f );
862 G4double scal = Scale();
863 G4double d_old = HowNear( hx->position( dist_along ) );
864 G4double rh = hx->GetRadius(); // radius of Helix
865 G4Vector3D prp = hx->getPerp(); // perpendicular vector
866 G4double prpmag = prp.Magnitude();
867 G4double rhp = rh / prpmag;
868 G4int max_iter = one_over_f * HELIX_MAX_TURNS;
869// Take up to a user-settable number of turns along the Helix,
870// groping for an intersection point.
871 for ( G4int k = 1; k < max_iter; k++ ) {
872 turn_angle = 2.0 * pi * k / one_over_f;
873 dist_along = turn_angle * std::fabs( rhp );
874 d_new = HowNear( hx->position( dist_along ) );
875 if ( ( d_old < 0.0 && d_new > 0.0 ) ||
876 ( d_old > 0.0 && d_new < 0.0 ) ) {
877 d_old = d_new;
878// Old and new points are on opposite sides of the G4CylindricalSurface, therefore
879// a solution lies in between, use a binary search to pin the point down
880// to the surface precision, but don't do more than 50 iterations.
881 G4int itr = 0;
882 while ( std::fabs( d_new / scal ) > SURFACE_PRECISION ) {
883 itr++;
884 if ( itr > 50 )
885 return turn_angle;
886 turn_angle -= fk * pi;
887 dist_along = turn_angle * std::fabs( rhp );
888 d_new = HowNear( hx->position( dist_along ) );
889 if ( ( d_old < 0.0 && d_new > 0.0 ) ||
890 ( d_old > 0.0 && d_new < 0.0 ) )
891 fk *= -0.5;
892 else
893 fk *= 0.5;
894 d_old = d_new;
895 } // end of while loop
896 return turn_angle; // this is the best solution
897 } // end of if condition
898 } // end of for loop
899// Get here only if no solution is found, so return -1.0 to indicate that.
900 return -1.0;
901}
902*/
Note: See TracBrowser for help on using the repository browser.