source: trunk/source/geometry/solids/BREPS/src/G4Ray.cc

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

tag geant4.9.4 beta 1 + modifs locales
File size: 7.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//
27// $Id: G4Ray.cc,v 1.12 2008/07/08 10:00:58 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
29//
30// ----------------------------------------------------------------------
31// GEANT 4 class source file
32//
33// G4Ray.cc
34//
35// ----------------------------------------------------------------------
36
37#include "G4Ray.hh"
38#include "G4PointRat.hh"
39
40G4Ray::G4Ray()
41{
42}
43
44G4Ray::G4Ray(const G4Point3D& start0, const G4Vector3D& dir0)
45{
46 Init(start0, dir0);
47}
48
49G4Ray::~G4Ray()
50{
51}
52
53
54const G4Plane& G4Ray::GetPlane(G4int number_of_plane) const
55{
56 if(number_of_plane==1)
57 { return plane2; }
58 else
59 { return plane1; }
60}
61
62
63void G4Ray::CreatePlanes()
64{
65 // Creates two orthogonal planes(plane1,plane2) the ray (rray)
66 // situated in the intersection of the planes. The planes are
67 // used to project the surface (nurb) in two dimensions.
68
69 G4Vector3D RayDir = dir;
70 G4Point3D RayOrigin = start;
71
72 G4Point3D p1, p2, p3, p4;
73 G4Vector3D dir1, dir2;
74 G4Vector3D invdir = G4Vector3D( PINFINITY );
75
76 if(!NearZero(RayDir.x(), SQRT_SMALL_FASTF))
77 { invdir.setX(1.0 / RayDir.x()); }
78
79 if(!NearZero(RayDir.y(), SQRT_SMALL_FASTF))
80 { invdir.setY(1.0 / RayDir.y()); }
81
82 if(!NearZero(RayDir.z(), SQRT_SMALL_FASTF))
83 { invdir.setZ(1.0 / RayDir.z()); }
84
85 MatVecOrtho(dir1, RayDir);
86
87 Vcross( dir2, RayDir, dir1);
88 Vmove(p1, RayOrigin);
89 Vadd2(p2, RayOrigin, RayDir);
90 Vadd2(p3, RayOrigin, dir1);
91 Vadd2(p4, RayOrigin, dir2);
92
93 CalcPlane3Pts( plane1, p1, p3, p2);
94 CalcPlane3Pts( plane2, p1, p2, p4);
95}
96
97
98void G4Ray::MatVecOrtho(register G4Vector3D &out,
99 register const G4Vector3D &in )
100{
101 register G4double f;
102 G4int i_Which;
103
104 if( NearZero(in.x(), 0.0001)
105 && NearZero(in.y(), 0.0001)
106 && NearZero(in.z(), 0.0001) )
107 {
108 Vsetall( out, 0 );
109 return;
110 }
111
112 // Find component closest to zero
113 f = std::fabs(in.x());
114 i_Which=0;
115
116 if( std::fabs(in.y()) < f )
117 {
118 f = std::fabs(in.y());
119 i_Which=1;
120 }
121
122 if( std::fabs(in.z()) < f )
123 {
124 i_Which=2;
125 }
126
127 if(!i_Which)
128 {
129 f = std::sqrt((in.y())*(in.y())+(in.z())*(in.z())); // hypot(in.y(),in.z())
130 }
131 else
132 {
133 if(i_Which==1)
134 {
135 f = std::sqrt((in.z())*(in.z())+(in.x())*(in.x())); // hypot(in.z(),in.x())
136 }
137 else
138 {
139 f = std::sqrt((in.x())*(in.x())+(in.y())*(in.y())); // hypot(in.x(),in.y())
140 }
141 }
142 if( NearZero( f, SMALL ) )
143 {
144 Vsetall( out, 0 );
145 return;
146 }
147
148 f = 1.0/f;
149
150 if(!i_Which)
151 {
152 out.setX(0.0);
153 out.setY(-in.z()*f);
154 out.setZ( in.y()*f);
155 }
156 else
157 {
158 if(i_Which==1)
159 {
160 out.setY(0.0);
161 out.setZ(-in.x()*f);
162 out.setX( in.y()*f);
163 }
164 else
165 {
166 out.setZ(0.0);
167 out.setX(-in.z()*f);
168 out.setY( in.y()*f);
169 }
170 }
171}
172
173
174// CALC_PLANE_3PTS
175//
176// Find the equation of a G4Plane that contains three points.
177// Note that Normal vector created is expected to point out (see vmath.h),
178// so the vector from A to C had better be counter-clockwise
179// (about the point A) from the vector from A to B.
180// This follows the outward-pointing Normal convention, and the
181// right-hand rule for cross products.
182//
183/*
184 C
185 *
186 |\
187 | \
188 ^ N | \
189 | \ | \
190 | \ | \
191 |C-A \ | \
192 | \ | \
193 | \ | \
194 \| \
195 *---------*
196 A B
197 ----->
198 B-A
199*/
200// If the points are given in the order A B C (eg, *counter*-clockwise),
201// then the outward pointing surface Normal N = (B-A) x (C-A).
202//
203// Explicit Return -
204// 0 OK
205// -1 Failure. At least two of the points were not distinct,
206// or all three were colinear.
207//
208// Implicit Return -
209// G4Plane The G4Plane equation is stored here.
210
211
212G4int G4Ray::CalcPlane3Pts(G4Plane &plane1,
213 const G4Point3D& a,
214 const G4Point3D& b,
215 const G4Point3D& c )
216{
217 // Creates the two orthogonal planes which are needed in projecting the
218 // surface into 2D.
219
220 G4Vector3D B_A;
221 G4Vector3D C_A;
222 G4Vector3D C_B;
223
224 register G4double mag;
225
226 Vsub2( B_A, b, a );
227 Vsub2( C_A, c, a );
228 Vsub2( C_B, c, b );
229
230 Vcross( plane1, B_A, C_A );
231
232 // Ensure unit length Normal
233 mag = Magnitude(plane1);
234 if( mag <= SQRT_SMALL_FASTF )
235 {
236 return(-1);// FAIL
237 }
238
239 mag = 1/mag;
240
241 G4Plane pl2(plane1);
242 Vscale( plane1, pl2, mag );
243
244 // Find distance from the origin to the G4Plane
245 plane1.d = Vdot( plane1, a );
246
247 return(0); //ok
248}
249
250
251void G4Ray::RayCheck()
252{
253 // Check that the ray has a G4Vector3D...
254 if (dir==G4Vector3D(0, 0, 0))
255 {
256 G4Exception("G4Ray::RayCheck()", "InvalidInput", FatalException,
257 "Invalid zero direction given !");
258 }
259
260 // Make sure that the vector is unit length
261 dir= dir.unit();
262 r_min = 0;
263 r_max = 0;
264}
265
266
267void G4Ray::Vcross(G4Plane &a, const G4Vector3D &b, const G4Vector3D &c)
268{
269 a.a = b.y() * c.z() - b.z() * c.y() ;
270 a.b = b.z() * c.x() - b.x() * c.z() ;
271 a.c = b.x() * c.y() - b.y() * c.x() ;
272}
273
274
275void G4Ray::Vcross(G4Vector3D &a, const G4Vector3D &b, const G4Vector3D &c)
276{
277 a.setX(b.y() * c.z() - b.z() * c.y()) ;
278 a.setY(b.z() * c.x() - b.x() * c.z()) ;
279 a.setZ(b.x() * c.y() - b.y() * c.x()) ;
280}
281
282
283void G4Ray::Vmove(G4Point3D &a, const G4Point3D &b)
284{
285 a.setX(b.x());
286 a.setY(b.y());
287 a.setZ(b.z());
288}
289
290
291void G4Ray::Vadd2(G4Point3D &a, const G4Point3D &b, const G4Vector3D &c)
292{
293 a.setX(b.x() + c.x()) ;
294 a.setY(b.y() + c.y()) ;
295 a.setZ(b.z() + c.z()) ;
296}
297
298
299void G4Ray::Vsub2(G4Vector3D &a, const G4Point3D &b, const G4Point3D &c)
300{
301 a.setX(b.x() - c.x());
302 a.setY(b.y() - c.y());
303 a.setZ(b.z() - c.z());
304}
305
306
307void G4Ray::Vscale(G4Plane& a, const G4Plane& b, G4double c)
308{
309 a.a = b.a * c;
310 a.b = b.b * c;
311 a.c = b.c * c;
312}
313
314
315G4double G4Ray::Vdot(const G4Plane &a, const G4Point3D &b)
316{
317 return (a.a * b.x() +
318 a.b * b.y() +
319 a.c * b.z());
320}
321
322
323G4double G4Ray::Magsq(const G4Plane &a)
324{
325 return ( a.a * a.a + a.b * a.b + a.c *a.c );
326}
327
328
329G4double G4Ray::Magnitude(const G4Plane &a)
330{
331 return (std::sqrt( Magsq( a )) );
332}
Note: See TracBrowser for help on using the repository browser.