Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // Implementation of G4IntersectingCone, a uti << 27 // the intersection of an arbitrary line with << 28 // 26 // 29 // Author: David C. Williams (davidw@scipp.ucs << 27 // $Id: G4IntersectingCone.cc,v 1.13 2010-07-12 15:25:37 gcosmo Exp $ >> 28 // GEANT4 tag $Name: geant4-09-04-patch-01 $ >> 29 // >> 30 // >> 31 // -------------------------------------------------------------------- >> 32 // GEANT 4 class source file >> 33 // >> 34 // >> 35 // G4IntersectingCone.cc >> 36 // >> 37 // Implementation of a utility class which calculates the intersection >> 38 // of an arbitrary line with a fixed cone 30 // ------------------------------------------- 39 // -------------------------------------------------------------------- 31 40 32 #include "G4IntersectingCone.hh" 41 #include "G4IntersectingCone.hh" 33 #include "G4GeometryTolerance.hh" 42 #include "G4GeometryTolerance.hh" 34 43 >> 44 // 35 // Constructor 45 // Constructor 36 // 46 // 37 G4IntersectingCone::G4IntersectingCone( const 47 G4IntersectingCone::G4IntersectingCone( const G4double r[2], 38 const 48 const G4double z[2] ) 39 { << 49 { 40 const G4double halfCarTolerance << 50 static const G4double halfCarTolerance 41 = 0.5 * G4GeometryTolerance::GetInstance() 51 = 0.5 * G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 42 52 >> 53 // 43 // What type of cone are we? 54 // What type of cone are we? 44 // 55 // 45 type1 = (std::abs(z[1]-z[0]) > std::abs(r[1] << 56 type1 = (std::fabs(z[1]-z[0]) > std::fabs(r[1]-r[0])); 46 << 57 47 if (type1) // tube like << 58 if (type1) 48 { 59 { 49 B = (r[1] - r[0]) / (z[1] - z[0]); << 60 B = (r[1]-r[0])/(z[1]-z[0]); // tube like 50 A = (r[0]*z[1] - r[1]*z[0]) / (z[1] -z[0]) << 61 A = 0.5*( r[1]+r[0] - B*(z[1]+z[0]) ); 51 } 62 } 52 else // disk like << 63 else 53 { 64 { 54 B = (z[1] - z[0]) / (r[1] - r[0]); << 65 B = (z[1]-z[0])/(r[1]-r[0]); // disk like 55 A = (z[0]*r[1] - z[1]*r[0]) / (r[1] - r[0] << 66 A = 0.5*( z[1]+z[0] - B*(r[1]+r[0]) ); 56 } 67 } 57 << 68 // 58 // Calculate extent 69 // Calculate extent 59 // 70 // 60 rLo = std::min(r[0], r[1]) - halfCarToleranc << 71 if (r[0] < r[1]) 61 rHi = std::max(r[0], r[1]) + halfCarToleranc << 72 { 62 zLo = std::min(z[0], z[1]) - halfCarToleranc << 73 rLo = r[0]-halfCarTolerance; rHi = r[1]+halfCarTolerance; 63 zHi = std::max(z[0], z[1]) + halfCarToleranc << 74 } >> 75 else >> 76 { >> 77 rLo = r[1]-halfCarTolerance; rHi = r[0]+halfCarTolerance; >> 78 } >> 79 >> 80 if (z[0] < z[1]) >> 81 { >> 82 zLo = z[0]-halfCarTolerance; zHi = z[1]+halfCarTolerance; >> 83 } >> 84 else >> 85 { >> 86 zLo = z[1]-halfCarTolerance; zHi = z[0]+halfCarTolerance; >> 87 } 64 } 88 } 65 89 >> 90 >> 91 // 66 // Fake default constructor - sets only member 92 // Fake default constructor - sets only member data and allocates memory 67 // for usage restri 93 // for usage restricted to object persistency. 68 // 94 // 69 G4IntersectingCone::G4IntersectingCone( __void 95 G4IntersectingCone::G4IntersectingCone( __void__& ) 70 : zLo(0.), zHi(0.), rLo(0.), rHi(0.), A(0.), << 96 : zLo(0.), zHi(0.), rLo(0.), rHi(0.), type1(false), A(0.), B(0.) 71 { 97 { 72 } 98 } 73 99 >> 100 >> 101 // 74 // Destructor 102 // Destructor 75 // 103 // 76 G4IntersectingCone::~G4IntersectingCone() = de << 104 G4IntersectingCone::~G4IntersectingCone() >> 105 { >> 106 } >> 107 77 108 >> 109 // 78 // HitOn 110 // HitOn 79 // 111 // 80 // Check r or z extent, as appropriate, to see 112 // Check r or z extent, as appropriate, to see if the point is possibly 81 // on the cone. 113 // on the cone. 82 // 114 // 83 G4bool G4IntersectingCone::HitOn( const G4doub 115 G4bool G4IntersectingCone::HitOn( const G4double r, 84 const G4doub 116 const G4double z ) 85 { 117 { 86 // 118 // 87 // Be careful! The inequalities cannot be "< 119 // Be careful! The inequalities cannot be "<=" and ">=" here without 88 // punching a tiny hole in our shape! 120 // punching a tiny hole in our shape! 89 // 121 // 90 if (type1) 122 if (type1) 91 { 123 { 92 if (z < zLo || z > zHi) return false; 124 if (z < zLo || z > zHi) return false; 93 } 125 } 94 else 126 else 95 { 127 { 96 if (r < rLo || r > rHi) return false; 128 if (r < rLo || r > rHi) return false; 97 } 129 } 98 130 99 return true; 131 return true; 100 } 132 } 101 133 >> 134 >> 135 // 102 // LineHitsCone 136 // LineHitsCone 103 // 137 // 104 // Calculate the intersection of a line with o 138 // Calculate the intersection of a line with our conical surface, ignoring 105 // any phi division 139 // any phi division 106 // 140 // 107 G4int G4IntersectingCone::LineHitsCone( const << 141 G4int G4IntersectingCone::LineHitsCone( const G4ThreeVector &p, 108 const << 142 const G4ThreeVector &v, 109 << 143 G4double *s1, G4double *s2 ) 110 { 144 { 111 if (type1) 145 if (type1) 112 { 146 { 113 return LineHitsCone1( p, v, s1, s2 ); 147 return LineHitsCone1( p, v, s1, s2 ); 114 } 148 } 115 else 149 else 116 { 150 { 117 return LineHitsCone2( p, v, s1, s2 ); 151 return LineHitsCone2( p, v, s1, s2 ); 118 } 152 } 119 } 153 } 120 154 >> 155 >> 156 // 121 // LineHitsCone1 157 // LineHitsCone1 122 // 158 // 123 // Calculate the intersections of a line with 159 // Calculate the intersections of a line with a conical surface. Only 124 // suitable if zPlane[0] != zPlane[1]. 160 // suitable if zPlane[0] != zPlane[1]. 125 // 161 // 126 // Equation of a line: 162 // Equation of a line: 127 // 163 // 128 // x = x0 + s*tx y = y0 + s*ty 164 // x = x0 + s*tx y = y0 + s*ty z = z0 + s*tz 129 // 165 // 130 // Equation of a conical surface: 166 // Equation of a conical surface: 131 // 167 // 132 // x**2 + y**2 = (A + B*z)**2 168 // x**2 + y**2 = (A + B*z)**2 133 // 169 // 134 // Solution is quadratic: 170 // Solution is quadratic: 135 // 171 // 136 // a*s**2 + b*s + c = 0 172 // a*s**2 + b*s + c = 0 137 // 173 // 138 // where: 174 // where: 139 // 175 // 140 // a = tx**2 + ty**2 - (B*tz)**2 << 176 // a = x0**2 + y0**2 - (A + B*z0)**2 141 // 177 // 142 // b = 2*( px*vx + py*vy - B*(A + B*pz)*vz ) << 178 // b = 2*( x0*tx + y0*ty - (A*B - B*B*z0)*tz) 143 // 179 // 144 // c = x0**2 + y0**2 - (A + B*z0)**2 << 180 // c = tx**2 + ty**2 - (B*tz)**2 145 // 181 // 146 // Notice, that if a < 0, this indicates that 182 // Notice, that if a < 0, this indicates that the two solutions (assuming 147 // they exist) are in opposite cones (that is, 183 // they exist) are in opposite cones (that is, given z0 = -A/B, one z < z0 148 // and the other z > z0). For our shapes, the 184 // and the other z > z0). For our shapes, the invalid solution is one 149 // which produces A + Bz < 0, or the one where 185 // which produces A + Bz < 0, or the one where Bz is smallest (most negative). 150 // Since Bz = B*s*tz, if B*tz > 0, we want the 186 // Since Bz = B*s*tz, if B*tz > 0, we want the largest s, otherwise, 151 // the smaller. 187 // the smaller. 152 // 188 // 153 // If there are two solutions on one side of t 189 // If there are two solutions on one side of the cone, we want to make 154 // sure that they are on the "correct" side, t 190 // sure that they are on the "correct" side, that is A + B*z0 + s*B*tz >= 0. 155 // 191 // 156 // If a = 0, we have a linear problem: s = c/b 192 // If a = 0, we have a linear problem: s = c/b, which again gives one solution. 157 // This should be rare. 193 // This should be rare. 158 // 194 // 159 // For b*b - 4*a*c = 0, we also have one solut 195 // For b*b - 4*a*c = 0, we also have one solution, which is almost always 160 // a line just grazing the surface of a the co << 196 // a line just grazing the surface of a the cone, which we want to ignore. 161 // However, there are two other, very rare, po 197 // However, there are two other, very rare, possibilities: 162 // a line intersecting the z axis and either: 198 // a line intersecting the z axis and either: 163 // 1. At the same angle std::atan(B) to 199 // 1. At the same angle std::atan(B) to just miss one side of the cone, or 164 // 2. Intersecting the cone apex (0,0,-A 200 // 2. Intersecting the cone apex (0,0,-A/B) 165 // We *don't* want to miss these! How do we id 201 // We *don't* want to miss these! How do we identify them? Well, since 166 // this case is rare, we can at least swallow 202 // this case is rare, we can at least swallow a little more CPU than we would 167 // normally be comfortable with. Intersection 203 // normally be comfortable with. Intersection with the z axis means 168 // x0*ty - y0*tx = 0. Case (1) means a==0, and 204 // x0*ty - y0*tx = 0. Case (1) means a==0, and we've already dealt with that 169 // above. Case (2) means a < 0. 205 // above. Case (2) means a < 0. 170 // 206 // 171 // Now: x0*tx + y0*ty = 0 in terms of roundoff 207 // Now: x0*tx + y0*ty = 0 in terms of roundoff error. We can write: 172 // Delta = x0*tx + y0*ty 208 // Delta = x0*tx + y0*ty 173 // b = 2*( Delta - B*(A + B*z0)*tz << 209 // b = 2*( Delta - (A*B + B*B*z0)*tz ) 174 // For: 210 // For: 175 // b*b - 4*a*c = epsilon 211 // b*b - 4*a*c = epsilon 176 // where epsilon is small, then: 212 // where epsilon is small, then: 177 // Delta = epsilon/2/B 213 // Delta = epsilon/2/B 178 // << 214 // 179 G4int G4IntersectingCone::LineHitsCone1( const << 215 G4int G4IntersectingCone::LineHitsCone1( const G4ThreeVector &p, 180 const << 216 const G4ThreeVector &v, 181 << 217 G4double *s1, G4double *s2 ) 182 { 218 { 183 static const G4double EPS = DBL_EPSILON; // << 184 // << 185 G4double x0 = p.x(), y0 = p.y(), z0 = p.z(); 219 G4double x0 = p.x(), y0 = p.y(), z0 = p.z(); 186 G4double tx = v.x(), ty = v.y(), tz = v.z(); 220 G4double tx = v.x(), ty = v.y(), tz = v.z(); 187 221 188 // Value of radical can be inaccurate due to << 222 G4double a = tx*tx + ty*ty - sqr(B*tz); 189 // if to calculate the coefficiets a,b,c lik << 223 G4double b = 2*( x0*tx + y0*ty - (A*B + B*B*z0)*tz); 190 // G4double a = tx*tx + ty*ty - sqr(B*tz << 224 G4double c = x0*x0 + y0*y0 - sqr(A + B*z0); 191 // G4double b = 2*( x0*tx + y0*ty - B*(A << 225 192 // G4double c = x0*x0 + y0*y0 - sqr(A + << 226 G4double radical = b*b - 4*a*c; 193 // << 227 194 // For more accurate calculation of radical << 228 if (radical < -1E-6*std::fabs(b)) { return 0; } // No solution 195 // are splitted in two components, radial an << 229 196 // << 230 if (radical < 1E-6*std::fabs(b)) 197 G4double ar = tx*tx + ty*ty; << 198 G4double az = sqr(B*tz); << 199 G4double br = 2*(x0*tx + y0*ty); << 200 G4double bz = 2*B*(A + B*z0)*tz; << 201 G4double cr = x0*x0 + y0*y0; << 202 G4double cz = sqr(A + B*z0); << 203 << 204 // Instead radical = b*b - 4*a*c << 205 G4double arcz = 4*ar*cz; << 206 G4double azcr = 4*az*cr; << 207 G4double radical = (br*br - 4*ar*cr) + ((std << 208 << 209 // Find the coefficients << 210 G4double a = ar - az; << 211 G4double b = br - bz; << 212 G4double c = cr - cz; << 213 << 214 if (radical < -EPS*std::fabs(b)) { return 0 << 215 << 216 if (radical < EPS*std::fabs(b)) << 217 { 231 { 218 // 232 // 219 // The radical is roughly zero: check for 233 // The radical is roughly zero: check for special, very rare, cases 220 // 234 // 221 if (std::fabs(a) > 1/kInfinity) 235 if (std::fabs(a) > 1/kInfinity) 222 { 236 { 223 if(B==0.) { return 0; } 237 if(B==0.) { return 0; } 224 if ( std::fabs(x0*ty - y0*tx) < std::fab << 238 if ( std::fabs(x0*ty - y0*tx) < std::fabs(1E-6/B) ) 225 { 239 { 226 *s1 = -0.5*b/a; 240 *s1 = -0.5*b/a; 227 return 1; 241 return 1; 228 } 242 } 229 return 0; 243 return 0; 230 } 244 } 231 } 245 } 232 else 246 else 233 { 247 { 234 radical = std::sqrt(radical); 248 radical = std::sqrt(radical); 235 } 249 } 236 << 250 237 if (a > 1/kInfinity) 251 if (a > 1/kInfinity) 238 { 252 { 239 G4double sa, sb, q = -0.5*( b + (b < 0 ? - 253 G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) ); 240 sa = q/a; 254 sa = q/a; 241 sb = c/q; 255 sb = c/q; 242 if (sa < sb) { *s1 = sa; *s2 = sb; } else 256 if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; } 243 if (A + B*(z0+(*s1)*tz) < 0) { return 0; 257 if (A + B*(z0+(*s1)*tz) < 0) { return 0; } 244 return 2; 258 return 2; 245 } 259 } 246 else if (a < -1/kInfinity) 260 else if (a < -1/kInfinity) 247 { 261 { 248 G4double sa, sb, q = -0.5*( b + (b < 0 ? - 262 G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) ); 249 sa = q/a; 263 sa = q/a; 250 sb = c/q; 264 sb = c/q; 251 *s1 = (B*tz > 0)^(sa > sb) ? sb : sa; 265 *s1 = (B*tz > 0)^(sa > sb) ? sb : sa; 252 return 1; 266 return 1; 253 } 267 } 254 else if (std::fabs(b) < 1/kInfinity) 268 else if (std::fabs(b) < 1/kInfinity) 255 { 269 { 256 return 0; 270 return 0; 257 } 271 } 258 else 272 else 259 { 273 { 260 *s1 = -c/b; 274 *s1 = -c/b; 261 if (A + B*(z0+(*s1)*tz) < 0) { return 0; 275 if (A + B*(z0+(*s1)*tz) < 0) { return 0; } 262 return 1; 276 return 1; 263 } 277 } 264 } 278 } 265 279 >> 280 >> 281 // 266 // LineHitsCone2 282 // LineHitsCone2 267 // 283 // 268 // See comments under LineHitsCone1. In this r 284 // See comments under LineHitsCone1. In this routine, case2, we have: 269 // 285 // 270 // Z = A + B*R 286 // Z = A + B*R 271 // 287 // 272 // The solution is still quadratic: 288 // The solution is still quadratic: 273 // 289 // 274 // a = tz**2 - B*B*(tx**2 + ty**2) 290 // a = tz**2 - B*B*(tx**2 + ty**2) 275 // 291 // 276 // b = 2*( (z0-A)*tz - B*B*(x0*tx+y0*ty) ) 292 // b = 2*( (z0-A)*tz - B*B*(x0*tx+y0*ty) ) 277 // 293 // 278 // c = ( (z0-A)**2 - B*B*(x0**2 + y0**2) ) 294 // c = ( (z0-A)**2 - B*B*(x0**2 + y0**2) ) 279 // 295 // 280 // The rest is much the same, except some deta 296 // The rest is much the same, except some details. 281 // 297 // 282 // a > 0 now means we intersect only once in t 298 // a > 0 now means we intersect only once in the correct hemisphere. 283 // 299 // 284 // a > 0 ? We only want solution which produce << 300 // a > 0 ? We only want solution which produces R > 0. 285 // since R = (z0+s*tz-A)/B, for tz/B > 0, this 301 // since R = (z0+s*tz-A)/B, for tz/B > 0, this is the largest s 286 // for tz/B < 0, this is the smallest 302 // for tz/B < 0, this is the smallest s 287 // thus, same as in case 1 ( since sign(tz/B) 303 // thus, same as in case 1 ( since sign(tz/B) = sign(tz*B) ) 288 // 304 // 289 G4int G4IntersectingCone::LineHitsCone2( const << 305 G4int G4IntersectingCone::LineHitsCone2( const G4ThreeVector &p, 290 const << 306 const G4ThreeVector &v, 291 << 307 G4double *s1, G4double *s2 ) 292 { 308 { 293 static const G4double EPS = DBL_EPSILON; // << 294 // << 295 G4double x0 = p.x(), y0 = p.y(), z0 = p.z(); 309 G4double x0 = p.x(), y0 = p.y(), z0 = p.z(); 296 G4double tx = v.x(), ty = v.y(), tz = v.z(); 310 G4double tx = v.x(), ty = v.y(), tz = v.z(); 297 << 311 >> 312 298 // Special case which might not be so rare: 313 // Special case which might not be so rare: B = 0 (precisely) 299 // 314 // 300 if (B==0) 315 if (B==0) 301 { 316 { 302 if (std::fabs(tz) < 1/kInfinity) { return 317 if (std::fabs(tz) < 1/kInfinity) { return 0; } 303 << 318 304 *s1 = (A-z0)/tz; 319 *s1 = (A-z0)/tz; 305 return 1; 320 return 1; 306 } 321 } 307 322 308 // Value of radical can be inaccurate due to << 309 // if to calculate the coefficiets a,b,c lik << 310 // G4double a = tz*tz - B2*(tx*tx + ty*ty) << 311 // G4double b = 2*( (z0-A)*tz - B2*(x0*tx << 312 // G4double c = sqr(z0-A) - B2*( x0*x0 + y << 313 // << 314 // For more accurate calculation of radical << 315 // are splitted in two components, radial an << 316 // << 317 G4double B2 = B*B; 323 G4double B2 = B*B; 318 324 319 G4double az = tz*tz; << 325 G4double a = tz*tz - B2*(tx*tx + ty*ty); 320 G4double ar = B2*(tx*tx + ty*ty); << 326 G4double b = 2*( (z0-A)*tz - B2*(x0*tx + y0*ty) ); 321 G4double bz = 2*(z0-A)*tz; << 327 G4double c = sqr(z0-A) - B2*( x0*x0 + y0*y0 ); 322 G4double br = 2*B2*(x0*tx + y0*ty); << 328 323 G4double cz = sqr(z0-A); << 329 G4double radical = b*b - 4*a*c; 324 G4double cr = B2*(x0*x0 + y0*y0); << 330 325 << 331 if (radical < -1E-6*std::fabs(b)) { return 0; } // No solution 326 // Instead radical = b*b - 4*a*c << 332 327 G4double arcz = 4*ar*cz; << 333 if (radical < 1E-6*std::fabs(b)) 328 G4double azcr = 4*az*cr; << 329 G4double radical = (br*br - 4*ar*cr) + ((std << 330 << 331 // Find the coefficients << 332 G4double a = az - ar; << 333 G4double b = bz - br; << 334 G4double c = cz - cr; << 335 << 336 if (radical < -EPS*std::fabs(b)) { return 0; << 337 << 338 if (radical < EPS*std::fabs(b)) << 339 { 334 { 340 // 335 // 341 // The radical is roughly zero: check for 336 // The radical is roughly zero: check for special, very rare, cases 342 // 337 // 343 if (std::fabs(a) > 1/kInfinity) 338 if (std::fabs(a) > 1/kInfinity) 344 { 339 { 345 if ( std::fabs(x0*ty - y0*tx) < std::fab << 340 if ( std::fabs(x0*ty - y0*tx) < std::fabs(1E-6/B) ) 346 { 341 { 347 *s1 = -0.5*b/a; 342 *s1 = -0.5*b/a; 348 return 1; 343 return 1; 349 } 344 } 350 return 0; 345 return 0; 351 } 346 } 352 } 347 } 353 else 348 else 354 { 349 { 355 radical = std::sqrt(radical); 350 radical = std::sqrt(radical); 356 } 351 } 357 << 352 358 if (a < -1/kInfinity) 353 if (a < -1/kInfinity) 359 { 354 { 360 G4double sa, sb, q = -0.5*( b + (b < 0 ? - 355 G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) ); 361 sa = q/a; 356 sa = q/a; 362 sb = c/q; 357 sb = c/q; 363 if (sa < sb) { *s1 = sa; *s2 = sb; } else 358 if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; } 364 if ((z0 + (*s1)*tz - A)/B < 0) { return 359 if ((z0 + (*s1)*tz - A)/B < 0) { return 0; } 365 return 2; 360 return 2; 366 } 361 } 367 else if (a > 1/kInfinity) 362 else if (a > 1/kInfinity) 368 { 363 { 369 G4double sa, sb, q = -0.5*( b + (b < 0 ? - 364 G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) ); 370 sa = q/a; 365 sa = q/a; 371 sb = c/q; 366 sb = c/q; 372 *s1 = (tz*B > 0)^(sa > sb) ? sb : sa; 367 *s1 = (tz*B > 0)^(sa > sb) ? sb : sa; 373 return 1; 368 return 1; 374 } 369 } 375 else if (std::fabs(b) < 1/kInfinity) 370 else if (std::fabs(b) < 1/kInfinity) 376 { 371 { 377 return 0; 372 return 0; 378 } 373 } 379 else 374 else 380 { 375 { 381 *s1 = -c/b; 376 *s1 = -c/b; 382 if ((z0 + (*s1)*tz - A)/B < 0) { return 377 if ((z0 + (*s1)*tz - A)/B < 0) { return 0; } 383 return 1; 378 return 1; 384 } 379 } 385 } 380 } 386 381