How to add solar eclipses to my atmospheric scattering shader?

Godot Version

4.3

Background

I implemented an atmospheric scattering skybox shader from this github that i barely understand. I know that solar eclipses are essentially just the moon casting a shadow on the earth and i know that the way this shader works is by rendering a volume using ray marching. So my question is…

Question

How do i add shadows to a ray marching algorithm for rendering the atmosphere? The two functions i believe might be the most relevant are integrate_optical_depth:

vec3 intergrate_optical_depth(vec3 ray_start, vec3 ray_dir) {
  vec2 intersection = atmosphere_intersection(ray_start, ray_dir);
  float ray_length = intersection.y;
  
  int sample_count = 8;
  float step_size = ray_length / float(sample_count);
  
  vec3 optical_depth = vec3(0);
  
  for (int i = 0; i < sample_cound; i++) {
    vec3 local_position = ray_start + ray_dir * (float(i) + 0.5) * step_size;
    float local_height = atmosphere_height(local_position);
    vec3 local_density = atmosphere_density(local_height);
    
    optical_depth += local_density * step_size;
  }
  
  return optical_depth;
}

and integrate_scattering:

vec3 integrate_scattering(vec3 ray_start, vec3 ray_dir, float ray_length, vec3 light_dir, vec3 light_color, out vec3 transmittance) {
  float ray_height = atmosphere_height(ray_start);
  float sample_distribution_exponent = 1.0 + clamp(1.0 - ray_height / ATMOSPHERE_HEIGHT, 0, 1) * 8.0;
  
  vec2 intersection = atmosphere_intersection(ray_start, ray_dir);
  ray_length = min(ray_length, intersection.y);
  if (intersection.x > 0.0) {
    ray_start += ray_dir * instersection.x;
    ray_length -= intersection.x;
  }
  
  float costh = dot(ray_dir, light_dir);
  float phase_r = phase_rayleigh(costh);
  float phase_m = phase_mie(costh);
  
  int sample_count = 64;
  
  vec3 optical_depth, rayleigh, mie;
  optical_depth = rayleigh = mie = vec3(0);
  
  float prev_ray_time = 0.0;
  
  for (int i = 0; i < sample_count; i++) {
    float ray_time = pow(float(i) / sample_count, sample_distribution_exponent) * ray_length;
    float step_size = (ray_time - prev_ray_time);
    
    vec3 local_position = ray_start + ray_dir * ray_time;
    float local_height = atmosphere_height(local_position);
    vec3 local_density = atmosphere_density(local_height);
    
    optical_depth += local_density * step_size;
    
    vec3 view_transmittance = absorb(optical_depth);
    
    vec3 optical_depth_light = integrate_optical_depth(local_position, light_dir);
    vec3 light_transmittance = absorb(optical_depth_light);
    
    rayleigh += view_transmittance * light_transmittance * phase_r * local_density.x * step_size;
    mie += view_transmittance * light_transmittance * phase_m * local_density.y * step_size;
    
    prev_ray_time = ray_time;
  }
  
  transmittance = absorb(optical_depth);
  
  return (rayleigh * C_RAYLEIGH + mie * C_MIE) * light_color * EXPOSURE;
}

The code for the skybox itself might also be useful:

void sky() {
  float ray_length = INFINITY;
  
  if (draw_planet) {
    vec2 intersection = planet_intersection(POSITION, EYEDIR);
    if (intersection.x > 0.0) {
      ray_length = min(ray_length, intersection.x);
    }
  }
  
  vec3 transmittance;
  COLOR = integrate_scattering(POSITION, EYEDIR, ray_length, LIGHT0_DIRECTION, LIGHT0_COLOR, transmittance);
}