source: src/vector.cpp@ eddea2

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since eddea2 was eddea2, checked in by Frederik Heber <heber@…>, 15 years ago

Added case '-T' (periodic translation) to testsuite.

  • BUGFIX: Vector::WrapPeriodically() added half of the box on top which is not as intended.
    • this broke testcase '-b', which relied on the half of the box adding.
    • new function molecule::DetermineCenterOfBox()
    • is used instead to additionally add this center in molecule::CenterInBox()
  • BUGFIX: molecule::TranslatePeriodically() subtracted vector instead of adding.
  • Property mode set to 100644
File size: 18.7 KB
Line 
1/** \file vector.cpp
2 *
3 * Function implementations for the class vector.
4 *
5 */
6
7
8#include "vector.hpp"
9#include "verbose.hpp"
10#include "World.hpp"
11#include "Helpers/Assert.hpp"
12#include "Helpers/fast_functions.hpp"
13
14#include <iostream>
15
16using namespace std;
17
18
19/************************************ Functions for class vector ************************************/
20
21/** Constructor of class vector.
22 */
23Vector::Vector()
24{
25 x[0] = x[1] = x[2] = 0.;
26};
27
28/**
29 * Copy constructor
30 */
31
32Vector::Vector(const Vector& src)
33{
34 x[0] = src[0];
35 x[1] = src[1];
36 x[2] = src[2];
37}
38
39/** Constructor of class vector.
40 */
41Vector::Vector(const double x1, const double x2, const double x3)
42{
43 x[0] = x1;
44 x[1] = x2;
45 x[2] = x3;
46};
47
48/**
49 * Assignment operator
50 */
51Vector& Vector::operator=(const Vector& src){
52 // check for self assignment
53 if(&src!=this){
54 x[0] = src[0];
55 x[1] = src[1];
56 x[2] = src[2];
57 }
58 return *this;
59}
60
61/** Desctructor of class vector.
62 */
63Vector::~Vector() {};
64
65/** Calculates square of distance between this and another vector.
66 * \param *y array to second vector
67 * \return \f$| x - y |^2\f$
68 */
69double Vector::DistanceSquared(const Vector &y) const
70{
71 double res = 0.;
72 for (int i=NDIM;i--;)
73 res += (x[i]-y[i])*(x[i]-y[i]);
74 return (res);
75};
76
77/** Calculates distance between this and another vector.
78 * \param *y array to second vector
79 * \return \f$| x - y |\f$
80 */
81double Vector::distance(const Vector &y) const
82{
83 return (sqrt(DistanceSquared(y)));
84};
85
86Vector Vector::getClosestPoint(const Vector &point) const{
87 // the closest point to a single point space is always the single point itself
88 return *this;
89}
90
91/** Calculates distance between this and another vector in a periodic cell.
92 * \param *y array to second vector
93 * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
94 * \return \f$| x - y |\f$
95 */
96double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
97{
98 double res = distance(y), tmp, matrix[NDIM*NDIM];
99 Vector Shiftedy, TranslationVector;
100 int N[NDIM];
101 matrix[0] = cell_size[0];
102 matrix[1] = cell_size[1];
103 matrix[2] = cell_size[3];
104 matrix[3] = cell_size[1];
105 matrix[4] = cell_size[2];
106 matrix[5] = cell_size[4];
107 matrix[6] = cell_size[3];
108 matrix[7] = cell_size[4];
109 matrix[8] = cell_size[5];
110 // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
111 for (N[0]=-1;N[0]<=1;N[0]++)
112 for (N[1]=-1;N[1]<=1;N[1]++)
113 for (N[2]=-1;N[2]<=1;N[2]++) {
114 // create the translation vector
115 TranslationVector.Zero();
116 for (int i=NDIM;i--;)
117 TranslationVector[i] = (double)N[i];
118 TranslationVector.MatrixMultiplication(matrix);
119 // add onto the original vector to compare with
120 Shiftedy = y + TranslationVector;
121 // get distance and compare with minimum so far
122 tmp = distance(Shiftedy);
123 if (tmp < res) res = tmp;
124 }
125 return (res);
126};
127
128/** Calculates distance between this and another vector in a periodic cell.
129 * \param *y array to second vector
130 * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
131 * \return \f$| x - y |^2\f$
132 */
133double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
134{
135 double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
136 Vector Shiftedy, TranslationVector;
137 int N[NDIM];
138 matrix[0] = cell_size[0];
139 matrix[1] = cell_size[1];
140 matrix[2] = cell_size[3];
141 matrix[3] = cell_size[1];
142 matrix[4] = cell_size[2];
143 matrix[5] = cell_size[4];
144 matrix[6] = cell_size[3];
145 matrix[7] = cell_size[4];
146 matrix[8] = cell_size[5];
147 // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
148 for (N[0]=-1;N[0]<=1;N[0]++)
149 for (N[1]=-1;N[1]<=1;N[1]++)
150 for (N[2]=-1;N[2]<=1;N[2]++) {
151 // create the translation vector
152 TranslationVector.Zero();
153 for (int i=NDIM;i--;)
154 TranslationVector[i] = (double)N[i];
155 TranslationVector.MatrixMultiplication(matrix);
156 // add onto the original vector to compare with
157 Shiftedy = y + TranslationVector;
158 // get distance and compare with minimum so far
159 tmp = DistanceSquared(Shiftedy);
160 if (tmp < res) res = tmp;
161 }
162 return (res);
163};
164
165/** Keeps the vector in a periodic cell, defined by the symmetric \a *matrix.
166 * \param *out ofstream for debugging messages
167 * Tries to translate a vector into each adjacent neighbouring cell.
168 */
169void Vector::KeepPeriodic(const double * const matrix)
170{
171 // int N[NDIM];
172 // bool flag = false;
173 //vector Shifted, TranslationVector;
174 // Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
175 // Log() << Verbose(2) << "Vector is: ";
176 // Output(out);
177 // Log() << Verbose(0) << endl;
178 InverseMatrixMultiplication(matrix);
179 for(int i=NDIM;i--;) { // correct periodically
180 if (at(i) < 0) { // get every coefficient into the interval [0,1)
181 at(i) += ceil(at(i));
182 } else {
183 at(i) -= floor(at(i));
184 }
185 }
186 MatrixMultiplication(matrix);
187 // Log() << Verbose(2) << "New corrected vector is: ";
188 // Output(out);
189 // Log() << Verbose(0) << endl;
190 // Log() << Verbose(1) << "End of KeepPeriodic." << endl;
191};
192
193/** Calculates scalar product between this and another vector.
194 * \param *y array to second vector
195 * \return \f$\langle x, y \rangle\f$
196 */
197double Vector::ScalarProduct(const Vector &y) const
198{
199 double res = 0.;
200 for (int i=NDIM;i--;)
201 res += x[i]*y[i];
202 return (res);
203};
204
205
206/** Calculates VectorProduct between this and another vector.
207 * -# returns the Product in place of vector from which it was initiated
208 * -# ATTENTION: Only three dim.
209 * \param *y array to vector with which to calculate crossproduct
210 * \return \f$ x \times y \f&
211 */
212void Vector::VectorProduct(const Vector &y)
213{
214 Vector tmp;
215 tmp[0] = x[1]* (y[2]) - x[2]* (y[1]);
216 tmp[1] = x[2]* (y[0]) - x[0]* (y[2]);
217 tmp[2] = x[0]* (y[1]) - x[1]* (y[0]);
218 (*this) = tmp;
219};
220
221
222/** projects this vector onto plane defined by \a *y.
223 * \param *y normal vector of plane
224 * \return \f$\langle x, y \rangle\f$
225 */
226void Vector::ProjectOntoPlane(const Vector &y)
227{
228 Vector tmp;
229 tmp = y;
230 tmp.Normalize();
231 tmp.Scale(ScalarProduct(tmp));
232 *this -= tmp;
233};
234
235/** Calculates the minimum distance vector of this vector to the plane.
236 * \param *out output stream for debugging
237 * \param *PlaneNormal normal of plane
238 * \param *PlaneOffset offset of plane
239 * \return distance to plane
240 * \return distance vector onto to plane
241 */
242Vector Vector::GetDistanceVectorToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
243{
244 Vector temp = (*this) - PlaneOffset;
245 temp.MakeNormalTo(PlaneNormal);
246 temp.Scale(-1.);
247 // then add connecting vector from plane to point
248 temp += (*this)-PlaneOffset;
249 double sign = temp.ScalarProduct(PlaneNormal);
250 if (fabs(sign) > MYEPSILON)
251 sign /= fabs(sign);
252 else
253 sign = 0.;
254
255 temp.Normalize();
256 temp.Scale(sign);
257 return temp;
258};
259
260
261/** Calculates the minimum distance of this vector to the plane.
262 * \sa Vector::GetDistanceVectorToPlane()
263 * \param *out output stream for debugging
264 * \param *PlaneNormal normal of plane
265 * \param *PlaneOffset offset of plane
266 * \return distance to plane
267 */
268double Vector::DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
269{
270 return GetDistanceVectorToPlane(PlaneNormal,PlaneOffset).Norm();
271};
272
273/** Calculates the projection of a vector onto another \a *y.
274 * \param *y array to second vector
275 */
276void Vector::ProjectIt(const Vector &y)
277{
278 (*this) += (-ScalarProduct(y))*y;
279};
280
281/** Calculates the projection of a vector onto another \a *y.
282 * \param *y array to second vector
283 * \return Vector
284 */
285Vector Vector::Projection(const Vector &y) const
286{
287 Vector helper = y;
288 helper.Scale((ScalarProduct(y)/y.NormSquared()));
289
290 return helper;
291};
292
293/** Calculates norm of this vector.
294 * \return \f$|x|\f$
295 */
296double Vector::Norm() const
297{
298 return (sqrt(NormSquared()));
299};
300
301/** Calculates squared norm of this vector.
302 * \return \f$|x|^2\f$
303 */
304double Vector::NormSquared() const
305{
306 return (ScalarProduct(*this));
307};
308
309/** Normalizes this vector.
310 */
311void Vector::Normalize()
312{
313 double factor = Norm();
314 (*this) *= 1/factor;
315};
316
317/** Zeros all components of this vector.
318 */
319void Vector::Zero()
320{
321 at(0)=at(1)=at(2)=0;
322};
323
324/** Zeros all components of this vector.
325 */
326void Vector::One(const double one)
327{
328 at(0)=at(1)=at(2)=one;
329};
330
331/** Checks whether vector has all components zero.
332 * @return true - vector is zero, false - vector is not
333 */
334bool Vector::IsZero() const
335{
336 return (fabs(x[0])+fabs(x[1])+fabs(x[2]) < MYEPSILON);
337};
338
339/** Checks whether vector has length of 1.
340 * @return true - vector is normalized, false - vector is not
341 */
342bool Vector::IsOne() const
343{
344 return (fabs(Norm() - 1.) < MYEPSILON);
345};
346
347/** Checks whether vector is normal to \a *normal.
348 * @return true - vector is normalized, false - vector is not
349 */
350bool Vector::IsNormalTo(const Vector &normal) const
351{
352 if (ScalarProduct(normal) < MYEPSILON)
353 return true;
354 else
355 return false;
356};
357
358/** Checks whether vector is normal to \a *normal.
359 * @return true - vector is normalized, false - vector is not
360 */
361bool Vector::IsEqualTo(const Vector &a) const
362{
363 bool status = true;
364 for (int i=0;i<NDIM;i++) {
365 if (fabs(x[i] - a[i]) > MYEPSILON)
366 status = false;
367 }
368 return status;
369};
370
371/** Calculates the angle between this and another vector.
372 * \param *y array to second vector
373 * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
374 */
375double Vector::Angle(const Vector &y) const
376{
377 double norm1 = Norm(), norm2 = y.Norm();
378 double angle = -1;
379 if ((fabs(norm1) > MYEPSILON) && (fabs(norm2) > MYEPSILON))
380 angle = this->ScalarProduct(y)/norm1/norm2;
381 // -1-MYEPSILON occured due to numerical imprecision, catch ...
382 //Log() << Verbose(2) << "INFO: acos(-1) = " << acos(-1) << ", acos(-1+MYEPSILON) = " << acos(-1+MYEPSILON) << ", acos(-1-MYEPSILON) = " << acos(-1-MYEPSILON) << "." << endl;
383 if (angle < -1)
384 angle = -1;
385 if (angle > 1)
386 angle = 1;
387 return acos(angle);
388};
389
390
391double& Vector::operator[](size_t i){
392 ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
393 return x[i];
394}
395
396const double& Vector::operator[](size_t i) const{
397 ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
398 return x[i];
399}
400
401double& Vector::at(size_t i){
402 return (*this)[i];
403}
404
405const double& Vector::at(size_t i) const{
406 return (*this)[i];
407}
408
409double* Vector::get(){
410 return x;
411}
412
413/** Compares vector \a to vector \a b component-wise.
414 * \param a base vector
415 * \param b vector components to add
416 * \return a == b
417 */
418bool Vector::operator==(const Vector& b) const
419{
420 return IsEqualTo(b);
421};
422
423bool Vector::operator!=(const Vector& b) const
424{
425 return !IsEqualTo(b);
426}
427
428/** Sums vector \a to this lhs component-wise.
429 * \param a base vector
430 * \param b vector components to add
431 * \return lhs + a
432 */
433const Vector& Vector::operator+=(const Vector& b)
434{
435 this->AddVector(b);
436 return *this;
437};
438
439/** Subtracts vector \a from this lhs component-wise.
440 * \param a base vector
441 * \param b vector components to add
442 * \return lhs - a
443 */
444const Vector& Vector::operator-=(const Vector& b)
445{
446 this->SubtractVector(b);
447 return *this;
448};
449
450/** factor each component of \a a times a double \a m.
451 * \param a base vector
452 * \param m factor
453 * \return lhs.x[i] * m
454 */
455const Vector& operator*=(Vector& a, const double m)
456{
457 a.Scale(m);
458 return a;
459};
460
461/** Sums two vectors \a and \b component-wise.
462 * \param a first vector
463 * \param b second vector
464 * \return a + b
465 */
466Vector const Vector::operator+(const Vector& b) const
467{
468 Vector x = *this;
469 x.AddVector(b);
470 return x;
471};
472
473/** Subtracts vector \a from \b component-wise.
474 * \param a first vector
475 * \param b second vector
476 * \return a - b
477 */
478Vector const Vector::operator-(const Vector& b) const
479{
480 Vector x = *this;
481 x.SubtractVector(b);
482 return x;
483};
484
485/** Factors given vector \a a times \a m.
486 * \param a vector
487 * \param m factor
488 * \return m * a
489 */
490Vector const operator*(const Vector& a, const double m)
491{
492 Vector x(a);
493 x.Scale(m);
494 return x;
495};
496
497/** Factors given vector \a a times \a m.
498 * \param m factor
499 * \param a vector
500 * \return m * a
501 */
502Vector const operator*(const double m, const Vector& a )
503{
504 Vector x(a);
505 x.Scale(m);
506 return x;
507};
508
509ostream& operator<<(ostream& ost, const Vector& m)
510{
511 ost << "(";
512 for (int i=0;i<NDIM;i++) {
513 ost << m[i];
514 if (i != 2)
515 ost << ",";
516 }
517 ost << ")";
518 return ost;
519};
520
521
522void Vector::ScaleAll(const double *factor)
523{
524 for (int i=NDIM;i--;)
525 x[i] *= factor[i];
526};
527
528
529
530void Vector::Scale(const double factor)
531{
532 for (int i=NDIM;i--;)
533 x[i] *= factor;
534};
535
536/** Given a box by its matrix \a *M and its inverse *Minv the vector is made to point within that box.
537 * \param *M matrix of box
538 * \param *Minv inverse matrix
539 */
540void Vector::WrapPeriodically(const double * const M, const double * const Minv)
541{
542 MatrixMultiplication(Minv);
543 // truncate to [0,1] for each axis
544 for (int i=0;i<NDIM;i++) {
545 //x[i] += 0.5; // set to center of box
546 while (x[i] >= 1.)
547 x[i] -= 1.;
548 while (x[i] < 0.)
549 x[i] += 1.;
550 }
551 MatrixMultiplication(M);
552};
553
554/** Do a matrix multiplication.
555 * \param *matrix NDIM_NDIM array
556 */
557void Vector::MatrixMultiplication(const double * const M)
558{
559 // do the matrix multiplication
560 at(0) = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
561 at(1) = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
562 at(2) = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
563};
564
565/** Do a matrix multiplication with the \a *A' inverse.
566 * \param *matrix NDIM_NDIM array
567 */
568bool Vector::InverseMatrixMultiplication(const double * const A)
569{
570 double B[NDIM*NDIM];
571 double detA = RDET3(A);
572 double detAReci;
573
574 // calculate the inverse B
575 if (fabs(detA) > MYEPSILON) {; // RDET3(A) yields precisely zero if A irregular
576 detAReci = 1./detA;
577 B[0] = detAReci*RDET2(A[4],A[5],A[7],A[8]); // A_11
578 B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]); // A_12
579 B[2] = detAReci*RDET2(A[1],A[2],A[4],A[5]); // A_13
580 B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]); // A_21
581 B[4] = detAReci*RDET2(A[0],A[2],A[6],A[8]); // A_22
582 B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]); // A_23
583 B[6] = detAReci*RDET2(A[3],A[4],A[6],A[7]); // A_31
584 B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]); // A_32
585 B[8] = detAReci*RDET2(A[0],A[1],A[3],A[4]); // A_33
586
587 // do the matrix multiplication
588 at(0) = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
589 at(1) = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
590 at(2) = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
591
592 return true;
593 } else {
594 return false;
595 }
596};
597
598
599/** Creates this vector as the b y *factors' components scaled linear combination of the given three.
600 * this vector = x1*factors[0] + x2* factors[1] + x3*factors[2]
601 * \param *x1 first vector
602 * \param *x2 second vector
603 * \param *x3 third vector
604 * \param *factors three-component vector with the factor for each given vector
605 */
606void Vector::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
607{
608 (*this) = (factors[0]*x1) +
609 (factors[1]*x2) +
610 (factors[2]*x3);
611};
612
613/** Mirrors atom against a given plane.
614 * \param n[] normal vector of mirror plane.
615 */
616void Vector::Mirror(const Vector &n)
617{
618 double projection;
619 projection = ScalarProduct(n)/n.NormSquared(); // remove constancy from n (keep as logical one)
620 // withdraw projected vector twice from original one
621 for (int i=NDIM;i--;)
622 at(i) -= 2.*projection*n[i];
623};
624
625/** Calculates orthonormal vector to one given vectors.
626 * Just subtracts the projection onto the given vector from this vector.
627 * The removed part of the vector is Vector::Projection()
628 * \param *x1 vector
629 * \return true - success, false - vector is zero
630 */
631bool Vector::MakeNormalTo(const Vector &y1)
632{
633 bool result = false;
634 double factor = y1.ScalarProduct(*this)/y1.NormSquared();
635 Vector x1;
636 x1 = factor * y1;
637 SubtractVector(x1);
638 for (int i=NDIM;i--;)
639 result = result || (fabs(x[i]) > MYEPSILON);
640
641 return result;
642};
643
644/** Creates this vector as one of the possible orthonormal ones to the given one.
645 * Just scan how many components of given *vector are unequal to zero and
646 * try to get the skp of both to be zero accordingly.
647 * \param *vector given vector
648 * \return true - success, false - failure (null vector given)
649 */
650bool Vector::GetOneNormalVector(const Vector &GivenVector)
651{
652 int Components[NDIM]; // contains indices of non-zero components
653 int Last = 0; // count the number of non-zero entries in vector
654 int j; // loop variables
655 double norm;
656
657 for (j=NDIM;j--;)
658 Components[j] = -1;
659
660 // in two component-systems we need to find the one position that is zero
661 int zeroPos = -1;
662 // find two components != 0
663 for (j=0;j<NDIM;j++){
664 if (fabs(GivenVector[j]) > MYEPSILON)
665 Components[Last++] = j;
666 else
667 // this our zero Position
668 zeroPos = j;
669 }
670
671 switch(Last) {
672 case 3: // threecomponent system
673 // the position of the zero is arbitrary in three component systems
674 zeroPos = Components[2];
675 case 2: // two component system
676 norm = sqrt(1./(GivenVector[Components[1]]*GivenVector[Components[1]]) + 1./(GivenVector[Components[0]]*GivenVector[Components[0]]));
677 at(zeroPos) = 0.;
678 // in skp both remaining parts shall become zero but with opposite sign and third is zero
679 at(Components[1]) = -1./GivenVector[Components[1]] / norm;
680 at(Components[0]) = 1./GivenVector[Components[0]] / norm;
681 return true;
682 break;
683 case 1: // one component system
684 // set sole non-zero component to 0, and one of the other zero component pendants to 1
685 at((Components[0]+2)%NDIM) = 0.;
686 at((Components[0]+1)%NDIM) = 1.;
687 at(Components[0]) = 0.;
688 return true;
689 break;
690 default:
691 return false;
692 }
693};
694
695/** Adds vector \a *y componentwise.
696 * \param *y vector
697 */
698void Vector::AddVector(const Vector &y)
699{
700 for(int i=NDIM;i--;)
701 x[i] += y[i];
702}
703
704/** Adds vector \a *y componentwise.
705 * \param *y vector
706 */
707void Vector::SubtractVector(const Vector &y)
708{
709 for(int i=NDIM;i--;)
710 x[i] -= y[i];
711}
712
713/**
714 * Checks whether this vector is within the parallelepiped defined by the given three vectors and
715 * their offset.
716 *
717 * @param offest for the origin of the parallelepiped
718 * @param three vectors forming the matrix that defines the shape of the parallelpiped
719 */
720bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
721{
722 Vector a = (*this)-offset;
723 a.InverseMatrixMultiplication(parallelepiped);
724 bool isInside = true;
725
726 for (int i=NDIM;i--;)
727 isInside = isInside && ((a[i] <= 1) && (a[i] >= 0));
728
729 return isInside;
730}
731
732
733// some comonly used vectors
734const Vector zeroVec(0,0,0);
735const Vector e1(1,0,0);
736const Vector e2(0,1,0);
737const Vector e3(0,0,1);
Note: See TracBrowser for help on using the repository browser.