simulations – Sphere – AABOX Edge detection

I am trying to implement a particle-AABOX edge collision, the below images represent, two timesteps(dt) where the spherical particle is accelerated by gravity.

enter image description here
enter image description here

I have the particle center C and its radius r, along with AABOX coordinates as B_min and Bmax.

I have to check if the particle collides with any edges, and then bounce back, if a collision happens.

If I use Arvo’s AABB vs Sphere collision algorithm, then I can find if there is an intersection within the radius, but I can’t get the hitpoint and hitdistance.

And after getting the hitpoint and hitdistance, how to reflect it back? (for every timestep, I will have the current position C, and also the previous position P).

I tried to calculate it by implementing,

BOX_N(0) = make_vector(-1, 0, 0);// x-min;
BOX_P(0) = make_point(bb_min.x, bb_min.y + (bb_max.y - bb_min.y) / 2.0f , bb_min.z + (bb_max.z - bb_min.z) / 2.0f);

BOX_N(1) = make_vector(1, 0, 0); // x-max
BOX_P(1) = make_point(bb_max.x, bb_min.y + (bb_max.y - bb_min.y) / 2.0f , bb_min.z + (bb_max.z - bb_min.z) / 2.0f);

BOX_N(2) = make_vector(0, -1, 0);// y-min;
BOX_P(2) = make_point(bb_min.x + (bb_max.x - bb_min.x) / 2.0f, bb_min.y, bb_min.z + (bb_max.z - bb_min.z) / 2.0f);

BOX_N(3) = make_vector(0, 1, 0); // y-max
BOX_P(3) = make_point(bb_min.x + (bb_max.x - bb_min.x) / 2.0f, bb_max.y, bb_min.z + (bb_max.z - bb_min.z) / 2.0f);

BOX_N(4) = make_vector(0, 0, -1);// z-min;
BOX_P(4) = make_point(bb_min.x + (bb_max.x - bb_min.x) / 2.0f, bb_min.y + (bb_max.y - bb_min.y) / 2.0f , bb_min.z);

BOX_N(5) = make_vector(0, 0, 1); // z-max
BOX_P(5) = make_point(bb_min.x + (bb_max.x - bb_min.x) / 2.0f, bb_min.y + (bb_max.y - bb_min.y) / 2.0f , bb_max.z);



__device__ int box_collision(const Point& previous_position, const Point& current_position, const Vector& direction, const float& radius, Point& hitPoint, Vector& normal, float& hit_distance) {
    Point boxPoint;
    int index;

         if (current_position.x < (bb_min.x + radius)) { index = 1; normal = BOX_N(0); boxPoint = BOX_P(0); boxPoint.x += radius;}
    else if (current_position.x > (bb_max.x - radius)) { index = 1; normal = BOX_N(1); boxPoint = BOX_P(1); boxPoint.x -= radius;}

    else if (current_position.y < (bb_min.y + radius)) { index = 2; normal = BOX_N(2); boxPoint = BOX_P(2); boxPoint.y += radius;}
    else if (current_position.y > (bb_max.y - radius)) { index = 2; normal = BOX_N(3); boxPoint = BOX_P(3); boxPoint.y -= radius;}

    else if (current_position.z < (bb_min.z + radius)) { index = 3; normal = BOX_N(4); boxPoint = BOX_P(4); boxPoint.z += radius;}
    else if (current_position.z > (bb_max.z - radius)) { index = 3; normal = BOX_N(5); boxPoint = BOX_P(5); boxPoint.z -= radius;}
    else return 0;

    auto denom = vdot(direction, normal);
    if (denom < 1e-6) return 0;

    hit_distance = vdot(normal, boxPoint - previous_position) / denom;
    if (hit_distance < 0) return 0;

    hitPoint = previous_position + hit_distance * direction;
    
    return index;
}

inline __device__ bool compute_box_collision(Point& previous_position, Point& current_position, const float& radius) {
    Point hitPoint;
    Vector normal;
    float hit_distance;

    Vector direction = vnormalize(current_position - previous_position);

    int collision = box_collision(previous_position, current_position, direction, radius, hitPoint, normal, hit_distance);

    if (!collision) return true;

    Vector damping{1, 1, 1};
    if (collision == 1) damping.x = bounce_factor;
    else if (collision == 2) damping.y = bounce_factor;
    else if (collision == 3) damping.z = bounce_factor;
        
    float relection_length = vlength(hitPoint - current_position);
    
    Vector R = vnormalize(direction - 2.0f * vdot(normal, direction) * normal) * damping;
            
    current_position = hitPoint + R * relection_length;
    previous_position = hitPoint - R * hit_distance;

    return false;
}

by calling the 2nd function to check for collision,

inline __device__ bool compute_box_collision(Point& previous_position, Point& current_position, const float& radius)

In the above code, I used the ray-plane interesection, and then use reflection of the direction to get the bounce direction, but after 500-550 timesteps, the solution deviates and becomes numerically unstable(possibly they settle at the bottom and collides continously). Is there something, I am missing here? Or am I doing something wrong?

Note: I am using simiulating 10k particles, so the solution deviation compounds by a large amount.