i have a 3d fluid grid i am stuck at advecttion:

its described like this:
Here is how it works. For each voxel of the velocity grid we
trace its midpoint backwards through the velocity field over a time step dt. This point will end up somewhere else in
the grid. We then linearly interpolate the velocity at that point from the neighboring voxels and transfer this velocity
back to the departure voxel.

code attached to the pdf looks like:

dt0 = dt*N;
for ( i=1 ; i<=N ; i++ ) {
for ( j=1 ; j<=N ; j++ ) {

x = i-dt0*u[IX(i,j)]; y = j-dt0*v[IX(i,j)]; //destination point

if (x<0.5) x=0.5; if (x>N+0.5) x=N+ 0.5; i0=(int)x; i1=i0+1; //i0 and j0 are integers so we can calculate difference between integer and floating result
if (y<0.5) y=0.5; if (y>N+0.5) y=N+ 0.5; j0=(int)y; j1=j0+1;

s1 = x-i0; s0 = 1-s1; t1 = y-j0; t0 = 1-t1; //that is such difference lets say x = 10.3 so i0 is 10 result s1 = 0.3; then s0 = 0.7, if t1= 0.1 then t0 = 0.9

TargetDensityAt[i][j] = s0* (t0*dens[IX(i0,j0)]+t1*dens[IX(i0,j1)])+
s1*(t0*dens[IX(i1,j0)]+t1*dens[IX(i1,j1)]);

TargetDensityAt[i][j] = 0.7 (0.9*DestPointDens+0.1*destpointdens+1y)+
0.3*(0.9*destPointDens+1x+0.1*destpointdens+1x+1y);

}

well it doesnt make any sense for me thus DestPoint[x+1][y+1] shouldnt even be checked (in my opinion)

btw from explanation is should be interpolated from points on the grid like: [x-1][y] [x+1][y] [x][y-1] [x][y+1]...

but to the point i have a 3d grid of this and i have no idea how to do lerp of these 6 corresponding 'voxels' indices in grid

that is [x+1][y][z] [x-1][y][z] [x][y-1][z] [x][y+1][z] [x][y][z+1] [x][y][z-1]

In 1D, a linearly interpolated point value combines the 2 closest samples -- the ones immediately to the left/right of the point.
In 2D, it's the same, but doubling the procedure over the up/down axis, combining the (2x2=)4 closest pixels (up-left, up-right, down-left, down-right -- NOT up, down, left, right).
The 3D case is the same, but done on the layer of voxels under the point and the layer of voxels above the point, resulting in (4x2=)8 voxels being combined (not 6).

ok i wrote a code based on this: http://en.wikipedia.org/wiki/Trilinear_interpolation

it seems old code was working the same as the new one, so i guess my bug is somewhere else