I am trying to implement the concept of distributed ray tracing on Rayleigh scatterting to render sky a bit more realistic. The code I use is borrowed fromthis website. It takes 16 samples along the ray and 8 samples for sunlight.This works fine and I was able to generate this image using only 1 sample per pixel(not sure why the sun is so small, any idea?):I then tried to use Monte Carlo to do this. 1 sample uniformly selected along the ray and 1 sample for sunlight. The only changes I did are on four lines:line 7:  uint32_t numSamples = 1; //16line 8:  uint32_t numSamplesLight = 1; //8 //used drand48() instead of 0.5line 18: Vec3f samplePosition = orig + (tCurrent+                                segmentLength * drand48()) * dir; line 32: Vec3f samplePositionLight = samplePosition + (tCurrentLight +                                 segmentLengthLight * drand48) * sunDirection; The pdf is 1/segmentLength and the sample is taken uniformly in this distance.Here is my implementation in c# code:  public override void Render(int pxlX, int pxlY, Random rand)        {            Vector estimatedColor = new Vector();            Ray ray = Scene.Camera.GenerateRay(pxlX + rand.NextDouble(), pxlY + rand.NextDouble(), rand.NextDouble(), rand.NextDouble());            ISect isect = Scene.Trace(ray);            if (ISect.IsNull(isect))                estimatedColor = Atmosph.ComputeIncidentLight(ray, rand);            else                estimatedColor = isect.Thing.Material.Emission +                                       ComputeRadiance(isect, rand, 0);            ImageFilm.AddToPixel(pxlX, pxlY, estimatedColor/numPixelSamples);        }  public Vector ComputeIncidentLight(Ray r, Random rand)        {            r.Start.Z += earthRadius;            ISect isect = earth.Intersect(r);            if (ISect.IsNull(isect)) return Vector.ZERO;            int numSamples = 1;            int numSamplesLight = 1;            double segmentLength = isect.Dist / numSamples;            double tCurrent = 0;            Vector sumR = Vector.ZERO;            Vector sumM = Vector.ZERO;// mie and rayleigh contribution             double opticalDepthR = 0, opticalDepthM = 0;            double mu = r.Dir.Dot(sunDirection); // mu in the paper which is the cosine of the angle between the sun direction and the ray direction             double phaseR = 3d / (16d * Math.PI) * (1 + mu * mu);            double g = 0.76f;            double phaseM = 3d / (8d * Math.PI) * ((1d - g * g) * (1d + mu * mu)) / ((2d + g * g) * Math.Pow(1d + g * g - 2d * g * mu, 1.5d));            for (int i = 0; i < numSamples; ++i)            {                Vector samplePosition = r.Start + r.Dir * (tCurrent + segmentLength * rand.NextDouble());                double height = samplePosition.Length() - earthRadius;                // compute optical depth for light                double hr = Math.Exp(-height / Hr) * segmentLength;                double hm = Math.Exp(-height / Hm) * segmentLength;                opticalDepthR += hr;                opticalDepthM += hm;                // light optical depth                ISect isectSun = earth.Intersect(new Ray(samplePosition, sunDirection));                double segmentLengthLight = isectSun.Dist/ numSamplesLight, tCurrentLight = 0;                double opticalDepthLightR = 0, opticalDepthLightM = 0;                int j;                for (j = 0; j < numSamplesLight; ++j)                {                    Vector samplePositionLight = samplePosition + sunDirection * (tCurrentLight + segmentLengthLight * rand.NextDouble());                    double heightLight = samplePositionLight.Length() - earthRadius;                    if (heightLight < 0) break;                    opticalDepthLightR += Math.Exp(-heightLight / Hr) * segmentLengthLight;                    opticalDepthLightM += Math.Exp(-heightLight / Hm) * segmentLengthLight;                    tCurrentLight += segmentLengthLight;                }                if (j == numSamplesLight)                {                    Vector tau = betaR * (opticalDepthR + opticalDepthLightR) + betaM * 1.1f * (opticalDepthM + opticalDepthLightM);                    Vector attenuation = new Vector(Math.Exp(-tau.X), Math.Exp(-tau.Y), Math.Exp(-tau.Z));                    sumR += attenuation * hr;                    sumM += attenuation * hm;                }                tCurrent += segmentLength;            }            // We use a magic number here for the intensity of the sun (20). We will make it more            // scientific in a future revision of this lesson/code            return (sumR * betaR * phaseR + sumM * betaM * phaseM) * sunIrradiance;        }The below image rendered with 400 samples per pixel, however is different from the above:Could anyone shed some light where I've got it wrong?Thanks,
Using Monte carlo on Rayleigh scattering
pathtracing;global illumination;distribution;volumetric scattering