builtin_data/_sky.glsl

   1 struct AtmosphericEvents
   2 {
   3         vec3 rayleigh_scatter;
   4         vec3 mie_scatter;
   5         vec3 mie_absorb;
   6         vec3 ozone_absorb;
   7 };
   8
   9 uniform Atmosphere
  10 {
  11         AtmosphericEvents events;
  12         float rayleigh_density_decay;
  13         float mie_density_decay;
  14         float ozone_band_center;
  15         float ozone_band_extent;
  16         float planet_radius;
  17         float atmosphere_thickness;
  18         vec3 ground_albedo;
  19         int n_steps;
  20 };
  21
  22 uniform View
  23 {
  24         float view_height;
  25         vec4 light_color;
  26         vec3 light_dir;
  27 };
  28
  29 struct OpticalPathInfo
  30 {
  31         vec3 optical_depth;
  32         vec3 luminance;
  33 };
  34
  35 const float pi = 3.1415926535;
  36 const float mie_asymmetry = 0.8;
  37
  38 uniform sampler2D transmittance_lookup;
  39
  40 vec3 rayleigh_density(vec3 base, float height)
  41 {
  42         return base*exp(height/rayleigh_density_decay);
  43 }
  44
  45 float rayleigh_phase(float cos_theta)
  46 {
  47         return 3.0*(1.0+cos_theta*cos_theta)/(16.0*pi);
  48 }
  49
  50 vec3 mie_density(vec3 base, float height)
  51 {
  52         return base*exp(height/mie_density_decay);
  53 }
  54
  55 float mie_phase(float cos_theta)
  56 {
  57         float g = mie_asymmetry;
  58         float num = (1.0-g*g)*(1.0+cos_theta*cos_theta);
  59         float denom = (2.0+g*g)*pow(1.0+g*g-2.0*g*cos_theta, 1.5);
  60         return 3.0/(8.0*pi)*num/denom;
  61 }
  62
  63 vec3 ozone_density(vec3 base, float height)
  64 {
  65         return base*max(1.0-abs(height-ozone_band_center)/ozone_band_extent, 0.0);
  66 }
  67
  68 AtmosphericEvents calculate_events(float height)
  69 {
  70         AtmosphericEvents ev;
  71         ev.rayleigh_scatter = rayleigh_density(events.rayleigh_scatter, height);
  72         ev.mie_scatter = mie_density(events.mie_scatter, height);
  73         ev.mie_absorb = mie_density(events.mie_absorb, height);
  74         ev.ozone_absorb = ozone_density(events.ozone_absorb, height);
  75         return ev;
  76 }
  77
  78 vec3 total_extinction(AtmosphericEvents ev)
  79 {
  80         return ev.rayleigh_scatter+ev.mie_scatter+ev.mie_absorb+ev.ozone_absorb;
  81 }
  82
  83 float ray_sphere_intersect(vec3 ray_start, vec3 ray_dir, vec3 sphere_center, float sphere_radius)
  84 {
  85         float t = dot(sphere_center-ray_start, ray_dir);
  86         vec3 nearest = ray_start+t*ray_dir-sphere_center;
  87         float d_sq = dot(nearest, nearest);
  88         float r_sq = sphere_radius*sphere_radius;
  89         if(d_sq>r_sq)
  90                 return -1.0;
  91
  92         float offset = sqrt(r_sq-d_sq);
  93         if(offset<t)
  94                 return t-offset;
  95         else if(offset>-t)
  96                 return t+offset;
  97         else
  98                 return -1.0;
  99 }
 100
 101 #pragma MSP stage(fragment)
 102 OpticalPathInfo raymarch_path(float start_height, vec3 look_dir)
 103 {
 104         float cos_theta = dot(look_dir, light_dir);
 105         float p_rayleigh = rayleigh_phase(cos_theta);
 106         float p_mie = mie_phase(cos_theta);
 107
 108         vec3 planet_center = vec3(0.0, 0.0, -planet_radius);
 109         vec3 pos = vec3(0.0, 0.0, start_height);
 110         vec3 path_luminance = vec3(0.0);
 111         vec3 path_extinction = vec3(0.0);
 112
 113         float ground_t = ray_sphere_intersect(pos, look_dir, planet_center, planet_radius);
 114         float space_t = ray_sphere_intersect(pos, look_dir, planet_center, planet_radius+atmosphere_thickness);
 115         float ray_length = (ground_t>0.0 ? ground_t : space_t);
 116         float step_size = ray_length/n_steps;
 117
 118         for(int i=0; i<=n_steps; ++i)
 119         {
 120                 vec3 from_center = pos-planet_center;
 121                 float height = length(from_center);
 122                 float light_z = dot(from_center/height, light_dir);
 123                 height -= planet_radius;
 124
 125                 AtmosphericEvents ev = calculate_events(height);
 126                 vec3 transmittance = exp(-path_extinction);
 127                 vec3 in_transmittance = texture(transmittance_lookup, vec2(sqrt(height/atmosphere_thickness), light_z)).rgb;
 128                 vec3 in_luminance = (ev.rayleigh_scatter*p_rayleigh+ev.mie_scatter*p_mie)*step_size;
 129                 if(i==n_steps && ground_t>0.0)
 130                         in_luminance += ground_albedo*light_z/pi;
 131                 path_luminance += transmittance*in_transmittance*in_luminance;
 132
 133                 path_extinction += total_extinction(ev)*step_size;
 134                 pos += look_dir*step_size;
 135         }
 136
 137         return OpticalPathInfo(path_extinction, path_luminance);
 138 }