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 "G4FConicalSurface.hh"
00037 #include "G4PhysicalConstants.hh"
00038 #include "G4Sort.hh"
00039 #include "G4CircularCurve.hh"
00040
00041
00042 G4FConicalSurface::G4FConicalSurface()
00043 {
00044 length = 1.0;
00045 small_radius = 0.0;
00046 large_radius = 1.0;
00047 tan_angle = (large_radius-small_radius)/length;
00048 }
00049
00050 G4FConicalSurface::~G4FConicalSurface()
00051 {
00052 }
00053
00054 G4FConicalSurface::G4FConicalSurface(const G4Point3D& o,
00055 const G4Vector3D& a,
00056 G4double l,
00057 G4double srad,
00058 G4double lr
00059 )
00060 {
00061
00062
00063
00064
00065
00066
00067
00068 G4Vector3D dir(1,1,1);
00069 Position.Init(dir, a, o);
00070 origin = o;
00071
00072
00073 if (l >=0)
00074 length = l;
00075 else
00076 {
00077 std::ostringstream message;
00078 message << "Negative length." << G4endl
00079 << "Default length of 0.0 is used.";
00080 G4Exception("G4FConicalSurface::G4FConicalSurface()",
00081 "GeomSolids1001", JustWarning, message);
00082
00083 length = 0.0;
00084 }
00085
00086
00087 if ( srad >= 0.0 )
00088 small_radius = srad;
00089 else
00090 {
00091 std::ostringstream message;
00092 message << "Negative small radius." << G4endl
00093 << "Default value of 0.0 is used.";
00094 G4Exception("G4FConicalSurface::G4FConicalSurface()",
00095 "GeomSolids1001", JustWarning, message);
00096
00097 small_radius = 0.0;
00098 }
00099
00100
00101 if ( lr > small_radius )
00102 large_radius = lr;
00103 else
00104 {
00105 std::ostringstream message;
00106 message << "Large radius must exceed small radius" << G4endl
00107 << "Default value of small radius +1 is used.";
00108 G4Exception("G4FConicalSurface::G4FConicalSurface()",
00109 "GeomSolids1001", JustWarning, message);
00110
00111 large_radius = small_radius + 1.0;
00112 }
00113
00114
00115 tan_angle = ( large_radius - small_radius ) / length ;
00116 }
00117
00118
00119 const char* G4FConicalSurface::Name() const
00120 {
00121 return "G4FConicalSurface";
00122 }
00123
00124
00125 void G4FConicalSurface::CalcBBox()
00126 {
00127 G4Point3D Max = G4Point3D(-PINFINITY);
00128 G4Point3D Min = G4Point3D( PINFINITY);
00129 G4Point3D Tmp;
00130
00131 G4Point3D Origin = Position.GetLocation();
00132 G4Point3D EndOrigin = G4Point3D( Origin + (length * Position.GetAxis()) );
00133
00134 G4double radius = large_radius;
00135 G4Point3D Radius(radius, radius, 0);
00136
00137
00138 G4Point3D Tolerance(kCarTolerance, kCarTolerance, kCarTolerance);
00139 G4Point3D BoxMin(Origin-Tolerance);
00140 G4Point3D BoxMax(Origin+Tolerance);
00141
00142 bbox = new G4BoundingBox3D();
00143 bbox->Init(BoxMin, BoxMax);
00144
00145 Tmp = (Origin - Radius);
00146 bbox->Extend(Tmp);
00147
00148 Tmp = Origin + Radius;
00149 bbox->Extend(Tmp);
00150
00151 Tmp = EndOrigin - Radius;
00152 bbox->Extend(Tmp);
00153
00154 Tmp = EndOrigin + Radius;
00155 bbox->Extend(Tmp);
00156 }
00157
00158
00159 void G4FConicalSurface::PrintOn( std::ostream& os ) const
00160 {
00161
00162 os << "G4FConicalSurface with origin: " << origin << "\t"
00163 << "and axis: " << Position.GetAxis() << "\n"
00164 << "\t small radius: " << small_radius
00165 << "\t large radius: " << large_radius
00166 << "\t and length: " << length << "\n";
00167 }
00168
00169
00170 G4int G4FConicalSurface::operator==( const G4FConicalSurface& c ) const
00171 {
00172 return ( origin == c.origin &&
00173 Position.GetAxis() == c.Position.GetAxis() &&
00174 small_radius == c.small_radius &&
00175 large_radius == c.large_radius &&
00176 length == c.length &&
00177 tan_angle == c.tan_angle );
00178 }
00179
00180
00181 G4int G4FConicalSurface::WithinBoundary( const G4Vector3D& x ) const
00182 {
00183
00184
00185 G4Vector3D q = G4Vector3D( x - origin );
00186
00187 G4double qmag = q.mag();
00188 G4double ss = std::sin( std::atan2(large_radius-small_radius, length) );
00189 G4double ls = small_radius / ss;
00190 G4double ll = large_radius / ss;
00191
00192 if ( ( qmag >= ls ) && ( qmag <= ll ) )
00193 return 1;
00194 else
00195 return 0;
00196 }
00197
00198
00199 G4double G4FConicalSurface::Scale() const
00200 {
00201
00202
00203
00204 if ( small_radius == 0.0 )
00205 return large_radius;
00206 else
00207 return small_radius;
00208 }
00209
00210
00211 G4double G4FConicalSurface::Area() const
00212 {
00213
00214 G4double rdif = large_radius - small_radius;
00215
00216 return ( pi * ( small_radius + large_radius ) *
00217 std::sqrt( length * length + rdif * rdif ) );
00218 }
00219
00220
00221 void G4FConicalSurface::resize( G4double l, G4double srad, G4double lr )
00222 {
00223
00224
00225
00226
00227
00228
00229 if ( l >= 0.0 )
00230 length = l;
00231 else
00232 {
00233 std::ostringstream message;
00234 message << "Negative length." << G4endl
00235 << "Original value of " << length << " is retained.";
00236 G4Exception("G4FConicalSurface::resize()",
00237 "GeomSolids1001", JustWarning, message);
00238 }
00239
00240
00241 if ( srad >= 0.0 )
00242 small_radius = srad;
00243 else
00244 {
00245 std::ostringstream message;
00246 message << "Negative small radius." << G4endl
00247 << "Original value of " << small_radius << " is retained.";
00248 G4Exception("G4FConicalSurface::resize()",
00249 "GeomSolids1001", JustWarning, message);
00250 }
00251
00252
00253 if ( lr > small_radius )
00254 large_radius = lr;
00255 else
00256 {
00257 G4double r = small_radius + 1.0;
00258 lr = ( large_radius <= small_radius ) ? r : large_radius;
00259 large_radius = lr;
00260
00261 std::ostringstream message;
00262 message << "Large radius must exceed small radius." << G4endl
00263 << "Default value of " << large_radius << " is used.";
00264 G4Exception("G4FConicalSurface::resize()",
00265 "GeomSolids1001", JustWarning, message);
00266 }
00267
00268
00269 tan_angle = ( large_radius - small_radius ) / length ;
00270
00271 }
00272
00273
00274 G4int G4FConicalSurface::Intersect(const G4Ray& ry )
00275 {
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287 distance = kInfinity;
00288 closest_hit = PINFINITY;
00289
00290
00291 G4Point3D x = ry.GetStart();
00292 G4Vector3D dhat = ry.GetDir();
00293
00294
00295 G4double ta = tan_angle;
00296 G4Vector3D ahat = Position.GetAxis();
00297
00298
00299 G4double sol[2];
00300 sol[0]=-1.0;
00301 sol[1]=-1.0;
00302
00303
00304 G4Vector3D gamma = G4Vector3D( x - Position.GetLocation() );
00305
00306 G4double t = 1 + ta * ta;
00307 G4double ga = gamma * ahat;
00308 G4double da = dhat * ahat;
00309
00310 G4double A = t * da * da - dhat * dhat;
00311 G4double B = 2 * ( -gamma * dhat + t * ga * da - large_radius * ta * da);
00312 G4double C = ( -gamma * gamma + t * ga * ga
00313 - 2 * large_radius * ta * ga
00314 + large_radius * large_radius );
00315
00316 G4double radical = B * B - 4.0 * A * C;
00317
00318 if ( radical < 0.0 )
00319
00320 return 0;
00321 else
00322 {
00323 G4double root = std::sqrt( radical );
00324 sol[0] = ( - B + root ) / ( 2. * A );
00325 sol[1] = ( - B - root ) / ( 2. * A );
00326 }
00327
00328
00329
00330 G4Point3D p0 = G4Point3D( x + sol[0]*dhat );
00331 G4Point3D p1 = G4Point3D( x + sol[1]*dhat );
00332
00333 if( !GetBBox()->Inside(p0) )
00334 sol[0] = kInfinity;
00335
00336 if( !GetBBox()->Inside(p1) )
00337 sol[1] = kInfinity;
00338
00339
00340
00341 G4int nbinter = 0;
00342 distance = kInfinity;
00343
00344 for ( G4int i = 0; i < 2; i++ )
00345 {
00346 if(sol[i] < kInfinity) {
00347 if ( (sol[i] > kCarTolerance*0.5) ) {
00348 nbinter++;
00349 if ( distance > (sol[i]*sol[i]) ) {
00350 distance = sol[i]*sol[i];
00351 }
00352 }
00353 }
00354 }
00355
00356 return nbinter;
00357 }
00358
00359
00360 G4double G4FConicalSurface::HowNear( const G4Vector3D& x ) const
00361 {
00362
00363
00364
00365
00366
00367
00368 G4double hownear ;
00369
00370 G4Vector3D upcorner = G4Vector3D ( small_radius, 0 , origin.z()+Position.GetAxis().z()*length);
00371 G4Vector3D downcorner = G4Vector3D ( large_radius, 0 , origin.z());
00372 G4Vector3D xd;
00373
00374 xd = G4Vector3D ( std::sqrt ( x.x()*x.x() + x.y()*x.y() ) , 0 , x.z() );
00375
00376 G4double r = (upcorner.z() - downcorner.z()) / (upcorner.x() - downcorner.x());
00377 G4double q = (downcorner.z()*upcorner.x() - upcorner.z()*downcorner.x()) /
00378 (upcorner.x() - downcorner.x());
00379
00380 G4double Zinter = (xd.z()*r*r + xd.x()*r +q)/(1+r*r) ;
00381
00382 if ( ((Zinter >= downcorner.z()) && (Zinter <=upcorner.z())) ||
00383 ((Zinter >= upcorner.z()) && (Zinter <=downcorner.z())) ) {
00384 hownear = std::fabs(r*xd.x()-xd.z()+q)/std::sqrt(1+r*r);
00385 return hownear;
00386 } else {
00387 hownear = std::min ( (xd-upcorner).mag() , (xd-downcorner).mag() );
00388 return hownear;
00389 }
00390 }
00391
00392
00393 G4Vector3D G4FConicalSurface::SurfaceNormal( const G4Point3D& p ) const
00394 {
00395
00396
00397 G4Vector3D ss = G4Vector3D( p - origin );
00398 G4double da = ss * Position.GetAxis();
00399 G4double r = std::sqrt( ss*ss - da*da);
00400 G4double z = tan_angle * r;
00401
00402 if (Position.GetAxis().z() < 0)
00403 z = -z;
00404
00405 G4Vector3D n(p.x(), p.y(), z);
00406 n = n.unit();
00407
00408 if( !sameSense )
00409 n = -n;
00410
00411 return n;
00412 }
00413
00414 G4int G4FConicalSurface::Inside ( const G4Vector3D& x ) const
00415 {
00416
00417 if ( HowNear( x ) >= -0.5*kCarTolerance )
00418 return 1;
00419 else
00420 return 0;
00421 }
00422