source: trunk/source/geometry/solids/specific/src/G4TwistTubsSide.cc@ 1036

Last change on this file since 1036 was 850, checked in by garnier, 17 years ago

geant4.8.2 beta
File size: 34.9 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4TwistTubsSide.cc,v 1.5 2006/06/29 18:49:18 gunter Exp $
28// GEANT4 tag $Name: HEAD $
29//
30//
31// --------------------------------------------------------------------
32// GEANT 4 class source file
33//
34//
35// G4TwistTubsSide.cc
36//
37// Author:
38// 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
39//
40// History:
41// 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
42// from original version in Jupiter-2.5.02 application.
43// 29-Apr-2004 - O.Link. Bug fixed in GetAreaCode
44// --------------------------------------------------------------------
45
46#include "G4TwistTubsSide.hh"
47
48//=====================================================================
49//* constructors ------------------------------------------------------
50
51G4TwistTubsSide::G4TwistTubsSide(const G4String &name,
52 const G4RotationMatrix &rot,
53 const G4ThreeVector &tlate,
54 G4int handedness,
55 const G4double kappa,
56 const EAxis axis0,
57 const EAxis axis1,
58 G4double axis0min,
59 G4double axis1min,
60 G4double axis0max,
61 G4double axis1max)
62 : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
63 axis0min, axis1min, axis0max, axis1max),
64 fKappa(kappa)
65{
66 if (axis0 == kZAxis && axis1 == kXAxis) {
67 G4Exception("G4TwistTubsSide::G4TwistTubsSide()", "InvalidSetup",
68 FatalException, "Should swap axis0 and axis1!");
69 }
70 fIsValidNorm = false;
71 SetCorners();
72 SetBoundaries();
73}
74
75G4TwistTubsSide::G4TwistTubsSide(const G4String &name,
76 G4double EndInnerRadius[2],
77 G4double EndOuterRadius[2],
78 G4double DPhi,
79 G4double EndPhi[2],
80 G4double EndZ[2],
81 G4double InnerRadius,
82 G4double OuterRadius,
83 G4double Kappa,
84 G4int handedness)
85 : G4VTwistSurface(name)
86{
87 fHandedness = handedness; // +z = +ve, -z = -ve
88 fAxis[0] = kXAxis; // in local coordinate system
89 fAxis[1] = kZAxis;
90 fAxisMin[0] = InnerRadius; // Inner-hype radius at z=0
91 fAxisMax[0] = OuterRadius; // Outer-hype radius at z=0
92 fAxisMin[1] = EndZ[0];
93 fAxisMax[1] = EndZ[1];
94
95 fKappa = Kappa;
96 fRot.rotateZ( fHandedness > 0
97 ? -0.5*DPhi
98 : 0.5*DPhi );
99 fTrans.set(0, 0, 0);
100 fIsValidNorm = false;
101
102 SetCorners( EndInnerRadius, EndOuterRadius, EndPhi, EndZ) ;
103 SetBoundaries();
104}
105
106
107//=====================================================================
108//* Fake default constructor ------------------------------------------
109
110G4TwistTubsSide::G4TwistTubsSide( __void__& a )
111 : G4VTwistSurface(a)
112{
113}
114
115
116//=====================================================================
117//* destructor --------------------------------------------------------
118
119G4TwistTubsSide::~G4TwistTubsSide()
120{
121}
122
123//=====================================================================
124//* GetNormal ---------------------------------------------------------
125
126G4ThreeVector G4TwistTubsSide::GetNormal(const G4ThreeVector &tmpxx,
127 G4bool isGlobal)
128{
129 // GetNormal returns a normal vector at a surface (or very close
130 // to surface) point at tmpxx.
131 // If isGlobal=true, it returns the normal in global coordinate.
132 //
133 G4ThreeVector xx;
134 if (isGlobal) {
135 xx = ComputeLocalPoint(tmpxx);
136 if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) {
137 return ComputeGlobalDirection(fCurrentNormal.normal);
138 }
139 } else {
140 xx = tmpxx;
141 if (xx == fCurrentNormal.p) {
142 return fCurrentNormal.normal;
143 }
144 }
145
146 G4ThreeVector er(1, fKappa * xx.z(), 0);
147 G4ThreeVector ez(0, fKappa * xx.x(), 1);
148 G4ThreeVector normal = fHandedness*(er.cross(ez));
149
150 if (isGlobal) {
151 fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
152 } else {
153 fCurrentNormal.normal = normal.unit();
154 }
155 return fCurrentNormal.normal;
156}
157
158//=====================================================================
159//* DistanceToSurface -------------------------------------------------
160
161G4int G4TwistTubsSide::DistanceToSurface(const G4ThreeVector &gp,
162 const G4ThreeVector &gv,
163 G4ThreeVector gxx[],
164 G4double distance[],
165 G4int areacode[],
166 G4bool isvalid[],
167 EValidate validate)
168{
169 // Coordinate system:
170 //
171 // The coordinate system is so chosen that the intersection of
172 // the twisted surface with the z=0 plane coincides with the
173 // x-axis.
174 // Rotation matrix from this coordinate system (local system)
175 // to global system is saved in fRot field.
176 // So the (global) particle position and (global) velocity vectors,
177 // p and v, should be rotated fRot.inverse() in order to convert
178 // to local vectors.
179 //
180 // Equation of a twisted surface:
181 //
182 // x(rho(z=0), z) = rho(z=0)
183 // y(rho(z=0), z) = rho(z=0)*K*z
184 // z(rho(z=0), z) = z
185 // with
186 // K = std::tan(fPhiTwist/2)/fZHalfLen
187 //
188 // Equation of a line:
189 //
190 // gxx = p + t*v
191 // with
192 // p = fRot.inverse()*gp
193 // v = fRot.inverse()*gv
194 //
195 // Solution for intersection:
196 //
197 // Required time for crossing is given by solving the
198 // following quadratic equation:
199 //
200 // a*t^2 + b*t + c = 0
201 //
202 // where
203 //
204 // a = K*v_x*v_z
205 // b = K*(v_x*p_z + v_z*p_x) - v_y
206 // c = K*p_x*p_z - p_y
207 //
208 // Out of the possible two solutions you must choose
209 // the one that gives a positive rho(z=0).
210 //
211 //
212
213 fCurStatWithV.ResetfDone(validate, &gp, &gv);
214
215 if (fCurStatWithV.IsDone()) {
216 G4int i;
217 for (i=0; i<fCurStatWithV.GetNXX(); i++) {
218 gxx[i] = fCurStatWithV.GetXX(i);
219 distance[i] = fCurStatWithV.GetDistance(i);
220 areacode[i] = fCurStatWithV.GetAreacode(i);
221 isvalid[i] = fCurStatWithV.IsValid(i);
222 }
223 return fCurStatWithV.GetNXX();
224 } else {
225 // initialize
226 G4int i;
227 for (i=0; i<2; i++) {
228 distance[i] = kInfinity;
229 areacode[i] = sOutside;
230 isvalid[i] = false;
231 gxx[i].set(kInfinity, kInfinity, kInfinity);
232 }
233 }
234
235 G4ThreeVector p = ComputeLocalPoint(gp);
236 G4ThreeVector v = ComputeLocalDirection(gv);
237 G4ThreeVector xx[2];
238
239
240 //
241 // special case!
242 // p is origin or
243 //
244
245 G4double absvz = std::fabs(v.z());
246
247 if ((absvz < DBL_MIN) && (std::fabs(p.x() * v.y() - p.y() * v.x()) < DBL_MIN)) {
248 // no intersection
249
250 isvalid[0] = false;
251 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
252 isvalid[0], 0, validate, &gp, &gv);
253 return 0;
254 }
255
256 //
257 // special case end
258 //
259
260
261 G4double a = fKappa * v.x() * v.z();
262 G4double b = fKappa * (v.x() * p.z() + v.z() * p.x()) - v.y();
263 G4double c = fKappa * p.x() * p.z() - p.y();
264 G4double D = b * b - 4 * a * c; // discriminant
265
266 if (std::fabs(a) < DBL_MIN) {
267 if (std::fabs(b) > DBL_MIN) {
268
269 // single solution
270
271 distance[0] = - c / b;
272 xx[0] = p + distance[0]*v;
273 gxx[0] = ComputeGlobalPoint(xx[0]);
274
275 if (validate == kValidateWithTol) {
276 areacode[0] = GetAreaCode(xx[0]);
277 if (!IsOutside(areacode[0])) {
278 if (distance[0] >= 0) isvalid[0] = true;
279 }
280 } else if (validate == kValidateWithoutTol) {
281 areacode[0] = GetAreaCode(xx[0], false);
282 if (IsInside(areacode[0])) {
283 if (distance[0] >= 0) isvalid[0] = true;
284 }
285 } else { // kDontValidate
286 // we must omit x(rho,z) = rho(z=0) < 0
287 if (xx[0].x() > 0) {
288 areacode[0] = sInside;
289 if (distance[0] >= 0) isvalid[0] = true;
290 } else {
291 distance[0] = kInfinity;
292 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
293 areacode[0], isvalid[0],
294 0, validate, &gp, &gv);
295 return 0;
296 }
297 }
298
299 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
300 isvalid[0], 1, validate, &gp, &gv);
301 return 1;
302
303 } else {
304 // if a=b=0 , v.y=0 and (v.x=0 && p.x=0) or (v.z=0 && p.z=0) .
305 // if v.x=0 && p.x=0, no intersection unless p is on z-axis
306 // (in that case, v is paralell to surface).
307 // if v.z=0 && p.z=0, no intersection unless p is on x-axis
308 // (in that case, v is paralell to surface).
309 // return distance = infinity.
310
311 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
312 isvalid[0], 0, validate, &gp, &gv);
313
314 return 0;
315 }
316
317 } else if (D > DBL_MIN) {
318
319 // double solutions
320
321 D = std::sqrt(D);
322 G4double factor = 0.5/a;
323 G4double tmpdist[2] = {kInfinity, kInfinity};
324 G4ThreeVector tmpxx[2];
325 G4int tmpareacode[2] = {sOutside, sOutside};
326 G4bool tmpisvalid[2] = {false, false};
327 G4int i;
328
329 for (i=0; i<2; i++) {
330 G4double bminusD = - b - D;
331
332 // protection against round off error
333 //G4double protection = 1.0e-6;
334 G4double protection = 0;
335 if ( b * D < 0 && std::fabs(bminusD / D) < protection ) {
336 G4double acovbb = (a*c)/(b*b);
337 tmpdist[i] = - c/b * ( 1 - acovbb * (1 + 2*acovbb));
338 } else {
339 tmpdist[i] = factor * bminusD;
340 }
341
342 D = -D;
343 tmpxx[i] = p + tmpdist[i]*v;
344
345 if (validate == kValidateWithTol) {
346 tmpareacode[i] = GetAreaCode(tmpxx[i]);
347 if (!IsOutside(tmpareacode[i])) {
348 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
349 continue;
350 }
351 } else if (validate == kValidateWithoutTol) {
352 tmpareacode[i] = GetAreaCode(tmpxx[i], false);
353 if (IsInside(tmpareacode[i])) {
354 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
355 continue;
356 }
357 } else { // kDontValidate
358 // we must choose x(rho,z) = rho(z=0) > 0
359 if (tmpxx[i].x() > 0) {
360 tmpareacode[i] = sInside;
361 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
362 continue;
363 } else {
364 tmpdist[i] = kInfinity;
365 continue;
366 }
367 }
368 }
369
370 if (tmpdist[0] <= tmpdist[1]) {
371 distance[0] = tmpdist[0];
372 distance[1] = tmpdist[1];
373 xx[0] = tmpxx[0];
374 xx[1] = tmpxx[1];
375 gxx[0] = ComputeGlobalPoint(tmpxx[0]);
376 gxx[1] = ComputeGlobalPoint(tmpxx[1]);
377 areacode[0] = tmpareacode[0];
378 areacode[1] = tmpareacode[1];
379 isvalid[0] = tmpisvalid[0];
380 isvalid[1] = tmpisvalid[1];
381 } else {
382 distance[0] = tmpdist[1];
383 distance[1] = tmpdist[0];
384 xx[0] = tmpxx[1];
385 xx[1] = tmpxx[0];
386 gxx[0] = ComputeGlobalPoint(tmpxx[1]);
387 gxx[1] = ComputeGlobalPoint(tmpxx[0]);
388 areacode[0] = tmpareacode[1];
389 areacode[1] = tmpareacode[0];
390 isvalid[0] = tmpisvalid[1];
391 isvalid[1] = tmpisvalid[0];
392 }
393
394 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
395 isvalid[0], 2, validate, &gp, &gv);
396 fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
397 isvalid[1], 2, validate, &gp, &gv);
398
399 // protection against roundoff error
400
401 for (G4int k=0; k<2; k++) {
402 if (!isvalid[k]) continue;
403
404 G4ThreeVector xxonsurface(xx[k].x(), fKappa * std::fabs(xx[k].x())
405 * xx[k].z() , xx[k].z());
406 G4double deltaY = (xx[k] - xxonsurface).mag();
407
408 if ( deltaY > 0.5*kCarTolerance ) {
409
410 G4int maxcount = 10;
411 G4int l;
412 G4double lastdeltaY = deltaY;
413 G4ThreeVector last = deltaY;
414 for (l=0; l<maxcount; l++) {
415 G4ThreeVector surfacenormal = GetNormal(xxonsurface);
416 distance[k] = DistanceToPlaneWithV(p, v, xxonsurface,
417 surfacenormal, xx[k]);
418 deltaY = (xx[k] - xxonsurface).mag();
419 if (deltaY > lastdeltaY) {
420
421 }
422 gxx[k] = ComputeGlobalPoint(xx[k]);
423
424 if (deltaY <= 0.5*kCarTolerance) {
425
426 break;
427 }
428 xxonsurface.set(xx[k].x(),
429 fKappa * std::fabs(xx[k].x()) * xx[k].z(),
430 xx[k].z());
431 }
432 if (l == maxcount) {
433 G4cerr << "ERROR - G4TwistTubsSide::DistanceToSurface(p,v)"
434 << G4endl
435 << " maxloop count " << maxcount << G4endl;
436 G4Exception("G4TwistTubsFlatSide::DistanceToSurface(p,v)",
437 "InvalidSetup", FatalException,
438 "Exceeded maxloop count!");
439 }
440 }
441
442 }
443 return 2;
444 } else {
445 // if D<0, no solution
446 // if D=0, just grazing the surfaces, return kInfinity
447
448 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
449 isvalid[0], 0, validate, &gp, &gv);
450
451 return 0;
452 }
453 G4Exception("G4TwistTubsSide::DistanceToSurface(p,v)",
454 "InvalidCondition", FatalException, "Illegal operation !");
455 return 1;
456}
457
458//=====================================================================
459//* DistanceToSurface -------------------------------------------------
460
461G4int G4TwistTubsSide::DistanceToSurface(const G4ThreeVector &gp,
462 G4ThreeVector gxx[],
463 G4double distance[],
464 G4int areacode[])
465{
466 fCurStat.ResetfDone(kDontValidate, &gp);
467 G4int i = 0;
468 if (fCurStat.IsDone()) {
469 for (i=0; i<fCurStat.GetNXX(); i++) {
470 gxx[i] = fCurStat.GetXX(i);
471 distance[i] = fCurStat.GetDistance(i);
472 areacode[i] = fCurStat.GetAreacode(i);
473 }
474 return fCurStat.GetNXX();
475 } else {
476 // initialize
477 for (i=0; i<2; i++) {
478 distance[i] = kInfinity;
479 areacode[i] = sOutside;
480 gxx[i].set(kInfinity, kInfinity, kInfinity);
481 }
482 }
483
484 static const G4double halftol = 0.5 * kCarTolerance;
485
486 G4ThreeVector p = ComputeLocalPoint(gp);
487 G4ThreeVector xx;
488 G4int parity = (fKappa >= 0 ? 1 : -1);
489
490 //
491 // special case!
492 // If p is on surface, or
493 // p is on z-axis,
494 // return here immediatery.
495 //
496
497 G4ThreeVector lastgxx[2];
498 G4double distfromlast[2];
499 for (i=0; i<2; i++) {
500 lastgxx[i] = fCurStatWithV.GetXX(i);
501 distfromlast[i] = (gp - lastgxx[i]).mag();
502 }
503
504 if ((gp - lastgxx[0]).mag() < halftol
505 || (gp - lastgxx[1]).mag() < halftol) {
506 // last winner, or last poststep point is on the surface.
507 xx = p;
508 distance[0] = 0;
509 gxx[0] = gp;
510
511 G4bool isvalid = true;
512 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
513 isvalid, 1, kDontValidate, &gp);
514 return 1;
515 }
516
517 if (p.getRho() == 0) {
518 // p is on z-axis. Namely, p is on twisted surface (invalid area).
519 // We must return here, however, returning distance to x-minimum
520 // boundary is better than return 0-distance.
521 //
522 G4bool isvalid = true;
523 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
524 distance[0] = DistanceToBoundary(sAxis0 & sAxisMin, xx, p);
525 areacode[0] = sInside;
526 } else {
527 distance[0] = 0;
528 xx.set(0., 0., 0.);
529 }
530 gxx[0] = ComputeGlobalPoint(xx);
531 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
532 isvalid, 0, kDontValidate, &gp);
533 return 1;
534 }
535
536 //
537 // special case end
538 //
539
540 // set corner points of quadrangle try area ...
541
542 G4ThreeVector A; // foot of normal from p to boundary of sAxis0 & sAxisMin
543 G4ThreeVector C; // foot of normal from p to boundary of sAxis0 & sAxisMax
544 G4ThreeVector B; // point on boundary sAxis0 & sAxisMax at z = A.z()
545 G4ThreeVector D; // point on boundary sAxis0 & sAxisMin at z = C.z()
546 G4double distToA; // distance from p to A
547 G4double distToC; // distance from p to C
548
549 distToA = DistanceToBoundary(sAxis0 & sAxisMin, A, p);
550 distToC = DistanceToBoundary(sAxis0 & sAxisMax, C, p);
551
552 // is p.z between a.z and c.z?
553 // p.z must be bracketed a.z and c.z.
554 if (A.z() > C.z()) {
555 if (p.z() > A.z()) {
556 A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
557 } else if (p.z() < C.z()) {
558 C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
559 }
560 } else {
561 if (p.z() > C.z()) {
562 C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
563 } else if (p.z() < A.z()) {
564 A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
565 }
566 }
567
568 G4ThreeVector d[2]; // direction vectors of boundary
569 G4ThreeVector x0[2]; // foot of normal from line to p
570 G4int btype[2]; // boundary type
571
572 for (i=0; i<2; i++) {
573 if (i == 0) {
574 GetBoundaryParameters((sAxis0 & sAxisMax), d[i], x0[i], btype[i]);
575 B = x0[i] + ((A.z() - x0[i].z()) / d[i].z()) * d[i];
576 // x0 + t*d , d is direction unit vector.
577 } else {
578 GetBoundaryParameters((sAxis0 & sAxisMin), d[i], x0[i], btype[i]);
579 D = x0[i] + ((C.z() - x0[i].z()) / d[i].z()) * d[i];
580 }
581 }
582
583 // In order to set correct diagonal, swap A and D, C and B if needed.
584 G4ThreeVector pt(p.x(), p.y(), 0.);
585 G4double rc = std::fabs(p.x());
586 G4ThreeVector surfacevector(rc, rc * fKappa * p.z(), 0.);
587 G4int pside = AmIOnLeftSide(pt, surfacevector);
588 G4double test = (A.z() - C.z()) * parity * pside;
589
590 if (test == 0) {
591 if (pside == 0) {
592 // p is on surface.
593 xx = p;
594 distance[0] = 0;
595 gxx[0] = gp;
596
597 G4bool isvalid = true;
598 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
599 isvalid, 1, kDontValidate, &gp);
600 return 1;
601 } else {
602 // A.z = C.z(). return distance to line.
603 d[0] = C - A;
604 distance[0] = DistanceToLine(p, A, d[0], xx);
605 areacode[0] = sInside;
606 gxx[0] = ComputeGlobalPoint(xx);
607 G4bool isvalid = true;
608 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
609 isvalid, 1, kDontValidate, &gp);
610 return 1;
611 }
612
613 } else if (test < 0) {
614
615 // wrong diagonal. vector AC is crossing the surface!
616 // swap A and D, C and B
617 G4ThreeVector tmp;
618 tmp = A;
619 A = D;
620 D = tmp;
621 tmp = C;
622 C = B;
623 B = tmp;
624
625 } else {
626 // correct diagonal. nothing to do.
627 }
628
629 // Now, we chose correct diaglnal.
630 // First try. divide quadrangle into double triangle by diagonal and
631 // calculate distance to both surfaces.
632
633 G4ThreeVector xxacb; // foot of normal from plane ACB to p
634 G4ThreeVector nacb; // normal of plane ACD
635 G4ThreeVector xxcad; // foot of normal from plane CAD to p
636 G4ThreeVector ncad; // normal of plane CAD
637 G4ThreeVector AB(A.x(), A.y(), 0);
638 G4ThreeVector DC(C.x(), C.y(), 0);
639
640 G4double distToACB = G4VTwistSurface::DistanceToPlane(p, A, C-A, AB, xxacb, nacb) * parity;
641 G4double distToCAD = G4VTwistSurface::DistanceToPlane(p, C, C-A, DC, xxcad, ncad) * parity;
642
643 // if calculated distance = 0, return
644
645 if (std::fabs(distToACB) <= halftol || std::fabs(distToCAD) <= halftol) {
646 xx = (std::fabs(distToACB) < std::fabs(distToCAD) ? xxacb : xxcad);
647 areacode[0] = sInside;
648 gxx[0] = ComputeGlobalPoint(xx);
649 distance[0] = 0;
650 G4bool isvalid = true;
651 fCurStat.SetCurrentStatus(0, gxx[0], distance[0] , areacode[0],
652 isvalid, 1, kDontValidate, &gp);
653 return 1;
654 }
655
656 if (distToACB * distToCAD > 0 && distToACB < 0) {
657 // both distToACB and distToCAD are negative.
658 // divide quadrangle into double triangle by diagonal
659 G4ThreeVector normal;
660 distance[0] = DistanceToPlane(p, A, B, C, D, parity, xx, normal);
661 } else {
662 if (distToACB * distToCAD > 0) {
663 // both distToACB and distToCAD are positive.
664 // Take smaller one.
665 if (distToACB <= distToCAD) {
666 distance[0] = distToACB;
667 xx = xxacb;
668 } else {
669 distance[0] = distToCAD;
670 xx = xxcad;
671 }
672 } else {
673 // distToACB * distToCAD is negative.
674 // take positive one
675 if (distToACB > 0) {
676 distance[0] = distToACB;
677 xx = xxacb;
678 } else {
679 distance[0] = distToCAD;
680 xx = xxcad;
681 }
682 }
683
684 }
685 areacode[0] = sInside;
686 gxx[0] = ComputeGlobalPoint(xx);
687 G4bool isvalid = true;
688 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
689 isvalid, 1, kDontValidate, &gp);
690 return 1;
691}
692
693//=====================================================================
694//* DistanceToPlane ---------------------------------------------------
695
696G4double G4TwistTubsSide::DistanceToPlane(const G4ThreeVector &p,
697 const G4ThreeVector &A,
698 const G4ThreeVector &B,
699 const G4ThreeVector &C,
700 const G4ThreeVector &D,
701 const G4int parity,
702 G4ThreeVector &xx,
703 G4ThreeVector &n)
704{
705 static const G4double halftol = 0.5 * kCarTolerance;
706
707 G4ThreeVector M = 0.5*(A + B);
708 G4ThreeVector N = 0.5*(C + D);
709 G4ThreeVector xxanm; // foot of normal from p to plane ANM
710 G4ThreeVector nanm; // normal of plane ANM
711 G4ThreeVector xxcmn; // foot of normal from p to plane CMN
712 G4ThreeVector ncmn; // normal of plane CMN
713
714 G4double distToanm = G4VTwistSurface::DistanceToPlane(p, A, (N - A), (M - A), xxanm, nanm) * parity;
715 G4double distTocmn = G4VTwistSurface::DistanceToPlane(p, C, (M - C), (N - C), xxcmn, ncmn) * parity;
716
717 // if p is behind of both surfaces, abort.
718 if (distToanm * distTocmn > 0 && distToanm < 0) {
719 G4Exception("G4TwistTubsSide::DistanceToPlane()",
720 "InvalidCondition", FatalException,
721 "Point p is behind the surfaces.");
722 }
723
724 // if p is on surface, return 0.
725 if (std::fabs(distToanm) <= halftol) {
726 xx = xxanm;
727 n = nanm * parity;
728 return 0;
729 } else if (std::fabs(distTocmn) <= halftol) {
730 xx = xxcmn;
731 n = ncmn * parity;
732 return 0;
733 }
734
735 if (distToanm <= distTocmn) {
736 if (distToanm > 0) {
737 // both distanses are positive. take smaller one.
738 xx = xxanm;
739 n = nanm * parity;
740 return distToanm;
741 } else {
742 // take -ve distance and call the function recursively.
743 return DistanceToPlane(p, A, M, N, D, parity, xx, n);
744 }
745 } else {
746 if (distTocmn > 0) {
747 // both distanses are positive. take smaller one.
748 xx = xxcmn;
749 n = ncmn * parity;
750 return distTocmn;
751 } else {
752 // take -ve distance and call the function recursively.
753 return DistanceToPlane(p, C, N, M, B, parity, xx, n);
754 }
755 }
756}
757
758//=====================================================================
759//* GetAreaCode -------------------------------------------------------
760
761G4int G4TwistTubsSide::GetAreaCode(const G4ThreeVector &xx,
762 G4bool withTol)
763{
764 // We must use the function in local coordinate system.
765 // See the description of DistanceToSurface(p,v).
766
767 static const G4double ctol = 0.5 * kCarTolerance;
768 G4int areacode = sInside;
769
770 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
771 G4int xaxis = 0;
772 G4int zaxis = 1;
773
774 if (withTol) {
775
776 G4bool isoutside = false;
777
778 // test boundary of xaxis
779
780 if (xx.x() < fAxisMin[xaxis] + ctol) {
781 areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
782 if (xx.x() <= fAxisMin[xaxis] - ctol) isoutside = true;
783
784 } else if (xx.x() > fAxisMax[xaxis] - ctol) {
785 areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
786 if (xx.x() >= fAxisMax[xaxis] + ctol) isoutside = true;
787 }
788
789 // test boundary of z-axis
790
791 if (xx.z() < fAxisMin[zaxis] + ctol) {
792 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
793
794 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
795 else areacode |= sBoundary;
796 if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
797
798 } else if (xx.z() > fAxisMax[zaxis] - ctol) {
799 areacode |= (sAxis1 & (sAxisZ | sAxisMax));
800
801 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
802 else areacode |= sBoundary;
803 if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
804 }
805
806 // if isoutside = true, clear inside bit.
807 // if not on boundary, add axis information.
808
809 if (isoutside) {
810 G4int tmpareacode = areacode & (~sInside);
811 areacode = tmpareacode;
812 } else if ((areacode & sBoundary) != sBoundary) {
813 areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
814 }
815
816 } else {
817
818 // boundary of x-axis
819
820 if (xx.x() < fAxisMin[xaxis] ) {
821 areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
822 } else if (xx.x() > fAxisMax[xaxis]) {
823 areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
824 }
825
826 // boundary of z-axis
827
828 if (xx.z() < fAxisMin[zaxis]) {
829 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
830 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
831 else areacode |= sBoundary;
832
833 } else if (xx.z() > fAxisMax[zaxis]) {
834 areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ;
835 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
836 else areacode |= sBoundary;
837 }
838
839 if ((areacode & sBoundary) != sBoundary) {
840 areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
841 }
842 }
843 return areacode;
844 } else {
845 G4Exception("G4TwistTubsSide::GetAreaCode()",
846 "NotImplemented", FatalException,
847 "Feature NOT implemented !");
848 }
849 return areacode;
850}
851
852//=====================================================================
853//* SetCorners( arglist ) -------------------------------------------------
854
855void G4TwistTubsSide::SetCorners(
856 G4double endInnerRad[2],
857 G4double endOuterRad[2],
858 G4double endPhi[2],
859 G4double endZ[2])
860{
861 // Set Corner points in local coodinate.
862
863 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
864
865 G4int zmin = 0 ; // at -ve z
866 G4int zmax = 1 ; // at +ve z
867
868 G4double x, y, z;
869
870 // corner of Axis0min and Axis1min
871 x = endInnerRad[zmin]*std::cos(endPhi[zmin]);
872 y = endInnerRad[zmin]*std::sin(endPhi[zmin]);
873 z = endZ[zmin];
874 SetCorner(sC0Min1Min, x, y, z);
875
876 // corner of Axis0max and Axis1min
877 x = endOuterRad[zmin]*std::cos(endPhi[zmin]);
878 y = endOuterRad[zmin]*std::sin(endPhi[zmin]);
879 z = endZ[zmin];
880 SetCorner(sC0Max1Min, x, y, z);
881
882 // corner of Axis0max and Axis1max
883 x = endOuterRad[zmax]*std::cos(endPhi[zmax]);
884 y = endOuterRad[zmax]*std::sin(endPhi[zmax]);
885 z = endZ[zmax];
886 SetCorner(sC0Max1Max, x, y, z);
887
888 // corner of Axis0min and Axis1max
889 x = endInnerRad[zmax]*std::cos(endPhi[zmax]);
890 y = endInnerRad[zmax]*std::sin(endPhi[zmax]);
891 z = endZ[zmax];
892 SetCorner(sC0Min1Max, x, y, z);
893
894 } else {
895 G4cerr << "ERROR - G4TwistTubsFlatSide::SetCorners()" << G4endl
896 << " fAxis[0] = " << fAxis[0] << G4endl
897 << " fAxis[1] = " << fAxis[1] << G4endl;
898 G4Exception("G4TwistTubsSide::SetCorners()",
899 "NotImplemented", FatalException,
900 "Feature NOT implemented !");
901 }
902}
903
904//=====================================================================
905//* SetCorners() ------------------------------------------------------
906
907void G4TwistTubsSide::SetCorners()
908{
909 G4Exception("G4TwistTubsSide::SetCorners()",
910 "NotImplemented", FatalException,
911 "Method NOT implemented !");
912}
913
914//=====================================================================
915//* SetBoundaries() ---------------------------------------------------
916
917void G4TwistTubsSide::SetBoundaries()
918{
919 // Set direction-unit vector of boundary-lines in local coodinate.
920 //
921 G4ThreeVector direction;
922
923 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
924
925 // sAxis0 & sAxisMin
926 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
927 direction = direction.unit();
928 SetBoundary(sAxis0 & (sAxisX | sAxisMin), direction,
929 GetCorner(sC0Min1Min), sAxisZ) ;
930
931 // sAxis0 & sAxisMax
932 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
933 direction = direction.unit();
934 SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction,
935 GetCorner(sC0Max1Min), sAxisZ);
936
937 // sAxis1 & sAxisMin
938 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
939 direction = direction.unit();
940 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
941 GetCorner(sC0Min1Min), sAxisX);
942
943 // sAxis1 & sAxisMax
944 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
945 direction = direction.unit();
946 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
947 GetCorner(sC0Min1Max), sAxisX);
948
949 } else {
950 G4cerr << "ERROR - G4TwistTubsFlatSide::SetBoundaries()" << G4endl
951 << " fAxis[0] = " << fAxis[0] << G4endl
952 << " fAxis[1] = " << fAxis[1] << G4endl;
953 G4Exception("G4TwistTubsSide::SetCorners()",
954 "NotImplemented", FatalException,
955 "Feature NOT implemented !");
956 }
957}
958
959//=====================================================================
960//* GetFacets() -------------------------------------------------------
961
962void G4TwistTubsSide::GetFacets( G4int m, G4int n, G4double xyz[][3],
963 G4int faces[][4], G4int iside )
964{
965
966 G4double z ; // the two parameters for the surface equation
967 G4double x,xmin,xmax ;
968
969 G4ThreeVector p ; // a point on the surface, given by (z,u)
970
971 G4int nnode ;
972 G4int nface ;
973
974 // calculate the (n-1)*(m-1) vertices
975
976 G4int i,j ;
977
978 for ( i = 0 ; i<n ; i++ )
979 {
980
981 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
982
983 for ( j = 0 ; j<m ; j++ ) {
984
985 nnode = GetNode(i,j,m,n,iside) ;
986
987 xmin = GetBoundaryMin(z) ;
988 xmax = GetBoundaryMax(z) ;
989
990 if (fHandedness < 0) {
991 x = xmin + j*(xmax-xmin)/(m-1) ;
992 } else {
993 x = xmax - j*(xmax-xmin)/(m-1) ;
994 }
995
996 p = SurfacePoint(x,z,true) ; // surface point in global coord.system
997
998 xyz[nnode][0] = p.x() ;
999 xyz[nnode][1] = p.y() ;
1000 xyz[nnode][2] = p.z() ;
1001
1002 if ( i<n-1 && j<m-1 ) { // clock wise filling
1003
1004 nface = GetFace(i,j,m,n,iside) ;
1005
1006 faces[nface][0] = GetEdgeVisibility(i,j,m,n,0,1) * ( GetNode(i ,j ,m,n,iside)+1) ;
1007 faces[nface][1] = GetEdgeVisibility(i,j,m,n,1,1) * ( GetNode(i+1,j ,m,n,iside)+1) ;
1008 faces[nface][2] = GetEdgeVisibility(i,j,m,n,2,1) * ( GetNode(i+1,j+1,m,n,iside)+1) ;
1009 faces[nface][3] = GetEdgeVisibility(i,j,m,n,3,1) * ( GetNode(i ,j+1,m,n,iside)+1) ;
1010
1011 }
1012 }
1013 }
1014}
Note: See TracBrowser for help on using the repository browser.