source: trunk/source/geometry/solids/specific/src/G4Ellipsoid.cc@ 1358

Last change on this file since 1358 was 1337, checked in by garnier, 15 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 31.9 KB
Line 
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// $Id: G4Ellipsoid.cc,v 1.24 2009/09/24 15:51:02 gcosmo Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-01 $
28//
29// class G4Ellipsoid
30//
31// Implementation for G4Ellipsoid class
32//
33// History:
34//
35// 10.11.99 G.Horton-Smith -- first writing, based on G4Sphere class
36// 25.02.05 G.Guerrieri -- Modified for future Geant4 release
37//
38// --------------------------------------------------------------------
39
40#include "globals.hh"
41
42#include "G4Ellipsoid.hh"
43
44#include "G4VoxelLimits.hh"
45#include "G4AffineTransform.hh"
46#include "G4GeometryTolerance.hh"
47
48#include "meshdefs.hh"
49
50#include "Randomize.hh"
51
52#include "G4VGraphicsScene.hh"
53#include "G4Polyhedron.hh"
54#include "G4NURBS.hh"
55#include "G4NURBSbox.hh"
56#include "G4VisExtent.hh"
57
58using namespace CLHEP;
59
60///////////////////////////////////////////////////////////////////////////////
61//
62// constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
63// - note if pDPhi>2PI then reset to 2PI
64
65G4Ellipsoid::G4Ellipsoid(const G4String& pName,
66 G4double pxSemiAxis,
67 G4double pySemiAxis,
68 G4double pzSemiAxis,
69 G4double pzBottomCut,
70 G4double pzTopCut)
71 : G4VSolid(pName), fpPolyhedron(0), fCubicVolume(0.), fSurfaceArea(0.),
72 zBottomCut(0.), zTopCut(0.)
73{
74 // note: for users that want to use the full ellipsoid it is useful
75 // to include a default for the cuts
76
77 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
78
79 // Check Semi-Axis
80 if ( (pxSemiAxis>0.) && (pySemiAxis>0.) && (pzSemiAxis>0.) )
81 {
82 SetSemiAxis(pxSemiAxis, pySemiAxis, pzSemiAxis);
83 }
84 else
85 {
86 G4cerr << "ERROR - G4Ellipsoid::G4Ellipsoid(): " << GetName() << G4endl
87 << " Invalid semi-axis !"
88 << G4endl;
89 G4Exception("G4Ellipsoid::G4Ellipsoid()", "InvalidSetup",
90 FatalException, "Invalid semi-axis.");
91 }
92
93 if ( pzBottomCut == 0 && pzTopCut == 0 )
94 {
95 SetZCuts(-pzSemiAxis, pzSemiAxis);
96 }
97 else if ( (pzBottomCut < pzSemiAxis) && (pzTopCut > -pzSemiAxis)
98 && (pzBottomCut < pzTopCut) )
99 {
100 SetZCuts(pzBottomCut, pzTopCut);
101 }
102 else
103 {
104 G4cerr << "ERROR - G4Ellipsoid::G4Ellipsoid(): " << GetName() << G4endl
105 << " Invalid z-coordinate for cutting plane !"
106 << G4endl;
107 G4Exception("G4Ellipsoid::G4Ellipsoid()", "InvalidSetup",
108 FatalException, "Invalid z-coordinate for cutting plane.");
109 }
110}
111
112///////////////////////////////////////////////////////////////////////////////
113//
114// Fake default constructor - sets only member data and allocates memory
115// for usage restricted to object persistency.
116//
117G4Ellipsoid::G4Ellipsoid( __void__& a )
118 : G4VSolid(a), fpPolyhedron(0), fCubicVolume(0.), fSurfaceArea(0.)
119{
120}
121
122///////////////////////////////////////////////////////////////////////////////
123//
124// Destructor
125
126G4Ellipsoid::~G4Ellipsoid()
127{
128}
129
130///////////////////////////////////////////////////////////////////////////////
131//
132// Calculate extent under transform and specified limit
133
134G4bool
135G4Ellipsoid::CalculateExtent(const EAxis pAxis,
136 const G4VoxelLimits& pVoxelLimit,
137 const G4AffineTransform& pTransform,
138 G4double& pMin, G4double& pMax) const
139{
140 if (!pTransform.IsRotated())
141 {
142 // Special case handling for unrotated solid ellipsoid
143 // Compute x/y/z mins and maxs for bounding box respecting limits,
144 // with early returns if outside limits. Then switch() on pAxis,
145 // and compute exact x and y limit for x/y case
146
147 G4double xoffset,xMin,xMax;
148 G4double yoffset,yMin,yMax;
149 G4double zoffset,zMin,zMax;
150
151 G4double maxDiff,newMin,newMax;
152 G4double xoff,yoff;
153
154 xoffset=pTransform.NetTranslation().x();
155 xMin=xoffset - xSemiAxis;
156 xMax=xoffset + xSemiAxis;
157 if (pVoxelLimit.IsXLimited())
158 {
159 if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
160 || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
161 {
162 return false;
163 }
164 else
165 {
166 if (xMin<pVoxelLimit.GetMinXExtent())
167 {
168 xMin=pVoxelLimit.GetMinXExtent();
169 }
170 if (xMax>pVoxelLimit.GetMaxXExtent())
171 {
172 xMax=pVoxelLimit.GetMaxXExtent();
173 }
174 }
175 }
176
177 yoffset=pTransform.NetTranslation().y();
178 yMin=yoffset - ySemiAxis;
179 yMax=yoffset + ySemiAxis;
180 if (pVoxelLimit.IsYLimited())
181 {
182 if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
183 || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
184 {
185 return false;
186 }
187 else
188 {
189 if (yMin<pVoxelLimit.GetMinYExtent())
190 {
191 yMin=pVoxelLimit.GetMinYExtent();
192 }
193 if (yMax>pVoxelLimit.GetMaxYExtent())
194 {
195 yMax=pVoxelLimit.GetMaxYExtent();
196 }
197 }
198 }
199
200 zoffset=pTransform.NetTranslation().z();
201 zMin=zoffset + (-zSemiAxis > zBottomCut ? -zSemiAxis : zBottomCut);
202 zMax=zoffset + ( zSemiAxis < zTopCut ? zSemiAxis : zTopCut);
203 if (pVoxelLimit.IsZLimited())
204 {
205 if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
206 || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
207 {
208 return false;
209 }
210 else
211 {
212 if (zMin<pVoxelLimit.GetMinZExtent())
213 {
214 zMin=pVoxelLimit.GetMinZExtent();
215 }
216 if (zMax>pVoxelLimit.GetMaxZExtent())
217 {
218 zMax=pVoxelLimit.GetMaxZExtent();
219 }
220 }
221 }
222
223 // if here, then known to cut bounding box around ellipsoid
224 //
225 xoff = (xoffset < xMin) ? (xMin-xoffset)
226 : (xoffset > xMax) ? (xoffset-xMax) : 0.0;
227 yoff = (yoffset < yMin) ? (yMin-yoffset)
228 : (yoffset > yMax) ? (yoffset-yMax) : 0.0;
229
230 // detailed calculations
231 // NOTE: does not use X or Y offsets to adjust Z range,
232 // and does not use Z offset to adjust X or Y range,
233 // which is consistent with G4Sphere::CalculateExtent behavior
234 //
235 switch (pAxis)
236 {
237 case kXAxis:
238 if (yoff==0.)
239 {
240 // YZ limits cross max/min x => no change
241 //
242 pMin=xMin;
243 pMax=xMax;
244 }
245 else
246 {
247 // YZ limits don't cross max/min x => compute max delta x,
248 // hence new mins/maxs
249 //
250 maxDiff= 1.0-sqr(yoff/ySemiAxis);
251 if (maxDiff < 0.0) { return false; }
252 maxDiff= xSemiAxis * std::sqrt(maxDiff);
253 newMin=xoffset-maxDiff;
254 newMax=xoffset+maxDiff;
255 pMin=(newMin<xMin) ? xMin : newMin;
256 pMax=(newMax>xMax) ? xMax : newMax;
257 }
258 break;
259 case kYAxis:
260 if (xoff==0.)
261 {
262 // XZ limits cross max/min y => no change
263 //
264 pMin=yMin;
265 pMax=yMax;
266 }
267 else
268 {
269 // XZ limits don't cross max/min y => compute max delta y,
270 // hence new mins/maxs
271 //
272 maxDiff= 1.0-sqr(xoff/xSemiAxis);
273 if (maxDiff < 0.0) { return false; }
274 maxDiff= ySemiAxis * std::sqrt(maxDiff);
275 newMin=yoffset-maxDiff;
276 newMax=yoffset+maxDiff;
277 pMin=(newMin<yMin) ? yMin : newMin;
278 pMax=(newMax>yMax) ? yMax : newMax;
279 }
280 break;
281 case kZAxis:
282 pMin=zMin;
283 pMax=zMax;
284 break;
285 default:
286 break;
287 }
288
289 pMin-=kCarTolerance;
290 pMax+=kCarTolerance;
291 return true;
292 }
293 else // not rotated
294 {
295 G4int i,j,noEntries,noBetweenSections;
296 G4bool existsAfterClip=false;
297
298 // Calculate rotated vertex coordinates
299
300 G4int noPolygonVertices=0;
301 G4ThreeVectorList* vertices =
302 CreateRotatedVertices(pTransform,noPolygonVertices);
303
304 pMin=+kInfinity;
305 pMax=-kInfinity;
306
307 noEntries=vertices->size(); // noPolygonVertices*noPhiCrossSections
308 noBetweenSections=noEntries-noPolygonVertices;
309
310 G4ThreeVectorList ThetaPolygon;
311 for (i=0;i<noEntries;i+=noPolygonVertices)
312 {
313 for(j=0;j<(noPolygonVertices/2)-1;j++)
314 {
315 ThetaPolygon.push_back((*vertices)[i+j]);
316 ThetaPolygon.push_back((*vertices)[i+j+1]);
317 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]);
318 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]);
319 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
320 ThetaPolygon.clear();
321 }
322 }
323 for (i=0;i<noBetweenSections;i+=noPolygonVertices)
324 {
325 for(j=0;j<noPolygonVertices-1;j++)
326 {
327 ThetaPolygon.push_back((*vertices)[i+j]);
328 ThetaPolygon.push_back((*vertices)[i+j+1]);
329 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]);
330 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]);
331 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
332 ThetaPolygon.clear();
333 }
334 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]);
335 ThetaPolygon.push_back((*vertices)[i]);
336 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]);
337 ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]);
338 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
339 ThetaPolygon.clear();
340 }
341 if ( (pMin!=kInfinity) || (pMax!=-kInfinity) )
342 {
343 existsAfterClip=true;
344
345 // Add 2*tolerance to avoid precision troubles
346 //
347 pMin-=kCarTolerance;
348 pMax+=kCarTolerance;
349
350 }
351 else
352 {
353 // Check for case where completely enveloping clipping volume
354 // If point inside then we are confident that the solid completely
355 // envelopes the clipping volume. Hence set min/max extents according
356 // to clipping volume extents along the specified axis.
357 //
358 G4ThreeVector
359 clipCentre((pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
360 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
361 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
362
363 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
364 {
365 existsAfterClip=true;
366 pMin=pVoxelLimit.GetMinExtent(pAxis);
367 pMax=pVoxelLimit.GetMaxExtent(pAxis);
368 }
369 }
370 delete vertices;
371 return existsAfterClip;
372 }
373}
374
375///////////////////////////////////////////////////////////////////////////////
376//
377// Return whether point inside/outside/on surface
378// Split into radius, phi, theta checks
379// Each check modifies `in', or returns as approprate
380
381EInside G4Ellipsoid::Inside(const G4ThreeVector& p) const
382{
383 G4double rad2oo, // outside surface outer tolerance
384 rad2oi; // outside surface inner tolerance
385 EInside in;
386
387 static const G4double halfRadTolerance=kRadTolerance*0.5;
388
389 // check this side of z cut first, because that's fast
390 //
391 if (p.z() < zBottomCut-halfRadTolerance) { return in=kOutside; }
392 if (p.z() > zTopCut+halfRadTolerance) { return in=kOutside; }
393
394 rad2oo= sqr(p.x()/(xSemiAxis+halfRadTolerance))
395 + sqr(p.y()/(ySemiAxis+halfRadTolerance))
396 + sqr(p.z()/(zSemiAxis+halfRadTolerance));
397
398 if (rad2oo > 1.0) { return in=kOutside; }
399
400 rad2oi= sqr(p.x()*(1.0+halfRadTolerance/xSemiAxis)/xSemiAxis)
401 + sqr(p.y()*(1.0+halfRadTolerance/ySemiAxis)/ySemiAxis)
402 + sqr(p.z()*(1.0+halfRadTolerance/zSemiAxis)/zSemiAxis);
403
404 // Check radial surfaces
405 // sets `in' (already checked for rad2oo > 1.0)
406 //
407 if (rad2oi < 1.0)
408 {
409 in = ( (p.z() < zBottomCut+halfRadTolerance)
410 || (p.z() > zTopCut-halfRadTolerance) ) ? kSurface : kInside;
411 if ( rad2oi > 1.0-halfRadTolerance ) { in=kSurface; }
412 }
413 else
414 {
415 in = kSurface;
416 }
417 return in;
418
419}
420
421///////////////////////////////////////////////////////////////////////////////
422//
423// Return unit normal of surface closest to p not protected against p=0
424
425G4ThreeVector G4Ellipsoid::SurfaceNormal( const G4ThreeVector& p) const
426{
427 G4double distR, distZBottom, distZTop;
428
429 // normal vector with special magnitude: parallel to normal, units 1/length
430 // norm*p == 1.0 if on surface, >1.0 if outside, <1.0 if inside
431 //
432 G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
433 p.y()/(ySemiAxis*ySemiAxis),
434 p.z()/(zSemiAxis*zSemiAxis));
435 G4double radius = 1.0/norm.mag();
436
437 // approximate distance to curved surface
438 //
439 distR = std::fabs( (p*norm - 1.0) * radius ) / 2.0;
440
441 // Distance to z-cut plane
442 //
443 distZBottom = std::fabs( p.z() - zBottomCut );
444 distZTop = std::fabs( p.z() - zTopCut );
445
446 if ( (distZBottom < distR) || (distZTop < distR) )
447 {
448 return G4ThreeVector(0.,0.,(distZBottom < distZTop) ? -1.0 : 1.0);
449 }
450 return ( norm *= radius );
451}
452
453///////////////////////////////////////////////////////////////////////////////
454//
455// Calculate distance to shape from outside, along normalised vector
456// - return kInfinity if no intersection, or intersection distance <= tolerance
457//
458
459G4double G4Ellipsoid::DistanceToIn( const G4ThreeVector& p,
460 const G4ThreeVector& v ) const
461{
462 static const G4double halfCarTolerance=kCarTolerance*0.5;
463 static const G4double halfRadTolerance=kRadTolerance*0.5;
464
465 G4double distMin = std::min(xSemiAxis,ySemiAxis);
466 const G4double dRmax = 100.*std::min(distMin,zSemiAxis);
467 distMin= kInfinity;
468
469 // check to see if Z plane is relevant
470 if (p.z() <= zBottomCut+halfCarTolerance)
471 {
472 if (v.z() <= 0.0) { return distMin; }
473 G4double distZ = (zBottomCut - p.z()) / v.z();
474
475 if ( (distZ > -halfRadTolerance) && (Inside(p+distZ*v) != kOutside) )
476 {
477 // early exit since can't intercept curved surface if we reach here
478 if ( std::abs(distZ) < halfRadTolerance ) { distZ=0.; }
479 return distMin= distZ;
480 }
481 }
482 if (p.z() >= zTopCut-halfCarTolerance)
483 {
484 if (v.z() >= 0.0) { return distMin;}
485 G4double distZ = (zTopCut - p.z()) / v.z();
486 if ( (distZ > -halfRadTolerance) && (Inside(p+distZ*v) != kOutside) )
487 {
488 // early exit since can't intercept curved surface if we reach here
489 if ( std::abs(distZ) < halfRadTolerance ) { distZ=0.; }
490 return distMin= distZ;
491 }
492 }
493 // if fZCut1 <= p.z() <= fZCut2, then must hit curved surface
494
495 // now check curved surface intercept
496 G4double A,B,C;
497
498 A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis);
499 C= sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) + sqr(p.z()/zSemiAxis) - 1.0;
500 B= 2.0 * ( p.x()*v.x()/(xSemiAxis*xSemiAxis)
501 + p.y()*v.y()/(ySemiAxis*ySemiAxis)
502 + p.z()*v.z()/(zSemiAxis*zSemiAxis) );
503
504 C= B*B - 4.0*A*C;
505 if (C > 0.0)
506 {
507 G4double distR= (-B - std::sqrt(C)) / (2.0*A);
508 G4double intZ = p.z()+distR*v.z();
509 if ( (distR > halfRadTolerance)
510 && (intZ >= zBottomCut-halfRadTolerance)
511 && (intZ <= zTopCut+halfRadTolerance) )
512 {
513 distMin = distR;
514 }
515 else if( (distR >- halfRadTolerance)
516 && (intZ >= zBottomCut-halfRadTolerance)
517 && (intZ <= zTopCut+halfRadTolerance) )
518 {
519 // p is on the curved surface, DistanceToIn returns 0 or kInfinity:
520 // DistanceToIn returns 0, if second root is positive (means going inside)
521 // If second root is negative, DistanceToIn returns kInfinity (outside)
522 //
523 distR = (-B + std::sqrt(C) ) / (2.0*A);
524 if(distR>0.) { distMin=0.; }
525 }
526 else
527 {
528 distR= (-B + std::sqrt(C)) / (2.0*A);
529 intZ = p.z()+distR*v.z();
530 if ( (distR > halfRadTolerance)
531 && (intZ >= zBottomCut-halfRadTolerance)
532 && (intZ <= zTopCut+halfRadTolerance) )
533 {
534 G4ThreeVector norm=SurfaceNormal(p);
535 if (norm.dot(v)<0.) { distMin = distR; }
536 }
537 }
538 if ( (distMin!=kInfinity) && (distMin>dRmax) )
539 { // Avoid rounding errors due to precision issues on
540 // 64 bits systems. Split long distances and recompute
541 G4double fTerm = distMin-std::fmod(distMin,dRmax);
542 distMin = fTerm + DistanceToIn(p+fTerm*v,v);
543 }
544 }
545
546 if (std::abs(distMin)<halfRadTolerance) { distMin=0.; }
547 return distMin;
548}
549
550///////////////////////////////////////////////////////////////////////////////
551//
552// Calculate distance (<= actual) to closest surface of shape from outside
553// - Return 0 if point inside
554
555G4double G4Ellipsoid::DistanceToIn(const G4ThreeVector& p) const
556{
557 G4double distR, distZ;
558
559 // normal vector: parallel to normal, magnitude 1/(characteristic radius)
560 //
561 G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
562 p.y()/(ySemiAxis*ySemiAxis),
563 p.z()/(zSemiAxis*zSemiAxis));
564 G4double radius= 1.0/norm.mag();
565
566 // approximate distance to curved surface ( <= actual distance )
567 //
568 distR= (p*norm - 1.0) * radius / 2.0;
569
570 // Distance to z-cut plane
571 //
572 distZ= zBottomCut - p.z();
573 if (distZ < 0.0)
574 {
575 distZ = p.z() - zTopCut;
576 }
577
578 // Distance to closest surface from outside
579 //
580 if (distZ < 0.0)
581 {
582 return (distR < 0.0) ? 0.0 : distR;
583 }
584 else if (distR < 0.0)
585 {
586 return distZ;
587 }
588 else
589 {
590 return (distZ < distR) ? distZ : distR;
591 }
592}
593
594///////////////////////////////////////////////////////////////////////////////
595//
596// Calculate distance to surface of shape from `inside', allowing for tolerance
597
598G4double G4Ellipsoid::DistanceToOut(const G4ThreeVector& p,
599 const G4ThreeVector& v,
600 const G4bool calcNorm,
601 G4bool *validNorm,
602 G4ThreeVector *n ) const
603{
604 G4double distMin;
605 enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
606
607 distMin= kInfinity;
608 surface= kNoSurf;
609
610 // check to see if Z plane is relevant
611 //
612 if (v.z() < 0.0)
613 {
614 G4double distZ = (zBottomCut - p.z()) / v.z();
615 if (distZ < 0.0)
616 {
617 distZ= 0.0;
618 if (!calcNorm) {return 0.0;}
619 }
620 distMin= distZ;
621 surface= kPlaneSurf;
622 }
623 if (v.z() > 0.0)
624 {
625 G4double distZ = (zTopCut - p.z()) / v.z();
626 if (distZ < 0.0)
627 {
628 distZ= 0.0;
629 if (!calcNorm) {return 0.0;}
630 }
631 distMin= distZ;
632 surface= kPlaneSurf;
633 }
634
635 // normal vector: parallel to normal, magnitude 1/(characteristic radius)
636 //
637 G4ThreeVector nearnorm(p.x()/(xSemiAxis*xSemiAxis),
638 p.y()/(ySemiAxis*ySemiAxis),
639 p.z()/(zSemiAxis*zSemiAxis));
640
641 // now check curved surface intercept
642 //
643 G4double A,B,C;
644
645 A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis);
646 C= (p * nearnorm) - 1.0;
647 B= 2.0 * (v * nearnorm);
648
649 C= B*B - 4.0*A*C;
650 if (C > 0.0)
651 {
652 G4double distR= (-B + std::sqrt(C) ) / (2.0*A);
653 if (distR < 0.0)
654 {
655 distR= 0.0;
656 if (!calcNorm) {return 0.0;}
657 }
658 if (distR < distMin)
659 {
660 distMin= distR;
661 surface= kCurvedSurf;
662 }
663 }
664
665 // set normal if requested
666 //
667 if (calcNorm)
668 {
669 if (surface == kNoSurf)
670 {
671 *validNorm = false;
672 }
673 else
674 {
675 *validNorm = true;
676 switch (surface)
677 {
678 case kPlaneSurf:
679 *n= G4ThreeVector(0.,0.,(v.z() > 0.0 ? 1. : -1.));
680 break;
681 case kCurvedSurf:
682 {
683 G4ThreeVector pexit= p + distMin*v;
684 G4ThreeVector truenorm(pexit.x()/(xSemiAxis*xSemiAxis),
685 pexit.y()/(ySemiAxis*ySemiAxis),
686 pexit.z()/(zSemiAxis*zSemiAxis));
687 truenorm *= 1.0/truenorm.mag();
688 *n= truenorm;
689 } break;
690 default:
691 G4cout.precision(16);
692 G4cout << G4endl;
693 DumpInfo();
694 G4cout << "Position:" << G4endl << G4endl;
695 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
696 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
697 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
698 G4cout << "Direction:" << G4endl << G4endl;
699 G4cout << "v.x() = " << v.x() << G4endl;
700 G4cout << "v.y() = " << v.y() << G4endl;
701 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
702 G4cout << "Proposed distance :" << G4endl << G4endl;
703 G4cout << "distMin = " << distMin/mm << " mm" << G4endl << G4endl;
704 G4Exception("G4Ellipsoid::DistanceToOut(p,v,..)",
705 "Notification", JustWarning,
706 "Undefined side for valid surface normal to solid.");
707 break;
708 }
709 }
710 }
711
712 return distMin;
713}
714
715///////////////////////////////////////////////////////////////////////////////
716//
717// Calculate distance (<=actual) to closest surface of shape from inside
718
719G4double G4Ellipsoid::DistanceToOut(const G4ThreeVector& p) const
720{
721 G4double distR, distZ;
722
723#ifdef G4SPECSDEBUG
724 if( Inside(p) == kOutside )
725 {
726 G4cout.precision(16) ;
727 G4cout << G4endl ;
728 DumpInfo();
729 G4cout << "Position:" << G4endl << G4endl ;
730 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
731 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
732 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
733 G4Exception("G4Ellipsoid::DistanceToOut(p)", "Notification", JustWarning,
734 "Point p is outside !?" );
735 }
736#endif
737
738 // Normal vector: parallel to normal, magnitude 1/(characteristic radius)
739 //
740 G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
741 p.y()/(ySemiAxis*ySemiAxis),
742 p.z()/(zSemiAxis*zSemiAxis));
743
744 // the following is a safe inlined "radius= min(1.0/norm.mag(),p.mag())
745 //
746 G4double radius= p.mag();
747 G4double tmp= norm.mag();
748 if ( (tmp > 0.0) && (1.0 < radius*tmp) ) {radius = 1.0/tmp;}
749
750 // Approximate distance to curved surface ( <= actual distance )
751 //
752 distR = (1.0 - p*norm) * radius / 2.0;
753
754 // Distance to z-cut plane
755 //
756 distZ = p.z() - zBottomCut;
757 if (distZ < 0.0) {distZ= zTopCut - p.z();}
758
759 // Distance to closest surface from inside
760 //
761 if ( (distZ < 0.0) || (distR < 0.0) )
762 {
763 return 0.0;
764 }
765 else
766 {
767 return (distZ < distR) ? distZ : distR;
768 }
769}
770
771///////////////////////////////////////////////////////////////////////////////
772//
773// Create a List containing the transformed vertices
774// Ordering [0-3] -fDz cross section
775// [4-7] +fDz cross section such that [0] is below [4],
776// [1] below [5] etc.
777// Note:
778// Caller has deletion resposibility
779// Potential improvement: For last slice, use actual ending angle
780// to avoid rounding error problems.
781
782G4ThreeVectorList*
783G4Ellipsoid::CreateRotatedVertices(const G4AffineTransform& pTransform,
784 G4int& noPolygonVertices) const
785{
786 G4ThreeVectorList *vertices;
787 G4ThreeVector vertex;
788 G4double meshAnglePhi, meshRMaxFactor,
789 crossAnglePhi, coscrossAnglePhi, sincrossAnglePhi, sAnglePhi;
790 G4double meshTheta, crossTheta, startTheta;
791 G4double rMaxX, rMaxY, rMaxZ, rMaxMax, rx, ry, rz;
792 G4int crossSectionPhi, noPhiCrossSections, crossSectionTheta, noThetaSections;
793
794 // Phi cross sections
795 //
796 noPhiCrossSections=G4int (twopi/kMeshAngleDefault)+1;
797
798 if (noPhiCrossSections<kMinMeshSections)
799 {
800 noPhiCrossSections=kMinMeshSections;
801 }
802 else if (noPhiCrossSections>kMaxMeshSections)
803 {
804 noPhiCrossSections=kMaxMeshSections;
805 }
806 meshAnglePhi=twopi/(noPhiCrossSections-1);
807
808 // Set start angle such that mesh will be at fRMax
809 // on the x axis. Will give better extent calculations when not rotated.
810
811 sAnglePhi = -meshAnglePhi*0.5;
812
813 // Theta cross sections
814
815 noThetaSections = G4int(pi/kMeshAngleDefault)+3;
816
817 if (noThetaSections<kMinMeshSections)
818 {
819 noThetaSections=kMinMeshSections;
820 }
821 else if (noThetaSections>kMaxMeshSections)
822 {
823 noThetaSections=kMaxMeshSections;
824 }
825 meshTheta= pi/(noThetaSections-2);
826
827 // Set start angle such that mesh will be at fRMax
828 // on the z axis. Will give better extent calculations when not rotated.
829
830 startTheta = -meshTheta*0.5;
831
832 meshRMaxFactor = 1.0/std::cos(0.5*
833 std::sqrt(meshAnglePhi*meshAnglePhi+meshTheta*meshTheta));
834 rMaxMax= (xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis);
835 if (zSemiAxis > rMaxMax) rMaxMax= zSemiAxis;
836 rMaxX= xSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
837 rMaxY= ySemiAxis + rMaxMax*(meshRMaxFactor-1.0);
838 rMaxZ= zSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
839 G4double* cosCrossTheta = new G4double[noThetaSections];
840 G4double* sinCrossTheta = new G4double[noThetaSections];
841 vertices=new G4ThreeVectorList(noPhiCrossSections*noThetaSections);
842 if (vertices && cosCrossTheta && sinCrossTheta)
843 {
844 for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
845 crossSectionTheta++)
846 {
847 // Compute sine and cosine table (for historical reasons)
848 //
849 crossTheta=startTheta+crossSectionTheta*meshTheta;
850 cosCrossTheta[crossSectionTheta]=std::cos(crossTheta);
851 sinCrossTheta[crossSectionTheta]=std::sin(crossTheta);
852 }
853 for (crossSectionPhi=0; crossSectionPhi<noPhiCrossSections;
854 crossSectionPhi++)
855 {
856 crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
857 coscrossAnglePhi=std::cos(crossAnglePhi);
858 sincrossAnglePhi=std::sin(crossAnglePhi);
859 for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
860 crossSectionTheta++)
861 {
862 // Compute coordinates of cross section at section crossSectionPhi
863 //
864 rx= sinCrossTheta[crossSectionTheta]*coscrossAnglePhi*rMaxX;
865 ry= sinCrossTheta[crossSectionTheta]*sincrossAnglePhi*rMaxY;
866 rz= cosCrossTheta[crossSectionTheta]*rMaxZ;
867 if (rz < zBottomCut)
868 { rz= zBottomCut; }
869 if (rz > zTopCut)
870 { rz= zTopCut; }
871 vertex= G4ThreeVector(rx,ry,rz);
872 vertices->push_back(pTransform.TransformPoint(vertex));
873 } // Theta forward
874 } // Phi
875 noPolygonVertices = noThetaSections ;
876 }
877 else
878 {
879 DumpInfo();
880 G4Exception("G4Ellipsoid::CreateRotatedVertices()",
881 "FatalError", FatalException,
882 "Error in allocation of vertices. Out of memory !");
883 }
884
885 delete[] cosCrossTheta;
886 delete[] sinCrossTheta;
887
888 return vertices;
889}
890
891//////////////////////////////////////////////////////////////////////////
892//
893// G4EntityType
894
895G4GeometryType G4Ellipsoid::GetEntityType() const
896{
897 return G4String("G4Ellipsoid");
898}
899
900//////////////////////////////////////////////////////////////////////////
901//
902// Stream object contents to an output stream
903
904std::ostream& G4Ellipsoid::StreamInfo( std::ostream& os ) const
905{
906 os << "-----------------------------------------------------------\n"
907 << " *** Dump for solid - " << GetName() << " ***\n"
908 << " ===================================================\n"
909 << " Solid type: G4Ellipsoid\n"
910 << " Parameters: \n"
911
912 << " semi-axis x: " << xSemiAxis/mm << " mm \n"
913 << " semi-axis y: " << ySemiAxis/mm << " mm \n"
914 << " semi-axis z: " << zSemiAxis/mm << " mm \n"
915 << " max semi-axis: " << semiAxisMax/mm << " mm \n"
916 << " lower cut plane level z: " << zBottomCut/mm << " mm \n"
917 << " upper cut plane level z: " << zTopCut/mm << " mm \n"
918 << "-----------------------------------------------------------\n";
919
920 return os;
921}
922
923////////////////////////////////////////////////////////////////////
924//
925// GetPointOnSurface
926
927G4ThreeVector G4Ellipsoid::GetPointOnSurface() const
928{
929 G4double aTop, aBottom, aCurved, chose, xRand, yRand, zRand, phi, theta;
930 G4double cosphi, sinphi, costheta, sintheta, alpha, beta, max1, max2, max3;
931
932 max1 = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis;
933 max1 = max1 > zSemiAxis ? max1 : zSemiAxis;
934 if (max1 == xSemiAxis) { max2 = ySemiAxis; max3 = zSemiAxis; }
935 else if (max1 == ySemiAxis) { max2 = xSemiAxis; max3 = zSemiAxis; }
936 else { max2 = xSemiAxis; max3 = ySemiAxis; }
937
938 phi = RandFlat::shoot(0.,twopi);
939 theta = RandFlat::shoot(0.,pi);
940
941 cosphi = std::cos(phi); sinphi = std::sin(phi);
942 costheta = RandFlat::shoot(zBottomCut,zTopCut)/zSemiAxis;
943 sintheta = std::sqrt(1.-sqr(costheta));
944
945 alpha = 1.-sqr(max2/max1); beta = 1.-sqr(max3/max1);
946
947 aTop = pi*xSemiAxis*ySemiAxis*(1 - sqr(zTopCut/zSemiAxis));
948 aBottom = pi*xSemiAxis*ySemiAxis*(1 - sqr(zBottomCut/zSemiAxis));
949
950 // approximation
951 // from:" http://www.citr.auckland.ac.nz/techreports/2004/CITR-TR-139.pdf"
952 aCurved = 4.*pi*max1*max2*(1.-1./6.*(alpha+beta)-
953 1./120.*(3.*sqr(alpha)+2.*alpha*beta+3.*sqr(beta)));
954
955 aCurved *= 0.5*(1.2*zTopCut/zSemiAxis - 1.2*zBottomCut/zSemiAxis);
956
957 if( ( zTopCut >= zSemiAxis && zBottomCut <= -1.*zSemiAxis )
958 || ( zTopCut == 0 && zBottomCut ==0 ) )
959 {
960 aTop = 0; aBottom = 0;
961 }
962
963 chose = RandFlat::shoot(0.,aTop + aBottom + aCurved);
964
965 if(chose < aCurved)
966 {
967 xRand = xSemiAxis*sintheta*cosphi;
968 yRand = ySemiAxis*sintheta*sinphi;
969 zRand = zSemiAxis*costheta;
970 return G4ThreeVector (xRand,yRand,zRand);
971 }
972 else if(chose >= aCurved && chose < aCurved + aTop)
973 {
974 xRand = RandFlat::shoot(-1.,1.)*xSemiAxis
975 * std::sqrt(1-sqr(zTopCut/zSemiAxis));
976 yRand = RandFlat::shoot(-1.,1.)*ySemiAxis
977 * std::sqrt(1.-sqr(zTopCut/zSemiAxis)-sqr(xRand/xSemiAxis));
978 zRand = zTopCut;
979 return G4ThreeVector (xRand,yRand,zRand);
980 }
981 else
982 {
983 xRand = RandFlat::shoot(-1.,1.)*xSemiAxis
984 * std::sqrt(1-sqr(zBottomCut/zSemiAxis));
985 yRand = RandFlat::shoot(-1.,1.)*ySemiAxis
986 * std::sqrt(1.-sqr(zBottomCut/zSemiAxis)-sqr(xRand/xSemiAxis));
987 zRand = zBottomCut;
988 return G4ThreeVector (xRand,yRand,zRand);
989 }
990}
991
992/////////////////////////////////////////////////////////////////////////////
993//
994// Methods for visualisation
995
996void G4Ellipsoid::DescribeYourselfTo (G4VGraphicsScene& scene) const
997{
998 scene.AddSolid(*this);
999}
1000
1001G4VisExtent G4Ellipsoid::GetExtent() const
1002{
1003 // Define the sides of the box into which the G4Ellipsoid instance would fit.
1004 //
1005 return G4VisExtent (-semiAxisMax, semiAxisMax,
1006 -semiAxisMax, semiAxisMax,
1007 -semiAxisMax, semiAxisMax);
1008}
1009
1010G4NURBS* G4Ellipsoid::CreateNURBS () const
1011{
1012 // Box for now!!!
1013 //
1014 return new G4NURBSbox(semiAxisMax, semiAxisMax, semiAxisMax);
1015}
1016
1017G4Polyhedron* G4Ellipsoid::CreatePolyhedron () const
1018{
1019 return new G4PolyhedronEllipsoid(xSemiAxis, ySemiAxis, zSemiAxis,
1020 zBottomCut, zTopCut);
1021}
1022
1023G4Polyhedron* G4Ellipsoid::GetPolyhedron () const
1024{
1025 if (!fpPolyhedron ||
1026 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1027 fpPolyhedron->GetNumberOfRotationSteps())
1028 {
1029 delete fpPolyhedron;
1030 fpPolyhedron = CreatePolyhedron();
1031 }
1032 return fpPolyhedron;
1033}
Note: See TracBrowser for help on using the repository browser.