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 "G4Ellipse.hh"
00037
00038 #include "G4PhysicalConstants.hh"
00039 #include "G4SystemOfUnits.hh"
00040 #include "G4GeometryTolerance.hh"
00041
00042 G4Ellipse::G4Ellipse()
00043 : semiAxis1(0.), semiAxis2(0.), ratioAxis2Axis1(0.), forTangent(.0)
00044 {
00045 }
00046
00047 G4Ellipse::~G4Ellipse()
00048 {
00049 }
00050
00051 G4Ellipse::G4Ellipse(const G4Ellipse& right)
00052 : G4Conic(), semiAxis1(right.semiAxis1), semiAxis2(right.semiAxis2),
00053 ratioAxis2Axis1(right.ratioAxis2Axis1), toUnitCircle(right.toUnitCircle),
00054 forTangent(right.forTangent)
00055 {
00056 pShift = right.pShift;
00057 position = right.position;
00058 bBox = right.bBox;
00059 start = right.start;
00060 end = right.end;
00061 pStart = right.pStart;
00062 pEnd = right.pEnd;
00063 pRange = right.pRange;
00064 bounded = right.bounded;
00065 sameSense = right.sameSense;
00066 }
00067
00068 G4Ellipse& G4Ellipse::operator=(const G4Ellipse& right)
00069 {
00070 if (&right == this) return *this;
00071
00072 semiAxis1 = right.semiAxis1;
00073 semiAxis2 = right.semiAxis2;
00074 ratioAxis2Axis1 = right.ratioAxis2Axis1;
00075 toUnitCircle = right.toUnitCircle;
00076 forTangent = right.forTangent;
00077 pShift = right.pShift;
00078 position = right.position;
00079 bBox = right.bBox;
00080 start = right.start;
00081 end = right.end;
00082 pStart = right.pStart;
00083 pEnd = right.pEnd;
00084 pRange = right.pRange;
00085 bounded = right.bounded;
00086 sameSense = right.sameSense;
00087
00088 return *this;
00089 }
00090
00091 G4Curve* G4Ellipse::Project(const G4Transform3D& tr)
00092 {
00093 G4Point3D newLocation = tr*position.GetLocation();
00094 newLocation.setZ(0);
00095 G4double axisZ = ( tr*position.GetPZ() ).unit().z();
00096
00097 if (std::abs(axisZ)<G4GeometryTolerance::GetInstance()->GetAngularTolerance())
00098 { return 0; }
00099
00100 G4Vector3D newAxis(0, 0, axisZ>0? +1: -1);
00101
00102
00103
00104 G4Vector3D xPrime= tr*position.GetPX();
00105 xPrime.setZ(0);
00106 G4Vector3D yPrime= tr*position.GetPY();
00107 yPrime.setZ(0);
00108
00109 G4Vector3D a = G4Vector3D( semiAxis1*xPrime );
00110 G4Vector3D b = G4Vector3D( semiAxis2*yPrime );
00111
00112 G4double u;
00113 G4double abmag = a.mag2()-b.mag2();
00114 G4double prod = 2*a*b;
00115
00116 if ((abmag > FLT_MAX) && (prod < -FLT_MAX))
00117 u = -pi/8;
00118 else if ((abmag < -FLT_MAX) && (prod > FLT_MAX))
00119 u = 3*pi/8;
00120 else if ((abmag < -FLT_MAX) && (prod < -FLT_MAX))
00121 u = -3*pi/8;
00122 else if ((std::abs(abmag) < perMillion) && (std::abs(prod) < perMillion))
00123 u = 0.;
00124 else
00125 u = std::atan2(prod,abmag) / 2;
00126
00127
00128 G4Vector3D sAxis1 = G4Vector3D( a*std::cos(u)+b*std::sin(u) );
00129 G4Vector3D sAxis2 = G4Vector3D( a*std::cos(u+pi/2)+b*std::sin(u+pi/2) );
00130 G4double newSemiAxis1 = sAxis1.mag();
00131 G4double newSemiAxis2 = sAxis2.mag();
00132 G4Vector3D newRefDirection = sAxis1;
00133
00134
00135 G4Axis2Placement3D newPosition;
00136 newPosition.Init(newRefDirection, newAxis, newLocation);
00137 G4Ellipse* r= new G4Ellipse;
00138 r->Init(newPosition, newSemiAxis1, newSemiAxis2);
00139
00140
00141
00142 r->SetPShift(u);
00143
00144
00145 if (IsBounded())
00146 r->SetBounds(GetPStart(), GetPEnd());
00147
00148
00149
00150 r->SetSameSense(GetSameSense());
00151
00152 return r;
00153 }
00154
00155 void G4Ellipse::InitBounded()
00156 {
00157
00158
00159
00160
00161
00162
00163
00164
00165 bBox.Init(GetStart(), GetEnd());
00166
00167
00168
00169 for (G4int i=0; i<3; i++)
00170 {
00171 G4double u= std::atan2(position.GetPY()(i)*semiAxis2,
00172 position.GetPX()(i)*semiAxis1);
00173 if (IsPOn(u))
00174 bBox.Extend(GetPoint(u));
00175
00176 if (IsPOn(u+pi))
00177 bBox.Extend(GetPoint(u+pi));
00178 }
00179 }
00180
00181
00182 G4bool G4Ellipse::Tangent(G4CurvePoint& cp, G4Vector3D& v)
00183 {
00184
00185
00186
00187
00188
00189 const G4Axis2Placement3D& pos = *(GetPosition());
00190 G4Point3D p= pos.GetToPlacementCoordinates() * cp.GetPoint();
00191
00192 v=forTangent*p.y()*pos.GetPX() + p.x()*pos.GetPY();
00193 if(GetSameSense())
00194 v = -v;
00195
00196 return true;
00197 }
00198