Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F72758627
vector_3.h
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Jul 16, 22:03
Size
18 KB
Mime Type
text/x-c++
Expires
Thu, Jul 18, 22:03 (2 d)
Engine
blob
Format
Raw Data
Handle
19092351
Attached To
rLAMMPS lammps
vector_3.h
View Options
This document is not UTF8. It was detected as ISO-8859-1 (Latin 1) and converted to UTF8 for display.
/*s***************************************************************************
*
* Copyright (c), Ilya Valuev 2005 All Rights Reserved.
*
* Author : Ilya Valuev, MIPT, Moscow, Russia
*
* Project : GridMD, ivutils
*
*****************************************************************************/
/*s****************************************************************************
* $Log: vector_3.h,v $
* Revision 1.2 2011/06/11 16:53:55 valuev
* sync with LAMMPS
*
* Revision 1.1 2011/06/10 17:15:07 morozov
* First Windows project with the correct directory structure
*
* Revision 1.22 2010/11/17 02:13:32 valuev
* Added analysis phase, fixed some memory leaks
*
* Revision 1.21 2009/09/09 14:43:35 valuev
* fixed trajreader
*
* Revision 1.20 2009/07/24 05:08:46 valuev
* Sync with FDTD, added molecule setup
*
* Revision 1.33 2009/06/01 13:01:42 valuev
* Added ShearBox
*
* Revision 1.32 2009/05/28 07:49:00 valuev
* updated chemosensor
*
* Revision 1.31 2009/05/22 08:53:52 valuev
* new uiExperiment interface
*
* Revision 1.30 2009/02/22 10:49:56 lesha
* Vector_Nt constructors are modified
*
* Revision 1.29 2009/02/04 14:23:31 valuev
* fixed bug in maxabscoord and minabscoord functions
*
* Revision 1.28 2009/02/04 10:56:30 lesha
* SINGLE_PRECISION is recovered
*
* Revision 1.27 2009/02/03 00:47:35 lesha
* *** empty log message ***
*
* Revision 1.26 2009/01/30 13:54:05 valuev
* restructured as a library
*
* Revision 1.16 2008/08/18 21:40:09 valuev
* added Gurski-Krasko potential
*
* Revision 1.15 2008/05/05 17:27:43 morozov
* cvector_3.h is the new header for class cVector_3. Old one is moved to cvector_3old.h
* class hmatrix is added to pairhash.h
* wavepackets.h contains class WavePacket
*
* Revision 1.14 2008/04/22 12:44:17 valuev
* made gcc 4.12 compilable
*
* Revision 1.13 2008/04/21 23:13:44 valuev
* made gcc 4.12 compilable
*
* Revision 1.12 2008/04/21 22:42:30 valuev
* *** empty log message ***
*
* Revision 1.11 2008/04/15 13:11:41 valuev
* Added antisymmetrized wave packets
*
*
*******************************************************************************/
/*r @file vector_3.h @brief ðàáîòà ñ òðåõìåðíûìè âåêòîðàìè
*/
# ifndef VECTOR_3_H
# define VECTOR_3_H
# define _USE_MATH_DEFINES
# include <iostream>
# include <cmath>
# include <limits>
# include <stdlib.h>
// some compilers don't define PI!
# ifndef M_PI
# define M_PI 3.1415926535897932385
# endif
# ifndef fmod
//r äåëåíèå ïî ìîäóëþ ÷èñåë ñ ïëàâàþùåé òî÷êîé
# define fmod(a,b) ((a)-((long)((a)/(b))*(b)))
# endif
using
namespace
std
;
#ifndef SINGLE_PRECISION
typedef
double
vec_type
;
#else
typedef
float
vec_type
;
#endif
//e "infinitely" large number in Vector_3 sense
//r "áåñêîíå÷íî áîëüøîå" ÷èñëî â ñìûñëå Vector_3
//# ifndef SINGLE_PRECISION
//# define VEC_INFTY 1.e20
//# else
//# define VEC_INFTY 1.e20f
//# endif
#define VEC_INFTY numeric_limits<vec_type>::max()
//e "infinitely" small number in Vector_3 sense (used for vector comparisons)
//r "áåñêîíå÷íî ìàëîå" ÷èñëî â ñìûñëå Vector_3 (èñïîëüçóåòñÿ äëÿ ñðàâíåíèé)
//# ifndef SINGLE_PRECISION
//# define VEC_ZERO 1.e-20
//# else
//# define VEC_ZERO 1.e-20f
//# endif
#define VEC_ZERO 512*numeric_limits<vec_type>::epsilon()
//r N-ìåðíûé âåêòîð òèïà T, ñ íåêîòîðûìè ïîëåçíûìè îïåðàöèÿìè
template
<
class
T
,
const
int
N
=
3
>
struct
Vector_Nt
{
typedef
T
value_type
;
T
v
[
N
];
Vector_Nt
(
const
T
&
a
=
0
){
for
(
int
i
=
0
;
i
<
N
;
i
++
)
v
[
i
]
=
a
;
}
explicit
Vector_Nt
(
const
T
&
a1
,
const
T
&
a2
)
{
if
(
N
>
0
)
v
[
0
]
=
a1
;
if
(
N
>
1
)
v
[
1
]
=
a2
;
for
(
int
i
=
2
;
i
<
N
;
i
++
)
v
[
i
]
=
0
;
}
explicit
Vector_Nt
(
const
T
&
s1
,
const
T
&
s2
,
const
T
&
s3
)
{
if
(
N
>
0
)
v
[
0
]
=
s1
;
if
(
N
>
1
)
v
[
1
]
=
s2
;
if
(
N
>
2
)
v
[
2
]
=
s3
;
for
(
int
i
=
3
;
i
<
N
;
i
++
)
v
[
i
]
=
0
;
}
//e construct from input iterator (or array)
template
<
class
A
>
Vector_Nt
(
const
A
*
beg
)
{
for
(
int
i
=
0
;
i
<
N
;
i
++
,
++
beg
)
v
[
i
]
=*
beg
;
};
//e construct from another vector
template
<
class
A
>
Vector_Nt
(
const
Vector_Nt
<
A
,
N
>&
v1
)
{
for
(
int
i
=
0
;
i
<
N
;
i
++
)
v
[
i
]
=
v1
[
i
];
};
//r Êîïèðóåò ñîäåðæèìîå âåêòîðà â it
template
<
class
A
>
void
copy_to
(
A
*
beg
)
const
{
for
(
int
i
=
0
;
i
<
N
;
i
++
,
++
beg
)
*
beg
=
v
[
i
];
};
//r ïîëó÷åíèå ýëåìåíòà
inline
T
&
operator
[](
int
i
)
const
{
return
(
T
&
)
v
[
i
];
};
inline
Vector_Nt
&
operator
=
(
const
Vector_Nt
&
vect
){
for
(
int
i
=
0
;
i
<
N
;
i
++
)
v
[
i
]
=
vect
.
v
[
i
];
return
*
this
;
}
//r ïðèñâàèâàåò âñåì êîìïîíåíòàì çíà÷åíèå a.
inline
Vector_Nt
&
operator
=
(
const
T
&
a
){
for
(
int
i
=
0
;
i
<
N
;
i
++
)
v
[
i
]
=
a
;
return
*
this
;
};
//r ñðàâíåíèå. Ïðè îòëè÷èè ìåíüøå ÷åì íà VEC_ZERO êîìïîíåíòû ñ÷èòàþòñÿ îäèíàêîâûìè
inline
bool
operator
==
(
const
Vector_Nt
&
vect
)
const
{
for
(
int
i
=
0
;
i
<
N
;
i
++
)
if
(
fabs
(
v
[
i
]
-
vect
.
v
[
i
])
>
VEC_ZERO
)
return
false
;
return
true
;
};
inline
bool
operator
!=
(
const
Vector_Nt
&
vect
)
const
{
return
(
!
(
this
->
operator
==
(
vect
)));
};
inline
Vector_Nt
operator
+
(
const
Vector_Nt
&
vect
)
const
{
Vector_Nt
result
;
for
(
int
i
=
0
;
i
<
N
;
i
++
)
result
.
v
[
i
]
=
v
[
i
]
+
vect
.
v
[
i
];
return
result
;
}
inline
Vector_Nt
operator
-
(
const
Vector_Nt
&
vect
)
const
{
Vector_Nt
result
;
for
(
int
i
=
0
;
i
<
N
;
i
++
)
result
.
v
[
i
]
=
v
[
i
]
-
vect
.
v
[
i
];
return
result
;
}
//r Ñêàëÿðíîå ïðîèçâåäåíèå âåêòîðîâ
inline
T
operator
*
(
const
Vector_Nt
&
vect
)
const
{
T
result
=
0
;
for
(
int
i
=
0
;
i
<
N
;
i
++
)
result
+=
v
[
i
]
*
vect
.
v
[
i
];
return
result
;
}
//r Ïîêîìïîíåíòíîå óìíîæåíèå íà êîýôôèöèåíò
inline
Vector_Nt
operator
*
(
const
T
&
coeff
)
const
{
Vector_Nt
result
;
for
(
int
i
=
0
;
i
<
N
;
i
++
)
result
[
i
]
=
coeff
*
v
[
i
];
return
result
;
}
//e vector multiplication (N=3 only)
//r Âåêòîðíîå ïðîèçâåäåíèå
inline
Vector_Nt
operator
%
(
const
Vector_Nt
&
r
)
const
{
//reserved for N specializations
if
(
N
==
3
){
return
Vector_Nt
(
v
[
1
]
*
r
.
v
[
2
]
-
v
[
2
]
*
r
.
v
[
1
],
v
[
2
]
*
r
.
v
[
0
]
-
v
[
0
]
*
r
.
v
[
2
],
v
[
0
]
*
r
.
v
[
1
]
-
v
[
1
]
*
r
.
v
[
0
]);
}
return
*
this
;
}
//r Óìíîæåíèå ÷èñëà íà âåêòîð (ïåðåñòàâëåíû ìåñòàìè ìíîæèòåëè).
// friend Vector_Nt operator*(T coeff,const Vector_Nt& vec);
//r Ïîêîìïîíåíòíîå äåëåíèå íà êîýôôèöèåíò
inline
Vector_Nt
operator
/
(
const
T
&
coeff
)
const
{
Vector_Nt
result
;
for
(
int
i
=
0
;
i
<
N
;
i
++
)
result
[
i
]
=
v
[
i
]
/
coeff
;
return
result
;
}
//r Óìíîæåíèå âåêòîðà íà -1
inline
Vector_Nt
operator
-
()
const
{
Vector_Nt
r
;
for
(
int
i
=
0
;
i
<
N
;
i
++
)
r
.
v
[
i
]
=-
v
[
i
];
return
r
;
}
//r Ñëîæåíèå ñ ïðèñâàèâàíèåì
inline
Vector_Nt
&
operator
+=
(
const
Vector_Nt
&
vect
){
for
(
int
i
=
0
;
i
<
N
;
i
++
)
v
[
i
]
+=
vect
.
v
[
i
];
return
*
this
;
}
//r Âû÷èòàíèå ñ ïðèñâàèâàíèåì
inline
Vector_Nt
&
operator
-=
(
const
Vector_Nt
&
vect
){
for
(
int
i
=
0
;
i
<
N
;
i
++
)
v
[
i
]
-=
vect
.
v
[
i
];
return
*
this
;
}
//r Óìíîæåíèå íà êîýôôèöèåíò ñ ïðèñâàèâàíèåì
inline
Vector_Nt
&
operator
*=
(
const
T
&
coeff
){
for
(
int
i
=
0
;
i
<
N
;
i
++
)
v
[
i
]
*=
coeff
;
return
*
this
;
}
//r Äåëåíèå íà ñêàëÿð ñ ïðèñâàèâàíèåì
inline
Vector_Nt
&
operator
/=
(
const
T
&
coeff
){
for
(
int
i
=
0
;
i
<
N
;
i
++
)
v
[
i
]
/=
coeff
;
return
*
this
;
}
//r Êâàäðàò íîðìû âåêòîðà
T
norm2
()
const
{
T
result
=
0
;
for
(
int
i
=
0
;
i
<
N
;
i
++
)
result
+=
v
[
i
]
*
v
[
i
];
return
result
;
}
//r Íîðìà âåêòîðà
T
norm
()
const
{
return
sqrt
(
norm2
());
}
//r Âîçâðàùàåò íîðìó è íîðìàëèçóåò âåêòîð (ïîñëå ýòîãî åãî norm() âåðíåò newnorm).
T
normalize
(
T
newnorm
=
1.
){
T
norm
=
this
->
norm
();
if
(
norm
>=
VEC_ZERO
){
T
c
=
newnorm
/
norm
;
for
(
int
i
=
0
;
i
<
N
;
i
++
)
v
[
i
]
*=
c
;
}
return
norm
;
}
//e Normalizes a vector that may have infinite components
T
inormalize
(
T
newnorm
=
1.
){
T
result
=
0
;
for
(
int
i
=
0
;
i
<
N
;
i
++
){
if
(
fabs
(
v
[
i
])
>=
VEC_INFTY
){
if
(
result
>=
0
)
result
=
0.
;
result
-=
1
;
v
[
i
]
=
v
[
i
]
>
0
?
1
:
-
1
;
}
else
if
(
result
>=
0
)
//still summing the normal components
result
+=
v
[
i
]
*
v
[
i
];
else
v
[
i
]
=
0.
;
}
if
(
fabs
(
result
)
<
VEC_ZERO
)
return
0.
;
newnorm
/=
sqrt
(
fabs
(
result
));
for
(
int
i
=
0
;
i
<
N
;
i
++
)
v
[
i
]
*=
newnorm
;
return
result
<
0
?
VEC_INFTY
:
result
;
}
//e nearest image distance within rectangular cell (FOR DISTANCE MEASUREMENTS)
//e assumes that each coordinate absolute value is in the range [0,cell[i])
//e returned vector is in the range [-cell[i]/2,cell[i]/2)
//e flags indicate the periodicity in specific directions: 0x1 for X, 0x2 for Y, 0x4 for Z
//r Áëèæàéøèé îáðàç â ïðÿìîóãîëüíîé ÿ÷åéêå
/*r Ñ÷èòàåì, ÷òî âñå ïðîñòðàíñòâî ðàçäåëåíî íà ÿ÷åéêè - ïàðàëëåëåïèïåäû ñ ðåáðàìè, ïàðàëëåëüíûìè
îñÿì êîîðäèíàò è äèàãîíàëüþ, çàäàííîé âåêòîðîì rcell, ïðè÷åì íà÷àëî êîîðäèíàò ÿâëÿåòñÿ
öåíòðîì îäíîé èç ÿ÷ååê. Åñëè *this íàõîäèòñÿ öåíòðàëüíîé ÿ÷åéêå, âîçâðàùàåòñÿ êîïèÿ *this.\n
Èíà÷å, åñëè *this íàõîäèòñÿ â êóáå 3*3 ÿ÷åéêè ñ öåíòðîì â íà÷àëå êîîðäèíàò, òî ñîçäàåò îáðàç
*this â öåíòðàëüíîé ÿ÷åéêå.\n
Èíà÷å, âîçâðàùàåò íåîïðåäåëåííîå çíà÷åíèå.
*/
Vector_Nt
rcell1
(
const
Vector_Nt
&
cell
,
int
flags
=
0xffff
)
const
{
Vector_Nt
ret
(
*
this
);
int
i
;
for
(
i
=
0
;
i
<
N
;
i
++
){
if
(
flags
&
(
0x1
<<
i
)){
if
(
v
[
i
]
>
cell
[
i
]
/
2
){
ret
[
i
]
-=
cell
[
i
];
}
else
if
(
v
[
i
]
<-
cell
[
i
]
/
2
){
ret
[
i
]
+=
cell
[
i
];
}
}
}
return
ret
;
}
//e @return a vector projection on a given axis
Vector_Nt
prj
(
int
i
)
const
{
Vector_Nt
res
;
res
[
i
]
=
v
[
i
];
return
res
;
}
//e @return a vector projection on a given vector (k*v)*k
Vector_Nt
prj
(
const
Vector_Nt
&
k
)
const
{
return
k
*
(
k
*
(
*
this
));
}
//e @return a vector of length l parallel to this
Vector_Nt
unit
(
T
l
=
1
)
const
{
Vector_Nt
res
(
*
this
);
res
.
normalize
(
l
);
return
res
;
}
//e reduction to elementary cell [0, cell[i]) (FOR REDUCTION TO ELEMENTARY CELL)
//e flags indicate the periodicity in specific directions: 0x1 for X, 0x2 for Y, 0x4 for Z
//r Ïî÷òè òî æå, ÷òî è rcell1, íî áåç îãðàíè÷åíèÿ íà ïîëîæåíèå *this è ñ äðóãîé ñèñòåìîé ÿ÷ååê.
/*r  íà÷àëå êîîðäèíàò íàõîäèòñÿ íå öåíòð ÿ÷åéêè, à åå óãîë. Ìîæåò ðàáîòàòü ìåäëåííåå èç-çà íàëè÷èÿ
îïåðàöèè äåëåíèÿ ïî ìîäóëþ ñ ïëàâàþùåé òî÷êîé */
Vector_Nt
rcell
(
const
Vector_Nt
&
cell
,
int
flags
=
0xffff
)
const
{
Vector_Nt
ret
(
*
this
);
for
(
int
i
=
0
,
flag
=
1
;
i
<
N
;
i
++
,
flag
<<=
1
)
{
if
(
flags
&
flag
){
ret
.
v
[
i
]
=
fmod
(
v
[
i
],
cell
[
i
]);
if
(
ret
.
v
[
i
]
<
0
)
ret
.
v
[
i
]
+=
cell
[
i
];
// if(ret.v[i]<0 || ret.v[i]>cell[i])
// printf("!");
}
}
return
ret
;
}
Vector_Nt
rpcell
(
const
Vector_Nt
&
p1
,
const
Vector_Nt
&
cell
,
int
flags
=
0xfff
)
const
{
Vector_Nt
ret
(
*
this
);
for
(
int
i
=
0
,
flag
=
1
;
i
<
N
;
i
++
,
flag
<<=
1
)
{
if
(
flags
&
flag
){
if
(
ret
.
v
[
i
]
<
p1
[
i
]
||
ret
.
v
[
i
]
>
p1
[
i
]
+
cell
[
i
])
{
ret
.
v
[
i
]
=
fmod
(
v
[
i
]
-
p1
[
i
],
cell
[
i
])
+
p1
[
i
];
if
(
ret
.
v
[
i
]
<
p1
[
i
])
ret
.
v
[
i
]
+=
cell
[
i
];
}
}
}
return
ret
;
}
//r Âîçâðàùàåò ìàêñèìàëüíóþ êîìïîíåíòó âåêòîðà è åå èíäåêñ â ind
T
maxcoord
(
int
*
ind
=
NULL
)
const
{
int
im
=
0
;
T
vv
=
v
[
0
];
for
(
int
i
=
1
;
i
<
N
;
i
++
)
{
if
(
v
[
i
]
>
vv
){
im
=
i
;
vv
=
v
[
i
];
}
}
if
(
ind
)
*
ind
=
im
;
return
vv
;
}
//e returns the corrd having maximal absolute value
T
maxabscoord
(
int
*
ind
=
NULL
)
const
{
int
im
=
0
;
T
vv
=
fabs
(
v
[
0
]);
for
(
int
i
=
1
;
i
<
N
;
i
++
)
{
if
(
fabs
(
v
[
i
])
>
vv
){
im
=
i
;
vv
=
fabs
(
v
[
i
]);
}
}
if
(
ind
)
*
ind
=
im
;
return
v
[
im
];
}
//e returns the corrd having minimal absolute value
T
minabscoord
(
int
*
ind
=
NULL
)
const
{
int
im
=
0
;
T
vv
=
fabs
(
v
[
0
]);
for
(
int
i
=
1
;
i
<
N
;
i
++
)
{
if
(
fabs
(
v
[
i
])
<
vv
){
im
=
i
;
vv
=
fabs
(
v
[
i
]);
}
}
if
(
ind
)
*
ind
=
im
;
return
v
[
im
];
}
//r Âîçâðàùàåò ìèíèìàëüíóþ êîìïîíåíòó âåêòîðà è åå èíäåêñ â ind
T
mincoord
(
int
*
ind
=
NULL
)
const
{
int
im
=
0
;
T
vv
=
v
[
0
];
for
(
int
i
=
1
;
i
<
N
;
i
++
)
{
if
(
v
[
i
]
<
vv
){
im
=
i
;
vv
=
v
[
i
];
}
}
if
(
ind
)
*
ind
=
im
;
return
vv
;
}
//r Âûâîäèò âåêòîð â ïîòîê âûâîäà ïî óìîë÷àíèþ, â ôîðìàòå (x,y,z)\\n
void
print
()
const
{
cout
<<
"("
;
for
(
int
i
=
0
;
i
<
N
;
i
++
){
cout
<<
v
[
i
];
if
(
i
!=
N
-
1
)
cout
<<
", "
;
}
cout
<<
")
\n
"
;
}
//e returns true if the vector has infinite components
bool
infinite
()
const
{
for
(
int
i
=
0
;
i
<
N
;
i
++
){
if
(
fabs
(
v
[
i
])
>=
VEC_INFTY
)
return
true
;
}
return
false
;
}
};
template
<
class
T
,
int
N
>
Vector_Nt
<
T
,
N
>
operator
*
(
const
T
&
coeff
,
const
Vector_Nt
<
T
,
N
>
&
vec
){
return
vec
*
coeff
;
}
template
<
class
T
,
int
N
>
Vector_Nt
<
T
,
N
>
operator
*
(
int
coeff
,
const
Vector_Nt
<
T
,
N
>
&
vec
){
return
vec
*
coeff
;
}
//template <> Vector_Nt<double, 3> operator*(const double &coeff,const Vector_Nt<double,3> &vec);
// old Vector_3 compatibility typedefs and functions
typedef
Vector_Nt
<
vec_type
,
3
>
Vector_3
;
typedef
Vector_3
*
Vector_3P
;
typedef
Vector_Nt
<
vec_type
,
2
>
Vector_2
;
typedef
Vector_2
*
Vector_2P
;
template
<
int
N
>
class
Vector_N
:
public
Vector_Nt
<
vec_type
,
N
>
{
};
//Vector_3 operator*(int coeff,const Vector_3 &vec){
// return vec*((vec_type)coeff);
//}
//e finds the maximum distance between vector pairs
//r Íàõîäèò ìàêñèìàëüíîå ðàññòîÿíèå ìåæäó âåêòîðàìè va1[i], va2[i], i=1..n
/*r @param va1 - ìàññèâ Vector_3[n]
@param n - äëèíà ìàññèâîâ va1 è va2
*/
vec_type
dist_max
(
Vector_3
*
va1
,
Vector_3
*
va2
,
int
n
);
//e finds average distance between vector pairs
//r Íàõîäèò ñðåäíåå ðàññòîÿíèå ìåæäó âåêòîðàìè va1[i], va2[i], i=1..n
vec_type
dist_av
(
Vector_3
*
va1
,
Vector_3
*
va2
,
int
n
);
//e finds the average difference norm between two vector sets of the same length
/*e optionally gives the indexes for maximal and minimal difference
va2 can be NULL, then the norm of va1 is used */
//r Íàõîäèò ñðåäíåå ðàññòîÿíèå ìåæäó va1[i] è va2[i], à òàêæå, ïî æåëàíèþ, èíäåêñû, íà êîòîðûõ äîñòèãàåòñÿ min è max ðàññòîÿíèå
vec_type
diff_av
(
Vector_3
*
va1
,
Vector_3
*
va2
,
int
n
,
int
*
minind
=
0
,
int
*
maxind
=
0
);
//e finds suitable perpendicular to a vector
//r Íàõîäèò ïåðïåíäèêóëÿð ê âåêòîðó vAB
Vector_3
FindPerp
(
const
Vector_3
&
vAB
);
//e Returns the average (center) vector of the vector array
//e and cooordinates of a minimal box in which it is contained
Vector_3
GetScope
(
const
Vector_3
*
varr
,
long
n
,
Vector_3
*
box_min
,
Vector_3
*
box_max
);
//e the same with long index array
Vector_3
GetIScope
(
const
Vector_3
*
varr
,
long
*
indarr
,
long
n
,
Vector_3
*
box_min
,
Vector_3
*
box_max
);
//e the same with int index array
Vector_3
GetIScopei
(
const
Vector_3
*
varr
,
int
*
indarr
,
int
n
,
Vector_3
*
box_min
,
Vector_3
*
box_max
);
// neue Funktionen
//e clears vector array with optional integer index
//r Î÷èñòêà ìàññèâà âåêòîðîâ, ñ ïîääåðæêîé èíäåêñèðîâàíèÿ
/*r
 äàííîì Vector_3 vec[] îáíóëÿåò n êîîðäèíàò. Åñëè ind==NULL, òî
î÷èùàåò ïåðâûå n ýëåìåíòîâ. Åñëè ind!=NULL, òî äëÿ i=0..n-1
î÷èùàåò vec[ind[i]]
Ñì. @ref indexed_calculations.
*/
void
clear_vecarri
(
int
n
,
Vector_3
*
vec
,
int
*
ind
=
0
);
//e reflects the vector ini+dir*t+0.5*force*t^2 to be inside a box limited by 0 and box sizes
//e changes dir according to the final state
//e fills crossed dir with bit flags corresponding directions along which the walls were crossed
Vector_3
Reflect
(
Vector_3
&
ini
,
double
t
,
Vector_3
&
dir
,
double
*
box
,
int
flag
=
0x7
,
const
Vector_3
&
force
=
Vector_3
());
inline
vec_type
vec_area
(
const
Vector_2
&
vect1
,
const
Vector_2
&
vect2
)
{
return
fabs
(
vect1
[
0
]
*
vect2
[
1
]
-
vect1
[
1
]
*
vect2
[
0
])
/
2
;
};
inline
vec_type
vec_area
(
const
Vector_3
&
vect1
,
const
Vector_3
&
vect2
)
{
return
(
vect1
%
vect2
).
norm
()
/
2
;
};
// remake for vec_type!
template
<
class
num_t
,
class
denum_t
,
class
val_t
>
denum_t
acccomp
(
const
num_t
x
,
const
denum_t
y
,
const
val_t
val
)
{
num_t
eps_num
=
numeric_limits
<
num_t
>::
epsilon
();
denum_t
eps_denum
=
numeric_limits
<
denum_t
>::
epsilon
();
return
(
eps_num
>=
eps_denum
)
?
acccomp
(
x
,
(
num_t
)
y
,
(
num_t
)
val
)
:
acccomp
((
denum_t
)
x
,
y
,
(
denum_t
)
val
);
}
template
<
class
num_t
,
class
denum_t
>
denum_t
acccomp
(
const
num_t
x
,
const
denum_t
y
)
{
return
acccomp
(
x
,
y
,
num_t
(
1
));
}
template
<
class
T
>
bool
acccomp
(
const
T
x
,
const
T
y
,
const
T
val
)
{
T
eps
=
numeric_limits
<
T
>::
epsilon
();
int
iexp
;
frexp
(
val
,
&
iexp
);
T
mult
=
ldexp
(
.5
,
iexp
+
1
);
T
err
=
T
(
256
)
*
mult
*
eps
;
return
fabs
(
x
-
y
)
<=
err
;
}
template
<
class
T
>
bool
acccomp
(
const
T
x
,
const
T
y
)
{
return
acccomp
(
x
,
y
,
T
(
1
));
}
template
<
class
num_t
,
class
denum_t
>
denum_t
accdiv
(
const
num_t
x
,
const
denum_t
y
)
{
num_t
eps_num
=
numeric_limits
<
num_t
>::
epsilon
();
denum_t
eps_denum
=
numeric_limits
<
denum_t
>::
epsilon
();
return
(
eps_num
>=
eps_denum
)
?
accdiv1
(
x
,
(
num_t
)
y
)
:
accdiv1
((
denum_t
)
x
,
y
);
}
template
<
class
T
>
T
accdiv1
(
T
x
,
T
y
)
{
T
result
;
T
eps
=
numeric_limits
<
T
>::
epsilon
();
T
fr
=
x
/
y
;
T
ifr
=
floor
(
fr
+
T
(
.5
));
// T err=64*eps*(ifr!=0 ? fabs(ifr) : 1);
int
mult
;
if
(
ifr
<=
512
)
mult
=
512
;
else
{
int
iexp
;
frexp
(
ifr
,
&
iexp
);
mult
=
int
(
ldexp
(
.5
,
iexp
+
1
));
}
T
err
=
mult
*
eps
;
if
(
fabs
(
fr
-
ifr
)
<=
err
)
result
=
ifr
;
else
result
=
fr
;
return
result
;
}
template
<
class
num_t
,
class
denum_t
,
class
P
>
denum_t
accdiv_rmd
(
const
num_t
x
,
const
denum_t
y
,
P
*
remainder
)
{
num_t
eps_num
=
numeric_limits
<
num_t
>::
epsilon
();
denum_t
eps_denum
=
numeric_limits
<
denum_t
>::
epsilon
();
return
(
eps_num
>=
eps_denum
)
?
accdiv_rmd
(
x
,
(
num_t
)
y
,
remainder
)
:
accdiv_rmd
((
denum_t
)
x
,
y
,
remainder
);
}
template
<
class
T
,
class
P
>
T
accdiv_rmd
(
const
T
x
,
const
T
y
,
P
*
remainder
)
{
T
result
=
accdiv
(
x
,
y
);
T
iresult
=
floor
(
result
);
if
(
result
-
iresult
>
0
)
*
remainder
=
x
-
iresult
*
y
;
else
*
remainder
=
0
;
return
result
;
}
//e returns random unit vector uniformely distributed in space (?? check this)
inline
Vector_3
randdir
(){
vec_type
xi1
=
2
*
((
vec_type
)
rand
())
/
RAND_MAX
-
1.
;
vec_type
xi2
=
((
vec_type
)
rand
())
/
RAND_MAX
;
vec_type
r1
=
sqrt
(
1.
-
xi1
*
xi1
);
return
Vector_3
(
r1
*
cos
(
2
*
M_PI
*
xi2
),
r1
*
sin
(
2
*
M_PI
*
xi2
),
xi1
);
}
///\en Calculates extent of the vector container.
/// \return the center of the vector set, optionally
/// (if arguments are not NULL) fills the bounding box in \a box_min, \a box_max.
template
<
class
vec_inp_it
>
Vector_3
get_extent
(
vec_inp_it
beg
,
vec_inp_it
end
,
Vector_3
*
box_min
=
NULL
,
Vector_3
*
box_max
=
NULL
){
if
(
beg
==
end
)
return
Vector_3
();
Vector_3
center
(
*
beg
++
);
Vector_3
cube1
(
center
),
cube2
(
center
);
size_t
n
=
1
;
for
(;
beg
!=
end
;
++
beg
){
Vector_3
vec
=*
beg
;
center
+=
vec
;
for
(
size_t
j
=
0
;
j
<
3
;
j
++
){
if
(
cube1
[
j
]
>
vec
[
j
])
cube1
[
j
]
=
vec
[
j
];
if
(
cube2
[
j
]
<
vec
[
j
])
cube2
[
j
]
=
vec
[
j
];
}
n
++
;
}
if
(
box_min
)
*
box_min
=
cube1
;
if
(
box_max
)
*
box_max
=
cube2
;
return
center
/
n
;
}
///\en Returns \a true if absolute value of \a a is at least \a factor times less than
/// the absolute value of \a b.
template
<
class
T
>
bool
much_less
(
const
T
&
a
,
const
T
&
b
,
const
T
&
factor
){
if
(
fabs
(
a
)
<
fabs
(
b
))
return
fabs
(
a
*
factor
/
b
)
<
1
?
true
:
false
;
else
return
false
;
}
# endif
// __VECTOR_3_H
Event Timeline
Log In to Comment