///
/// Conversion from Luminance Factor to Munsell Value with High Accuracy
/// Author: Mikhail Semenov. 2016-04-03
///

#include <iostream>
#include <cmath>

double MunsellValue(double Y)
{    
    double a = log(Y + 1) * 0.25;
    double p;
    if (Y < 1.18)
    {
        p = (((((16191949.3 * a - 7941654.8) * a + 510539.3) * a + 91655.6) * a + 89334.3) * a + 33567.0) * a;
    }
    else
    {
        p = (((((((((393050 * a - 2768832.5) * a + 8676160.0) * a - 15985982.5) * a + 19187963.5) * a
            - 15642658.5) * a + 8737970.2) * a - 3237397.0) * a + 745292.7) * a - 30440.7) * a + 2488.1;
    }
    return p * 0.0001;
}

double Luminance(double V)
{
    return (((((0.0008404*V) - 0.021009)*V+ 0.23951)*V - 0.23111)*V+ 1.2219)*V;        
}

double Luminance100(double V)
{
    return (((((0.00081939*V) - 0.020484)*V + 0.23352)*V - 0.22533)*V + 1.1914)*V;    
}

int main()
{
    double max_error = 0;
    double max_Y = 0;
    double max_V = 0;
    double V = 0.0;
    double h = 0.001;
    
    for (int i = 0; i <= 10000; i++)
    {
        double Y = Luminance100(V);
        double V2 = MunsellValue(Y);
        double error = fabs(V - V2);
        if (error > max_error)
        {
            max_error = error;
            max_Y = Y;
            max_V = V;
        }
        V += h;
    }

    std::cout << "Y100 to V. max error: " << max_error << " V: " << max_V << " Y: " << max_Y << std::endl;
    max_error = 0;
    max_Y = 0;
    max_V = 0;
    V = 0.0;
    h = 0.001;
    
    for (int i = 0; i <= 10000; i++)    
    {
        double Y = Luminance(V);
        double V2 = MunsellValue(Y/1.02568 /*0.974968*/);
        double error = fabs(V - V2);
        if (error > max_error)
        {
            max_error = error;
            max_Y = Y;
            max_V = V;
        }
        V += h;
    }

    std::cout << "Y to V. max error: " << max_error << " V: " << max_V << " Y: " << max_Y << std::endl;    
}