Context Navigation

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

Last change on this file since 1315 was 1228, checked in by garnier, 14 years ago
update geant4.9.3 tag
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.6 2009/11/11 12:23:37 gcosmo Exp $
28	// GEANT4 tag $Name: geant4-09-03 $
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
51	G4TwistTubsSide::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
75	G4TwistTubsSide::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
110	G4TwistTubsSide::G4TwistTubsSide( __void__& a )
111	: G4VTwistSurface(a)
112	{
113	}
114
115
116	//=====================================================================
117	//* destructor --------------------------------------------------------
118
119	G4TwistTubsSide::~G4TwistTubsSide()
120	{
121	}
122
123	//=====================================================================
124	//* GetNormal ---------------------------------------------------------
125
126	G4ThreeVector 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
161	G4int 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)Kz
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	// at^2 + bt + c = 0
201	//
202	// where
203	//
204	// a = Kv_xv_z
205	// b = K(v_xp_z + v_z*p_x) - v_y
206	// c = Kp_xp_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 = (ac)/(bb);
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	for (l=0; l<maxcount; l++) {
414	G4ThreeVector surfacenormal = GetNormal(xxonsurface);
415	distance[k] = DistanceToPlaneWithV(p, v, xxonsurface,
416	surfacenormal, xx[k]);
417	deltaY = (xx[k] - xxonsurface).mag();
418	if (deltaY > lastdeltaY) {
419
420	}
421	gxx[k] = ComputeGlobalPoint(xx[k]);
422
423	if (deltaY <= 0.5*kCarTolerance) {
424
425	break;
426	}
427	xxonsurface.set(xx[k].x(),
428	fKappa * std::fabs(xx[k].x()) * xx[k].z(),
429	xx[k].z());
430	}
431	if (l == maxcount) {
432	G4cerr << "ERROR - G4TwistTubsSide::DistanceToSurface(p,v)"
433	<< G4endl
434	<< " maxloop count " << maxcount << G4endl;
435	G4Exception("G4TwistTubsFlatSide::DistanceToSurface(p,v)",
436	"InvalidSetup", FatalException,
437	"Exceeded maxloop count!");
438	}
439	}
440
441	}
442	return 2;
443	} else {
444	// if D<0, no solution
445	// if D=0, just grazing the surfaces, return kInfinity
446
447	fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
448	isvalid[0], 0, validate, &gp, &gv);
449
450	return 0;
451	}
452	G4Exception("G4TwistTubsSide::DistanceToSurface(p,v)",
453	"InvalidCondition", FatalException, "Illegal operation !");
454	return 1;
455	}
456
457	//=====================================================================
458	//* DistanceToSurface -------------------------------------------------
459
460	G4int G4TwistTubsSide::DistanceToSurface(const G4ThreeVector &gp,
461	G4ThreeVector gxx[],
462	G4double distance[],
463	G4int areacode[])
464	{
465	fCurStat.ResetfDone(kDontValidate, &gp);
466	G4int i = 0;
467	if (fCurStat.IsDone()) {
468	for (i=0; i<fCurStat.GetNXX(); i++) {
469	gxx[i] = fCurStat.GetXX(i);
470	distance[i] = fCurStat.GetDistance(i);
471	areacode[i] = fCurStat.GetAreacode(i);
472	}
473	return fCurStat.GetNXX();
474	} else {
475	// initialize
476	for (i=0; i<2; i++) {
477	distance[i] = kInfinity;
478	areacode[i] = sOutside;
479	gxx[i].set(kInfinity, kInfinity, kInfinity);
480	}
481	}
482
483	static const G4double halftol = 0.5 * kCarTolerance;
484
485	G4ThreeVector p = ComputeLocalPoint(gp);
486	G4ThreeVector xx;
487	G4int parity = (fKappa >= 0 ? 1 : -1);
488
489	//
490	// special case!
491	// If p is on surface, or
492	// p is on z-axis,
493	// return here immediatery.
494	//
495
496	G4ThreeVector lastgxx[2];
497	G4double distfromlast[2];
498	for (i=0; i<2; i++) {
499	lastgxx[i] = fCurStatWithV.GetXX(i);
500	distfromlast[i] = (gp - lastgxx[i]).mag();
501	}
502
503	if ((gp - lastgxx[0]).mag() < halftol
504	\|\| (gp - lastgxx[1]).mag() < halftol) {
505	// last winner, or last poststep point is on the surface.
506	xx = p;
507	distance[0] = 0;
508	gxx[0] = gp;
509
510	G4bool isvalid = true;
511	fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
512	isvalid, 1, kDontValidate, &gp);
513	return 1;
514	}
515
516	if (p.getRho() == 0) {
517	// p is on z-axis. Namely, p is on twisted surface (invalid area).
518	// We must return here, however, returning distance to x-minimum
519	// boundary is better than return 0-distance.
520	//
521	G4bool isvalid = true;
522	if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
523	distance[0] = DistanceToBoundary(sAxis0 & sAxisMin, xx, p);
524	areacode[0] = sInside;
525	} else {
526	distance[0] = 0;
527	xx.set(0., 0., 0.);
528	}
529	gxx[0] = ComputeGlobalPoint(xx);
530	fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
531	isvalid, 0, kDontValidate, &gp);
532	return 1;
533	}
534
535	//
536	// special case end
537	//
538
539	// set corner points of quadrangle try area ...
540
541	G4ThreeVector A; // foot of normal from p to boundary of sAxis0 & sAxisMin
542	G4ThreeVector C; // foot of normal from p to boundary of sAxis0 & sAxisMax
543	G4ThreeVector B; // point on boundary sAxis0 & sAxisMax at z = A.z()
544	G4ThreeVector D; // point on boundary sAxis0 & sAxisMin at z = C.z()
545	G4double distToA; // distance from p to A
546	G4double distToC; // distance from p to C
547
548	distToA = DistanceToBoundary(sAxis0 & sAxisMin, A, p);
549	distToC = DistanceToBoundary(sAxis0 & sAxisMax, C, p);
550
551	// is p.z between a.z and c.z?
552	// p.z must be bracketed a.z and c.z.
553	if (A.z() > C.z()) {
554	if (p.z() > A.z()) {
555	A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
556	} else if (p.z() < C.z()) {
557	C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
558	}
559	} else {
560	if (p.z() > C.z()) {
561	C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
562	} else if (p.z() < A.z()) {
563	A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
564	}
565	}
566
567	G4ThreeVector d[2]; // direction vectors of boundary
568	G4ThreeVector x0[2]; // foot of normal from line to p
569	G4int btype[2]; // boundary type
570
571	for (i=0; i<2; i++) {
572	if (i == 0) {
573	GetBoundaryParameters((sAxis0 & sAxisMax), d[i], x0[i], btype[i]);
574	B = x0[i] + ((A.z() - x0[i].z()) / d[i].z()) * d[i];
575	// x0 + t*d , d is direction unit vector.
576	} else {
577	GetBoundaryParameters((sAxis0 & sAxisMin), d[i], x0[i], btype[i]);
578	D = x0[i] + ((C.z() - x0[i].z()) / d[i].z()) * d[i];
579	}
580	}
581
582	// In order to set correct diagonal, swap A and D, C and B if needed.
583	G4ThreeVector pt(p.x(), p.y(), 0.);
584	G4double rc = std::fabs(p.x());
585	G4ThreeVector surfacevector(rc, rc * fKappa * p.z(), 0.);
586	G4int pside = AmIOnLeftSide(pt, surfacevector);
587	G4double test = (A.z() - C.z()) * parity * pside;
588
589	if (test == 0) {
590	if (pside == 0) {
591	// p is on surface.
592	xx = p;
593	distance[0] = 0;
594	gxx[0] = gp;
595
596	G4bool isvalid = true;
597	fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
598	isvalid, 1, kDontValidate, &gp);
599	return 1;
600	} else {
601	// A.z = C.z(). return distance to line.
602	d[0] = C - A;
603	distance[0] = DistanceToLine(p, A, d[0], xx);
604	areacode[0] = sInside;
605	gxx[0] = ComputeGlobalPoint(xx);
606	G4bool isvalid = true;
607	fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
608	isvalid, 1, kDontValidate, &gp);
609	return 1;
610	}
611
612	} else if (test < 0) {
613
614	// wrong diagonal. vector AC is crossing the surface!
615	// swap A and D, C and B
616	G4ThreeVector tmp;
617	tmp = A;
618	A = D;
619	D = tmp;
620	tmp = C;
621	C = B;
622	B = tmp;
623
624	} else {
625	// correct diagonal. nothing to do.
626	}
627
628	// Now, we chose correct diaglnal.
629	// First try. divide quadrangle into double triangle by diagonal and
630	// calculate distance to both surfaces.
631
632	G4ThreeVector xxacb; // foot of normal from plane ACB to p
633	G4ThreeVector nacb; // normal of plane ACD
634	G4ThreeVector xxcad; // foot of normal from plane CAD to p
635	G4ThreeVector ncad; // normal of plane CAD
636	G4ThreeVector AB(A.x(), A.y(), 0);
637	G4ThreeVector DC(C.x(), C.y(), 0);
638
639	G4double distToACB = G4VTwistSurface::DistanceToPlane(p, A, C-A, AB, xxacb, nacb) * parity;
640	G4double distToCAD = G4VTwistSurface::DistanceToPlane(p, C, C-A, DC, xxcad, ncad) * parity;
641
642	// if calculated distance = 0, return
643
644	if (std::fabs(distToACB) <= halftol \|\| std::fabs(distToCAD) <= halftol) {
645	xx = (std::fabs(distToACB) < std::fabs(distToCAD) ? xxacb : xxcad);
646	areacode[0] = sInside;
647	gxx[0] = ComputeGlobalPoint(xx);
648	distance[0] = 0;
649	G4bool isvalid = true;
650	fCurStat.SetCurrentStatus(0, gxx[0], distance[0] , areacode[0],
651	isvalid, 1, kDontValidate, &gp);
652	return 1;
653	}
654
655	if (distToACB * distToCAD > 0 && distToACB < 0) {
656	// both distToACB and distToCAD are negative.
657	// divide quadrangle into double triangle by diagonal
658	G4ThreeVector normal;
659	distance[0] = DistanceToPlane(p, A, B, C, D, parity, xx, normal);
660	} else {
661	if (distToACB * distToCAD > 0) {
662	// both distToACB and distToCAD are positive.
663	// Take smaller one.
664	if (distToACB <= distToCAD) {
665	distance[0] = distToACB;
666	xx = xxacb;
667	} else {
668	distance[0] = distToCAD;
669	xx = xxcad;
670	}
671	} else {
672	// distToACB * distToCAD is negative.
673	// take positive one
674	if (distToACB > 0) {
675	distance[0] = distToACB;
676	xx = xxacb;
677	} else {
678	distance[0] = distToCAD;
679	xx = xxcad;
680	}
681	}
682
683	}
684	areacode[0] = sInside;
685	gxx[0] = ComputeGlobalPoint(xx);
686	G4bool isvalid = true;
687	fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
688	isvalid, 1, kDontValidate, &gp);
689	return 1;
690	}
691
692	//=====================================================================
693	//* DistanceToPlane ---------------------------------------------------
694
695	G4double G4TwistTubsSide::DistanceToPlane(const G4ThreeVector &p,
696	const G4ThreeVector &A,
697	const G4ThreeVector &B,
698	const G4ThreeVector &C,
699	const G4ThreeVector &D,
700	const G4int parity,
701	G4ThreeVector &xx,
702	G4ThreeVector &n)
703	{
704	static const G4double halftol = 0.5 * kCarTolerance;
705
706	G4ThreeVector M = 0.5*(A + B);
707	G4ThreeVector N = 0.5*(C + D);
708	G4ThreeVector xxanm; // foot of normal from p to plane ANM
709	G4ThreeVector nanm; // normal of plane ANM
710	G4ThreeVector xxcmn; // foot of normal from p to plane CMN
711	G4ThreeVector ncmn; // normal of plane CMN
712
713	G4double distToanm = G4VTwistSurface::DistanceToPlane(p, A, (N - A), (M - A), xxanm, nanm) * parity;
714	G4double distTocmn = G4VTwistSurface::DistanceToPlane(p, C, (M - C), (N - C), xxcmn, ncmn) * parity;
715
716	// if p is behind of both surfaces, abort.
717	if (distToanm * distTocmn > 0 && distToanm < 0) {
718	G4Exception("G4TwistTubsSide::DistanceToPlane()",
719	"InvalidCondition", FatalException,
720	"Point p is behind the surfaces.");
721	}
722
723	// if p is on surface, return 0.
724	if (std::fabs(distToanm) <= halftol) {
725	xx = xxanm;
726	n = nanm * parity;
727	return 0;
728	} else if (std::fabs(distTocmn) <= halftol) {
729	xx = xxcmn;
730	n = ncmn * parity;
731	return 0;
732	}
733
734	if (distToanm <= distTocmn) {
735	if (distToanm > 0) {
736	// both distanses are positive. take smaller one.
737	xx = xxanm;
738	n = nanm * parity;
739	return distToanm;
740	} else {
741	// take -ve distance and call the function recursively.
742	return DistanceToPlane(p, A, M, N, D, parity, xx, n);
743	}
744	} else {
745	if (distTocmn > 0) {
746	// both distanses are positive. take smaller one.
747	xx = xxcmn;
748	n = ncmn * parity;
749	return distTocmn;
750	} else {
751	// take -ve distance and call the function recursively.
752	return DistanceToPlane(p, C, N, M, B, parity, xx, n);
753	}
754	}
755	}
756
757	//=====================================================================
758	//* GetAreaCode -------------------------------------------------------
759
760	G4int G4TwistTubsSide::GetAreaCode(const G4ThreeVector &xx,
761	G4bool withTol)
762	{
763	// We must use the function in local coordinate system.
764	// See the description of DistanceToSurface(p,v).
765
766	static const G4double ctol = 0.5 * kCarTolerance;
767	G4int areacode = sInside;
768
769	if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
770	G4int xaxis = 0;
771	G4int zaxis = 1;
772
773	if (withTol) {
774
775	G4bool isoutside = false;
776
777	// test boundary of xaxis
778
779	if (xx.x() < fAxisMin[xaxis] + ctol) {
780	areacode \|= (sAxis0 & (sAxisX \| sAxisMin)) \| sBoundary;
781	if (xx.x() <= fAxisMin[xaxis] - ctol) isoutside = true;
782
783	} else if (xx.x() > fAxisMax[xaxis] - ctol) {
784	areacode \|= (sAxis0 & (sAxisX \| sAxisMax)) \| sBoundary;
785	if (xx.x() >= fAxisMax[xaxis] + ctol) isoutside = true;
786	}
787
788	// test boundary of z-axis
789
790	if (xx.z() < fAxisMin[zaxis] + ctol) {
791	areacode \|= (sAxis1 & (sAxisZ \| sAxisMin));
792
793	if (areacode & sBoundary) areacode \|= sCorner; // xx is on the corner.
794	else areacode \|= sBoundary;
795	if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
796
797	} else if (xx.z() > fAxisMax[zaxis] - ctol) {
798	areacode \|= (sAxis1 & (sAxisZ \| sAxisMax));
799
800	if (areacode & sBoundary) areacode \|= sCorner; // xx is on the corner.
801	else areacode \|= sBoundary;
802	if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
803	}
804
805	// if isoutside = true, clear inside bit.
806	// if not on boundary, add axis information.
807
808	if (isoutside) {
809	G4int tmpareacode = areacode & (~sInside);
810	areacode = tmpareacode;
811	} else if ((areacode & sBoundary) != sBoundary) {
812	areacode \|= (sAxis0 & sAxisX) \| (sAxis1 & sAxisZ);
813	}
814
815	} else {
816
817	// boundary of x-axis
818
819	if (xx.x() < fAxisMin[xaxis] ) {
820	areacode \|= (sAxis0 & (sAxisX \| sAxisMin)) \| sBoundary;
821	} else if (xx.x() > fAxisMax[xaxis]) {
822	areacode \|= (sAxis0 & (sAxisX \| sAxisMax)) \| sBoundary;
823	}
824
825	// boundary of z-axis
826
827	if (xx.z() < fAxisMin[zaxis]) {
828	areacode \|= (sAxis1 & (sAxisZ \| sAxisMin));
829	if (areacode & sBoundary) areacode \|= sCorner; // xx is on the corner.
830	else areacode \|= sBoundary;
831
832	} else if (xx.z() > fAxisMax[zaxis]) {
833	areacode \|= (sAxis1 & (sAxisZ \| sAxisMax)) ;
834	if (areacode & sBoundary) areacode \|= sCorner; // xx is on the corner.
835	else areacode \|= sBoundary;
836	}
837
838	if ((areacode & sBoundary) != sBoundary) {
839	areacode \|= (sAxis0 & sAxisX) \| (sAxis1 & sAxisZ);
840	}
841	}
842	return areacode;
843	} else {
844	G4Exception("G4TwistTubsSide::GetAreaCode()",
845	"NotImplemented", FatalException,
846	"Feature NOT implemented !");
847	}
848	return areacode;
849	}
850
851	//=====================================================================
852	//* SetCorners( arglist ) -------------------------------------------------
853
854	void G4TwistTubsSide::SetCorners(
855	G4double endInnerRad[2],
856	G4double endOuterRad[2],
857	G4double endPhi[2],
858	G4double endZ[2])
859	{
860	// Set Corner points in local coodinate.
861
862	if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
863
864	G4int zmin = 0 ; // at -ve z
865	G4int zmax = 1 ; // at +ve z
866
867	G4double x, y, z;
868
869	// corner of Axis0min and Axis1min
870	x = endInnerRad[zmin]*std::cos(endPhi[zmin]);
871	y = endInnerRad[zmin]*std::sin(endPhi[zmin]);
872	z = endZ[zmin];
873	SetCorner(sC0Min1Min, x, y, z);
874
875	// corner of Axis0max and Axis1min
876	x = endOuterRad[zmin]*std::cos(endPhi[zmin]);
877	y = endOuterRad[zmin]*std::sin(endPhi[zmin]);
878	z = endZ[zmin];
879	SetCorner(sC0Max1Min, x, y, z);
880
881	// corner of Axis0max and Axis1max
882	x = endOuterRad[zmax]*std::cos(endPhi[zmax]);
883	y = endOuterRad[zmax]*std::sin(endPhi[zmax]);
884	z = endZ[zmax];
885	SetCorner(sC0Max1Max, x, y, z);
886
887	// corner of Axis0min and Axis1max
888	x = endInnerRad[zmax]*std::cos(endPhi[zmax]);
889	y = endInnerRad[zmax]*std::sin(endPhi[zmax]);
890	z = endZ[zmax];
891	SetCorner(sC0Min1Max, x, y, z);
892
893	} else {
894	G4cerr << "ERROR - G4TwistTubsFlatSide::SetCorners()" << G4endl
895	<< " fAxis[0] = " << fAxis[0] << G4endl
896	<< " fAxis[1] = " << fAxis[1] << G4endl;
897	G4Exception("G4TwistTubsSide::SetCorners()",
898	"NotImplemented", FatalException,
899	"Feature NOT implemented !");
900	}
901	}
902
903	//=====================================================================
904	//* SetCorners() ------------------------------------------------------
905
906	void G4TwistTubsSide::SetCorners()
907	{
908	G4Exception("G4TwistTubsSide::SetCorners()",
909	"NotImplemented", FatalException,
910	"Method NOT implemented !");
911	}
912
913	//=====================================================================
914	//* SetBoundaries() ---------------------------------------------------
915
916	void G4TwistTubsSide::SetBoundaries()
917	{
918	// Set direction-unit vector of boundary-lines in local coodinate.
919	//
920	G4ThreeVector direction;
921
922	if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
923
924	// sAxis0 & sAxisMin
925	direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
926	direction = direction.unit();
927	SetBoundary(sAxis0 & (sAxisX \| sAxisMin), direction,
928	GetCorner(sC0Min1Min), sAxisZ) ;
929
930	// sAxis0 & sAxisMax
931	direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
932	direction = direction.unit();
933	SetBoundary(sAxis0 & (sAxisX \| sAxisMax), direction,
934	GetCorner(sC0Max1Min), sAxisZ);
935
936	// sAxis1 & sAxisMin
937	direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
938	direction = direction.unit();
939	SetBoundary(sAxis1 & (sAxisZ \| sAxisMin), direction,
940	GetCorner(sC0Min1Min), sAxisX);
941
942	// sAxis1 & sAxisMax
943	direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
944	direction = direction.unit();
945	SetBoundary(sAxis1 & (sAxisZ \| sAxisMax), direction,
946	GetCorner(sC0Min1Max), sAxisX);
947
948	} else {
949	G4cerr << "ERROR - G4TwistTubsFlatSide::SetBoundaries()" << G4endl
950	<< " fAxis[0] = " << fAxis[0] << G4endl
951	<< " fAxis[1] = " << fAxis[1] << G4endl;
952	G4Exception("G4TwistTubsSide::SetCorners()",
953	"NotImplemented", FatalException,
954	"Feature NOT implemented !");
955	}
956	}
957
958	//=====================================================================
959	//* GetFacets() -------------------------------------------------------
960
961	void G4TwistTubsSide::GetFacets( G4int m, G4int n, G4double xyz[][3],
962	G4int faces[][4], G4int iside )
963	{
964
965	G4double z ; // the two parameters for the surface equation
966	G4double x,xmin,xmax ;
967
968	G4ThreeVector p ; // a point on the surface, given by (z,u)
969
970	G4int nnode ;
971	G4int nface ;
972
973	// calculate the (n-1)*(m-1) vertices
974
975	G4int i,j ;
976
977	for ( i = 0 ; i<n ; i++ )
978	{
979
980	z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
981
982	for ( j = 0 ; j<m ; j++ ) {
983
984	nnode = GetNode(i,j,m,n,iside) ;
985
986	xmin = GetBoundaryMin(z) ;
987	xmax = GetBoundaryMax(z) ;
988
989	if (fHandedness < 0) {
990	x = xmin + j*(xmax-xmin)/(m-1) ;
991	} else {
992	x = xmax - j*(xmax-xmin)/(m-1) ;
993	}
994
995	p = SurfacePoint(x,z,true) ; // surface point in global coord.system
996
997	xyz[nnode][0] = p.x() ;
998	xyz[nnode][1] = p.y() ;
999	xyz[nnode][2] = p.z() ;
1000
1001	if ( i<n-1 && j<m-1 ) { // clock wise filling
1002
1003	nface = GetFace(i,j,m,n,iside) ;
1004
1005	faces[nface][0] = GetEdgeVisibility(i,j,m,n,0,1) * ( GetNode(i ,j ,m,n,iside)+1) ;
1006	faces[nface][1] = GetEdgeVisibility(i,j,m,n,1,1) * ( GetNode(i+1,j ,m,n,iside)+1) ;
1007	faces[nface][2] = GetEdgeVisibility(i,j,m,n,2,1) * ( GetNode(i+1,j+1,m,n,iside)+1) ;
1008	faces[nface][3] = GetEdgeVisibility(i,j,m,n,3,1) * ( GetNode(i ,j+1,m,n,iside)+1) ;
1009
1010	}
1011	}
1012	}
1013	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: