Lagrangian Interpolation

I want to write a code in C++ for the Lagrangian/Polynomial interpolation for the given equation:

void Function(p_struct particle) {
 
for (int X = 0; X < Nx; ++X) {
   for (int Y = 0; Y < Ny; ++Y) {
       Ax[X][Y] = 0;
       Ay[X][Y] = 0;
   }
}
 
for (int n = 0; n < particle.num_nodes; ++n) {
   int xStart = static_cast <int> (particle.node[n].x - 3.0);
   int xEnd = static_cast <int> (particle.node[n].x + 3.0);
   int yStart = static_cast <int> (particle.node[n].y - 3.0);
   int yEnd = static_cast <int> (particle.node[n].y + 3.0);
   for (int X = xStart; X < xEnd; ++X) {
       for (int Y = yStart; Y <= yEnd; ++Y) {
          **// here I want to use Lagrangian Interpolation**
            const double xDistance = X - particle.node[n].x;
            const double yDistance = Y - particle.node[n].y;
            const double delta = dirac_4(xDistance, yDistance);
            Ax[X][Y] += (particle.node[n].Ax * delta);
            Ay[X][Y] += (particle.node[n].Ay * delta);
        }
    }
}
return;
}

How can I get interpolated values " Ax[X][Y] and Ay[X][Y]" by changing my code.

Best Regards