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

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

file release beta
File size: 30.5 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.14 2007/05/18 07:39:56 gcosmo Exp $
27// GEANT4 tag $Name: geant4-09-02-ref-02 $
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 // check this side of z cut first, because that's fast
388 //
389 if (p.z() < zBottomCut-kRadTolerance/2.0)
390 { return in=kOutside; }
391 if (p.z() > zTopCut+kRadTolerance/2.0)
392 { return in=kOutside; }
393
394 rad2oo= sqr(p.x()/(xSemiAxis+kRadTolerance/2.))
395 + sqr(p.y()/(ySemiAxis+kRadTolerance/2.))
396 + sqr(p.z()/(zSemiAxis+kRadTolerance/2.));
397
398 if (rad2oo > 1.0)
399 { return in=kOutside; }
400
401 rad2oi= sqr(p.x()*(1.0+kRadTolerance/2./xSemiAxis)/xSemiAxis)
402 + sqr(p.y()*(1.0+kRadTolerance/2./ySemiAxis)/ySemiAxis)
403 + sqr(p.z()*(1.0+kRadTolerance/2./zSemiAxis)/zSemiAxis);
404
405 // Check radial surfaces
406 // sets `in' (already checked for rad2oo > 1.0)
407 //
408 if (rad2oi < 1.0)
409 {
410 in = ( (p.z() < zBottomCut+kRadTolerance/2.0)
411 || (p.z() > zTopCut-kRadTolerance/2.0) ) ? kSurface : kInside;
412 }
413 else
414 {
415 in = kSurface;
416 }
417
418 return in;
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 G4double distMin;
463
464 distMin= kInfinity;
465
466 // check to see if Z plane is relevant
467 if (p.z() < zBottomCut) {
468 if (v.z() <= 0.0)
469 return distMin;
470 G4double distZ = (zBottomCut - p.z()) / v.z();
471 if (distZ > kRadTolerance/2.0 && Inside(p+distZ*v) != kOutside )
472 {
473 // early exit since can't intercept curved surface if we reach here
474 return distMin= distZ;
475 }
476 }
477 if (p.z() > zTopCut) {
478 if (v.z() >= 0.0)
479 return distMin;
480 G4double distZ = (zTopCut - p.z()) / v.z();
481 if (distZ > kRadTolerance/2.0 && Inside(p+distZ*v) != kOutside )
482 {
483 // early exit since can't intercept curved surface if we reach here
484 return distMin= distZ;
485 }
486 }
487 // if fZCut1 <= p.z() <= fZCut2, then must hit curved surface
488
489 // now check curved surface intercept
490 G4double A,B,C;
491
492 A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis);
493 C= sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) + sqr(p.z()/zSemiAxis) - 1.0;
494 B= 2.0 * ( p.x()*v.x()/(xSemiAxis*xSemiAxis) + p.y()*v.y()/(ySemiAxis*ySemiAxis)
495 + p.z()*v.z()/(zSemiAxis*zSemiAxis) );
496
497 C= B*B - 4.0*A*C;
498 if (C > 0.0)
499 {
500 G4double distR= (-B - std::sqrt(C) ) / (2.0*A);
501 G4double intZ= p.z()+distR*v.z();
502 if (distR > kRadTolerance/2.0
503 && intZ >= zBottomCut-kRadTolerance/2.0
504 && intZ <= zTopCut+kRadTolerance/2.0)
505 {
506 distMin = distR;
507 }
508 else
509 {
510 distR= (-B + std::sqrt(C) ) / (2.0*A);
511 intZ= p.z()+distR*v.z();
512 if (distR > kRadTolerance/2.0
513 && intZ >= zBottomCut-kRadTolerance/2.0
514 && intZ <= zTopCut+kRadTolerance/2.0)
515 {
516 distMin = distR;
517 }
518 }
519 }
520
521 return distMin;
522}
523
524///////////////////////////////////////////////////////////////////////////////
525//
526// Calculate distance (<= actual) to closest surface of shape from outside
527// - Return 0 if point inside
528
529G4double G4Ellipsoid::DistanceToIn(const G4ThreeVector& p) const
530{
531 G4double distR, distZ;
532
533 // normal vector: parallel to normal, magnitude 1/(characteristic radius)
534 //
535 G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
536 p.y()/(ySemiAxis*ySemiAxis),
537 p.z()/(zSemiAxis*zSemiAxis));
538 G4double radius= 1.0/norm.mag();
539
540 // approximate distance to curved surface ( <= actual distance )
541 //
542 distR= (p*norm - 1.0) * radius / 2.0;
543
544 // Distance to z-cut plane
545 //
546 distZ= zBottomCut - p.z();
547 if (distZ < 0.0)
548 {
549 distZ = p.z() - zTopCut;
550 }
551
552 // Distance to closest surface from outside
553 //
554 if (distZ < 0.0)
555 {
556 return (distR < 0.0) ? 0.0 : distR;
557 }
558 else if (distR < 0.0)
559 {
560 return distZ;
561 }
562 else
563 {
564 return (distZ < distR) ? distZ : distR;
565 }
566}
567
568///////////////////////////////////////////////////////////////////////////////
569//
570// Calculate distance to surface of shape from `inside', allowing for tolerance
571
572G4double G4Ellipsoid::DistanceToOut(const G4ThreeVector& p,
573 const G4ThreeVector& v,
574 const G4bool calcNorm,
575 G4bool *validNorm,
576 G4ThreeVector *n ) const
577{
578 G4double distMin;
579 enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
580
581 distMin= kInfinity;
582 surface= kNoSurf;
583
584 // check to see if Z plane is relevant
585 //
586 if (v.z() < 0.0)
587 {
588 G4double distZ = (zBottomCut - p.z()) / v.z();
589 if (distZ < 0.0)
590 {
591 distZ= 0.0;
592 if (!calcNorm) {return 0.0;}
593 }
594 distMin= distZ;
595 surface= kPlaneSurf;
596 }
597 if (v.z() > 0.0)
598 {
599 G4double distZ = (zTopCut - p.z()) / v.z();
600 if (distZ < 0.0)
601 {
602 distZ= 0.0;
603 if (!calcNorm) {return 0.0;}
604 }
605 distMin= distZ;
606 surface= kPlaneSurf;
607 }
608
609 // normal vector: parallel to normal, magnitude 1/(characteristic radius)
610 //
611 G4ThreeVector nearnorm(p.x()/(xSemiAxis*xSemiAxis),
612 p.y()/(ySemiAxis*ySemiAxis),
613 p.z()/(zSemiAxis*zSemiAxis));
614
615 // now check curved surface intercept
616 //
617 G4double A,B,C;
618
619 A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis);
620 C= (p * nearnorm) - 1.0;
621 B= 2.0 * (v * nearnorm);
622
623 C= B*B - 4.0*A*C;
624 if (C > 0.0)
625 {
626 G4double distR= (-B + std::sqrt(C) ) / (2.0*A);
627 if (distR < 0.0)
628 {
629 distR= 0.0;
630 if (!calcNorm) {return 0.0;}
631 }
632 if (distR < distMin)
633 {
634 distMin= distR;
635 surface= kCurvedSurf;
636 }
637 }
638
639 // set normal if requested
640 //
641 if (calcNorm)
642 {
643 if (surface == kNoSurf)
644 {
645 *validNorm = false;
646 }
647 else
648 {
649 *validNorm = true;
650 switch (surface)
651 {
652 case kPlaneSurf:
653 *n= G4ThreeVector(0.,0.,(v.z() > 1.0 ? 1. : -1.));
654 break;
655 case kCurvedSurf:
656 {
657 G4ThreeVector pexit= p + distMin*v;
658 G4ThreeVector truenorm(pexit.x()/(xSemiAxis*xSemiAxis),
659 pexit.y()/(ySemiAxis*ySemiAxis),
660 pexit.z()/(zSemiAxis*zSemiAxis));
661 truenorm *= 1.0/truenorm.mag();
662 *n= truenorm;
663 } break;
664 default:
665 G4cout.precision(16);
666 G4cout << G4endl;
667 DumpInfo();
668 G4cout << "Position:" << G4endl << G4endl;
669 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
670 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
671 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
672 G4cout << "Direction:" << G4endl << G4endl;
673 G4cout << "v.x() = " << v.x() << G4endl;
674 G4cout << "v.y() = " << v.y() << G4endl;
675 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
676 G4cout << "Proposed distance :" << G4endl << G4endl;
677 G4cout << "distMin = " << distMin/mm << " mm" << G4endl << G4endl;
678 G4Exception("G4Ellipsoid::DistanceToOut(p,v,..)",
679 "Notification", JustWarning,
680 "Undefined side for valid surface normal to solid.");
681 break;
682 }
683 }
684 }
685 return distMin;
686}
687
688///////////////////////////////////////////////////////////////////////////////
689//
690// Calculate distance (<=actual) to closest surface of shape from inside
691
692G4double G4Ellipsoid::DistanceToOut(const G4ThreeVector& p) const
693{
694 G4double distR, distZ;
695
696#ifdef G4SPECSDEBUG
697 if( Inside(p) == kOutside )
698 {
699 G4cout.precision(16) ;
700 G4cout << G4endl ;
701 DumpInfo();
702 G4cout << "Position:" << G4endl << G4endl ;
703 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
704 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
705 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
706 G4Exception("G4Ellipsoid::DistanceToOut(p)", "Notification", JustWarning,
707 "Point p is outside !?" );
708 }
709#endif
710
711 // Normal vector: parallel to normal, magnitude 1/(characteristic radius)
712 //
713 G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
714 p.y()/(ySemiAxis*ySemiAxis),
715 p.z()/(zSemiAxis*zSemiAxis));
716
717 // the following is a safe inlined "radius= min(1.0/norm.mag(),p.mag())
718 //
719 G4double radius= p.mag();
720 G4double tmp= norm.mag();
721 if ( (tmp > 0.0) && (1.0 < radius*tmp) ) {radius = 1.0/tmp;}
722
723 // Approximate distance to curved surface ( <= actual distance )
724 //
725 distR = (1.0 - p*norm) * radius / 2.0;
726
727 // Distance to z-cut plane
728 //
729 distZ = p.z() - zBottomCut;
730 if (distZ < 0.0) {distZ= zTopCut - p.z();}
731
732 // Distance to closest surface from inside
733 //
734 if ( (distZ < 0.0) || (distR < 0.0) )
735 {
736 return 0.0;
737 }
738 else
739 {
740 return (distZ < distR) ? distZ : distR;
741 }
742}
743
744///////////////////////////////////////////////////////////////////////////////
745//
746// Create a List containing the transformed vertices
747// Ordering [0-3] -fDz cross section
748// [4-7] +fDz cross section such that [0] is below [4],
749// [1] below [5] etc.
750// Note:
751// Caller has deletion resposibility
752// Potential improvement: For last slice, use actual ending angle
753// to avoid rounding error problems.
754
755G4ThreeVectorList*
756G4Ellipsoid::CreateRotatedVertices(const G4AffineTransform& pTransform,
757 G4int& noPolygonVertices) const
758{
759 G4ThreeVectorList *vertices;
760 G4ThreeVector vertex;
761 G4double meshAnglePhi, meshRMaxFactor,
762 crossAnglePhi, coscrossAnglePhi, sincrossAnglePhi, sAnglePhi;
763 G4double meshTheta, crossTheta, startTheta;
764 G4double rMaxX, rMaxY, rMaxZ, rMaxMax, rx, ry, rz;
765 G4int crossSectionPhi, noPhiCrossSections, crossSectionTheta, noThetaSections;
766
767 // Phi cross sections
768 //
769 noPhiCrossSections=G4int (twopi/kMeshAngleDefault)+1;
770
771 if (noPhiCrossSections<kMinMeshSections)
772 {
773 noPhiCrossSections=kMinMeshSections;
774 }
775 else if (noPhiCrossSections>kMaxMeshSections)
776 {
777 noPhiCrossSections=kMaxMeshSections;
778 }
779 meshAnglePhi=twopi/(noPhiCrossSections-1);
780
781 // Set start angle such that mesh will be at fRMax
782 // on the x axis. Will give better extent calculations when not rotated.
783
784 sAnglePhi = -meshAnglePhi*0.5;
785
786 // Theta cross sections
787
788 noThetaSections = G4int(pi/kMeshAngleDefault)+3;
789
790 if (noThetaSections<kMinMeshSections)
791 {
792 noThetaSections=kMinMeshSections;
793 }
794 else if (noThetaSections>kMaxMeshSections)
795 {
796 noThetaSections=kMaxMeshSections;
797 }
798 meshTheta= pi/(noThetaSections-2);
799
800 // Set start angle such that mesh will be at fRMax
801 // on the z axis. Will give better extent calculations when not rotated.
802
803 startTheta = -meshTheta*0.5;
804
805 meshRMaxFactor = 1.0/std::cos(0.5*
806 std::sqrt(meshAnglePhi*meshAnglePhi+meshTheta*meshTheta));
807 rMaxMax= (xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis);
808 if (zSemiAxis > rMaxMax) rMaxMax= zSemiAxis;
809 rMaxX= xSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
810 rMaxY= ySemiAxis + rMaxMax*(meshRMaxFactor-1.0);
811 rMaxZ= zSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
812 G4double* cosCrossTheta = new G4double[noThetaSections];
813 G4double* sinCrossTheta = new G4double[noThetaSections];
814 vertices=new G4ThreeVectorList(noPhiCrossSections*noThetaSections);
815 if (vertices && cosCrossTheta && sinCrossTheta)
816 {
817 for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
818 crossSectionTheta++)
819 {
820 // Compute sine and cosine table (for historical reasons)
821 //
822 crossTheta=startTheta+crossSectionTheta*meshTheta;
823 cosCrossTheta[crossSectionTheta]=std::cos(crossTheta);
824 sinCrossTheta[crossSectionTheta]=std::sin(crossTheta);
825 }
826 for (crossSectionPhi=0; crossSectionPhi<noPhiCrossSections;
827 crossSectionPhi++)
828 {
829 crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
830 coscrossAnglePhi=std::cos(crossAnglePhi);
831 sincrossAnglePhi=std::sin(crossAnglePhi);
832 for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
833 crossSectionTheta++)
834 {
835 // Compute coordinates of cross section at section crossSectionPhi
836 //
837 rx= sinCrossTheta[crossSectionTheta]*coscrossAnglePhi*rMaxX;
838 ry= sinCrossTheta[crossSectionTheta]*sincrossAnglePhi*rMaxY;
839 rz= cosCrossTheta[crossSectionTheta]*rMaxZ;
840 if (rz < zBottomCut)
841 { rz= zBottomCut; }
842 if (rz > zTopCut)
843 { rz= zTopCut; }
844 vertex= G4ThreeVector(rx,ry,rz);
845 vertices->push_back(pTransform.TransformPoint(vertex));
846 } // Theta forward
847 } // Phi
848 noPolygonVertices = noThetaSections ;
849 }
850 else
851 {
852 DumpInfo();
853 G4Exception("G4Ellipsoid::CreateRotatedVertices()",
854 "FatalError", FatalException,
855 "Error in allocation of vertices. Out of memory !");
856 }
857
858 delete[] cosCrossTheta;
859 delete[] sinCrossTheta;
860
861 return vertices;
862}
863
864//////////////////////////////////////////////////////////////////////////
865//
866// G4EntityType
867
868G4GeometryType G4Ellipsoid::GetEntityType() const
869{
870 return G4String("G4Ellipsoid");
871}
872
873//////////////////////////////////////////////////////////////////////////
874//
875// Stream object contents to an output stream
876
877std::ostream& G4Ellipsoid::StreamInfo( std::ostream& os ) const
878{
879 os << "-----------------------------------------------------------\n"
880 << " *** Dump for solid - " << GetName() << " ***\n"
881 << " ===================================================\n"
882 << " Solid type: G4Ellipsoid\n"
883 << " Parameters: \n"
884
885 << " semi-axis x: " << xSemiAxis/mm << " mm \n"
886 << " semi-axis y: " << ySemiAxis/mm << " mm \n"
887 << " semi-axis z: " << zSemiAxis/mm << " mm \n"
888 << " max semi-axis: " << semiAxisMax/mm << " mm \n"
889 << " lower cut plane level z: " << zBottomCut/mm << " mm \n"
890 << " upper cut plane level z: " << zTopCut/mm << " mm \n"
891 << "-----------------------------------------------------------\n";
892
893 return os;
894}
895
896////////////////////////////////////////////////////////////////////
897//
898// GetPointOnSurface
899
900G4ThreeVector G4Ellipsoid::GetPointOnSurface() const
901{
902 G4double aTop, aBottom, aCurved, chose, xRand, yRand, zRand, phi, theta;
903 G4double cosphi, sinphi, costheta, sintheta, alpha, beta, max1, max2, max3;
904
905 max1 = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis;
906 max1 = max1 > zSemiAxis ? max1 : zSemiAxis;
907 if(max1 == xSemiAxis){max2 = ySemiAxis; max3 = zSemiAxis;}
908 else if(max1 == ySemiAxis){max2 = xSemiAxis; max3 = zSemiAxis;}
909 else {max2 = xSemiAxis; max3 = ySemiAxis; }
910
911 phi = RandFlat::shoot(0.,2.*pi);
912 theta = RandFlat::shoot(0.,pi);
913
914 cosphi = std::cos(phi); sinphi = std::sin(phi);
915 costheta = RandFlat::shoot(zBottomCut,zTopCut)/zSemiAxis;
916 sintheta = std::sqrt(1.-sqr(costheta));
917
918 alpha = 1.-sqr(max2/max1); beta = 1.-sqr(max3/max1);
919
920 aTop = pi*xSemiAxis*ySemiAxis*(1 - sqr(zTopCut/zSemiAxis));
921 aBottom = pi*xSemiAxis*ySemiAxis*(1 - sqr(zBottomCut/zSemiAxis));
922
923 // approximation
924 // from:" http://www.citr.auckland.ac.nz/techreports/2004/CITR-TR-139.pdf"
925 aCurved = 4.*pi*max1*max2*(1.-1./6.*(alpha+beta)-
926 1./120.*(3.*sqr(alpha)+2.*alpha*beta+3.*sqr(beta)));
927
928 aCurved *= 0.5*(1.2*zTopCut/zSemiAxis - 1.2*zBottomCut/zSemiAxis);
929
930 if( ( zTopCut >= zSemiAxis && zBottomCut <= -1.*zSemiAxis )
931 || ( zTopCut == 0 && zBottomCut ==0 ) )
932 {
933 aTop = 0; aBottom = 0;
934 }
935
936 chose = RandFlat::shoot(0.,aTop + aBottom + aCurved);
937
938 if(chose < aCurved)
939 {
940 xRand = xSemiAxis*sintheta*cosphi;
941 yRand = ySemiAxis*sintheta*sinphi;
942 zRand = zSemiAxis*costheta;
943 return G4ThreeVector (xRand,yRand,zRand);
944 }
945 else if(chose >= aCurved && chose < aCurved + aTop)
946 {
947 xRand = RandFlat::shoot(-1.,1.)*xSemiAxis
948 * std::sqrt(1-sqr(zTopCut/zSemiAxis));
949 yRand = RandFlat::shoot(-1.,1.)*ySemiAxis
950 * std::sqrt(1.-sqr(zTopCut/zSemiAxis)-sqr(xRand/xSemiAxis));
951 zRand = zTopCut;
952 return G4ThreeVector (xRand,yRand,zRand);
953 }
954 else
955 {
956 xRand = RandFlat::shoot(-1.,1.)*xSemiAxis
957 * std::sqrt(1-sqr(zBottomCut/zSemiAxis));
958 yRand = RandFlat::shoot(-1.,1.)*ySemiAxis
959 * std::sqrt(1.-sqr(zBottomCut/zSemiAxis)-sqr(xRand/xSemiAxis));
960 zRand = zBottomCut;
961 return G4ThreeVector (xRand,yRand,zRand);
962 }
963}
964
965/////////////////////////////////////////////////////////////////////////////
966//
967// Methods for visualisation
968
969void G4Ellipsoid::DescribeYourselfTo (G4VGraphicsScene& scene) const
970{
971 scene.AddSolid(*this);
972}
973
974G4VisExtent G4Ellipsoid::GetExtent() const
975{
976 // Define the sides of the box into which the G4Ellipsoid instance would fit.
977 //
978 return G4VisExtent (-semiAxisMax, semiAxisMax,
979 -semiAxisMax, semiAxisMax,
980 -semiAxisMax, semiAxisMax);
981}
982
983G4NURBS* G4Ellipsoid::CreateNURBS () const
984{
985 // Box for now!!!
986 //
987 return new G4NURBSbox(semiAxisMax, semiAxisMax, semiAxisMax);
988}
989
990G4Polyhedron* G4Ellipsoid::CreatePolyhedron () const
991{
992 return new G4PolyhedronEllipsoid(xSemiAxis, ySemiAxis, zSemiAxis,
993 zBottomCut, zTopCut);
994}
995
996G4Polyhedron* G4Ellipsoid::GetPolyhedron () const
997{
998 if (!fpPolyhedron ||
999 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1000 fpPolyhedron->GetNumberOfRotationSteps())
1001 {
1002 delete fpPolyhedron;
1003 fpPolyhedron = CreatePolyhedron();
1004 }
1005 return fpPolyhedron;
1006}
Note: See TracBrowser for help on using the repository browser.