21 return (make_pair(wpqs, wpqf)) ;
25 virtual double distance (
int id,
v1d loc) {
double res=0 ;
for (
int i=0 ; i<
d ; i++) res+=(loc[i]-
data->
pos[i][
id])*(loc[i]-
data->
pos[i][id]) ;
return sqrt(res) ; }
26 virtual double cutoff (
void) {
return 2.5*
w ;}
33 double Lucy (
double r) {
if (r>=
w)
return 0 ;
else {
double f=r/
w ;
return (
cst*(-3*f*f*f*f + 8*f*f*f - 6*f*f +1)) ; }}
48 double prev, cur, first ;
51 for (
int i=0 ; i<
d ;i++)
60 for (
int i=0 ; i<
Nsteps ; i++)
62 for (
int j=0; j<
d ; j++)
65 res += (cur+prev)/(2*
Nsteps) ;
70 return (make_pair(
window_avg(first,cur),res)) ;
72 double distancevec (
v1d l1,
v1d loc) {
double res=0 ;
for (
int i=0 ; i<
d ; i++) res+=(loc[i]-l1[i])*(loc[i]-l1[i]) ;
return sqrt(res) ;}
73 void set_integrationsteps (
int steps) {
if (steps<1) printf(
"Less than 1 step is not meaningful.") ;
if (steps>1000) printf(
"You've chosen a very large number of integration steps, you may want to reconsider") ;
Nsteps=steps ; }
85 double window (
double r) {
if (r>=
w)
return 0 ;
else return cst ;}
91 double window (
double r) {
if (r>=
w)
return 0 ;
else {
double a =1 ;
for (
int b=0 ; b<
d ; b++,
a*=
w) {}
return 1/
a ;}}
102 for (
int i=0 ; i<
d ; i++)
106 if (dst>
w+r)
return 0 ;
107 else if (dst<=fabs(
w-r))
112 return (1/(4./3. * M_PI * r*r*r)) ;
116 double vol = M_PI * (
w+r-dst) * (
w+r-dst) * (dst*dst + 2*dst*r - 3*r*r + 2*dst*
w + 6*r*
w - 3*
w*
w) / (12.*dst) ;
117 return (vol *
cst / (4./3.*M_PI*r*r*r)) ;
123 if (
data->
N==0)
return 2*
w ;
124 for (
int i=0 ; i<
data ->
N ; i++)
140 double normp=0, normlpq=0 ;
142 for (
int i=0 ; i<
d ;i++)
145 locp[i] =
data->
pos[i][p]-loc[i] ;
146 normp += locp[i]*locp[i] ;
147 normlpq +=
lpq[i]*
lpq[i] ;
148 b += locp[i]*
lpq[i] ;
152 double Delta = b*b-4*normlpq*(normp-
w*
w) ;
157 double alpha1 = (-b - sqrt(Delta))/(2*normlpq) ;
158 double alpha2 = (-b + sqrt(Delta))/(2*normlpq) ;
160 if (alpha1<0 && alpha2>0 && alpha2<1)
161 wpqf =
cst * alpha2 ;
162 else if (alpha1<0 && alpha2>1)
164 else if (alpha1>0 && alpha1<1 && alpha2>1)
165 wpqf =
cst * (1-alpha1) ;
166 else if (alpha1>0 && alpha1<1 && alpha2>0 && alpha2<1)
167 wpqf =
cst*(alpha2-alpha1) ;
172 return (make_pair(wpqs,wpqf)) ;
183 if (h<0 || h>2*r)
return 0 ;
187 return 1/2. *
nballvolume() * pow(r,
d) * boost::math::ibeta ((
d+1)/2., 1/2., (2*r*h-h*h)/(r*r)) ;
193 return nballvolume() * pow(r,
d) - 1/2. *
nballvolume() * pow(r,
d) * boost::math::ibeta((
d+1)/2., 1/2., (2*r*(2*r-h)-(2*r-h)*(2*r-h))/(r*r)) ;
201 for (
int i=0 ; i<
d ; i++)
206 if (r>r1+r2)
return 0 ;
207 else if (r<=fabs(r1-r2))
216 double h1 = (r2+r1-r)*(r2-r1+r)/(2*r) ;
217 double h2 = (r1+r2-r)*(r1-r2+r)/(2*r) ;
225 if (
data->
N==0)
return 2*
w ;
226 for (
int i=0 ; i<
data ->
N ; i++)
242 double normp=0, normlpq=0 ;
244 for (
int i=0 ; i<
d ;i++)
247 locp[i] =
data->
pos[i][p]-loc[i] ;
248 normp += locp[i]*locp[i] ;
249 normlpq +=
lpq[i]*
lpq[i] ;
250 b += locp[i]*
lpq[i] ;
254 double Delta = b*b-4*normlpq*(normp-
w*
w) ;
259 double alpha1 = (-b - sqrt(Delta))/(2*normlpq) ;
260 double alpha2 = (-b + sqrt(Delta))/(2*normlpq) ;
262 if (alpha1<0 && alpha2>0 && alpha2<1)
263 wpqf =
cst * alpha2 ;
264 else if (alpha1<0 && alpha2>1)
266 else if (alpha1>0 && alpha1<1 && alpha2>1)
267 wpqf =
cst * (1-alpha1) ;
268 else if (alpha1>0 && alpha1<1 && alpha2>0 && alpha2<1)
269 wpqf =
cst*(alpha2-alpha1) ;
274 return (make_pair(wpqs,wpqf)) ;
283 double window (
double r) {
if (r>=
w)
return 0 ;
else return cst*pow(cos(M_PI*r/(2.*
w)),2) ;}
291 double Vol=pow(M_PI,
d/2.)/(tgamma(
d/2.+1)) ;
292 scale = Vol *
d * (-3./(
d+4) + 8./(
d+3) - 6./(
d+2) + 1./
d) ;
295 printf(
"Lucy function scaling : %f \n",
scale) ;
298 double Lucy (
double r) {
if (r>=
w)
return 0 ;
else {
double f=r/
w ;
return (
scale*(-3*f*f*f*f + 8*f*f*f - 6*f*f +1)) ; }}
306 :
LibLucyND(D,ww,dd-__builtin_popcount(periodic))
308 printf(
"%g ",
scale) ;
310 for (
size_t i =0 ; i<boxes.size() ; i++)
311 if ((periodic&(1<<i)) && boxes[i]==1)
316 printf(
"NB: do not use periodic_atoms with LibLucyND_Periodic! %g \n",
scale) ;
322 for (
int i=0 ; i<
d ; i++)
Windows
Definition: WindowLibrary.h:2
A window base class that needs to be specialised to a specific CG window.
Definition: WindowLibrary.h:6
virtual std::pair< double, double > window_contact_weight(int p, int q, const v1d &loc)
Definition: WindowLibrary.h:16
virtual ~LibBase()
Definition: WindowLibrary.h:10
virtual double window(double)=0
virtual double distance(int id, v1d loc)
function for mixed particle id / vector informations.
Definition: WindowLibrary.h:25
virtual double window_int(double r1, double r2)
Definition: WindowLibrary.h:23
virtual double window_avg(double r1, double r2)
Definition: WindowLibrary.h:24
struct Data * data
Definition: WindowLibrary.h:12
double d
Definition: WindowLibrary.h:13
LibBase()
Definition: WindowLibrary.h:9
virtual double cutoff(void)
Definition: WindowLibrary.h:26
LibBase(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:8
double w
Definition: WindowLibrary.h:13
Definition: WindowLibrary.h:279
double window(double r)
Definition: WindowLibrary.h:283
LibHann3D(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:281
double cst
Definition: WindowLibrary.h:282
Definition: WindowLibrary.h:42
LibLucy3DFancyInt(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:44
int Nsteps
Definition: WindowLibrary.h:75
std::pair< double, double > window_contact_weight(int p, int q, const v1d &loc)
Definition: WindowLibrary.h:45
double distancevec(v1d l1, v1d loc)
Definition: WindowLibrary.h:72
void set_integrationsteps(int steps)
Definition: WindowLibrary.h:73
Definition: WindowLibrary.h:29
double window(double r)
Definition: WindowLibrary.h:34
double cutoff(void)
Definition: WindowLibrary.h:35
double cst
Definition: WindowLibrary.h:39
LibLucy3D(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:31
double Lucy(double r)
Definition: WindowLibrary.h:33
LibLucy3D()
Definition: WindowLibrary.h:32
Definition: WindowLibrary.h:303
LibLucyND_Periodic(struct Data *D, double ww, double dd, int periodic, vector< int >boxes, vector< double > deltas)
Definition: WindowLibrary.h:305
int maskperiodic
Definition: WindowLibrary.h:328
double distance(int id, v1d loc) override
function for mixed particle id / vector informations.
Definition: WindowLibrary.h:319
Definition: WindowLibrary.h:286
double scale
Definition: WindowLibrary.h:297
double Lucy(double r)
Definition: WindowLibrary.h:298
double window(double r) override
Definition: WindowLibrary.h:299
LibLucyND(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:288
Definition: WindowLibrary.h:81
double window(double r)
Definition: WindowLibrary.h:85
double cst
Definition: WindowLibrary.h:84
LibRect3D(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:83
Definition: WindowLibrary.h:88
LibRectND(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:90
double window(double r)
Definition: WindowLibrary.h:91
Definition: WindowLibrary.h:94
double cst
Definition: WindowLibrary.h:97
virtual double window_avg(double r1, double r2)
Definition: WindowLibrary.h:174
LibSphere3DIntersect(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:96
double windowreal(double r)
Definition: WindowLibrary.h:130
double window(double r)
Definition: WindowLibrary.h:129
double cutoff(void)
Definition: WindowLibrary.h:121
virtual std::pair< double, double > window_contact_weight(int p, int q, const v1d &loc)
Definition: WindowLibrary.h:131
double result
Definition: WindowLibrary.h:98
double distance(int id, v1d loc)
function for mixed particle id / vector informations.
Definition: WindowLibrary.h:99
Definition: WindowLibrary.h:177
double nballvolume()
Definition: WindowLibrary.h:180
double cutoff(void)
Definition: WindowLibrary.h:223
LibSphereNDIntersect(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:179
virtual std::pair< double, double > window_contact_weight(int p, int q, const v1d &loc)
Definition: WindowLibrary.h:233
double windowreal(double r)
Definition: WindowLibrary.h:232
virtual double window_avg(double r1, double r2)
Definition: WindowLibrary.h:276
double cst
Definition: WindowLibrary.h:196
double distance(int id, v1d loc)
function for mixed particle id / vector informations.
Definition: WindowLibrary.h:198
double window(double r)
Definition: WindowLibrary.h:231
double result
Definition: WindowLibrary.h:197
double ncapvolume(double r, double h)
Definition: WindowLibrary.h:181
double * radius
Particle radius.
Definition: Coarsing.h:97
int N
Number of particles.
Definition: Coarsing.h:96
vector< double * > pos
Particle positions.
Definition: Coarsing.h:100
@ lpq
Definition: Typedefs.h:19
@ radius
Definition: Typedefs.h:19
vector< double > v1d
Definition: Typedefs.h:9
const GenericPointer< typename T::ValueType > T2 T::AllocatorType & a
Definition: pointer.h:1181
Data structure handling point data and contact data.
Definition: Coarsing.h:93