/*
    Read

        http://z0b.kapsi.fi/snippets.php

    before using this code. Thank you.
*/

/*
Creates a "shade map". Takes a sun (or any light, actually) position in
latitude/longitude coordinates and the image size, then creates an 8-bit image
showing light/dark areas and smooth transitions between them.

Latitude (Y) is between -90.0 and +90.0 and longitude (X) is between -180.0 and +180.0.

The 'minShade' is the minimum light level. It can be used to prevent "too dark"
dark areas. It must be between 0 and 1, but values > 0.5 usually don't look good.

The output is an 8-bit transition image. It can be used directly to blend two
images, for example. 255 is fully lit, 0 is fully dark.

Returns false if fails, true if OK.
*/
static bool shadeMap(const float sunLatitude, const float sunLongitude, // sun (light) position
                     const uint32_t width, const uint32_t height,       // output size
                     std::vector<uint8_t> &output,                      // output
                     const float minShade = 0.0f)                       // darkest shade
{
    if (width < 1 || height < 1 || minShade < 0.0f || minShade > 1.0f)
        return false;

    // converts from degrees to radians (always multiply)
    static const float kToRad = M_PI / 180.0f;

    // This specifies how far away the light source is. If it's too small, the
    // transitions become too fuzzy. If it's too large, the transitions become
    // binary (dark or light, nothing between). We could use real Earth and Sun
    // diameters and distance between the two, then fudge the numbers, but why
    // bother? The distance varies anyway and the Sun isn't a point-light source.
    // Experiments have shown that values between 5 and 15 work well and
    // produces nice and reasonable results. Feel free to adjust.
    static const float kSunDistance = 10.0f;

    // use spherical mapping to project the 2D point into 3D space
    const float sunPos[3] =
    {
        sinf(sunLongitude * kToRad) * cosf(sunLatitude * kToRad) * kSunDistance,
        sinf( sunLatitude * kToRad)                              * kSunDistance,
        cosf(sunLongitude * kToRad) * cosf(sunLatitude * kToRad) * kSunDistance
    };

    // We'll project the 2D map points into 3D space in the same way as the sun
    // vector. Because the X is same for every column on every row, we can
    // precalculate sin(x) and cos(x) values. It saves us a LOT of repetitive
    // calculations.
    float sinX[width],
          cosX[width];

    for (uint32_t x = 0; x < width; x++) {
        // convert X to longitude, then to radians
        const float v = ((float)x / (float)width - 0.5f) * 360.0f * kToRad;

        sinX[x] = sinf(v);
        cosX[x] = cosf(v);
    }

    output.resize(width * height);
    uint8_t *out = &output[0];

    for (uint32_t y = 0; y < height; y++) {
        // Again, use spherical mapping to project the 2D points,
        // this gives us surface normals. We're doing a lot of dot
        // products and most of the values don't change, so why not
        // precompute them? We've precomputed the sin(x) and cos(x)
        // values already.

        // Convert Y to latitude, then to radians. Project the constants
        // (for this row) into 3D space and precalculate every dot
        // product factor we can.
        const float lat = ((float)y / (float)height - 0.5f) * -180.0f * kToRad,
                    latS = sinf(lat) * sunPos[1],
                    latC0 = cosf(lat) * sunPos[0],
                    latC2 = cosf(lat) * sunPos[2];

        for (uint32_t x = 0; x < width; x++) {
            // Dot product the surface normal with the sun vector.
            // It might not look like it, but the next line actually
            // combines a dot product and sunPos[] -style projection.
            // We just have precalculated most of the numbers already.
            const float dot = sinX[x] * latC0 + cosX[x] * latC2 + latS;

            // negative dot products mean the point is "behind" the
            // sphere surface and thus receives no light at all
            *out++ = (unsigned)(min(max(dot, minShade), 1.0f) * 255.0f) & 0xFF;
        }
    }

    return true;
}