Garfield++ 4.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Vector.c
Go to the documentation of this file.
1/*
2(c) 2005, Supratik Mukhopadhayay, Nayana Majumdar
3*/
4// Note that Opt = 1 yields a forward rotation while -1 yields a -ve rotation.
5#define DEFINE_VGLOBAL
6
7#include <math.h>
8#include <stdio.h>
9#include <stdlib.h>
10
11#include "Vector.h"
12
13#ifdef __cplusplus
14namespace neBEM {
15#endif
16
17void VectorRotate_Rect3D(double Xin, double Yin, double Zin, double RotX,
18 double RotY, double RotZ,
19 int Opt, // 1 forward, -1 backward
20 double *Xout, double *Yout, double *Zout) {
21 double R11, R12, R13, // Follow Numerical Methods Using Matlab
22 R21, R22, R23, // J.H.Mathews and K.D.Fink
23 R31, R32, R33; // 4th Edition, Prentice-Hall of India Pvt. Ltd.
24 double X_X, Y_X, Z_X, // New Delhi, 2004
25 X_Y, Y_Y, Z_Y; // p.115
26
27 if ((fabs(RotX) < 1.0e-12) // A most happy unrotated situation
28 && (fabs(RotY) < 1.0e-12) && (fabs(RotZ) < 1.0e-12)) {
29 *Xout = Xin;
30 *Yout = Yin;
31 *Zout = Zin;
32 return;
33 }
34
35 // Rotation in X
36 if (fabs(RotX) < 1.0e-12) {
37 X_X = Xin;
38 Y_X = Yin;
39 Z_X = Zin;
40 } else {
41 R11 = 1.0;
42 R12 = 0.0;
43 R13 = 0.0;
44 R21 = 0.0;
45 R22 = cos(Opt * RotX);
46 R23 = -sin(Opt * RotX);
47 R31 = 0.0;
48 R32 = sin(Opt * RotX);
49 R33 = cos(Opt * RotX);
50 X_X = R11 * Xin + R12 * Yin + R13 * Zin; // separate matrix multiplication
51 Y_X = R21 * Xin + R22 * Yin + R23 * Zin; // function avoided intentionally.
52 Z_X = R31 * Xin + R32 * Yin + R33 * Zin;
53 }
54
55 // Rotation in Y
56 if (fabs(RotY) < 1.0e-12) {
57 X_Y = X_X;
58 Y_Y = Y_X;
59 Z_Y = Z_X;
60 } else {
61 R11 = cos(Opt * RotY);
62 R12 = 0.0;
63 R13 = sin(Opt * RotY);
64 R21 = 0.0;
65 R22 = 1.0;
66 R23 = 0.0;
67 R31 = -sin(Opt * RotY);
68 R32 = 0.0;
69 R33 = cos(Opt * RotY);
70 X_Y = R11 * X_X + R12 * Y_X + R13 * Z_X; // separate matrix multiplication
71 Y_Y = R21 * X_X + R22 * Y_X + R23 * Z_X; // function avoided intentionally.
72 Z_Y = R31 * X_X + R32 * Y_X + R33 * Z_X;
73 }
74
75 // Rotation in Z - final rotation
76 if (fabs(RotZ) < 1.0e-12) {
77 *Xout = X_Y;
78 *Yout = Y_Y;
79 *Zout = Z_Y;
80 } else {
81 R11 = cos(Opt * RotZ);
82 R12 = -sin(Opt * RotZ);
83 R13 = 0.0;
84 R21 = sin(Opt * RotZ);
85 R22 = cos(Opt * RotZ);
86 R23 = 0.0;
87 R31 = 0.0;
88 R32 = 0.0;
89 R33 = 1.0;
90 *Xout =
91 R11 * X_Y + R12 * Y_Y + R13 * Z_Y; // separate matrix multiplication
92 *Yout =
93 R21 * X_Y + R22 * Y_Y + R23 * Z_Y; // function avoided intentionally.
94 *Zout = R31 * X_Y + R32 * Y_Y + R33 * Z_Y;
95 }
96}
97
98// Rotation of coordinate system is effectively oppsite to rotating the vector
99// keeping the coordinate system the same, as in the above routine.
100void CoordRotate_Rect3D(double Xin, double Yin, double Zin, double RotX,
101 double RotY, double RotZ,
102 int Opt, // 1 forward, -1 backward
103 double *Xout, double *Yout, double *Zout) {
104 double R11, R12, R13, // Follow Numerical Methods Using Matlab
105 R21, R22, R23, // J.H.Mathews and K.D.Fink
106 R31, R32, R33; // 4th Edition, Prentice-Hall of India Pvt. Ltd.
107 double X_X, Y_X, Z_X, // New Delhi, 2004
108 X_Y, Y_Y, Z_Y; // p.115
109
110 if ((fabs(RotX) < 1.0e-12) // A most happy unrotated situation
111 && (fabs(RotY) < 1.0e-12) && (fabs(RotZ) < 1.0e-12)) {
112 *Xout = Xin;
113 *Yout = Yin;
114 *Zout = Zin;
115 return;
116 }
117
118 // Rotation in X
119 if (fabs(RotX) < 1.0e-12) {
120 X_X = Xin;
121 Y_X = Yin;
122 Z_X = Zin;
123 } else {
124 R11 = 1.0;
125 R12 = 0.0;
126 R13 = 0.0;
127 R21 = 0.0;
128 R22 = cos(Opt * RotX);
129 R23 = sin(Opt * RotX);
130 R31 = 0.0;
131 R32 = -sin(Opt * RotX);
132 R33 = cos(Opt * RotX);
133 X_X = R11 * Xin + R12 * Yin + R13 * Zin; // separate matrix multiplication
134 Y_X = R21 * Xin + R22 * Yin + R23 * Zin; // function avoided intentionally.
135 Z_X = R31 * Xin + R32 * Yin + R33 * Zin;
136 }
137
138 // Rotation in Y
139 if (fabs(RotY) < 1.0e-12) {
140 X_Y = X_X;
141 Y_Y = Y_X;
142 Z_Y = Z_X;
143 } else {
144 R11 = cos(Opt * RotY);
145 R12 = 0.0;
146 R13 = -sin(Opt * RotY);
147 R21 = 0.0;
148 R22 = 1.0;
149 R23 = 0.0;
150 R31 = sin(Opt * RotY);
151 R32 = 0.0;
152 R33 = cos(Opt * RotY);
153 X_Y = R11 * X_X + R12 * Y_X + R13 * Z_X; // separate matrix multiplication
154 Y_Y = R21 * X_X + R22 * Y_X + R23 * Z_X; // function avoided intentionally.
155 Z_Y = R31 * X_X + R32 * Y_X + R33 * Z_X;
156 }
157
158 // Rotation in Z - final rotation
159 if (fabs(RotZ) < 1.0e-12) {
160 *Xout = X_Y;
161 *Yout = Y_Y;
162 *Zout = Z_Y;
163 } else {
164 R11 = cos(Opt * RotZ);
165 R12 = sin(Opt * RotZ);
166 R13 = 0.0;
167 R21 = -sin(Opt * RotZ);
168 R22 = cos(Opt * RotZ);
169 R23 = 0.0;
170 R31 = 0.0;
171 R32 = 0.0;
172 R33 = 1.0;
173 *Xout =
174 R11 * X_Y + R12 * Y_Y + R13 * Z_Y; // separate matrix multiplication
175 *Yout =
176 R21 * X_Y + R22 * Y_Y + R23 * Z_Y; // function avoided intentionally.
177 *Zout = R31 * X_Y + R32 * Y_Y + R33 * Z_Y;
178 }
179}
180
181// Create a 3D point
182Point3D CreatePoint3D(double x, double y, double z) {
183 Point3D p;
184 p.X = x;
185 p.Y = y;
186 p.Z = z;
187
188 return (p);
189}
190
191// Get distance between two 3D poionts
192double GetDistancePoint3D(Point3D *a, Point3D *b) {
193 return (sqrt((b->X - a->X) * (b->X - a->X) + (b->Y - a->Y) * (b->Y - a->Y) +
194 (b->Z - a->Z) * (b->Z - a->Z)));
195}
196
197// Create a 3D distance vector from point a to point b
198Vector3D CreateDistanceVector3D(Point3D *a, Point3D *b) {
199 Vector3D v;
200
201 v.X = b->X - a->X;
202 v.Y = b->Y - a->Y;
203 v.Z = b->Z - a->Z;
204
205 return (v);
206}
207
208// Computes magnitude of a 3D vector
209double MagVector3D(Vector3D *A) {
210 double mag;
211
212 mag = sqrt(A->X * A->X + A->Y * A->Y + A->Z * A->Z);
213
214 return (mag);
215}
216
217// Computes vector dot product
218double Vector3DDotProduct(Vector3D *A, Vector3D *B) {
219 double product;
220
221 product = A->X * B->X + A->Y * B->Y + A->Z * B->Z;
222
223 return (product);
224}
225
226// Computes unit vector for a given vector
227Vector3D UnitVector3D(Vector3D *v) {
228 Vector3D u;
229
230 double mag = MagVector3D(v);
231 if (fabs(mag) <= 1.0e-12) {
232 printf("UnitVector3D: magnitude smaller than 1.0e-12; no normalization.\n");
233 u.X = v->X;
234 u.Y = v->Y;
235 u.Z = v->Z;
236 } else {
237 u.X = v->X / mag;
238 u.Y = v->Y / mag;
239 u.Z = v->Z / mag;
240 }
241
242 return (u);
243}
244
245// Computes vector cross product
246Vector3D Vector3DCrossProduct(Vector3D *A, Vector3D *B) {
247 Vector3D product;
248
249 product.X = A->Y * B->Z - A->Z * B->Y;
250 product.Y = A->Z * B->X - A->X * B->Z;
251 product.Z = A->X * B->Y - A->Y * B->X;
252
253 return (product);
254}
255
256// prints out co-ordinates of a 3D point
257int PrintPoint3D(Point3D A) {
258 printf("%lg %lg %lg", A.X, A.Y, A.Z);
259 return (0);
260}
261
262// prints out components of a 3D vector
263int PrintVector3D(Vector3D A) {
264 printf("%lg %lg %lg", A.X, A.Y, A.Z);
265 return (0);
266}
267
268// prints out components of 3D direction cosines
269int PrintDirnCosn3D(DirnCosn3D A) {
270 printf("XUnit: ");
271 PrintVector3D(A.XUnit);
272 printf("\n");
273 printf("YUnit: ");
274 PrintVector3D(A.YUnit);
275 printf("\n");
276 printf("ZUnit: ");
277 PrintVector3D(A.ZUnit);
278 printf("\n");
279 return (0);
280}
281
282// Translates a point to a new origin specified in terms of the existing system
283// Translating the axes by (xT, yT, zT) is equivalent to translating the point
284// by (-xT, -yT, -zT)
285Point3D TranslatePoint3D(Point3D *A, Point3D *Origin, int Sense) {
286 double InitialVector[4];
287 double TranslationMatrix[4][4] = {{1.0, 0.0, 0.0, 0.0},
288 {0.0, 1.0, 0.0, 0.0},
289 {0.0, 0.0, 1.0, 0.0},
290 {0.0, 0.0, 0.0, 1.0}};
291 double FinalVector[4];
292 Point3D TranslatedPt;
293
294 InitialVector[0] = A->X;
295 InitialVector[1] = A->Y;
296 InitialVector[2] = A->Z;
297 InitialVector[3] = 1.0;
298
299 switch (Sense) {
300 case 1:
301 TranslationMatrix[0][3] = -Origin->X;
302 TranslationMatrix[1][3] = -Origin->Y;
303 TranslationMatrix[2][3] = -Origin->Z;
304 break;
305
306 case -1:
307 TranslationMatrix[0][3] = Origin->X;
308 TranslationMatrix[1][3] = Origin->Y;
309 TranslationMatrix[2][3] = Origin->Z;
310 break;
311
312 default:
313 printf("Only forward and inverse senses are allowed ...\n");
314 exit(-1);
315 }
316
317 for (int i = 0; i < 4; ++i) {
318 FinalVector[i] = 0.0;
319 for (int j = 0; j < 4; ++j) {
320 FinalVector[i] += TranslationMatrix[i][j] * InitialVector[j];
321 }
322 }
323
324 TranslatedPt.X = FinalVector[0];
325 TranslatedPt.Y = FinalVector[1];
326 TranslatedPt.Z = FinalVector[2];
327 return (TranslatedPt);
328}
329
330// Position of a point in a coord system with direction cosines specified
331// in terms of the existing system
332// Theory from Rotation representation (Wikipedia)
333// Sense +1 implies we are moving towards the new coordinate system whose
334// direction cosines we know in terms of the original coordinate system.
335// This option can be invoked by choosing 'global2local' defined in Vector.h
336// Sense -1 implies we are moving towards the original coordinate system in
337// which we know the direction cosines of the new coordinate system.
338// This option can be invoked by choosing 'local2global' defined in Vector.h
339Point3D RotatePoint3D(Point3D *A, DirnCosn3D *DC, int Sense) {
340
341 double TransformationMatrix[3][3] = {{0.0, 0.0, 0.0},
342 {0.0, 0.0, 0.0},
343 {0.0, 0.0, 0.0}};
344 switch (Sense) {
345 case 1:
346 TransformationMatrix[0][0] = DC->XUnit.X;
347 TransformationMatrix[0][1] = DC->XUnit.Y;
348 TransformationMatrix[0][2] = DC->XUnit.Z;
349 TransformationMatrix[1][0] = DC->YUnit.X;
350 TransformationMatrix[1][1] = DC->YUnit.Y;
351 TransformationMatrix[1][2] = DC->YUnit.Z;
352 TransformationMatrix[2][0] = DC->ZUnit.X;
353 TransformationMatrix[2][1] = DC->ZUnit.Y;
354 TransformationMatrix[2][2] = DC->ZUnit.Z;
355 break;
356
357 case -1:
358 TransformationMatrix[0][0] = DC->XUnit.X;
359 TransformationMatrix[0][1] = DC->YUnit.X;
360 TransformationMatrix[0][2] = DC->ZUnit.X;
361 TransformationMatrix[1][0] = DC->XUnit.Y;
362 TransformationMatrix[1][1] = DC->YUnit.Y;
363 TransformationMatrix[1][2] = DC->ZUnit.Y;
364 TransformationMatrix[2][0] = DC->XUnit.Z;
365 TransformationMatrix[2][1] = DC->YUnit.Z;
366 TransformationMatrix[2][2] = DC->ZUnit.Z;
367 break;
368
369 default:
370 printf("Only forward and inverse senses are allowed ...\n");
371 exit(-1);
372 }
373
374 double InitialVector[3] = {A->X, A->Y, A->Z};
375 double FinalVector[3] = {0., 0., 0.};
376 for (int i = 0; i < 3; ++i) {
377 for (int j = 0; j < 3; ++j) {
378 FinalVector[i] += TransformationMatrix[i][j] * InitialVector[j];
379 }
380 }
381 Point3D RotatedPt;
382 RotatedPt.X = FinalVector[0];
383 RotatedPt.Y = FinalVector[1];
384 RotatedPt.Z = FinalVector[2];
385 return (RotatedPt);
386}
387
388// Vector in a coord system with direction cosines specified
389// in terms of the existing system
390// Theory from Rotation representation (Wikipedia)
391// Sense +1 implies we are moving towards the new coordinate system whose
392// direction cosines we know in terms of the original coordinate system.
393// This option can be invoked by choosing 'global2local' defined in Vector.h
394// Sense -1 implies we are moving towards the original coordinate system in
395// which we know the direction cosines of the new coordinate system.
396// This option can be invoked by choosing 'local2global' defined in Vector.h
397Vector3D RotateVector3D(Vector3D *A, DirnCosn3D *DC, int Sense) {
398 double TransformationMatrix[4][4] = {{0.0, 0.0, 0.0, 0.0},
399 {0.0, 0.0, 0.0, 0.0},
400 {0.0, 0.0, 0.0, 0.0},
401 {0.0, 0.0, 0.0, 1.0}};
402 switch (Sense) {
403 case 1:
404 TransformationMatrix[0][0] = DC->XUnit.X;
405 TransformationMatrix[0][1] = DC->XUnit.Y;
406 TransformationMatrix[0][2] = DC->XUnit.Z;
407 TransformationMatrix[1][0] = DC->YUnit.X;
408 TransformationMatrix[1][1] = DC->YUnit.Y;
409 TransformationMatrix[1][2] = DC->YUnit.Z;
410 TransformationMatrix[2][0] = DC->ZUnit.X;
411 TransformationMatrix[2][1] = DC->ZUnit.Y;
412 TransformationMatrix[2][2] = DC->ZUnit.Z;
413 break;
414
415 case -1:
416 TransformationMatrix[0][0] = DC->XUnit.X;
417 TransformationMatrix[0][1] = DC->YUnit.X;
418 TransformationMatrix[0][2] = DC->ZUnit.X;
419 TransformationMatrix[1][0] = DC->XUnit.Y;
420 TransformationMatrix[1][1] = DC->YUnit.Y;
421 TransformationMatrix[1][2] = DC->ZUnit.Y;
422 TransformationMatrix[2][0] = DC->XUnit.Z;
423 TransformationMatrix[2][1] = DC->YUnit.Z;
424 TransformationMatrix[2][2] = DC->ZUnit.Z;
425 break;
426
427 default:
428 printf("Only forward and inverse senses are allowed ...\n");
429 exit(-1);
430 }
431
432 double InitialVector[3] = {A->X, A->Y, A->Z};
433 double FinalVector[3] = {0., 0., 0.};
434 for (int i = 0; i < 3; ++i) {
435 for (int j = 0; j < 3; ++j) {
436 FinalVector[i] += TransformationMatrix[i][j] * InitialVector[j];
437 }
438 }
439 Vector3D RotatedVector;
440 RotatedVector.X = FinalVector[0];
441 RotatedVector.Y = FinalVector[1];
442 RotatedVector.Z = FinalVector[2];
443 return (RotatedVector);
444}
445
446// Transform point: Get the new coordinates of a point described in the original
447// coordinate system. The new system of coordinates is also
448// described in terms of the original coordiate system, the
449// origin being NewOrigin and new direction cosines being
450// termed as NewDirns.
451// Translate-Rotate sequence
452// Inputs: point in the original coordinate system, origin of the new system
453// w.r.t the original coord system, direction cosines of the new
454// system w.r.t the original coord system
455Point3D TransformPoint3D(Point3D *initial, Point3D *NewOrigin,
456 DirnCosn3D *NewDirns) {
457 Point3D TmpPoint, final;
458
459 TmpPoint = TranslatePoint3D(initial, NewOrigin, 1);
460 final = RotatePoint3D(&TmpPoint, NewDirns, 1);
461 return (final);
462}
463
464// Reflect a 3D point (p1) on a mirror which is defined by a bi-vector (n)
465// perpendicular to the mirror
466// It is assumed that the mirror passes through the origin
467Point3D ReflectPoint3DByMirrorAtOrigin(Point3D *p1, Vector3D *n) {
468 double matrix[3][3];
469 matrix[0][0] = -n->X * n->X + n->Y * n->Y + n->Z * n->Z;
470 matrix[0][1] = -2.0 * n->X * n->Y;
471 matrix[0][2] = -2.0 * n->X * n->Z;
472 matrix[1][0] = -2.0 * n->X * n->Y;
473 matrix[1][1] = n->X * n->X - n->Y * n->Y + n->Z * n->Z;
474 matrix[1][2] = -2.0 * n->Y * n->Z;
475 matrix[2][0] = -2.0 * n->X * n->Z;
476 matrix[2][1] = -2.0 * n->Y * n->Z;
477 matrix[2][2] = n->X * n->X + n->Y * n->Y - n->Z * n->Z;
478
479 Point3D p2; // reflected point
480 p2.X = matrix[0][0] * p1->X + matrix[0][1] * p1->Y + matrix[0][2] * p1->Z;
481 p2.Y = matrix[1][0] * p1->X + matrix[1][1] * p1->Y + matrix[1][2] * p1->Z;
482 p2.Z = matrix[2][0] * p1->X + matrix[2][1] * p1->Y + matrix[2][2] * p1->Z;
483
484 return (p2);
485} // ReflectPoint3DByMirrorAtOrigin ends
486
487#ifdef __cplusplus
488} // namespace
489#endif
ReflectPoint3DByMirrorAtOrigin
Point3D ReflectPoint3DByMirrorAtOrigin(Point3D *p1, Vector3D *n)
Definition: Vector.c:467
CreateDistanceVector3D
Vector3D CreateDistanceVector3D(Point3D *a, Point3D *b)
Definition: Vector.c:198
VectorRotate_Rect3D
void VectorRotate_Rect3D(double Xin, double Yin, double Zin, double RotX, double RotY, double RotZ, int Opt, double *Xout, double *Yout, double *Zout)
Definition: Vector.c:17
RotatePoint3D
Point3D RotatePoint3D(Point3D *A, DirnCosn3D *DC, int Sense)
Definition: Vector.c:339
CoordRotate_Rect3D
void CoordRotate_Rect3D(double Xin, double Yin, double Zin, double RotX, double RotY, double RotZ, int Opt, double *Xout, double *Yout, double *Zout)
Definition: Vector.c:100
TranslatePoint3D
Point3D TranslatePoint3D(Point3D *A, Point3D *Origin, int Sense)
Definition: Vector.c:285
MagVector3D
double MagVector3D(Vector3D *A)
Definition: Vector.c:209
RotateVector3D
Vector3D RotateVector3D(Vector3D *A, DirnCosn3D *DC, int Sense)
Definition: Vector.c:397
Vector3DDotProduct
double Vector3DDotProduct(Vector3D *A, Vector3D *B)
Definition: Vector.c:218
TransformPoint3D
Point3D TransformPoint3D(Point3D *initial, Point3D *NewOrigin, DirnCosn3D *NewDirns)
Definition: Vector.c:455
PrintVector3D
int PrintVector3D(Vector3D A)
Definition: Vector.c:263
GetDistancePoint3D
double GetDistancePoint3D(Point3D *a, Point3D *b)
Definition: Vector.c:192
UnitVector3D
Vector3D UnitVector3D(Vector3D *v)
Definition: Vector.c:227
Vector3DCrossProduct
Vector3D Vector3DCrossProduct(Vector3D *A, Vector3D *B)
Definition: Vector.c:246
PrintDirnCosn3D
int PrintDirnCosn3D(DirnCosn3D A)
Definition: Vector.c:269
CreatePoint3D
Point3D CreatePoint3D(double x, double y, double z)
Definition: Vector.c:182
PrintPoint3D
int PrintPoint3D(Point3D A)
Definition: Vector.c:257
DirnCosn3D::ZUnit
Vector3D ZUnit
Definition: Vector.h:38
DirnCosn3D::YUnit
Vector3D YUnit
Definition: Vector.h:37
DirnCosn3D::XUnit
Vector3D XUnit
Definition: Vector.h:36
Point3D
Definition: Vector.h:21
Point3D::X
double X
Definition: Vector.h:22
Point3D::Z
double Z
Definition: Vector.h:24
Point3D::Y
double Y
Definition: Vector.h:23
Vector3D::Z
double Z
Definition: Vector.h:31
Vector3D::Y
double Y
Definition: Vector.h:30
Vector3D::X
double X
Definition: Vector.h:29