00001 // 00002 // ******************************************************************** 00003 // * License and Disclaimer * 00004 // * * 00005 // * The Geant4 software is copyright of the Copyright Holders of * 00006 // * the Geant4 Collaboration. It is provided under the terms and * 00007 // * conditions of the Geant4 Software License, included in the file * 00008 // * LICENSE and available at http://cern.ch/geant4/license . These * 00009 // * include a list of copyright holders. * 00010 // * * 00011 // * Neither the authors of this software system, nor their employing * 00012 // * institutes,nor the agencies providing financial support for this * 00013 // * work make any representation or warranty, express or implied, * 00014 // * regarding this software system or assume any liability for its * 00015 // * use. Please see the license in the file LICENSE and URL above * 00016 // * for the full disclaimer and the limitation of liability. * 00017 // * * 00018 // * This code implementation is the result of the scientific and * 00019 // * technical work of the GEANT4 collaboration. * 00020 // * By using, copying, modifying or distributing the software (or * 00021 // * any work based on the software) you agree to acknowledge its * 00022 // * use in resulting scientific publications, and indicate your * 00023 // * acceptance of all terms of the Geant4 Software license. * 00024 // ******************************************************************** 00025 // 00026 // 00027 // $Id$ 00028 // 00029 // -------------------------------------------------------------------- 00030 // GEANT 4 inline definitions file 00031 // 00032 // G4Ellipse.icc 00033 // 00034 // Implementation of inline methods of G4Ellipse 00035 // -------------------------------------------------------------------- 00036 00037 inline 00038 void G4Ellipse::Init(const G4Axis2Placement3D& position0, 00039 G4double semiAxis10, G4double semiAxis20) 00040 { 00041 position= position0; 00042 semiAxis1= semiAxis10; 00043 semiAxis2= semiAxis20; 00044 00045 ratioAxis2Axis1= semiAxis2/semiAxis1; 00046 00047 SetBounds(0, 0); 00048 00049 // needed only for 2D ellipses 00050 toUnitCircle = G4Scale3D(1/semiAxis1, 1/semiAxis2, 0) 00051 * position.GetToPlacementCoordinates(); 00052 00053 forTangent= -semiAxis1*semiAxis1/(semiAxis2*semiAxis2); 00054 } 00055 00056 inline 00057 G4double G4Ellipse::GetSemiAxis1() const 00058 { 00059 return semiAxis1; 00060 } 00061 00062 inline 00063 G4double G4Ellipse::GetSemiAxis2() const 00064 { 00065 return semiAxis2; 00066 } 00067 00069 00070 inline 00071 G4double G4Ellipse::GetPMax() const 00072 { 00073 return CLHEP::twopi; 00074 } 00075 00076 inline 00077 G4Point3D G4Ellipse::GetPoint(G4double param) const 00078 { 00079 param-= GetPShift(); 00080 return G4Point3D( position.GetLocation() 00081 + semiAxis1*std::cos(param)*position.GetPX() 00082 + semiAxis2*std::sin(param)*position.GetPY() ); 00083 } 00084 00085 inline 00086 G4double G4Ellipse::GetPPoint(const G4Point3D& pt) const 00087 { 00088 G4Point3D ptLocal= position.GetToPlacementCoordinates()*pt; 00089 G4double angle= std::atan2(ptLocal.y(), ptLocal.x()*ratioAxis2Axis1); 00090 G4double r= (angle<0)? angle+CLHEP::twopi: angle; 00091 return r+GetPShift(); 00092 } 00093 00095 00096 #include "G4CurveRayIntersection.hh" 00097 00098 /* 00099 inline 00100 void G4Ellipse::IntersectRay2D(const G4Ray& ray, 00101 G4CurveRayIntersection& is) 00102 { 00103 is.Init(*this, ray); 00104 00105 // transform s.t. the ellipse becomes the unit circle 00106 // with the center at the origin 00107 00108 // 2D operations would be faster 00109 G4Point3D s= toUnitCircle*ray.GetStart(); 00110 G4Vector3D d= toUnitCircle*ray.GetDir(); 00111 00112 // solve (s+i*t)^2 = 1 for i (the distance) 00113 00114 G4double sd= s*d; 00115 G4double dd= d.mag2(); // never 0 00116 G4double ss= s.mag2(); 00117 00118 G4double discr= sd*sd-dd*(ss-1); 00119 if (discr >= 0) { 00120 00121 // 2 intersections (maybe 1, but this case is rare) 00122 G4double sqrtdiscr= std::sqrt(discr); 00123 // find the smallest positive i 00124 G4double i= -sd-sqrtdiscr; 00125 if (i<kCarTolerance) { 00126 i= -sd+sqrtdiscr; 00127 if (i<kCarTolerance) { 00128 return; 00129 } 00130 } 00131 i/= dd; 00132 G4CurveRayIntersection isTmp(*this, ray); 00133 isTmp.ResetDistance(i); 00134 is.Update(isTmp); 00135 00136 } 00137 } 00138 */ 00139 00140 inline 00141 G4int G4Ellipse::IntersectRay2D(const G4Ray& ray) 00142 { 00143 // transform s.t. the ellipse becomes the unit circle 00144 // with the center at the origin 00145 00146 // 2D operations would be faster 00147 G4Point3D st= toUnitCircle*ray.GetStart(); 00148 G4Vector3D d= toUnitCircle*ray.GetDir(); 00149 00150 // solve (st+i*t)^2 = 1 for i (the distance) 00151 G4double sd= st*d; 00152 G4double dd= d.mag2(); // never 0 00153 G4double ss= st.mag2(); 00154 00155 G4double discr= sd*sd-dd*(ss-1); 00156 G4int nbinter = 0; 00157 G4double i; 00158 00159 if (discr > 0) 00160 { 00161 // 2 intersections 00162 G4double sqrtdiscr= std::sqrt(discr); 00163 00164 // if i is positive, we have an intersection 00165 i= -sd-sqrtdiscr; 00166 if (i > kCarTolerance) 00167 nbinter++; 00168 00169 i= -sd+sqrtdiscr; 00170 if (i > kCarTolerance) 00171 nbinter++; 00172 } 00173 00174 // if the ray is tangent on the circle 00175 if (discr == 0) 00176 nbinter = 1; 00177 00178 return nbinter; 00179 }