namespace BartenContrastSensitivityFunction 
{
    // Modell des Auges, abgetippt und durch!
    const double k      = 3.0;
    const double sigma0 = 0.5000/60;// arc grad
    const double Cab    = 0.0800/60;// arc grad/mm
    const double T      = 0.1;      // sec
    const double Xmax   = 12;       // arc grad
    const double Nmax   = 15;       // cycles
    const double eta    = 0.03;
    const double Phi0   = 3e-8;     // sec deg²
    const double u0     = 7;        // cycles/grad
    const double p      = 1.285e6;  // photonen / sec / grad² / Td

    const double sigma (const double d)
    {
        return sqrt (sigma0 * sigma0 + Cab * Cab * d * d);
    }
    const double d (const double L, const double X0)
    {
        return 5 - 3 * ::tanh (0.4 * log (L*X0*X0/40/40));
    }
    const double E (const double d, const double L)
    {
        return 3.1415926 * d*d / 4 * L * (1 - ::pow(d/9.7,2) - ::pow(d/12.4,4));
    }
    const double Mopt (const double u, const double d)
    {
        const double s = sigma (d);
        return ::exp (-2*3.1415926*3.1415926*s*s*u*u);
    }
    // Reziproke der min. sichtbare Modulation bei Helligkeit L und Ortsfrequenz u
    const double CSF (const double L, const double u) 
    {
        const double X0 = 13.0/u;
        const double _d = d(L, X0);
        const double Above = Mopt(u, _d) / k;
        const double Below = 2/T *(1/X0/X0 + 1/Xmax/Xmax + u*u/Nmax/Nmax) * (1/eta/p/E(_d,L) + Phi0 / (1-exp(-u*u/u0/u0)));
        return Above / ::sqrt(Below);       // Rückgabewert 100 = 1% Modulation sichtbar
    }
    // Reziproke der min. sichtbare Modulation bei Helligkeit L, optimale Ortsfrequenz
    const double CSFmax (const double L) 
    {
        double retmax = 0.0;
        for (int i = -2000; i <= 2000; i++)
        {
            const double C = CSF (L, ::pow(10, 0.001*i));
            if (C > retmax) retmax = C;
        }
        return retmax;
    }
    // SMPTE-definierte EOTF, N = 0...1 = 0...4095 bei 12 bit, Rückgabe 0...1 = 0...10000 cd/m²
    const double EOTF (const double N)
    {
        const double  m1  = 1305.0 / 8192;
        const double  m2  = 2523.0 /   32;
        const double  rm1 = 8192 / 1305.0;
        const double  rm2 =   32 / 2523.0;
        const double  c1  = 3424.0 / 4096;
        const double  c2  = 2413.0 /  128;
        const double  c3  =  299.0 /   16;

        return ::pow ( (::pow(N,rm2) - c1) / (c2 - c3*::pow(N,rm2) ), rm1);
    }

    void CalTable ()
    {
        double M = 10000;                   // maximale Helligkeit für den größten Wert in cd/m²
        double mod = 1;

        for (int i = 4095; i >= 0; i--)     // Helligkeitswerte 0...4095
        {
            for (int j = 0; j < 4; j++)     // bestimme iterativ die Modulationstiefe für den Wert zwischen L_i und L_i-1
            {
                const double x = 0.6621 / BartenContrastSensitivityFunction::CSFmax (M/sqrt(mod));
                mod = (1+x)/(1-x);   
            }
            ::printf ("%4u %16.9f %16.9f (%7.3f)\n", 
                i,                                                         // Helligkeitsindex
                M,                                                         // Helligkeit aus BartenContrastSensitivityFunction
                10000 * BartenContrastSensitivityFunction::EOTF(i/4095.0), // Helligkeit nach SMPTE ST-2084
                100*(mod-1));                                              // erkennbare Modulationstiefe 
            M /= mod;
        }
    }
}