if (!IS_GUID_VALID(comet) ||
comet >= pv_space_n_bodies_type_comet ||
!IS_BODY_COMET_VALID(pv_space_bodies_type_comet[comet]))
return (GUID_INVALID);
//printf("light_angle=%.3f\n", light_angle);
// FIXME: allow automatic upsampling on shadow boundaries? is it really worth it?
// FIXME2: pre-multiply shadow buffer with optical buffer? hum...
/* light_visibility = pv_atmosphere_sample_sbuffer(a, &to_vert, altitude); */
light_visibility = 1.0; // FIXME: temp debug
// if this ray received some light, compute its contribution...
if (light_visibility >= DELTA_LIGHT_VISIBILITY)
{
// get the distance the light ray travels before reaching the current sample point.
pv_atmosphere_sample_obuffer(a, &light_depth, light_angle, altitude);
// now, get the distance the inscattered light will have to travel to reach us.
pv_atmosphere_sample_obuffer(a, &sample_depth, sample_angle, altitude);
// optical depth from the lightsoure to the sample point + from the sample point to the viewpoint
rayleigh_depth = light_depth.RAYLEIGH_DEPTH + (camera_depth.RAYLEIGH_DEPTH - sample_depth.RAYLEIGH_DEPTH);
mie_depth = light_depth.MIE_DEPTH + (camera_depth.MIE_DEPTH - sample_depth.MIE_DEPTH);
mie_depth = 0.0; // FIXME: temp, ignore mie scattering...
// scattering equation:
//
// Lo = Esc * ((Els * Li) * phase);
// with:
// - Lo: Outcoming light
// - Esc: Extinction from the scattering point to the camera.
// - Els: Extinction from the lightsource to the scattering point.
// - Li: Incoming light
// - phase: phase function
//
// so basically, this gives us, for rayleigh scattering:
// Lo = exp(-(Re * Dsc)) * ((exp(-(Re * Dls)) * Li) * phase);
// and as we are applying this to a sample segment,
// we will have to multiply the whole scattered light at the sample point by the sample segment's optical depth.
// Lo = exp(-(Re * Dsc)) * ((exp(-(Re * Dls)) * Li) * (p(t) * Rs) * ODs);
//
// with:
// - Re: Rayleigh extinction cross-section
// - Dsc: optical depth from the sample point to the camera.
// - Dls: optical depth from the lightsource to the sample point
// - p(t): phase function, arg: angle between the viewing ray and the ray going from the lightsource to the sample point.
// - Rs: Rayleigh scattering cross-section
// - ODs: sample segment's optical depth.
//
// we can optimize this equation to:
//
// Lo = exp(-(Re * Dsc)) * exp(-(Re * Dls)) * Li * (p(t) * Rs) * ODs;
// Lo = exp(-(Re * Dsc + Re * Dls)) * Li * (p(t) * Rs) * ODs;
//
// apply the phase function...
// NOTE: if we assume the light rays to be parallel, we can "optimize" this, and move the phase
// function outside the loop, as it will be identical for all samples.
// however, this is totally inaccurate for high turbidities and specific high frequency Mie scattering patterns...
// (rainbows, coronas, glories, etc...)
// (but maybe we shouldn't do this here... anyway... keep it as general as possible atm..)
// allright, get rayleigh and mie phase function
en_atmosphere_sample_phase_lookup(a, phase_angle, &rayleigh_phase, &mie_phase);
rayleigh_phase *= light_visibility * sample_depth.RAYLEIGH_DENSITY * sample_length;
mie_phase *= light_visibility * sample_depth.MIE_DENSITY * sample_length;
mie_phase = 0.0; // FIXME: temp, ignore mie scattering...
// FIXME: optimize / vectorize these with SSE(2) when everything will work allright
for (s = 0; s < a->s_n_samples; s++)
{
double extinction;
double phase;
// compute the extinction from the lightsource to the viewpoint
// extinction = Rdepth * Rextinction + Mdepth * Mextinction
extinction = rayleigh_depth * a->s_sigma_wre[s] + mie_depth * a->s_sigma_wme[s];
// compute the phase function
// phase = Rp(t) * Rscattering + Mp(t) * Mscattering
phase = rayleigh_phase * a->s_sigma_wrs[s] + mie_phase * a->s_sigma_wms[s];
// accumulate the inscattered spectral radiance
a->s_spectrum[s] += exp(-extinction) * a->s_lightsource[s] * phase;
}
}
// jump to the center of the next sample ray
pos.x += sample_ray.x;
pos.y += sample_ray.y;
pos.z += sample_ray.z;
}