1 | /*
|
---|
2 | * Project: MoleCuilder
|
---|
3 | * Description: creates and alters molecular systems
|
---|
4 | * Copyright (C) 2010 University of Bonn. All rights reserved.
|
---|
5 | * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
|
---|
6 | */
|
---|
7 |
|
---|
8 | /** \file vector.cpp
|
---|
9 | *
|
---|
10 | * Function implementations for the class vector.
|
---|
11 | *
|
---|
12 | */
|
---|
13 |
|
---|
14 | // include config.h
|
---|
15 | #ifdef HAVE_CONFIG_H
|
---|
16 | #include <config.h>
|
---|
17 | #endif
|
---|
18 |
|
---|
19 | #include "CodePatterns/MemDebug.hpp"
|
---|
20 |
|
---|
21 | #include "Exceptions/MathException.hpp"
|
---|
22 | #include "LinearAlgebra/Vector.hpp"
|
---|
23 | #include "LinearAlgebra/VectorContent.hpp"
|
---|
24 | #include "CodePatterns/Assert.hpp"
|
---|
25 | #include "Helpers/defs.hpp"
|
---|
26 | #include "Helpers/fast_functions.hpp"
|
---|
27 | #include "CodePatterns/Verbose.hpp"
|
---|
28 |
|
---|
29 | #include <cmath>
|
---|
30 | #include <iostream>
|
---|
31 | #include <gsl/gsl_blas.h>
|
---|
32 | #include <gsl/gsl_vector.h>
|
---|
33 |
|
---|
34 |
|
---|
35 | using namespace std;
|
---|
36 |
|
---|
37 |
|
---|
38 | /************************************ Functions for class vector ************************************/
|
---|
39 |
|
---|
40 | /** Constructor of class vector.
|
---|
41 | */
|
---|
42 | Vector::Vector()
|
---|
43 | {
|
---|
44 | content = new VectorContent((size_t) NDIM);
|
---|
45 | };
|
---|
46 |
|
---|
47 | /** Copy constructor.
|
---|
48 | * \param &src source Vector reference
|
---|
49 | */
|
---|
50 | Vector::Vector(const Vector& src)
|
---|
51 | {
|
---|
52 | content = new VectorContent(*(src.content));
|
---|
53 | }
|
---|
54 |
|
---|
55 | /** Constructor of class vector.
|
---|
56 | * \param x1 first component
|
---|
57 | * \param x2 second component
|
---|
58 | * \param x3 third component
|
---|
59 | */
|
---|
60 | Vector::Vector(const double x1, const double x2, const double x3)
|
---|
61 | {
|
---|
62 | content = new VectorContent((size_t) NDIM);
|
---|
63 | content->at(0) = x1;
|
---|
64 | content->at(1) = x2;
|
---|
65 | content->at(2) = x3;
|
---|
66 | };
|
---|
67 |
|
---|
68 | /** Constructor of class vector.
|
---|
69 | * \param x[3] three values to initialize Vector with
|
---|
70 | */
|
---|
71 | Vector::Vector(const double x[3])
|
---|
72 | {
|
---|
73 | content = new VectorContent((size_t) NDIM);
|
---|
74 | for (size_t i = NDIM; i--; )
|
---|
75 | content->at(i) = x[i];
|
---|
76 | };
|
---|
77 |
|
---|
78 | /** Copy constructor of class vector from VectorContent.
|
---|
79 | * \note This is destructive, i.e. we take over _content.
|
---|
80 | */
|
---|
81 | Vector::Vector(VectorContent *&_content) :
|
---|
82 | content(_content)
|
---|
83 | {
|
---|
84 | _content = NULL;
|
---|
85 | }
|
---|
86 |
|
---|
87 | /** Copy constructor of class vector from VectorContent.
|
---|
88 | * \note This is non-destructive, i.e. _content is copied.
|
---|
89 | */
|
---|
90 | Vector::Vector(VectorContent &_content)
|
---|
91 | {
|
---|
92 | content = new VectorContent(_content);
|
---|
93 | }
|
---|
94 |
|
---|
95 | /** Assignment operator.
|
---|
96 | * \param &src source vector to assign \a *this to
|
---|
97 | * \return reference to \a *this
|
---|
98 | */
|
---|
99 | Vector& Vector::operator=(const Vector& src){
|
---|
100 | // check for self assignment
|
---|
101 | if(&src!=this){
|
---|
102 | *content = *(src.content);
|
---|
103 | }
|
---|
104 | return *this;
|
---|
105 | }
|
---|
106 |
|
---|
107 | /** Desctructor of class vector.
|
---|
108 | * Vector::content is deleted.
|
---|
109 | */
|
---|
110 | Vector::~Vector() {
|
---|
111 | delete content;
|
---|
112 | };
|
---|
113 |
|
---|
114 | /** Calculates square of distance between this and another vector.
|
---|
115 | * \param *y array to second vector
|
---|
116 | * \return \f$| x - y |^2\f$
|
---|
117 | */
|
---|
118 | double Vector::DistanceSquared(const Vector &y) const
|
---|
119 | {
|
---|
120 | double res = 0.;
|
---|
121 | for (int i=NDIM;i--;)
|
---|
122 | res += (at(i)-y[i])*(at(i)-y[i]);
|
---|
123 | return (res);
|
---|
124 | };
|
---|
125 |
|
---|
126 | /** Calculates distance between this and another vector.
|
---|
127 | * \param *y array to second vector
|
---|
128 | * \return \f$| x - y |\f$
|
---|
129 | */
|
---|
130 | double Vector::distance(const Vector &y) const
|
---|
131 | {
|
---|
132 | return (sqrt(DistanceSquared(y)));
|
---|
133 | };
|
---|
134 |
|
---|
135 | size_t Vector::GreatestComponent() const
|
---|
136 | {
|
---|
137 | int greatest = 0;
|
---|
138 | for (int i=1;i<NDIM;i++) {
|
---|
139 | if (at(i) > at(greatest))
|
---|
140 | greatest = i;
|
---|
141 | }
|
---|
142 | return greatest;
|
---|
143 | }
|
---|
144 |
|
---|
145 | size_t Vector::SmallestComponent() const
|
---|
146 | {
|
---|
147 | int smallest = 0;
|
---|
148 | for (int i=1;i<NDIM;i++) {
|
---|
149 | if (at(i) < at(smallest))
|
---|
150 | smallest = i;
|
---|
151 | }
|
---|
152 | return smallest;
|
---|
153 | }
|
---|
154 |
|
---|
155 |
|
---|
156 | Vector Vector::getClosestPoint(const Vector &point) const{
|
---|
157 | // the closest point to a single point space is always the single point itself
|
---|
158 | return *this;
|
---|
159 | }
|
---|
160 |
|
---|
161 | /** Calculates scalar product between this and another vector.
|
---|
162 | * \param *y array to second vector
|
---|
163 | * \return \f$\langle x, y \rangle\f$
|
---|
164 | */
|
---|
165 | double Vector::ScalarProduct(const Vector &y) const
|
---|
166 | {
|
---|
167 | double res = 0.;
|
---|
168 | gsl_blas_ddot(content->content, y.content->content, &res);
|
---|
169 | return (res);
|
---|
170 | };
|
---|
171 |
|
---|
172 |
|
---|
173 | /** Calculates VectorProduct between this and another vector.
|
---|
174 | * -# returns the Product in place of vector from which it was initiated
|
---|
175 | * -# ATTENTION: Only three dim.
|
---|
176 | * \param *y array to vector with which to calculate crossproduct
|
---|
177 | * \return \f$ x \times y \f&
|
---|
178 | */
|
---|
179 | void Vector::VectorProduct(const Vector &y)
|
---|
180 | {
|
---|
181 | Vector tmp;
|
---|
182 | for(int i=NDIM;i--;)
|
---|
183 | tmp[i] = at((i+1)%NDIM)*y[(i+2)%NDIM] - at((i+2)%NDIM)*y[(i+1)%NDIM];
|
---|
184 | (*this) = tmp;
|
---|
185 | };
|
---|
186 |
|
---|
187 |
|
---|
188 | /** projects this vector onto plane defined by \a *y.
|
---|
189 | * \param *y normal vector of plane
|
---|
190 | * \return \f$\langle x, y \rangle\f$
|
---|
191 | */
|
---|
192 | void Vector::ProjectOntoPlane(const Vector &y)
|
---|
193 | {
|
---|
194 | Vector tmp;
|
---|
195 | tmp = y;
|
---|
196 | tmp.Normalize();
|
---|
197 | tmp.Scale(ScalarProduct(tmp));
|
---|
198 | *this -= tmp;
|
---|
199 | };
|
---|
200 |
|
---|
201 | /** Calculates the minimum distance of this vector to the plane.
|
---|
202 | * \sa Vector::GetDistanceVectorToPlane()
|
---|
203 | * \param *out output stream for debugging
|
---|
204 | * \param *PlaneNormal normal of plane
|
---|
205 | * \param *PlaneOffset offset of plane
|
---|
206 | * \return distance to plane
|
---|
207 | */
|
---|
208 | double Vector::DistanceToSpace(const Space &space) const
|
---|
209 | {
|
---|
210 | return space.distance(*this);
|
---|
211 | };
|
---|
212 |
|
---|
213 | /** Calculates the projection of a vector onto another \a *y.
|
---|
214 | * \param *y array to second vector
|
---|
215 | */
|
---|
216 | void Vector::ProjectIt(const Vector &y)
|
---|
217 | {
|
---|
218 | (*this) += (-ScalarProduct(y))*y;
|
---|
219 | };
|
---|
220 |
|
---|
221 | /** Calculates the projection of a vector onto another \a *y.
|
---|
222 | * \param *y array to second vector
|
---|
223 | * \return Vector
|
---|
224 | */
|
---|
225 | Vector Vector::Projection(const Vector &y) const
|
---|
226 | {
|
---|
227 | Vector helper = y;
|
---|
228 | helper.Scale((ScalarProduct(y)/y.NormSquared()));
|
---|
229 |
|
---|
230 | return helper;
|
---|
231 | };
|
---|
232 |
|
---|
233 | /** Calculates norm of this vector.
|
---|
234 | * \return \f$|x|\f$
|
---|
235 | */
|
---|
236 | double Vector::Norm() const
|
---|
237 | {
|
---|
238 | return (content->Norm());
|
---|
239 | };
|
---|
240 |
|
---|
241 | /** Calculates squared norm of this vector.
|
---|
242 | * \return \f$|x|^2\f$
|
---|
243 | */
|
---|
244 | double Vector::NormSquared() const
|
---|
245 | {
|
---|
246 | return (content->NormSquared());
|
---|
247 | };
|
---|
248 |
|
---|
249 | /** Normalizes this vector.
|
---|
250 | */
|
---|
251 | void Vector::Normalize()
|
---|
252 | {
|
---|
253 | content->Normalize();
|
---|
254 | };
|
---|
255 |
|
---|
256 | Vector Vector::getNormalized() const{
|
---|
257 | Vector res= *this;
|
---|
258 | res.Normalize();
|
---|
259 | return res;
|
---|
260 | }
|
---|
261 |
|
---|
262 | /** Zeros all components of this vector.
|
---|
263 | */
|
---|
264 | void Vector::Zero()
|
---|
265 | {
|
---|
266 | at(0)=at(1)=at(2)=0;
|
---|
267 | };
|
---|
268 |
|
---|
269 | /** Zeros all components of this vector.
|
---|
270 | */
|
---|
271 | void Vector::One(const double one)
|
---|
272 | {
|
---|
273 | at(0)=at(1)=at(2)=one;
|
---|
274 | };
|
---|
275 |
|
---|
276 | /** Checks whether vector has all components zero.
|
---|
277 | * @return true - vector is zero, false - vector is not
|
---|
278 | */
|
---|
279 | bool Vector::IsZero() const
|
---|
280 | {
|
---|
281 | return (fabs(at(0))+fabs(at(1))+fabs(at(2)) < MYEPSILON);
|
---|
282 | };
|
---|
283 |
|
---|
284 | /** Checks whether vector has length of 1.
|
---|
285 | * @return true - vector is normalized, false - vector is not
|
---|
286 | */
|
---|
287 | bool Vector::IsOne() const
|
---|
288 | {
|
---|
289 | return (fabs(Norm() - 1.) < MYEPSILON);
|
---|
290 | };
|
---|
291 |
|
---|
292 | /** Checks whether vector is normal to \a *normal.
|
---|
293 | * @return true - vector is normalized, false - vector is not
|
---|
294 | */
|
---|
295 | bool Vector::IsNormalTo(const Vector &normal) const
|
---|
296 | {
|
---|
297 | if (ScalarProduct(normal) < MYEPSILON)
|
---|
298 | return true;
|
---|
299 | else
|
---|
300 | return false;
|
---|
301 | };
|
---|
302 |
|
---|
303 | /** Checks whether vector is normal to \a *normal.
|
---|
304 | * @return true - vector is normalized, false - vector is not
|
---|
305 | */
|
---|
306 | bool Vector::IsEqualTo(const Vector &a) const
|
---|
307 | {
|
---|
308 | bool status = true;
|
---|
309 | for (int i=0;i<NDIM;i++) {
|
---|
310 | if (fabs(at(i) - a[i]) > MYEPSILON)
|
---|
311 | status = false;
|
---|
312 | }
|
---|
313 | return status;
|
---|
314 | };
|
---|
315 |
|
---|
316 | /** Calculates the angle between this and another vector.
|
---|
317 | * \param *y array to second vector
|
---|
318 | * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
|
---|
319 | */
|
---|
320 | double Vector::Angle(const Vector &y) const
|
---|
321 | {
|
---|
322 | double norm1 = Norm(), norm2 = y.Norm();
|
---|
323 | double angle = -1;
|
---|
324 | if ((fabs(norm1) > MYEPSILON) && (fabs(norm2) > MYEPSILON))
|
---|
325 | angle = this->ScalarProduct(y)/norm1/norm2;
|
---|
326 | // -1-MYEPSILON occured due to numerical imprecision, catch ...
|
---|
327 | //Log() << Verbose(2) << "INFO: acos(-1) = " << acos(-1) << ", acos(-1+MYEPSILON) = " << acos(-1+MYEPSILON) << ", acos(-1-MYEPSILON) = " << acos(-1-MYEPSILON) << "." << endl;
|
---|
328 | if (angle < -1)
|
---|
329 | angle = -1;
|
---|
330 | if (angle > 1)
|
---|
331 | angle = 1;
|
---|
332 | return acos(angle);
|
---|
333 | };
|
---|
334 |
|
---|
335 |
|
---|
336 | double& Vector::operator[](size_t i){
|
---|
337 | ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
|
---|
338 | return *gsl_vector_ptr (content->content, i);
|
---|
339 | }
|
---|
340 |
|
---|
341 | const double& Vector::operator[](size_t i) const{
|
---|
342 | ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
|
---|
343 | return *gsl_vector_ptr (content->content, i);
|
---|
344 | }
|
---|
345 |
|
---|
346 | double& Vector::at(size_t i){
|
---|
347 | return (*this)[i];
|
---|
348 | }
|
---|
349 |
|
---|
350 | const double& Vector::at(size_t i) const{
|
---|
351 | return (*this)[i];
|
---|
352 | }
|
---|
353 |
|
---|
354 | VectorContent* Vector::get() const
|
---|
355 | {
|
---|
356 | return content;
|
---|
357 | }
|
---|
358 |
|
---|
359 | /** Compares vector \a to vector \a b component-wise.
|
---|
360 | * \param a base vector
|
---|
361 | * \param b vector components to add
|
---|
362 | * \return a == b
|
---|
363 | */
|
---|
364 | bool Vector::operator==(const Vector& b) const
|
---|
365 | {
|
---|
366 | return IsEqualTo(b);
|
---|
367 | };
|
---|
368 |
|
---|
369 | bool Vector::operator!=(const Vector& b) const
|
---|
370 | {
|
---|
371 | return !IsEqualTo(b);
|
---|
372 | }
|
---|
373 |
|
---|
374 | /** Sums vector \a to this lhs component-wise.
|
---|
375 | * \param a base vector
|
---|
376 | * \param b vector components to add
|
---|
377 | * \return lhs + a
|
---|
378 | */
|
---|
379 | const Vector& Vector::operator+=(const Vector& b)
|
---|
380 | {
|
---|
381 | this->AddVector(b);
|
---|
382 | return *this;
|
---|
383 | };
|
---|
384 |
|
---|
385 | /** Subtracts vector \a from this lhs component-wise.
|
---|
386 | * \param a base vector
|
---|
387 | * \param b vector components to add
|
---|
388 | * \return lhs - a
|
---|
389 | */
|
---|
390 | const Vector& Vector::operator-=(const Vector& b)
|
---|
391 | {
|
---|
392 | this->SubtractVector(b);
|
---|
393 | return *this;
|
---|
394 | };
|
---|
395 |
|
---|
396 | /** factor each component of \a *this times \a m.
|
---|
397 | * \param m factor
|
---|
398 | * \return \f$(\text{*this} \cdot m\f$
|
---|
399 | */
|
---|
400 | const Vector& Vector::operator*=(const double m)
|
---|
401 | {
|
---|
402 | Scale(m);
|
---|
403 | return *this;
|
---|
404 | };
|
---|
405 |
|
---|
406 | /** Sums two vectors \a and \b component-wise.
|
---|
407 | * \param a first vector
|
---|
408 | * \param b second vector
|
---|
409 | * \return a + b
|
---|
410 | */
|
---|
411 | Vector const Vector::operator+(const Vector& b) const
|
---|
412 | {
|
---|
413 | Vector x = *this;
|
---|
414 | x.AddVector(b);
|
---|
415 | return x;
|
---|
416 | };
|
---|
417 |
|
---|
418 | /** Subtracts vector \a from \b component-wise.
|
---|
419 | * \param a first vector
|
---|
420 | * \param b second vector
|
---|
421 | * \return a - b
|
---|
422 | */
|
---|
423 | Vector const Vector::operator-(const Vector& b) const
|
---|
424 | {
|
---|
425 | Vector x = *this;
|
---|
426 | x.SubtractVector(b);
|
---|
427 | return x;
|
---|
428 | };
|
---|
429 |
|
---|
430 | /** Factors given vector \a *this times \a m.
|
---|
431 | * \param m factor
|
---|
432 | * \return \f$(\text{*this} \cdot m)\f$
|
---|
433 | */
|
---|
434 | const Vector Vector::operator*(const double m) const
|
---|
435 | {
|
---|
436 | Vector x(*this);
|
---|
437 | x.Scale(m);
|
---|
438 | return x;
|
---|
439 | };
|
---|
440 |
|
---|
441 | /** Factors given vector \a a times \a m.
|
---|
442 | * \param m factor
|
---|
443 | * \param a vector
|
---|
444 | * \return m * a
|
---|
445 | */
|
---|
446 | Vector const operator*(const double m, const Vector& a )
|
---|
447 | {
|
---|
448 | Vector x(a);
|
---|
449 | x.Scale(m);
|
---|
450 | return x;
|
---|
451 | };
|
---|
452 |
|
---|
453 | ostream& operator<<(ostream& ost, const Vector& m)
|
---|
454 | {
|
---|
455 | ost << "(";
|
---|
456 | for (int i=0;i<NDIM;i++) {
|
---|
457 | ost << m[i];
|
---|
458 | if (i != 2)
|
---|
459 | ost << ",";
|
---|
460 | }
|
---|
461 | ost << ")";
|
---|
462 | return ost;
|
---|
463 | };
|
---|
464 |
|
---|
465 |
|
---|
466 | void Vector::ScaleAll(const double *factor)
|
---|
467 | {
|
---|
468 | for (int i=NDIM;i--;)
|
---|
469 | at(i) *= factor[i];
|
---|
470 | };
|
---|
471 |
|
---|
472 | void Vector::ScaleAll(const Vector &factor){
|
---|
473 | gsl_vector_mul(content->content, factor.content->content);
|
---|
474 | }
|
---|
475 |
|
---|
476 |
|
---|
477 | void Vector::Scale(const double factor)
|
---|
478 | {
|
---|
479 | gsl_vector_scale(content->content,factor);
|
---|
480 | };
|
---|
481 |
|
---|
482 | std::pair<Vector,Vector> Vector::partition(const Vector &rhs) const{
|
---|
483 | double factor = ScalarProduct(rhs)/rhs.NormSquared();
|
---|
484 | Vector res= factor * rhs;
|
---|
485 | return make_pair(res,(*this)-res);
|
---|
486 | }
|
---|
487 |
|
---|
488 | std::pair<pointset,Vector> Vector::partition(const pointset &points) const{
|
---|
489 | Vector helper = *this;
|
---|
490 | pointset res;
|
---|
491 | for(pointset::const_iterator iter=points.begin();iter!=points.end();++iter){
|
---|
492 | pair<Vector,Vector> currPart = helper.partition(*iter);
|
---|
493 | res.push_back(currPart.first);
|
---|
494 | helper = currPart.second;
|
---|
495 | }
|
---|
496 | return make_pair(res,helper);
|
---|
497 | }
|
---|
498 |
|
---|
499 | /** Creates this vector as the b y *factors' components scaled linear combination of the given three.
|
---|
500 | * this vector = x1*factors[0] + x2* factors[1] + x3*factors[2]
|
---|
501 | * \param *x1 first vector
|
---|
502 | * \param *x2 second vector
|
---|
503 | * \param *x3 third vector
|
---|
504 | * \param *factors three-component vector with the factor for each given vector
|
---|
505 | */
|
---|
506 | void Vector::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
|
---|
507 | {
|
---|
508 | (*this) = (factors[0]*x1) +
|
---|
509 | (factors[1]*x2) +
|
---|
510 | (factors[2]*x3);
|
---|
511 | };
|
---|
512 |
|
---|
513 | /** Calculates orthonormal vector to one given vectors.
|
---|
514 | * Just subtracts the projection onto the given vector from this vector.
|
---|
515 | * The removed part of the vector is Vector::Projection()
|
---|
516 | * \param *x1 vector
|
---|
517 | * \return true - success, false - vector is zero
|
---|
518 | */
|
---|
519 | bool Vector::MakeNormalTo(const Vector &y1)
|
---|
520 | {
|
---|
521 | bool result = false;
|
---|
522 | double factor = y1.ScalarProduct(*this)/y1.NormSquared();
|
---|
523 | Vector x1 = factor * y1;
|
---|
524 | SubtractVector(x1);
|
---|
525 | for (int i=NDIM;i--;)
|
---|
526 | result = result || (fabs(at(i)) > MYEPSILON);
|
---|
527 |
|
---|
528 | return result;
|
---|
529 | };
|
---|
530 |
|
---|
531 | /** Creates this vector as one of the possible orthonormal ones to the given one.
|
---|
532 | * Just scan how many components of given *vector are unequal to zero and
|
---|
533 | * try to get the skp of both to be zero accordingly.
|
---|
534 | * \param *vector given vector
|
---|
535 | * \return true - success, false - failure (null vector given)
|
---|
536 | */
|
---|
537 | bool Vector::GetOneNormalVector(const Vector &GivenVector)
|
---|
538 | {
|
---|
539 | int Components[NDIM]; // contains indices of non-zero components
|
---|
540 | int Last = 0; // count the number of non-zero entries in vector
|
---|
541 | int j; // loop variables
|
---|
542 | double norm;
|
---|
543 |
|
---|
544 | for (j=NDIM;j--;)
|
---|
545 | Components[j] = -1;
|
---|
546 |
|
---|
547 | // in two component-systems we need to find the one position that is zero
|
---|
548 | int zeroPos = -1;
|
---|
549 | // find two components != 0
|
---|
550 | for (j=0;j<NDIM;j++){
|
---|
551 | if (fabs(GivenVector[j]) > MYEPSILON)
|
---|
552 | Components[Last++] = j;
|
---|
553 | else
|
---|
554 | // this our zero Position
|
---|
555 | zeroPos = j;
|
---|
556 | }
|
---|
557 |
|
---|
558 | switch(Last) {
|
---|
559 | case 3: // threecomponent system
|
---|
560 | // the position of the zero is arbitrary in three component systems
|
---|
561 | zeroPos = Components[2];
|
---|
562 | case 2: // two component system
|
---|
563 | norm = sqrt(1./(GivenVector[Components[1]]*GivenVector[Components[1]]) + 1./(GivenVector[Components[0]]*GivenVector[Components[0]]));
|
---|
564 | at(zeroPos) = 0.;
|
---|
565 | // in skp both remaining parts shall become zero but with opposite sign and third is zero
|
---|
566 | at(Components[1]) = -1./GivenVector[Components[1]] / norm;
|
---|
567 | at(Components[0]) = 1./GivenVector[Components[0]] / norm;
|
---|
568 | return true;
|
---|
569 | break;
|
---|
570 | case 1: // one component system
|
---|
571 | // set sole non-zero component to 0, and one of the other zero component pendants to 1
|
---|
572 | at((Components[0]+2)%NDIM) = 0.;
|
---|
573 | at((Components[0]+1)%NDIM) = 1.;
|
---|
574 | at(Components[0]) = 0.;
|
---|
575 | return true;
|
---|
576 | break;
|
---|
577 | default:
|
---|
578 | return false;
|
---|
579 | }
|
---|
580 | };
|
---|
581 |
|
---|
582 | /** Adds vector \a *y componentwise.
|
---|
583 | * \param *y vector
|
---|
584 | */
|
---|
585 | void Vector::AddVector(const Vector &y)
|
---|
586 | {
|
---|
587 | gsl_vector_add(content->content, y.content->content);
|
---|
588 | }
|
---|
589 |
|
---|
590 | /** Adds vector \a *y componentwise.
|
---|
591 | * \param *y vector
|
---|
592 | */
|
---|
593 | void Vector::SubtractVector(const Vector &y)
|
---|
594 | {
|
---|
595 | gsl_vector_sub(content->content, y.content->content);
|
---|
596 | }
|
---|
597 |
|
---|
598 |
|
---|
599 | // some comonly used vectors
|
---|
600 | const Vector zeroVec(0,0,0);
|
---|
601 | const Vector unitVec[NDIM]={Vector(1,0,0),Vector(0,1,0),Vector(0,0,1)};
|
---|