source: trunk/source/geometry/solids/BREPS/src/G4FConicalSurface.cc@ 1310

Last change on this file since 1310 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 12.1 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//
27// $Id: G4FConicalSurface.cc,v 1.19 2006/06/29 18:42:12 gunter Exp $
28// GEANT4 tag $Name: geant4-09-03 $
29//
30// ----------------------------------------------------------------------
31// GEANT 4 class source file
32//
33// G4FConicalSurface.cc
34//
35// ----------------------------------------------------------------------
36
37#include "G4FConicalSurface.hh"
38#include "G4Sort.hh"
39#include "G4CircularCurve.hh"
40
41
42G4FConicalSurface::G4FConicalSurface()
43{
44 length = 1.0;
45 small_radius = 0.0;
46 large_radius = 1.0;
47 tan_angle = (large_radius-small_radius)/length;
48}
49
50G4FConicalSurface::~G4FConicalSurface()
51{
52}
53
54G4FConicalSurface::G4FConicalSurface(const G4Point3D& o,
55 const G4Vector3D& a,
56 G4double l,
57 G4double sr,
58 G4double lr
59 )
60{
61 // Make a G4FConicalSurface with origin o, axis a, length l, small radius
62 // sr, and large radius lr. The angle is calculated below and the SetAngle
63 // function of G4ConicalSurface is used to set it properly from the default
64 // value used above in the initialization.
65
66 // Create the position with origin o, axis a, and a direction
67
68 G4Vector3D dir(1,1,1);
69 Position.Init(dir, a, o);
70 origin = o;
71
72 // Require length to be nonnegative
73 if (l >=0)
74 length = l;
75 else
76 {
77 G4cerr << "Error in G4FConicalSurface::G4FConicalSurface"
78 << "--asked for negative length\n"
79 << "\tDefault length of 0.0 is used.\n";
80
81 length = 0.0;
82 }
83
84 // Require small radius to be non-negative (i.e., allow zero)
85 if ( sr >= 0.0 )
86 small_radius = sr;
87 else
88 {
89 G4cerr << "Error in G4FConicalSurface::G4FConicalSurface"
90 << "--asked for negative small radius\n"
91 << "\tDefault value of 0.0 is used.\n";
92
93 small_radius = 0.0;
94 }
95
96 // Require large radius to exceed small radius
97 if ( lr > small_radius )
98 large_radius = lr;
99 else
100 {
101 G4cerr << "Error in G4FConicalSurface::G4FConicalSurface"
102 << "--large radius must exceed small radius\n"
103 << "\tDefault value of small radius +1 is used.\n";
104
105 large_radius = small_radius + 1.0;
106 }
107
108 // Calculate the angle of the G4ConicalSurface from the length and radii
109 tan_angle = ( large_radius - small_radius ) / length ;
110}
111
112
113const char* G4FConicalSurface::Name() const
114{
115 return "G4FConicalSurface";
116}
117
118// Modified by L. Broglia (01/12/98)
119void G4FConicalSurface::CalcBBox()
120{
121 G4Point3D Max = G4Point3D(-PINFINITY);
122 G4Point3D Min = G4Point3D( PINFINITY);
123 G4Point3D Tmp;
124
125 G4Point3D Origin = Position.GetLocation();
126 G4Point3D EndOrigin = G4Point3D( Origin + (length * Position.GetAxis()) );
127
128 G4double radius = large_radius;
129 G4Point3D Radius(radius, radius, 0);
130
131 // Default BBox
132 G4Point3D Tolerance(kCarTolerance, kCarTolerance, kCarTolerance);
133 G4Point3D BoxMin(Origin-Tolerance);
134 G4Point3D BoxMax(Origin+Tolerance);
135
136 bbox = new G4BoundingBox3D();
137 bbox->Init(BoxMin, BoxMax);
138
139 Tmp = (Origin - Radius);
140 bbox->Extend(Tmp);
141
142 Tmp = Origin + Radius;
143 bbox->Extend(Tmp);
144
145 Tmp = EndOrigin - Radius;
146 bbox->Extend(Tmp);
147
148 Tmp = EndOrigin + Radius;
149 bbox->Extend(Tmp);
150}
151
152
153void G4FConicalSurface::PrintOn( std::ostream& os ) const
154{
155 // printing function using C++ std::ostream class
156 os << "G4FConicalSurface with origin: " << origin << "\t"
157 << "and axis: " << Position.GetAxis() << "\n"
158 << "\t small radius: " << small_radius
159 << "\t large radius: " << large_radius
160 << "\t and length: " << length << "\n";
161}
162
163
164G4int G4FConicalSurface::operator==( const G4FConicalSurface& c ) const
165{
166 return ( origin == c.origin &&
167 Position.GetAxis() == c.Position.GetAxis() &&
168 small_radius == c.small_radius &&
169 large_radius == c.large_radius &&
170 length == c.length &&
171 tan_angle == c.tan_angle );
172}
173
174
175G4int G4FConicalSurface::WithinBoundary( const G4Vector3D& x ) const
176{
177 // return 1 if point x is within the boundaries of the G4FConicalSurface
178 // return 0 otherwise (assume it is on the G4ConicalSurface)
179 G4Vector3D q = G4Vector3D( x - origin );
180
181 G4double qmag = q.mag();
182 G4double s = std::sin( std::atan2(large_radius-small_radius, length) );
183 G4double ls = small_radius / s;
184 G4double ll = large_radius / s;
185
186 if ( ( qmag >= ls ) && ( qmag <= ll ) )
187 return 1;
188 else
189 return 0;
190}
191
192
193G4double G4FConicalSurface::Scale() const
194{
195 // Returns the small radius of a G4FConicalSurface unless it is zero, in
196 // which case returns the large radius.
197 // Used for Scale-invariant tests of surface thickness.
198 if ( small_radius == 0.0 )
199 return large_radius;
200 else
201 return small_radius;
202}
203
204
205G4double G4FConicalSurface::Area() const
206{
207 // Returns the Area of a G4FConicalSurface
208 G4double rdif = large_radius - small_radius;
209
210 return ( pi * ( small_radius + large_radius ) *
211 std::sqrt( length * length + rdif * rdif ) );
212}
213
214
215void G4FConicalSurface::resize( G4double l, G4double sr, G4double lr )
216{
217 // Resize a G4FConicalSurface to a new length l, and new radii sr and lr.
218 // Must Reset angle of the G4ConicalSurface as well based on these new
219 // values.
220 // Require length to be non-negative
221
222 // if ( l > 0.0 )
223 if ( l >= 0.0 )
224 length = l;
225 else
226 {
227 G4cerr << "Error in G4FConicalSurface::resize"
228 << "--asked for negative length\n"
229 << "\tOriginal value of " << length << " is retained.\n";
230 }
231
232 // Require small radius to be non-negative (i.e., allow zero)
233 if ( sr >= 0.0 )
234 small_radius = sr;
235 else
236 {
237 G4cerr << "Error in G4FConicalSurface::resize"
238 << "--asked for negative small radius\n"
239 << "\tOriginal value of " << small_radius
240 << " is retained.\n";
241 }
242
243 // Require large radius to exceed small radius
244 if ( lr > small_radius )
245 large_radius = lr;
246 else
247 {
248 G4double r = small_radius + 1.0;
249 lr = ( large_radius <= small_radius ) ? r : large_radius;
250 large_radius = lr;
251
252 G4cerr << "Error in G4FConicalSurface::G4FConicalSurface"
253 << "--large radius must exceed small radius\n"
254 << "\tDefault value of " << large_radius << " is used.\n";
255 }
256
257 // Calculate the angle of the G4ConicalSurface from the length and radii
258 tan_angle = ( large_radius - small_radius ) / length ;
259
260}
261
262
263G4int G4FConicalSurface::Intersect(const G4Ray& ry )
264{
265 // This function count the number of intersections of a
266 // bounded conical surface by a ray.
267 // At first, calculates the intersections with the semi-infinite
268 // conical surfsace. After, count the intersections within the
269 // finite conical surface boundaries, and set "distance" to the
270 // closest distance from the start point to the nearest intersection
271 // If the point is on the surface it returns or the intersection with
272 // the opposite surface or kInfinity
273 // If no intersection is founded, set distance = kInfinity and
274 // return 0
275
276 distance = kInfinity;
277 closest_hit = PINFINITY;
278
279 // origin and direction of the ray
280 G4Point3D x = ry.GetStart();
281 G4Vector3D dhat = ry.GetDir();
282
283 // cone angle and axis
284 G4double ta = tan_angle;
285 G4Vector3D ahat = Position.GetAxis();
286
287 // array of solutions in distance along the ray
288 G4double s[2];
289 s[0]=-1.0;
290 s[1]=-1.0;
291
292 // calculate the two intersections (quadratic equation)
293 G4Vector3D gamma = G4Vector3D( x - Position.GetLocation() );
294
295 G4double t = 1 + ta * ta;
296 G4double ga = gamma * ahat;
297 G4double da = dhat * ahat;
298
299 G4double A = t * da * da - dhat * dhat;
300 G4double B = 2 * ( -gamma * dhat + t * ga * da - large_radius * ta * da);
301 G4double C = ( -gamma * gamma + t * ga * ga
302 - 2 * large_radius * ta * ga
303 + large_radius * large_radius );
304
305 G4double radical = B * B - 4.0 * A * C;
306
307 if ( radical < 0.0 )
308 // no intersection
309 return 0;
310 else
311 {
312 G4double root = std::sqrt( radical );
313 s[0] = ( - B + root ) / ( 2. * A );
314 s[1] = ( - B - root ) / ( 2. * A );
315 }
316
317 // validity of the solutions
318 // the hit point must be into the bounding box of the conical surface
319 G4Point3D p0 = G4Point3D( x + s[0]*dhat );
320 G4Point3D p1 = G4Point3D( x + s[1]*dhat );
321
322 if( !GetBBox()->Inside(p0) )
323 s[0] = kInfinity;
324
325 if( !GetBBox()->Inside(p1) )
326 s[1] = kInfinity;
327
328 // now loop over each positive solution, keeping the first one (smallest
329 // distance along the ray) which is within the boundary of the sub-shape
330 G4int nbinter = 0;
331 distance = kInfinity;
332
333 for ( G4int i = 0; i < 2; i++ )
334 {
335 if(s[i] < kInfinity) {
336 if ( (s[i] > kCarTolerance*0.5) ) {
337 nbinter++;
338 if ( distance > (s[i]*s[i]) ) {
339 distance = s[i]*s[i];
340 }
341 }
342 }
343 }
344
345 return nbinter;
346}
347
348
349G4double G4FConicalSurface::HowNear( const G4Vector3D& x ) const
350{
351// Shortest distance from the point x to the G4FConicalSurface.
352// The distance will be always positive
353// This function works only with Cone axis equal (0,0,1) or (0,0,-1), it project
354// the surface and the point on the x,z plane and compute the distance in analytical
355// way
356
357 G4double hownear ;
358
359 G4Vector3D upcorner = G4Vector3D ( small_radius, 0 , origin.z()+Position.GetAxis().z()*length);
360 G4Vector3D downcorner = G4Vector3D ( large_radius, 0 , origin.z());
361 G4Vector3D xd;
362
363 xd = G4Vector3D ( std::sqrt ( x.x()*x.x() + x.y()*x.y() ) , 0 , x.z() );
364
365 G4double m = (upcorner.z() - downcorner.z()) / (upcorner.x() - downcorner.x());
366 G4double q = (downcorner.z()*upcorner.x() - upcorner.z()*downcorner.x()) /
367 (upcorner.x() - downcorner.x());
368
369 G4double Zinter = (xd.z()*m*m + xd.x()*m +q)/(1+m*m) ;
370
371 if ( ((Zinter >= downcorner.z()) && (Zinter <=upcorner.z())) ||
372 ((Zinter >= upcorner.z()) && (Zinter <=downcorner.z())) ) {
373 hownear = std::fabs(m*xd.x()-xd.z()+q)/std::sqrt(1+m*m);
374 return hownear;
375 } else {
376 hownear = std::min ( (xd-upcorner).mag() , (xd-downcorner).mag() );
377 return hownear;
378 }
379
380
381}
382
383
384G4Vector3D G4FConicalSurface::SurfaceNormal( const G4Point3D& p ) const
385{
386 // return the Normal unit vector to the G4ConicalSurface at a point p
387 // on (or nearly on) the G4ConicalSurface
388 G4Vector3D s = G4Vector3D( p - origin );
389 G4double da = s * Position.GetAxis();
390 G4double r = std::sqrt( s*s - da*da);
391 G4double z = tan_angle * r;
392
393 if (Position.GetAxis().z() < 0)
394 z = -z;
395
396 G4Vector3D n(p.x(), p.y(), z);
397 n = n.unit();
398
399 if( !sameSense )
400 n = -n;
401
402 return n;
403}
404
405G4int G4FConicalSurface::Inside ( const G4Vector3D& x ) const
406{
407 // Return 0 if point x is outside G4ConicalSurface, 1 if Inside.
408 if ( HowNear( x ) >= -0.5*kCarTolerance )
409 return 1;
410 else
411 return 0;
412}
413
Note: See TracBrowser for help on using the repository browser.