Allink  v0.1
MatematicaVect.cpp
1 #include "../include/Matematica.h"
3  x = NULL;
4  if(N<0)return;
5  NDim = N;
6  //x = (double *)calloc(NDim,sizeof(double));
7  x = new double[NDim];
8 }
10  x = NULL;
11  NDim = 3;
12  //x = (double *)calloc(NDim,sizeof(double));
13  x = new double[NDim];
14 }
15 Vettore::Vettore(double *Pos,int N){
16  x = NULL;
17  NDim = N;
18  x = new double[NDim];
19  // x = (double *)calloc(NDim,sizeof(double));
20  for(int d=0;d<3;d++)
21  x[d] = Pos[d];
22 };
23 Vettore::Vettore(double xx,double yy){
24  x = NULL;
25  NDim = 2;
26  x = new double[NDim];
27  // x = (double *)calloc(NDim,sizeof(double));
28  x[0] = xx;
29  x[1] = yy;
30 }
31 Vettore::Vettore(double xx,double yy,double zz){
32  x = NULL;
33  NDim = 3;
34  x = new double[NDim];
35  // x = (double *)calloc(NDim,sizeof(double));
36  x[0] = xx;
37  x[1] = yy;
38  x[2] = zz;
39 }
41  if(x) delete [] x;
42  //if(x) free(x);
43 }
45  assert(NDim == u.NDim);
46  Vettore Resp(this->NDim);
47  for(int d=0;d<this->NDim;d++){
48  Resp.x[d] = this->x[d] + u.x[d];
49  }
50  return Resp;
51 }
53  assert(NDim == u.NDim);
54  Vettore Resp(this->NDim);
55  for(int d=0;d<this->NDim;d++){
56  Resp.x[d] = this->x[d] + u.x[d];
57  }
58  return Resp;
59 }
61  for(int d=0;d<this->NDim;d++){
62  this->x[d] += u.x[d];
63  }
64  return *this;
65 }
67  Vettore Resp(this->NDim);
68  for(int d=0;d<this->NDim;d++){
69  Resp.x[d] = this->x[d] - u.x[d];
70  }
71  return Resp;
72 }
74  Vettore Resp(this->NDim);
75  for(int d=0;d<this->NDim;d++){
76  Resp.x[d] = this->x[d] - u.x[d];
77  }
78  return Resp;
79 }
81  for(int d=0;d<this->NDim;d++){
82  this->x[d] -= u.x[d];
83  }
84  return *this;
85 }
87  assert(NDim == u.NDim);
88  Vettore Resp(this->NDim);
89  for(int d=0;d<this->NDim;d++){
90  Resp.x[d] = this->x[d] * u.x[d];
91  }
92  return Resp;
93 }
95  for(int d=0;d<this->NDim;d++){
96  this->x[d] *= u.x[d];
97  }
98  return *this;
99 }
100 Vettore Vettore::operator*(const double Fact){
101  Vettore Resp(this->NDim);
102  for(int d=0;d<this->NDim;d++){
103  Resp.x[d] = this->x[d]*Fact;
104  }
105  return Resp;
106 }
107 Vettore& Vettore::operator*=(const double Fact){
108  Vettore Resp(this->NDim);
109  for(int d=0;d<this->NDim;d++){
110  this->x[d] *= Fact;
111  }
112  return *this;
113 }
114 //inline
115 Vettore operator*(double Fact, const Vettore& vec) {
116  Vettore Resp(vec.NDim);
117  for(int d=0;d<vec.NDim;d++){
118  Resp.x[d] = vec.x[d]*Fact;
119  }
120  return Resp;
121 }
123  assert(NDim == u.NDim);
124  Vettore Resp(this->NDim);
125  for(int d=0;d<this->NDim;d++){
126  Resp.x[d] = this->x[d] / u.x[d];
127  }
128  return Resp;
129 }
131  for(int d=0;d<this->NDim;d++){
132  this->x[d] /= u.x[d];
133  }
134  return *this;
135 }
136 double Vettore::operator%(const Vettore &u){
137  assert(NDim == u.NDim);
138  double Resp = 0.;
139  for(int d=0;d<this->NDim;d++){
140  Resp += this->x[d] * u.x[d];
141  }
142  return Resp;
143 }
145  Vettore Resp(this->NDim);
146  for(int d=0;d<this->NDim;d++){
147  Resp.x[d] = u.x[d];
148  }
149  return Resp;
150 }
152  assert(NDim == u.NDim);
153  Vettore Resp(NDim);
154  for(int d=0;d<NDim;d++){
155  int NUno = (d+1)%NDim;
156  int NDue = (d+2)%NDim;
157  Resp.x[d] = u.x[NUno] * x[NDue] - u.x[NDue]*x[NUno];
158  }
159  return Resp;
160 }
162  for(int d=0;d<NDim;d++){
163  int NUno = (d+1)%NDim;
164  int NDue = (d+2)%NDim;
165  this->x[d] = this->x[NUno] * u.x[NDue] - this->x[NDue]*u.x[NUno];
166  }
167  return *this;
168 }
170  for(int d=0;d<NDim;d++)
171  printf("%d) %lf\n",d,x[d]);
172 }
173 double Vettore::Abs(){
174  double Resp = 0.;
175  for(int d=0;d<NDim;d++)
176  Resp += x[d];
177  return ABS(Resp);
178 }
179 double Vettore::Norm(){
180  double Resp=0.;
181  for(int d=0;d<NDim;d++)
182  Resp += QUAD((x[d]));
183  return sqrt(Resp);
184 }
186  double Div = Norm();
187  if(Div > 0.0)
188  for(int d=0;d<NDim;d++)
189  x[d] /= Div;
190  return Div;
191 }
192 void Vettore::Subs(const Vettore *u,const Vettore *v){
193  for(int d=0;d<NDim;d++)
194  x[d] = u->x[d] - v->x[d];
195 }
196 void Vettore::Mult(double Fact){
197  for(int d=0;d<NDim;d++)
198  x[d] = Fact*x[d];
199 }
200 double Vettore::ScalS(const Vettore *u,const Vettore *v){
201 #ifdef VECT_DEBUG
202  if(NDim != u->NDim || NDim != v->NDim){
203  printf("Incompatible vectors. Dim: %d %d\n",NDim,u->NDim);
204  return 0.;
205  }
206 #endif
207  double Resp=0.;
208  for(int d=0;d<NDim;d++)
209  Resp += u->x[d] * v->x[d];
210  return Resp;
211 }
212 double Vettore::Col(int N){
213 #ifdef VECT_DEBUG
214  if( N >= NDim || N < 0) return 0.;
215 #endif
216  return x[N];
217 }
218 void Vettore::Set(double Val,int N){
219 #ifdef VECT_DEBUG
220  if( N >= NDim || N < 0) return 0.;
221 #endif
222  x[N] = Val;
223 }
225  double Resp = ScalS(u,v);
226  Resp /= v->Norm()*u->Norm();
227  return Resp;
228 }
230 #ifdef VECT_DEBUG
231 #endif
232  double Resp = 0.;
233  for(int d=0;d<NDim;d++)
234  Resp += x[d] * u->x[d];
235  Resp /= Norm();
236  Resp /= u->Norm();
237  return Resp;
238 }
240 #ifdef VECT_DEBUG
241 #endif
242  Vettore w(u->NDim);
243  w.VetV(u,v);
244  double Resp = w.Norm()/(u->Norm()*v->Norm());
245  return Resp;
246 }
248 #ifdef VECT_DEBUG
249 #endif
250  double Resp = 0.;
251  for(int d=0;d<NDim;d++){
252  int NUno = (d+1)%NDim;
253  int NDue = (d+2)%NDim;
254  Resp += SQR(x[NUno] * u->x[NDue] - x[NDue]*u->x[NUno]);
255  }
256  Resp = sqrt(Resp)/(Norm()*u->Norm());
257  return Resp;
258 }
260  return acos(CosAngle(u,v));
261 }
263  return acos(CosAngle(u));
264 }
266 #ifdef VECT_DEBUG
267  if(NDim != u->NDim || NDim != v->NDim){
268  printf("Incompatible vectors. Dim: %d %d\n",NDim,u->NDim);
269  return ;
270  }
271 #endif
272  for(int d=0;d<NDim;d++) x[d] = .5*(u->x[d] + v->x[d]);
273 }
274 void Vettore::Normal(const Vettore *u,const Vettore *v){
275 #ifdef VECT_DEBUG
276  if(v->NDim != u->NDim){
277  printf("Incompatible vectors. Dim: %d %d\n",u->NDim,v->NDim);
278  return ;
279  }
280 #endif
281  VetV(u,v);
282  Normalize();
283 }
284 void Vettore::NormalSurf(const Vettore *u,const Vettore *v,const Vettore *w){
285 #ifdef VECT_DEBUG
286  if(v->NDim != u->NDim){
287  printf("Incompatible vectors. Dim: %d %d\n",u->NDim,v->NDim);
288  return ;
289  }
290 #endif
291  for(int d=0;d<NDim;d++){
292  int NUno = (d+1)%NDim;
293  int NDue = (d+2)%NDim;
294  double dUno = (u->x[NUno] - w->x[NUno])*(v->x[NDue] - w->x[NDue]);
295  double dDue = (u->x[NDue] - w->x[NDue])*(v->x[NUno] - w->x[NUno]);
296  x[d] = dUno - dDue;
297  //printf("%d -(%d %d)/%d %lf\n",d,NUno,NDue,NDim,x[d]);
298  }
299 }
300 void Vettore::ScalV(const Vettore *u,const Vettore *v){
301 #ifdef VECT_DEBUG
302  if(NDim != u->NDim || NDim != v->NDim){
303  printf("Incompatible vectors. Dim: %d %d\n",NDim,u->NDim);
304  return ;
305  }
306 #endif
307  for(int d=0;d<NDim;d++)
308  x[d] = u->x[d] * v->x[d];
309 }
310 double Vettore::VetV(const Vettore *u,const Vettore *v){
311  assert(NDim == u->NDim);
312  double Area = 0.;
313  for(int d=0;d<NDim;d++){
314  int NUno = (d+1)%NDim;
315  int NDue = (d+2)%NDim;
316  x[d] = u->x[NUno] * v->x[NDue] - u->x[NDue]*v->x[NUno];
317  Area += x[d];
318  }
319  return ABS(Area);
320 }
321 double Vettore::VetV(const Vettore *u){
322  assert(NDim == u->NDim);
323  double Area = 0.;
324  Vettore w(NDim);
325  for(int d=0;d<NDim;d++){
326  int NUno = (d+1)%NDim;
327  int NDue = (d+2)%NDim;
328  w.x[d] = x[NUno] * u->x[NDue] - x[NDue]*u->x[NUno];
329  Area += w.x[d];
330  }
331  for(int d=0;d<NDim;d++){
332  x[d] = w.x[d];
333  }
334  return ABS(Area);
335 }
336 double Vettore::VetV3(const Vettore *u,const Vettore *v){
337  double Area = 0.;
338  x[0] = u->x[1]*v->x[2] - u->x[2]*v->x[1];
339  x[1] = u->x[2]*v->x[0] - u->x[0]*v->x[2];
340  x[2] = u->x[0]*v->x[1] - u->x[1]*v->x[0];
341  for(int d=0;d<3;d++){
342  Area += x[d];
343  }
344  return ABS(Area);
345 }
346 void Vettore::Rescale(double Length){
347  double Fact = Length/Norm();
348  for(int d=0;d<NDim;d++)
349  x[d] *= Fact;
350 }
352  double Pre = 0.;
353  double Norma = 0.;
354  for(int d=0;d<NDim;d++){
355  Pre += x[d]*Axis->x[d];
356  Norma += SQR(Axis->x[d]);
357  }
358  double Length = 0.;
359  for(int d=0;d<NDim;d++)
360  Length += x[d] = Axis->x[d]*Pre/Norma;
361  return Length;
362 }
364  double Cos = CosAngle(Pos,Axis);
365  double InvNorma = 1./Axis->Norm();
366  return Cos*InvNorma;
367 }
369  for(int d=0;d<NDim;d++)
370  x[d] = o->x[d] - x[d];
371 }
373  for(int d=0;d<NDim;d++)
374  x[d] = c->x[d];
375 }
376 void Vettore::Export(double *xx){
377  for(int d=0;d<NDim;d++)
378  xx[d] = x[d];
379 }
381  double Sin = SinAngle(a);
382  double Norma = Norm();
383  Vettore w(NDim);
384  for(int d=0;d<NDim;d++){
385  int NUno = (d+1)%NDim;
386  int NDue = (d+2)%NDim;
387  w.x[d] = x[NUno]*a->x[NDue] - x[NDue]*a->x[NUno];
388  }
389  VetV(&w,a);
390  Normalize();
391  Rescale(Sin*Norma);
392 }
394  Vettore Perp1(Axis->NDim);
395  Vettore Perp2(Axis->NDim);
396  Perp1.VetV(Pos,Axis);
397  Perp2.VetV(&Perp1,Axis);
398  double Distance = Perp1.Norm()/Axis->Norm();
399  double NormaInv = 1./Perp2.Norm();
400  for(int d=0;d<Axis->NDim;d++)
401  x[d] = Perp2[d]*Distance*NormaInv;
402  return Distance;
403  // double Sin = SinAngle(Axis,Pos);
404  // double Norma = Norm();
405  // Vettore w(NDim);
406  // for(int d=0;d<NDim;d++){
407  // int NUno = (d+1)%NDim;
408  // int NDue = (d+2)%NDim;
409  // w.x[d] = Axis->x[NUno] * a.x[NDue] - x[NDue]*a->x[NUno];
410  // }
411  // VetV(w,a);
412  // Normalize();
413  // Rescale(Sin*Norma);
414 }
416  Vettore Perp1(3);
417  Vettore Perp2(3);
418  Perp1.VetV3(Pos,Axis);
419  Perp2.VetV3(&Perp1,Axis);
420  double Distance = Perp1.Norm()/Axis->Norm();
421  double NormaInv = 1./Perp2.Norm();
422  for(int d=0;d<3;d++)
423  x[d] = Perp2[d]*Distance*NormaInv;
424  return Distance;
425 }
427  double Known[3];
428  double UnKnown[3];
429  Vettore Normal(3);
430  Normal.NormalSurf(S1,S2,S3);
431  Normal.Normalize();
432  Matrice Mat(3,3);
433  for(int d=0;d<3;d++){
434  Known[d] = S2->Val(d) + P->Val(d);
435  Mat.Set(0,d,S1->Val(d)-S2->Val(d));
436  Mat.Set(1,d,S2->Val(d)-S3->Val(d));
437  Mat.Set(2,d,Normal.Val(d));
438  }
439  Mat.Solve(Known,UnKnown);
440  //double Norma = 0.;
441  for(int d=0;d<3;d++){
442  //Norma += SQR(UnKnown[d]);
443  x[d] = UnKnown[d];
444  }
445  return Norm();
446 }
void Export(double *x)
Export.
void Subs(const Vettore *u, const Vettore *v)
substruct two Vettore
void Copy(Vettore *o)
Copy the vector.
double operator%(const Vettore &Vet)
Scalar product.
int Solve(double *Known, double *UnKnown)
Solve a system A|b = y.
double Norm()
Return the norm of a Vettore.
void ApplyOn(Vettore *o)
Apply on a origin.
double Normalize()
Normlizes a Vettore.
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
~Vettore()
Frees.
Vettore operator^(const Vettore &vec)
Vectorial product.
Vettore operator=(const Vettore &Vet)
Assigns the entries of the rhs Vettore to the lhs.
void Print()
Prints the components.
double PerpTo3(Vettore *Pos, Vettore *Axis)
The vector perpendicolar in three dimension (faster)
void ScalV(const Vettore *u, const Vettore *v)
Scalar product between two Vettore.
Vettore operator*(const Vettore &Vet)
Moltiplication component by component.
double Val(int N)
Value of the N column.
double * x
Where the data are stored.
void Axis(Vettore *u, Vettore *v)
Calculates the axis formed by two Vettore.
double CosAngle(Vettore *u)
Computes the cosine with respect to.
double ScalS(const Vettore *u, const Vettore *v)
Multiplies the components of a Vettore for a scalar.
void Rescale(double NewLength)
Rescale the total length of the vector.
Vettore & operator-=(const Vettore &Vet)
Difference component by component.
double Col(int N)
Value of the N column.
Vettore()
Allocates.
Vettore operator-(const Vettore &Vet)
Difference component by component.
Vettore & operator+=(const Vettore &Vet)
Sum component by component.
double ProjOnAxis(Vettore *a)
Projects along the axis.
Vettore operator/(const Vettore &vec)
Division component by component.
double VetV(const Vettore *u, const Vettore *v)
Vectorial product between two Vettore returns the area.
Matrice computes the algebric operations on matrices.
bool Set(int row, int column, double Val)
Set a coefficient.
Vettore & operator*=(const Vettore &vec)
Scalar product.
int NDim
Dimension allocated.
double Abs()
Return the absolute value.
double Angle(Vettore *u, Vettore *v)
Computes the angle between two Vetttore.
Vettore & operator/=(const Vettore &vec)
Division component by component.
double ProjOnSurf(Vettore *S1, Vettore *S2, Vettore *S3, Vettore *P)
Project a point P on the point PS perpendicular to the surface described by the points S1...
void Mult(double Fact)
multiply by a scalar
double VetV3(const Vettore *u, const Vettore *v)
Vectorial product between two Vettore in three dimension (faster) returns the area.
void Normal(const Vettore *u, const Vettore *v)
Computes the normal with respect to the Vettore u and v.
Vettore operator+(const Vettore &Vet)
Sum component by component.
double SinAngle(Vettore *u, Vettore *v)
Computes the sine between two Vettore.
void Set(double Val, int Col)
Set the N column.
Vettore & operator^=(const Vettore &vec)
Vectorial product.
void NormalSurf(const Vettore *u, const Vettore *v, const Vettore *w)
Computes the normal to the plane described by u,v,w.
void PerpTo(Vettore *o)
The vector perpendicolar.