-
Notifications
You must be signed in to change notification settings - Fork 1
/
CloudInCell.c
109 lines (85 loc) · 3.34 KB
/
CloudInCell.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
int cic_assign(double x, double y, double z, double p){
// Assigns an object at (x, y, z) into a 3D density grid with weight p using Cloud-in-Cell assignment.
// e.g xCellSize is the cell size in the x direction.
int indx, indy, indz;
int iindx, iindy, iindz;
double dx, dy, dz;
double hx0, hy0, hz0, hx1, hy1, hz1;
double nuCellSize, newlim;
switch(fold){
case 0:
x /= xCellSize;
y /= yCellSize;
z /= zCellSize;
break;
case 1:
newlim = AxisLimsArray[1][2]/FOLDFACTOR;
nuCellSize = xCellSize/FOLDFACTOR;
x = fmod(x, newlim)/nuCellSize;
y = fmod(y, newlim)/nuCellSize;
z = fmod(z, newlim)/nuCellSize;
break;
case 2:
newlim = AxisLimsArray[1][2]/(FOLDFACTOR*FOLDFACTOR);
nuCellSize = xCellSize/(FOLDFACTOR*FOLDFACTOR);
x = fmod(x, newlim)/nuCellSize;
y = fmod(y, newlim)/nuCellSize;
z = fmod(z, newlim)/nuCellSize;
break;
}
// same as "box label" co-ordinates.
indx = (int) floor(x);
indy = (int) floor(y);
indz = (int) floor(z);
// distance from centre of cell.
dx = x - (double) indx - 0.5; // 0.5 for centre of cell labelled by (indx, indy, indz).
dy = y - (double) indy - 0.5;
dz = z - (double) indz - 0.5;
hx0 = 1.0 - fabs(dx);
hx1 = fabs(dx);
hy0 = 1.0 - fabs(dy);
hy1 = fabs(dy);
hz0 = 1.0 - fabs(dz);
hz1 = fabs(dz);
// periodic boundaries.
if(indx >= n2) indx -= n2;
if(indy >= n1) indy -= n1;
if(indz >= n0) indz -= n0;
// right of cell centre -> up a cell, otherwise down.
if(dx>0) iindx=indx+1; else iindx=indx-1;
if(dy>0) iindy=indy+1; else iindy=indy-1;
if(dz>0) iindz=indz+1; else iindz=indz-1;
// periodic boundaries.
if(iindx >= n2) iindx -= n2;
if(iindy >= n1) iindy -= n1;
if(iindz >= n0) iindz -= n0;
if(iindx < 0) iindx += n2;
if(iindy < 0) iindy += n1;
if(iindz < 0) iindz += n0;
// Assignment of 8 corners of cube in which the particle sits,
// Weights are ------
//
// ------
//
// W = product_i of W_i(x_i - x_cell) = 1. - |x_i - x_cell|
// Note: e.g. weight at bottom left = 1. - |dx|, at bottom right = 1. - (1. - dx)
/*
overdensity[iindx + n2*iindy + n1*n2*iindz][0] += hx1*hy1*hz1*p;
overdensity[indx + n2*iindy + n1*n2*iindz][0] += hx0*hy1*hz1*p;
overdensity[iindx + n2*indy + n1*n2*iindz][0] += hx1*hy0*hz1*p;
overdensity[indx + n2*indy + n1*n2*iindz][0] += hx0*hy0*hz1*p;
overdensity[iindx + n2*iindy + n1*n2*indz ][0] += hx1*hy1*hz0*p;
overdensity[indx + n2*iindy + n1*n2*indz ][0] += hx0*hy1*hz0*p;
overdensity[iindx + n2*indy + n1*n2*indz ][0] += hx1*hy0*hz0*p;
overdensity[indx + n2*indy + n1*n2*indz ][0] += hx0*hy0*hz0*p;
*/
overdensity[iindx + n2*iindy + n1*n2*iindz] += hx1*hy1*hz1*p;
overdensity[indx + n2*iindy + n1*n2*iindz] += hx0*hy1*hz1*p;
overdensity[iindx + n2*indy + n1*n2*iindz] += hx1*hy0*hz1*p;
overdensity[indx + n2*indy + n1*n2*iindz] += hx0*hy0*hz1*p;
overdensity[iindx + n2*iindy + n1*n2*indz ] += hx1*hy1*hz0*p;
overdensity[indx + n2*iindy + n1*n2*indz ] += hx0*hy1*hz0*p;
overdensity[iindx + n2*indy + n1*n2*indz ] += hx1*hy0*hz0*p;
overdensity[indx + n2*indy + n1*n2*indz ] += hx0*hy0*hz0*p;
return 0;
}