00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include "G4Ray.hh"
00037 #include "G4PointRat.hh"
00038
00039 G4Ray::G4Ray()
00040 : r_min(0.), r_max(0.)
00041 {
00042 }
00043
00044 G4Ray::G4Ray(const G4Point3D& start0, const G4Vector3D& dir0)
00045 {
00046 Init(start0, dir0);
00047 }
00048
00049 G4Ray::~G4Ray()
00050 {
00051 }
00052
00053
00054 const G4Plane& G4Ray::GetPlane(G4int number_of_plane) const
00055 {
00056 if(number_of_plane==1)
00057 { return plane2; }
00058 else
00059 { return plane1; }
00060 }
00061
00062
00063 void G4Ray::CreatePlanes()
00064 {
00065
00066
00067
00068
00069 G4Vector3D RayDir = dir;
00070 G4Point3D RayOrigin = start;
00071
00072 G4Point3D p1, p2, p3, p4;
00073 G4Vector3D dir1, dir2;
00074 G4Vector3D invdir = G4Vector3D( PINFINITY );
00075
00076 if(!NearZero(RayDir.x(), SQRT_SMALL_FASTF))
00077 { invdir.setX(1.0 / RayDir.x()); }
00078
00079 if(!NearZero(RayDir.y(), SQRT_SMALL_FASTF))
00080 { invdir.setY(1.0 / RayDir.y()); }
00081
00082 if(!NearZero(RayDir.z(), SQRT_SMALL_FASTF))
00083 { invdir.setZ(1.0 / RayDir.z()); }
00084
00085 MatVecOrtho(dir1, RayDir);
00086
00087 Vcross( dir2, RayDir, dir1);
00088 Vmove(p1, RayOrigin);
00089 Vadd2(p2, RayOrigin, RayDir);
00090 Vadd2(p3, RayOrigin, dir1);
00091 Vadd2(p4, RayOrigin, dir2);
00092
00093 CalcPlane3Pts( plane1, p1, p3, p2);
00094 CalcPlane3Pts( plane2, p1, p2, p4);
00095 }
00096
00097
00098 void G4Ray::MatVecOrtho(register G4Vector3D &out,
00099 register const G4Vector3D &in )
00100 {
00101 register G4double f;
00102 G4int i_Which;
00103
00104 if( NearZero(in.x(), 0.0001)
00105 && NearZero(in.y(), 0.0001)
00106 && NearZero(in.z(), 0.0001) )
00107 {
00108 Vsetall( out, 0 );
00109 return;
00110 }
00111
00112
00113 f = std::fabs(in.x());
00114 i_Which=0;
00115
00116 if( std::fabs(in.y()) < f )
00117 {
00118 f = std::fabs(in.y());
00119 i_Which=1;
00120 }
00121
00122 if( std::fabs(in.z()) < f )
00123 {
00124 i_Which=2;
00125 }
00126
00127 if(!i_Which)
00128 {
00129 f = std::sqrt((in.y())*(in.y())+(in.z())*(in.z()));
00130 }
00131 else
00132 {
00133 if(i_Which==1)
00134 {
00135 f = std::sqrt((in.z())*(in.z())+(in.x())*(in.x()));
00136 }
00137 else
00138 {
00139 f = std::sqrt((in.x())*(in.x())+(in.y())*(in.y()));
00140 }
00141 }
00142 if( NearZero( f, SMALL ) )
00143 {
00144 Vsetall( out, 0 );
00145 return;
00146 }
00147
00148 f = 1.0/f;
00149
00150 if(!i_Which)
00151 {
00152 out.setX(0.0);
00153 out.setY(-in.z()*f);
00154 out.setZ( in.y()*f);
00155 }
00156 else
00157 {
00158 if(i_Which==1)
00159 {
00160 out.setY(0.0);
00161 out.setZ(-in.x()*f);
00162 out.setX( in.y()*f);
00163 }
00164 else
00165 {
00166 out.setZ(0.0);
00167 out.setX(-in.z()*f);
00168 out.setY( in.y()*f);
00169 }
00170 }
00171 }
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212 G4int G4Ray::CalcPlane3Pts(G4Plane &plane1,
00213 const G4Point3D& a,
00214 const G4Point3D& b,
00215 const G4Point3D& c )
00216 {
00217
00218
00219
00220 G4Vector3D B_A;
00221 G4Vector3D C_A;
00222 G4Vector3D C_B;
00223
00224 register G4double mag;
00225
00226 Vsub2( B_A, b, a );
00227 Vsub2( C_A, c, a );
00228 Vsub2( C_B, c, b );
00229
00230 Vcross( plane1, B_A, C_A );
00231
00232
00233 mag = Magnitude(plane1);
00234 if( mag <= SQRT_SMALL_FASTF )
00235 {
00236 return(-1);
00237 }
00238
00239 mag = 1/mag;
00240
00241 G4Plane pl2(plane1);
00242 Vscale( plane1, pl2, mag );
00243
00244
00245 plane1.d = Vdot( plane1, a );
00246
00247 return(0);
00248 }
00249
00250
00251 void G4Ray::RayCheck()
00252 {
00253
00254 if (dir==G4Vector3D(0, 0, 0))
00255 {
00256 G4Exception("G4Ray::RayCheck()", "GeomSolids0002", FatalException,
00257 "Invalid zero direction given !");
00258 }
00259
00260
00261 dir= dir.unit();
00262 r_min = 0;
00263 r_max = 0;
00264 }
00265
00266
00267 void G4Ray::Vcross(G4Plane &a, const G4Vector3D &b, const G4Vector3D &c)
00268 {
00269 a.a = b.y() * c.z() - b.z() * c.y() ;
00270 a.b = b.z() * c.x() - b.x() * c.z() ;
00271 a.c = b.x() * c.y() - b.y() * c.x() ;
00272 }
00273
00274
00275 void G4Ray::Vcross(G4Vector3D &a, const G4Vector3D &b, const G4Vector3D &c)
00276 {
00277 a.setX(b.y() * c.z() - b.z() * c.y()) ;
00278 a.setY(b.z() * c.x() - b.x() * c.z()) ;
00279 a.setZ(b.x() * c.y() - b.y() * c.x()) ;
00280 }
00281
00282
00283 void G4Ray::Vmove(G4Point3D &a, const G4Point3D &b)
00284 {
00285 a.setX(b.x());
00286 a.setY(b.y());
00287 a.setZ(b.z());
00288 }
00289
00290
00291 void G4Ray::Vadd2(G4Point3D &a, const G4Point3D &b, const G4Vector3D &c)
00292 {
00293 a.setX(b.x() + c.x()) ;
00294 a.setY(b.y() + c.y()) ;
00295 a.setZ(b.z() + c.z()) ;
00296 }
00297
00298
00299 void G4Ray::Vsub2(G4Vector3D &a, const G4Point3D &b, const G4Point3D &c)
00300 {
00301 a.setX(b.x() - c.x());
00302 a.setY(b.y() - c.y());
00303 a.setZ(b.z() - c.z());
00304 }
00305
00306
00307 void G4Ray::Vscale(G4Plane& a, const G4Plane& b, G4double c)
00308 {
00309 a.a = b.a * c;
00310 a.b = b.b * c;
00311 a.c = b.c * c;
00312 }
00313
00314
00315 G4double G4Ray::Vdot(const G4Plane &a, const G4Point3D &b)
00316 {
00317 return (a.a * b.x() +
00318 a.b * b.y() +
00319 a.c * b.z());
00320 }
00321
00322
00323 G4double G4Ray::Magsq(const G4Plane &a)
00324 {
00325 return ( a.a * a.a + a.b * a.b + a.c *a.c );
00326 }
00327
00328
00329 G4double G4Ray::Magnitude(const G4Plane &a)
00330 {
00331 return (std::sqrt( Magsq( a )) );
00332 }