#include <fstream>
#include <vector>
#include <string>
#include <algorithm>
#include <stdexcept>
#include <complex>
#include <stdint.h>
#include <tr1/cmath>

class Function
{
    public:
        virtual ~Function() {}

        //virtual std::string name() const = 0;
        virtual std::complex<double> eval(double x, double y) const = 0;
};

class Identity : public Function
{
    std::complex<double> eval(double x, double y) const
    {
        return std::complex<double>(x, y);
    }
};

class Propagator : public Function
{
    public:
        Propagator(double mass/*, char const *name*/) : m(mass)/*, n(name)*/ {}
        //std::string name() const { return n/*"propagator"*/; }
        std::complex<double> eval(double x, double y) const
        {
            double s = y*y-x*x;  // y is time, x is a space dimension.
            //double const m = 0.0000000001;//.0001;
            if(s >= -1.e-20 && s <= 1.e-20)
                return 1.e20;
            if(s >= 0)
                return m/(8*M_PI*std::sqrt(s))*std::complex<double>(std::tr1::cyl_bessel_j(1, m*std::sqrt(s)), std::tr1::cyl_neumann(1, m*std::sqrt(s)));
            else
                return m/(4*M_PI*M_PI*std::sqrt(-s))*std::complex<double>(0, -std::tr1::cyl_bessel_k(1, m*std::sqrt(-s)));
            return std::complex<double>(x, y);
        }

    private:
        double m;
        //char const *n;
};

struct Pixel
{
    Pixel() { col[0] = 0xFF; col[1] = 0x00; col[2] = 0xFF; }
    Pixel(double r, double g, double b)
    {
        col[0] = std::max(0, std::min(255, int(0.5+255*r)));
        col[1] = std::max(0, std::min(255, int(0.5+255*g)));
        col[2] = std::max(0, std::min(255, int(0.5+255*b)));
    }
    //Pixel(std::complex<double> const &val) { col[0] = 0; col[1] = std::max(0, std::min(255, int(127.5+63.75*std::real(val)))); col[2] = std::max(0, std::min(255, int(127.5+63.75*std::imag(val)))); }

    uint8_t col[3];
};

Pixel colourComplex(std::complex<double> const &z)
{
    double r = std::abs(z);
    double theta = std::arg(z);

    double y = std::tanh(r/2);
    double m = y*(1-y);
    double s = m*std::sin(theta);
    double c = m*std::cos(theta);

    double const q = .5;
    double const w = std::sqrt(.75);

    return Pixel(y + c, y - q*c + w*s, y - q*c - w*s);
}

class Image
{
    public:
        Image(unsigned x, unsigned y, Function const &function_, double x1_, double y1_, double x2_, double y2_)
            : sx(x), sy(y), function(&function_), /*filename(function_.name() + ".ppm"),*/ x1(x1_), y1(y1_), x2(x2_), y2(y2_)
        {
            data.resize(x*y);
//#pragma omp parallel for
            for(int py = 0; py < int(y); ++py)
                for(unsigned px = 0; px != x; ++px)
                    data[px+py*sx] = colourComplex(function->eval(x1+(x2-x1)*(px/(sx-1.)), y1+(y2-y1)*(py/(sy-1.))));
        }
        void save(std::string const &name)
        {
            std::string filename = name+".ppm";
            std::ofstream f(filename.c_str(), std::ios::binary);
            if(!f)
                throw std::runtime_error("Couldn't open file \""+filename+"\" for writing.");
            f << "P6\n" << sx << " " << sy << "\n255\n";
            f.write(data, 3*sx*sy);
        }

    private:
        unsigned sx, sy;
        Function const *function;
        //std::string filename;
        double x1, y1, x2, y2;
        std::vector<Pixel> data;
};

int main()
{
    Image identity(1001, 1001, Identity(), -2, -2, 2, 2);
    identity.save("PropagatorColours");

    Image propagator1(1001, 1001, Propagator(.2), -2, -2, 2, 2);
    propagator1.save("FeynmanPropagatorWithMass0.2");
    Image propagator2(1001, 1001, Propagator(2), -2, -2, 2, 2);
    propagator2.save("FeynmanPropagatorWithMass2");
    Image propagator3(1001, 1001, Propagator(20), -2, -2, 2, 2);
    propagator3.save("FeynmanPropagatorWithMass20");
    Image propagator4(1001, 1001, Propagator(200), -2, -2, 2, 2);
    propagator4.save("FeynmanPropagatorWithMass200");

    return 0;
}