// Mesa is an open-source OpenGL implementation.
// OSMesa is its off-screen driver, used for rendering into a memory buffer
// rather than through a graphics card. This is great for making OpenGL work
// pretty much on anything, including DOSBox on its emulated basic VGA hardware.
//
// You can download my DJGPP-compiled version of OSMesa at:
// http://bisqwit.iki.fi/jutut/kuvat/programming_examples/djgpp_mesa.zip
#define PRERENDER   0
#define USE_FLOAT   0
#define USE_HDRI    0
#define LIGHTWEIGHT 0
#define REPLAY      0
#define TRUECOLOR_ONLY 0

#if PRERENDER < 2
#include <osmesa.h>
#include <GL/glu.h>
#else
typedef float GLfloat;
#include <cmath>
#include <cstdio>
#endif
#include <sys/time.h>
#include <unistd.h>

// Standard C++ includes:
#include <algorithm>
#include <vector>

#if PRERENDER < 2
// DJGPP-specific include files, for accessing the screen & keyboard etc.:
#include <pc.h> // For outportb
#include <go32.h> // For _dos_ds
#include <conio.h> // For kbhit, getch, textmode
#include <sys/farptr.h> // For _farsetsel and _farnspokeb
#include <dpmi.h> // For __dpmi_int (mouse access)

#if USE_FLOAT && USE_HDRI
#include "fattal02/tmo_fattal02.cpp"
#include "fattal02/pde.cpp"
#endif

#if USE_FLOAT
const double BrightnessSteps[] = {16.0, 4.0, 2.0, 1.0, 0.5, 0.1};//0.25, 0.5, 1.0, 2.0};
const unsigned NumBrightnessSteps = sizeof(BrightnessSteps) / sizeof(*BrightnessSteps);
const double TargetBrightness     = 2.0;
bool TexturesInstalled = false;
unsigned CurrentBrightnessStep;
#endif

namespace PC
{
#if LIGHTWEIGHT < 1
    const bool TemporalDithering = true;
#endif
    const unsigned DitheringBits = 6;
    const unsigned W = 320, H = 200;
    //const unsigned R = 6, G = 7, B = 6; // 6*7*6 regular palette (252 colors)
    //const unsigned R = 6, G = 8, B = 5; // 6*8*5 regular palette (240 colors)
    const unsigned R = 7, G = 9, B = 4; // 7*9*4 regular palette (252 colors)
    //const unsigned R = 8, G = 8, B = 4; // 8*8*4 regular palette (256 colors)
    //const unsigned R = 3, G = 3, B = 3;
    //const unsigned R = 2, G = 2, B = 2; // 2*2*2 regular palette (8 colors)
    unsigned ImageBuffer[W*H];
#if USE_FLOAT
    GLfloat  ImageBufferFloat[W*H*3], ImageSumFloat[W*H*3], ImageSumBrightness[W*H*3], ImageSample[W*H*3];
#endif
    const double PaletteGamma = 1.5;            // Apply this gamma to palette
    //const double PaletteGamma = 1.0;            // Apply this gamma to palette
    const double DitherGamma = 2.0/PaletteGamma;// Apply this gamma to dithering
    unsigned char ColorConvert[3][256][256], Dither8x8[8][8];

    void Init() // Initialize graphics
    {
        // Create bayer 8x8 dithering matrix.
        for(unsigned y=0; y<8; ++y)
            for(unsigned x=0; x<8; ++x)
                Dither8x8[y][x] =
                    ((x  ) & 4)/4u + ((x  ) & 2)*2u + ((x  ) & 1)*16u
                  + ((x^y) & 4)/2u + ((x^y) & 2)*4u + ((x^y) & 1)*32u;
        // Create gamma-corrected look-up tables for dithering.
        double dtab[256], ptab[256];
        for(unsigned n=0; n<256; ++n)
        {
            dtab[n] = (255.0/256.0) - std::pow(n/256.0, 1/DitherGamma);
            ptab[n] = std::pow( n/255.0, 1.0 / PaletteGamma);
         #if LIGHTWEIGHT >= 2
            dtab[n] = 0.5;
         #endif
        }
        for(unsigned n=0; n<256; ++n)
            for(unsigned d=0; d<256; ++d)
            {
                ColorConvert[0][n][d] =     std::min(B-1, (unsigned)(ptab[n]*(B-1) + dtab[d]));
                ColorConvert[1][n][d] =   B*std::min(G-1, (unsigned)(ptab[n]*(G-1) + dtab[d]));
                ColorConvert[2][n][d] = G*B*std::min(R-1, (unsigned)(ptab[n]*(R-1) + dtab[d]));
            }

        // Set VGA mode 13h (320x200, 256 colors)
        textmode(0x13); // Will set graphics mode despite the name.
        // Set up regular palette as configured earlier.
        // However, bias the colors towards darker ones in an exponential curve.
        outportb(0x3C8, 0);
        for(unsigned color=0; color< R*G*B; ++color)
        {
            outportb(0x3C9, std::pow(((color/(B*G))%R)*1./(R-1), PaletteGamma) *63);
            outportb(0x3C9, std::pow(((color/    B)%G)*1./(G-1), PaletteGamma) *63);
            outportb(0x3C9, std::pow(((color      )%B)*1./(B-1), PaletteGamma) *63);
        }

        #ifdef __DJGPP__
        __dpmi_regs regs = { };
        regs.x.ax = 0;
        __dpmi_int(0x33, &regs); // Initialize mouse
        #endif
    }

#if USE_FLOAT
    void UpdateFloatPixmap()
    {
        double brightness = 1.0 / BrightnessSteps[CurrentBrightnessStep];
        if(CurrentBrightnessStep == 0)
        {
            // Before
            for(unsigned p=0; p<W*H*3; ++p)
            {
                ImageSumFloat[p]      = 0.0;
                ImageSumBrightness[p] = 0.0;
            }
        }
        for(unsigned p=0; p<W*H*3; ++p)
        {
            float v = ImageBufferFloat[p];
            if(v > 0.0 && v < 1.0)
            {
                ImageSumFloat[p]      += v * brightness;
                ImageSumBrightness[p] += brightness;
            }
            ImageSample[p] = v*brightness;
        }
        if(CurrentBrightnessStep == NumBrightnessSteps-1)
        {
            // After
            for(unsigned p=0; p<W*H*3; ++p)
                if(ImageSumBrightness[p] == 0.0)
                    ImageSumFloat[p] = ImageSample[p];
                else
                    ImageSumFloat[p] *= TargetBrightness / ImageSumBrightness[p];
        }
    }
    void ConvertFloatPixmap()
    {
    #if USE_HDRI
        float Y[W*H], L[W*H];
        double X[W*H], Z[W*H];
        #pragma omp parallel for
        for(unsigned p=0; p<W*H; ++p)
        {
            double r=ImageSumFloat[p*3+0];//+1e-15;
            double g=ImageSumFloat[p*3+1];//+1e-15;
            double b=ImageSumFloat[p*3+2];//+1e-15;
            X[p]=0.4124240 *r+0.3575790 *g+0.1804640 *b;
            Y[p]=0.2126560f*r+0.7151580f*g+0.0721856f*b;
            Z[p]=0.0193324 *r+0.1191930 *g+0.9504440 *b;
        }
        float alpha=1.0f, beta=0.9f, sat=1.0f;
        tmo_fattal02(W,H, Y, L, alpha, beta, alpha*0.01);
    #endif
        #pragma omp parallel for
        for(unsigned p=0; p<W*H; ++p)
        {
        #if USE_HDRI
            float y = L[p]; double l = y;
            float x = std::pow(X[p]/Y[p], sat) * y;
            float z = std::pow(Z[p]/Y[p], sat) * y;
            double rgbf[3] = {
                ( 3.2404542*x-1.5371385*y-0.4985314*z + std::pow(ImageSumFloat[p*3+0],2.)*2.8) * 0.5,
                (-0.9692660*x+1.8760108*y+0.0415560*z + std::pow(ImageSumFloat[p*3+1],2.)*2.8) * 0.5,
                ( 0.0556434*x-0.2040259*y+1.0572252*z + std::pow(ImageSumFloat[p*3+2],2.)*2.8) * 0.5 };
        #else
            double rgbf[3] = {
                ImageSumFloat[p*3+0],
                ImageSumFloat[p*3+1],
                ImageSumFloat[p*3+2] };
            double l = rgbf[0]*0.299 + rgbf[1]*0.587 + rgbf[2]*0.114;
        #endif
            double s=1.0;
            if(l >= 1.0)      ImageBuffer[p] = 0xFFFFFF;
            else if(l <= 0.0) ImageBuffer[p] = 0x000000;
            else
            {
                for(unsigned n=0; n<3; ++n)
                    if(rgbf[n] > 1.0) s = std::min(s, (l-1.0) / (l-rgbf[n]));
                    else if(rgbf[n] < 0.0) s = std::min(s, l  / (l-rgbf[n]));
                if(s != 1.0)
                    for(unsigned n=0; n<3; ++n) rgbf[n] = (rgbf[n] - l) * s + l;
                unsigned r = std::max(0.0, std::min(255.0, rgbf[2]*255.0));
                unsigned g = std::max(0.0, std::min(255.0, rgbf[1]*255.0));
                unsigned b = std::max(0.0, std::min(255.0, rgbf[0]*255.0));
                ImageBuffer[p] = (r<<16)|(g<<8)|b;
            }
        }
    }
#endif
    void Render() // Update the displayed screen
    {
    #if !TRUECOLOR_ONLY
        static unsigned f=0; ++f; // Frame number
        _farsetsel(_dos_ds);
        #pragma omp parallel for
        for(unsigned y=0; y<H; ++y)
            for(unsigned p=y*W, x=0; x<W; ++x, ++p)
            {
                unsigned rgb = ImageBuffer[p], d = Dither8x8[y&7][x&7]; // 0..63
                d &= (0x3F - (0x3F >> DitheringBits));
#if LIGHTWEIGHT < 1
                if(!TemporalDithering)
#endif
                    d *= 4; // No temporal dithering
#if LIGHTWEIGHT < 1
                else // Do temporal dithering
                {
                    d += ((f^y^(x&1)*2u ^ (x&2)/2u) & 3) << 6;
                    // ^ This step is optional. It is a tradeoff:
                    // Pros: it improves the spatial color resolution.
                    // Cons: it causes flickering. The higher the framerate,
                    //       the less noticeable the flickering is.
                }
#endif
                _farnspokeb(0xA0000 + p,
                    ColorConvert[0][(rgb >> 0) & 0xFF][d]
                  + ColorConvert[1][(rgb >> 8) & 0xFF][d]
                  + ColorConvert[2][(rgb >>16) & 0xFF][d]);
            }
    #endif
    #ifndef __DJGPP__
        updatescreen(ImageBuffer);
    #endif
    }
    void Close() // End graphics
    {
        textmode(C80); // set textmode again
    }
}
#endif

#define drand(m) ((std::rand()%1000-500)*5e-2*m)

// Vector mathematics
#define vdot(a,b)   (a[0]*b[0] + a[1]*b[1] + a[2]*b[2])
#define vcross(a,b) {a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]}
#define vsub(a,b)   {a[0]-b[0], a[1]-b[1], a[2]-b[2]}
#define vlen(a)     std::sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2])
#define vaddm(t,b,m) for(unsigned a=0; a<3; ++a) t[a] += b[a] * (m)

// List of the walls that comprise the level map. In no particular order...
static const struct maptype
{
    GLfloat normal[3];  // Normal = where the wall "looks" towards.
    GLfloat p[4][3];    // corner points (used for collision tests)
} map[] =
{
    {{ 0, 0, 1}, {{ 6, 2,1},{ 6, 5,1},{ 1, 5,1},{ 1, 2,1}} }, // 0
    {{ 0, 0, 1}, {{ 8, 4,1},{ 8, 5,1},{ 6, 5,1},{ 6, 4,1}} }, // 1
    {{ 0, 0, 1}, {{24, 4,1},{24, 5,1},{16, 5,1},{16, 4,1}} }, // 2
    {{ 0, 0, 1}, {{24, 3,1},{24, 4,1},{14, 4,1},{14, 3,1}} }, // 3
    {{ 0, 0, 1}, {{24, 5,1},{24, 6,1},{17, 6,1},{17, 5,1}} }, // 4
    {{ 0, 0, 1}, {{24, 6,1},{24, 9,1},{13, 9,1},{13, 6,1}} }, // 5
    {{ 0,-1, 0}, {{24, 9,1},{24, 9,3},{16, 9,3},{16, 9,1}} }, // 6
    {{ 0,-1, 0}, {{24, 9,4},{24, 9,7},{16, 9,7},{16, 9,4}} }, // 7
    {{ 0,-1, 0}, {{24, 9,3},{24, 9,4},{17, 9,4},{17, 9,3}} }, // 8
    {{ 0,-1, 0}, {{24, 9,7},{24, 9,8},{17, 9,8},{17, 9,7}} }, // 9
    {{ 0,-1, 0}, {{16, 9,1},{16, 9,8},{13, 9,8},{13, 9,1}} }, // 10
    {{ 0, 0,-1}, {{16, 9,8},{16, 8,8},{13, 8,8},{13, 9,8}} }, // 11
    {{ 0, 0,-1}, {{16, 7,8},{16, 6,8},{13, 6,8},{13, 7,8}} }, // 12
    {{ 0, 0,-1}, {{16, 8,8},{16, 7,8},{11, 7,8},{11, 8,8}} }, // 13
    {{ 0, 0,-1}, {{16, 6,8},{16, 2,8},{ 1, 2,8},{ 1, 6,8}} }, // 14
    {{ 0, 0,-1}, {{11, 4,5},{11, 2,5},{ 8, 2,5},{ 8, 4,5}} }, // 15
    {{ 0,-1, 0}, {{11, 4,1},{11, 4,5},{ 8, 4,5},{ 8, 4,1}} }, // 16
    {{ 0, 1, 0}, {{ 6, 4,1},{ 6, 4,2},{ 7, 4,2},{ 7, 4,1}} }, // 17
    {{ 1, 0, 0}, {{ 7, 4,1},{ 7, 4,2},{ 7, 2,2},{ 7, 2,1}} }, // 18
    {{ 0, 1, 0}, {{ 7, 2,1},{ 7, 2,2},{14, 2,2},{14, 2,1}} }, // 19
    {{ 0, 1, 0}, {{ 8, 2,2},{ 8, 2,5},{14, 2,5},{14, 2,2}} }, // 20
    {{ 0, 1, 0}, {{ 8, 2,6},{ 8, 2,8},{16, 2,8},{16, 2,6}} }, // 21
    {{ 0, 1, 0}, {{ 1, 2,1},{ 1, 2,8},{ 6, 2,8},{ 6, 2,1}} }, // 22
    {{ 1, 0, 0}, {{ 1, 5,1},{ 1, 5,8},{ 1, 2,8},{ 1, 2,1}} }, // 23
    {{ 1, 0, 0}, {{ 1, 6,4},{ 1, 6,8},{ 1, 5,8},{ 1, 5,4}} }, // 24
    {{ 0,-1, 0}, {{ 7, 6,4},{ 7, 6,8},{ 1, 6,8},{ 1, 6,4}} }, // 25
    {{ 0,-1, 0}, {{ 8, 6,4},{ 8, 6,6},{ 7, 6,6},{ 7, 6,4}} }, // 26
    {{ 0, 1, 0}, {{16, 6,4},{16, 6,7},{17, 6,7},{17, 6,4}} }, // 27
    {{ 0, 1, 0}, {{16, 6,1},{16, 6,3},{17, 6,3},{17, 6,1}} }, // 28
    {{ 1, 0, 0}, {{17, 6,1},{17, 6,3},{17, 5,3},{17, 5,1}} }, // 29
    {{ 1, 0, 0}, {{17, 6,4},{17, 6,7},{17, 5,7},{17, 5,4}} }, // 30
    {{ 0,-1, 0}, {{17, 5,4},{17, 5,7},{16, 5,7},{16, 5,4}} }, // 31
    {{ 0,-1, 0}, {{17, 5,1},{17, 5,3},{16, 5,3},{16, 5,1}} }, // 32
    {{ 0, 0,-1}, {{17, 5,3},{17, 2,3},{16, 2,3},{16, 5,3}} }, // 33
    {{ 0, 0,-1}, {{17, 5,7},{17, 2,7},{16, 2,7},{16, 5,7}} }, // 34
    {{ 0, 0,-1}, {{17, 9,7},{17, 6,7},{16, 6,7},{16, 9,7}} }, // 35
    {{ 0, 0,-1}, {{17, 9,3},{17, 6,3},{16, 6,3},{16, 9,3}} }, // 36
    {{ 0, 0, 1}, {{17, 2,4},{17, 5,4},{16, 5,4},{16, 2,4}} }, // 37
    {{ 0, 0, 1}, {{17, 6,4},{17, 9,4},{16, 9,4},{16, 6,4}} }, // 38
    {{ 0, 0, 1}, {{15, 2,4},{15, 3,4},{14, 3,4},{14, 2,4}} }, // 39
    {{ 0, 1, 0}, {{15, 2,1},{15, 2,4},{16, 2,4},{16, 2,1}} }, // 40
    {{ 1, 0, 0}, {{15, 3,1},{15, 3,4},{15, 2,4},{15, 2,1}} }, // 41
    {{ 0, 1, 0}, {{14, 3,1},{14, 3,4},{15, 3,4},{15, 3,1}} }, // 42
    {{-1, 0, 0}, {{14, 2,1},{14, 2,4},{14, 3,4},{14, 3,1}} }, // 43
    {{-1, 0, 0}, {{16, 2,3},{16, 2,4},{16, 9,4},{16, 9,3}} }, // 44
    {{-1, 0, 0}, {{16, 2,7},{16, 2,8},{16, 9,8},{16, 9,7}} }, // 45
    {{-1, 0, 0}, {{16, 5,1},{16, 5,3},{16, 6,3},{16, 6,1}} }, // 46
    {{-1, 0, 0}, {{16, 5,4},{16, 5,7},{16, 6,7},{16, 6,4}} }, // 47
    {{-1, 0, 0}, {{ 8, 5,4},{ 8, 5,5},{ 8, 6,5},{ 8, 6,4}} }, // 48
    {{-1, 0, 0}, {{ 8, 2,5},{ 8, 2,6},{ 8, 6,6},{ 8, 6,5}} }, // 49
    {{-1, 0, 0}, {{ 8, 4,1},{ 8, 4,5},{ 8, 5,5},{ 8, 5,1}} }, // 50
    {{ 0,-1, 0}, {{ 8, 5,1},{ 8, 5,4},{ 1, 5,4},{ 1, 5,1}} }, // 51
    {{ 0, 0, 1}, {{ 8, 5,4},{ 8, 6,4},{ 1, 6,4},{ 1, 5,4}} }, // 52
    {{ 0, 0, 1}, {{11, 2,6},{11, 6,6},{ 8, 6,6},{ 8, 2,6}} }, // 53
    {{ 0, 0, 1}, {{11, 6,6},{11,18,6},{ 7,18,6},{ 7, 6,6}} }, // 54
    {{ 0, 0, 1}, {{13, 7,6},{13, 8,6},{11, 8,6},{11, 7,6}} }, // 55
    {{ 0,-1, 0}, {{13, 8,6},{13, 8,8},{11, 8,8},{11, 8,6}} }, // 56
    {{ 0,-1, 0}, {{13, 6,1},{13, 6,8},{11, 6,8},{11, 6,1}} }, // 57
    {{ 1, 0, 0}, {{13, 7,1},{13, 7,8},{13, 6,8},{13, 6,1}} }, // 58
    {{ 1, 0, 0}, {{13, 9,1},{13, 9,6},{13, 7,6},{13, 7,1}} }, // 59
    {{ 1, 0, 0}, {{13, 9,6},{13, 9,8},{13, 8,8},{13, 8,6}} }, // 60
    {{ 1, 0, 0}, {{17, 9,7},{17, 9,8},{17, 4,8},{17, 4,7}} }, // 61
    {{ 1, 0, 0}, {{17, 9,3},{17, 9,4},{17, 2,4},{17, 2,3}} }, // 62
    {{ 0, 1, 0}, {{17, 2,3},{17, 2,4},{20, 2,4},{20, 2,3}} }, // 63
    {{ 0, 1, 0}, {{11, 2,5},{11, 2,6},{20, 2,6},{20, 2,5}} }, // 64
    {{ 0, 1, 0}, {{14, 2,4},{14, 2,5},{20, 2,5},{20, 2,4}} }, // 65
    {{ 0, 1, 0}, {{16, 2,6},{16, 2,7},{20, 2,7},{20, 2,6}} }, // 66
    {{ 0, 1, 0}, {{16, 2,1},{16, 2,3},{23, 2,3},{23, 2,1}} }, // 67
    {{-1, 0, 0}, {{23, 2,1},{23, 2,3},{23, 3,3},{23, 3,1}} }, // 68
    {{ 0, 1, 0}, {{23, 3,1},{23, 3,3},{24, 3,3},{24, 3,1}} }, // 69
    {{ 0, 1, 0}, {{22, 3,3},{22, 3,6},{24, 3,6},{24, 3,3}} }, // 70
    {{ 0, 1, 0}, {{20, 3,3},{20, 3,7},{22, 3,7},{22, 3,3}} }, // 71
    {{-1, 0, 0}, {{20, 2,3},{20, 2,7},{20, 3,7},{20, 3,3}} }, // 72
    {{-1, 0, 0}, {{20, 4,7},{20, 4,8},{20, 5,8},{20, 5,7}} }, // 73
    {{ 0, 1, 0}, {{20, 5,7},{20, 5,8},{22, 5,8},{22, 5,7}} }, // 74
    {{-1, 0, 0}, {{22, 5,7},{22, 5,8},{22, 7,8},{22, 7,7}} }, // 75
    {{ 0, 1, 0}, {{22, 7,7},{22, 7,8},{24, 7,8},{24, 7,7}} }, // 76
    {{-1, 0, 0}, {{24, 7,7},{24, 7,8},{24, 9,8},{24, 9,7}} }, // 77
    {{-1, 0, 0}, {{24, 3,1},{24, 3,6},{24, 9,6},{24, 9,1}} }, // 78
    {{-1, 0, 0}, {{24, 6,6},{24, 6,7},{24, 9,7},{24, 9,6}} }, // 79
    {{ 0, 1, 0}, {{22, 6,6},{22, 6,7},{24, 6,7},{24, 6,6}} }, // 80
    {{-1, 0, 0}, {{22, 3,6},{22, 3,7},{22, 6,7},{22, 6,6}} }, // 81
    {{ 0, 0,-1}, {{24, 6,6},{24, 3,6},{22, 3,6},{22, 6,6}} }, // 82
    {{ 0, 0,-1}, {{24, 9,8},{24, 7,8},{22, 7,8},{22, 9,8}} }, // 83
    {{ 0, 0,-1}, {{22, 9,8},{22, 5,8},{20, 5,8},{20, 9,8}} }, // 84
    {{ 0, 0,-1}, {{20, 9,8},{20, 4,8},{17, 4,8},{17, 9,8}} }, // 85
    {{ 0, 0,-1}, {{20, 4,7},{20, 2,7},{17, 2,7},{17, 4,7}} }, // 86
    {{ 0, 0,-1}, {{22, 5,7},{22, 3,7},{20, 3,7},{20, 5,7}} }, // 87
    {{ 0, 0,-1}, {{24, 7,7},{24, 6,7},{22, 6,7},{22, 7,7}} }, // 88
    {{ 0, 0,-1}, {{23, 3,3},{23, 2,3},{20, 2,3},{20, 3,3}} }, // 89
    {{ 0, 0, 1}, {{23, 2,1},{23, 3,1},{15, 3,1},{15, 2,1}} }, // 90
    {{ 0, 0, 1}, {{16, 4,1},{16, 6,1},{11, 6,1},{11, 4,1}} }, // 91
    {{ 1, 0, 0}, {{11, 6,1},{11, 6,5},{11, 4,5},{11, 4,1}} }, // 92
    {{ 1, 0, 0}, {{11, 6,5},{11, 6,6},{11, 2,6},{11, 2,5}} }, // 93
    {{-1, 0, 0}, {{11, 6,6},{11, 6,8},{11, 7,8},{11, 7,6}} }, // 94
    {{-1, 0, 0}, {{11, 8,6},{11, 8,8},{11,18,8},{11,18,6}} }, // 95
    {{ 0,-1, 0}, {{11,18,6},{11,18,8},{ 7,18,8},{ 7,18,6}} }, // 96
    {{ 0, 0,-1}, {{11,18,8},{11, 6,8},{ 7, 6,8},{ 7,18,8}} }, // 97
    {{ 0, 0, 1}, {{14, 2,1},{14, 4,1},{ 7, 4,1},{ 7, 2,1}} }, // 98
    {{ 0, 0, 1}, {{ 7, 2,2},{ 7, 4,2},{ 6, 4,2},{ 6, 2,2}} }, // 99
    {{-1, 0, 0}, {{ 6, 2,1},{ 6, 2,2},{ 6, 4,2},{ 6, 4,1}} }, // 100
    {{ 0, 1, 0}, {{ 6, 2,2},{ 6, 2,8},{ 8, 2,8},{ 8, 2,2}} }, // 101
    {{ 0, 1, 0}, {{11, 7,6},{11, 7,8},{13, 7,8},{13, 7,6}} }, // 102
    {{ 0, 1, 0}, {{17, 4,7},{17, 4,8},{20, 4,8},{20, 4,7}} }, // 103
    {{ 1, 0, 0}, {{ 7,18,6},{ 7,18,8},{ 7, 6,8},{ 7, 6,6}} }, // 104
};
static struct lighttype { float pos[3], dif[3]; } lights[] =
{
//    { {  1.5 , 5.5, 7.5 },  {.1,.1,.2 } }, // in begin ceiling
//    { { 22.5 , 5.7, 7.5 },  { .1, .3, .7 } },
  { { 17.3 , 5.7, 7.5 },   { 1, .2, .2 } },  // blue at the end
 //   { { 17.5 , 5.7, 7.5 },  {.06, .3, .7 } },
 //   { { 17.5 , 5.7, 4.5 },  {.06, .3, .7 } },
 //   { { 15.5 , 4.5, 7.5 },  {.02,.04, .1 } }, // middle (reddish)
//    { { 22.0 , 1.2, 1.5 },  {.03*.10, .13*.10, .2*.10 } },
    { { 15.2 , 2.2, 1.5 },  {.1, 0.6, 1 } },
    { { 9.5  , 17.5,   7 },  { 200,200,200 } }, // huge white in the ceiling tunnel
    { { 9.5  , 3.9, 1.1 },  { .2, .4, 2 } }  // blue in tunnel
};

struct HitRec
{
    unsigned wallno;
    float    distance, hit[3], alpha, beta;
    HitRec() : wallno(~0u), distance(0), hit() { }
    bool set() const { return wallno != ~0u; }
};
HitRec IntersectRay(const double org[3], const double dir[3])
{
    // Ray-quadrilateral intersection test by Ares Lagae and Philip Dutré.
    // Because our walls are simple rectangles, the calculation of
    // alpha11 and beta11 (both of which are always 1.0) was skipped,
    // and many extra branches were eliminated, simpifying the code.
    // Our walls are also perpendicular to main axis, but it is a coincidence.
    HitRec result;

    for(unsigned wallno=0; wallno < sizeof(map)/sizeof(*map); ++wallno)
    {
        const maptype& q = map[wallno];

        const double eps = 1e-9;
        #define halftest(A,B,C) \
            double e##A##B[3] = vsub(q.p[B], q.p[A]); \
            double e##A##C[3] = vsub(q.p[C], q.p[A]); \
            double vp[3]   = vcross(dir, e##A##C); \
            double det     = vdot(e##A##B, vp); \
            if(std::abs(det) < eps) continue; \
            double vt[3]   = vsub(org, q.p[A]); \
            double inv_det = 1. / det; \
            double alpha   = vdot(vt, vp) * inv_det; \
            if(alpha < 0.f) continue; \
            double vq[3]   = vcross(vt, e##A##B); \
            double beta    = vdot(dir, vq) * inv_det; \
            if(beta < 0.f) continue

        halftest(0,1,3);

        if( (alpha + beta) > 1.) { halftest(2,3,1); }

        // Compute the ray parameter of the intersection point, and
        // reject the ray if it does not hit q. Choose nearest hit.
        double t = vdot(e03, vq) * inv_det;
        if(t >= 0. && (!result.set() || t < result.distance))
        {
            result.wallno   = wallno;
            result.distance = t;
            result.alpha    = alpha;
            result.beta     = beta;
        }
    }
    if(result.set())
        for(unsigned a=0; a<3; ++a)
            result.hit[a] = org[a] + dir[a] * result.distance;
    return result;
}

HitRec IntersectSphereSweep(const double org[3], const double dir[3])
{
    // Ray-quadrilateral intersection test by Ares Lagae and Philip Dutré.
    // Because our walls are simple rectangles, the calculation of
    // alpha11 and beta11 (both of which are always 1.0) was skipped,
    // and many extra branches were eliminated, simpifying the code.
    // Our walls are also perpendicular to main axis, but it is a coincidence.
    HitRec result;

    for(unsigned wallno=0; wallno < sizeof(map)/sizeof(*map); ++wallno)
    {
        const maptype& q = map[wallno];

        const double eps = 1e-9;
        #define halftest(A,B,C) \
            double e##A##B[3] = vsub(q.p[B], q.p[A]); \
            double e##A##C[3] = vsub(q.p[C], q.p[A]); \
            double vp[3]   = vcross(dir, e##A##C); \
            double det     = vdot(e##A##B, vp); \
            if(std::abs(det) < eps) continue; \
            double vt[3]   = vsub(org, q.p[A]); \
            double inv_det = 1. / det; \
            double alpha   = vdot(vt, vp) * inv_det; \
            if(alpha < 0.f) continue; \
            double vq[3]   = vcross(vt, e##A##B); \
            double beta    = vdot(dir, vq) * inv_det; \
            if(beta < 0.f) continue

        halftest(0,1,3);

        if( (alpha + beta) > 1.) { halftest(2,3,1); }

        // Compute the ray parameter of the intersection point, and
        // reject the ray if it does not hit q. Choose nearest hit.
        double t = vdot(e03, vq) * inv_det;
        if(t >= 0. && (!result.set() || t < result.distance))
        {
            result.wallno   = wallno;
            result.distance = t;
            result.alpha    = alpha;
            result.beta     = beta;
        }
    }
    if(result.set())
        for(unsigned a=0; a<3; ++a)
            result.hit[a] = org[a] + dir[a] * result.distance;
    return result;
}

#if PRERENDER > 0

#define vadd(t,b)   for(unsigned a=0; a<3; ++a) t[a] += b[a]
#define vmul(t,b)   for(unsigned a=0; a<3; ++a) t[a] *= (b)
#define vdiv(t,b)   for(unsigned a=0; a<3; ++a) t[a] /= (b)
static void Bilinear(
    const float* sourcemap, unsigned sw, unsigned sh,
    float color[3], double sx, double sy)
{
    unsigned intsx = (unsigned)sx; sx -= intsx;
    unsigned intsy = (unsigned)sy; sy -= intsy;
    unsigned nextsx = std::min(sw-1u, intsx+1u);
    unsigned nextsy = std::min(sh-1u, intsy+1u);
    const float* rabove = &sourcemap[ 3 * intsy  * sw ];
    const float* rbelow = &sourcemap[ 3 * nextsy * sw ];
    for(unsigned c=0; c<3; ++c)
    {
        float above = rabove[c+3*intsx] + (rabove[c+3*nextsx]-rabove[c+3*intsx])*sx;
        float below = rbelow[c+3*intsx] + (rbelow[c+3*nextsx]-rbelow[c+3*intsx])*sx;
        color[c] += (above + (below-above)*sy);
    }
}
static void MakeLightMaps()
{
    const unsigned narealightcomponents = 64;
    const double area_light_radius       = 0.2;
    const unsigned nrandomvectors        = 8192;
    const unsigned maxrounds             = 100;
    const double fade_distance_diffuse   = 2.0;
    const double fade_distance_radiosity = 2.0;
    const double radiomul = 1.0;

    const unsigned nwalls = sizeof(map) / sizeof(*map);
    struct wallmeta
    {
        unsigned lmW, lmH;
        unsigned rmW, rmH;
        // The total rendered outcome:
        std::vector<float> total_lmap;
        // The outcome of last rendering round, not including preceding rounds:
        std::vector<float> last_lmap;
        // Outcome of current round
        std::vector<float> new_lmap;

        void ReadLightmapAt(double rx,double ry, float color[3]) const
        {
            Bilinear(&last_lmap[0], lmW,lmH, color, rx*lmW, ry*lmH);
        }
    } walls[nwalls];

    for(unsigned wallno=0; wallno<nwalls; ++wallno)
    {
        const maptype& m = map[wallno];
        float v30[3] = vsub(m.p[3],m.p[0]), v10[3] = vsub(m.p[1],m.p[0]);
        walls[wallno].lmW = 32u * vlen(v30);
        walls[wallno].lmH = 32u * vlen(v10);
        walls[wallno].rmW = 16u * vlen(v30);
        walls[wallno].rmH = 16u * vlen(v10);
        walls[wallno].total_lmap.resize( walls[wallno].lmW * walls[wallno].lmH * 3);
        walls[wallno].last_lmap.resize( walls[wallno].lmW * walls[wallno].lmH * 3);
    }
    for(unsigned round=1; round<=maxrounds; ++round)
    {
        fprintf(stderr, "Lighting calculation, round %u...\n", round);

        struct { double s, c, d; } rvec[nrandomvectors];
        for(unsigned n=0; n<nrandomvectors; ++n)
        {
            double angle = std::rand()/1e6;
            double angle2 = std::pow((std::rand()%10000)/10000.0, 1.3);
            rvec[n].s = std::sin(angle);
            rvec[n].c = std::cos(angle);
            rvec[n].d = std::tan(angle2*M_PI/2); // 0..+inf
        }

        for(unsigned wallno=0; wallno<nwalls; ++wallno)
        {
            const maptype& m = map[wallno];
            wallmeta& meta = walls[wallno];
            float v30[3] = vsub(m.p[3],m.p[0]);
            float v10[3] = vsub(m.p[1],m.p[0]);
            std::vector<float>& lmap = meta.new_lmap;
            // Set up a workbench
            lmap.clear();
            lmap.resize( meta.lmW * meta.lmH * 3 );

            if(round == 1)
            {
                // A lightsource is represented by a spherical cloud
                // of smaller lightsources around the actual lightsource.
                // This achieves smooth edges for the shadows.
                double tvec[narealightcomponents][3];
                for(unsigned qa=0; qa<narealightcomponents; ++qa)
                {
                    double* offset = tvec[qa];
                    for(unsigned v=0; v<3; ++v) offset[v] = drand(100.0);
                    double factor = area_light_radius / vlen(offset);
                    vmul(offset, factor);
                }

                // Round 1: Check lightsources
                for(unsigned y=0; y<meta.lmH; ++y)
                {
                    fprintf(stderr, "- Wall %u/%u %u/%u diffuse light...\r", wallno, nwalls, y,meta.lmH);
                    #pragma omp parallel for
                    for(unsigned x=0; x<meta.lmW; ++x)
                    {
                        float* color = &lmap[ (x + y*meta.lmW)*3 ];
                        GLfloat p[3] = { m.p[0][0],  m.p[0][1],  m.p[0][2] };
                        vaddm(p, v30, .01 + .98 * x/(double)(meta.lmW-1));
                        vaddm(p, v10, .01 + .98 * y/(double)(meta.lmH-1));
                        // For each lightsource, check if there is an obstacle
                        // in between this vertex and the lightsource. Calculate
                        // the ambient light levels from the fact.
                        // This simulates diffuse light.
                        double begin[3];
                        for(unsigned a=0; a<3; ++a) begin[a] = p[a] + m.normal[a]*1e-2;
                        for(unsigned qq=0; qq<sizeof(lights)/sizeof(*lights); ++qq)
                        {
                            const lighttype& l = lights[qq];
                            for(unsigned qa=0; qa<narealightcomponents; ++qa)
                            {
                                double towards[3] = vsub(l.pos, p);
                                vadd(towards, tvec[qa]);
                                double len = vlen(towards); vdiv(towards, len);
                                double cosine = vdot( m.normal, towards );
                                const double fade_distance = fade_distance_diffuse;
                                double power = cosine / (1. + std::pow(len/fade_distance, 2.0));
                                power /= (double) narealightcomponents;
                                if(power > 1e-7)
                                {
                                    HitRec tmp = IntersectRay(begin, towards);
                                    if(!tmp.set() || tmp.distance > len)
                                        vaddm(color, l.dif, power);
                }   }   }   }   }
            }
            else
            {
                // Round 2+: Radiosity
                float v10abs[3] = { v10[0],v10[1],v10[2] }, v10len = vlen(v10abs);
                float v30abs[3] = { v30[0],v30[1],v30[2] }, v30len = vlen(v30abs);
                vdiv(v10abs, v10len);
                vdiv(v30abs, v30len);

                // Apply the set of random vectors to this surface.
                // This produces a set of vectors all pointing away
                // from the wall to random directions.
                double tvec[nrandomvectors][4];
                for(unsigned qq=0; qq<nrandomvectors; ++qq)
                {
                    double* towards = tvec[qq];
                    for(unsigned v=0; v<3; ++v)
                        towards[v] =
                            m.normal[v] * rvec[qq].d
                          + v10abs[v] * rvec[qq].s
                          + v30abs[v] * rvec[qq].c;
                    double len = vlen(towards); vdiv(towards, len);
                    towards[3] = vdot(m.normal, towards); // precalculate the cosine
                }

                // Calculate radiosity in decreased resolution
                std::vector<float> rmap(meta.rmW * meta.rmH * 3);

                for(unsigned y=0; y<meta.rmH; ++y)
                {
                    fprintf(stderr, "- Wall %u/%u %u/%u radiosity...\r", wallno, nwalls, y,meta.rmH);
                    #pragma omp parallel for
                    for(unsigned x=0; x<meta.rmW; ++x)
                    {
                        float* color = &rmap[ (x + y*meta.rmW)*3 ];
                        GLfloat p[3] = { m.p[0][0],  m.p[0][1],  m.p[0][2] };
                        vaddm(p, v30, .01 + .98 * x/(double)(meta.rmW-1));
                        vaddm(p, v10, .01 + .98 * y/(double)(meta.rmH-1));
                        // Shoot rays to each random direction and see what it hits.
                        // Take the last round's light value from that location.
                        double begin[3];
                        for(unsigned a=0; a<3; ++a) begin[a] = p[a] + m.normal[a]*1e-2;
                        double basepower = radiomul / nrandomvectors;
                        for(unsigned qq=0; qq<nrandomvectors; ++qq)
                        {
                            const double* towards = tvec[qq];
                            double cosine = towards[3];

                            HitRec hit = IntersectRay(begin, towards);
                            if(hit.set())
                            {
                                double len = hit.distance;
                                float tmp[3] = {0,0,0};
                                walls[hit.wallno].ReadLightmapAt(hit.beta, hit.alpha, tmp);
                                const double fade_distance = fade_distance_radiosity;
                                double value = basepower/(1. + std::pow(len/fade_distance,2.0));
                                double power = value * cosine;
                                vaddm(color, tmp, power);
                }   }   }   }
                for(unsigned y=0; y<meta.lmH; ++y)
                {
                    fprintf(stderr, "- Wall %u/%u %u/%u scaling...\r", wallno, nwalls, y,meta.lmH);
                    #pragma omp parallel for
                    for(unsigned x=0; x<meta.lmW; ++x)
                    {
                        float* color = &lmap[ (x + y*meta.lmW)*3 ];
                        // Add the previously calculated radiosity. Because the
                        // radiosity bitmap is smaller, bilinearly interpolate it.
                        Bilinear(&rmap[0], meta.rmW, meta.rmH,
                            color,
                            x * (meta.rmW-1.0) / (meta.lmW-1),
                            y * (meta.rmH-1.0) / (meta.lmH-1) );
                }   }
            }
            fprintf(stderr, "\n");
        }

        fprintf(stderr, "Accumulating and normalizing...\n");
        float min = 1e9, max = 1e-9, avg=0, sum = 0;
        size_t avgcount=0;
        for(unsigned wallno=0; wallno<nwalls; ++wallno)
        {
            std::vector<float>& map = walls[wallno].new_lmap;
            std::vector<float>& map2 = walls[wallno].total_lmap;
            for(size_t a=0; a<map.size(); ++a)
            {
                map2[a] += map[a];
                if(map2[a] < min) min = map2[a];
                if(map2[a] > max) max = map2[a];
                avg += map2[a];
                sum += map[a];
            }
            avgcount += map.size();
            walls[wallno].last_lmap.swap( walls[wallno].new_lmap );
        }
        avg /= avgcount;
        fprintf(stderr, "Saving new lightmaps... Sum of change this round %.4f, min=%.4f,avg=%.4f,max=%.4f\n",
            sum,min,avg,max);
        // Install the newly calculated lightmaps
        for(unsigned wallno=0; wallno<nwalls; ++wallno)
        {
            wallmeta& meta = walls[wallno];
            std::vector<float>& map = meta.total_lmap;

            std::vector<float> out(map);
            for(size_t a=0; a<out.size(); ++a)
                out[a] =
                    (((out[a] - min) / (std::min(1.5f,avg) - min))
                     * ( 3.4 * 254.0/255) + 1/255.0);

            char Buf[64]; std::sprintf(Buf, "lmap%u.raw", wallno);
            FILE* fp = std::fopen(Buf, "wb");
             std::fwrite(&out[0], out.size(), sizeof(float), fp);
            std::fclose(fp);
        }
        if(sum < 1.0) break;
    }
}
#endif

#if PRERENDER < 2
static void InstallTexture(const void* data,int w,int h, int txno,
                           int type1, int type2, int filter, int wrap)
{
    glBindTexture(GL_TEXTURE_2D, txno);
    glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, wrap);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, wrap);
    #if LIGHTWEIGHT >= 1
    //filter = GL_NEAREST_MIPMAP_NEAREST;
    filter = GL_NEAREST;
    #endif
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, filter);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, filter);
    glTexEnvi (GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
    if(filter != GL_LINEAR && filter != GL_NEAREST)
        gluBuild2DMipmaps(GL_TEXTURE_2D, type1, w,h, type1, type2, data);
    else
        glTexImage2D(GL_TEXTURE_2D, 0, type1, w,h, 0, type1, type2, data);
}

// This function converts the level map into OpenGL quad primitives. This code
// is not particularly optimized (in particular, everything is always rendered).
static void ExtractLevelMap()
{
    // Walls are all created using this one texture.
    glActiveTextureARB(GL_TEXTURE0_ARB); glEnable(GL_TEXTURE_2D);
    glBindTexture(GL_TEXTURE_2D, 1u); glColor3f(1,1,1);
  #if LIGHTWEIGHT < 2
    glActiveTextureARB(GL_TEXTURE1_ARB); glEnable(GL_TEXTURE_2D);
    glShadeModel(GL_SMOOTH);
  #endif
    for(unsigned wallno=0; wallno < sizeof(map)/sizeof(*map); ++wallno)
    {
        const maptype& m = map[wallno];
        float v10[3] = vsub(m.p[1], m.p[0]);
        float v30[3] = vsub(m.p[3], m.p[0]);
        int width  = vlen(v30); // Number of times the texture
        int height = vlen(v10); // is repeated across the surface.

    #if LIGHTWEIGHT < 2
      #if USE_FLOAT
        if(!TexturesInstalled)
        {
      #endif
        // Create a diffuse + radiosity lightmap.
        unsigned lmW = width * 32, lmH = height * 32;
        std::vector<float> map(  lmW*lmH*3 );
        // Load it from a file I have created beforehand. The code to generate
        // the map adds about 150 lines to this program's size, and there's not
        // enough time in Youtube's limit for me to include that part here.
        // Simply put, it calls IntersectRay many times to see which lightsources
        // are visible from each spot, and does the same recursively to random
        // directions to catch indirect light ("global illumination").
        char Buf[64]; std::sprintf(Buf, "lmap%u.raw", wallno);
        FILE* fp = std::fopen(Buf, "rb");
        std::fread(&map[0], map.size(), sizeof(float), fp);
        for(size_t a=0; a<map.size(); ++a) map[a] *= 4.0;
        std::fclose(fp);
      #if USE_FLOAT
        for(unsigned step = 0; step < NumBrightnessSteps; ++step)
        {
            std::vector<float> tmp( map );
            for(size_t a=0; a<map.size(); ++a) tmp[a] *= BrightnessSteps[step];
            InstallTexture(&tmp[0],lmW,lmH,
                wallno*NumBrightnessSteps + 2+step,
                GL_RGB, GL_FLOAT, GL_LINEAR, GL_CLAMP_TO_EDGE);
            }
      #else
        #if LIGHTWEIGHT == 0
        InstallTexture(&map[0],lmW,lmH, 2+wallno, GL_RGB, GL_FLOAT, GL_LINEAR, GL_CLAMP_TO_EDGE);
        #else
        unsigned lW = lmW / 8, lH = lmH / 8;
        std::vector<float> smaller_map(lW * lH * 3);
        for(unsigned y=0; y<lH; ++y)
          for(unsigned x=0; x<lW; ++x)
            for(unsigned c=0; c<3; ++c)
              smaller_map[c + 3*(x+ y*lW)] = map[c + 3*((x*lmW/lW) + (y*lmH/lH)*lmW)];
        InstallTexture(&smaller_map[0],lW,lH, 2+wallno, GL_RGB, GL_FLOAT, GL_LINEAR, GL_CLAMP_TO_EDGE);
        #endif
      #endif
      #if USE_FLOAT
        } // TexturesInstalled
        glBindTexture(GL_TEXTURE_2D, wallno*NumBrightnessSteps + 2+CurrentBrightnessStep);
      #endif
    #endif
        glNormal3fv( m.normal );
        glBegin(GL_QUADS);
        for(unsigned e=0; e<4; ++e)
        {
            glMultiTexCoord2fARB(GL_TEXTURE0_ARB, width * !((e+2)&2), height * !((e+3)&2) );
            glMultiTexCoord2fARB(GL_TEXTURE1_ARB,     1 * !((e+2)&2),      1 * !((e+3)&2) );
            glVertex3fv( m.p[e] );
        }
        glEnd();
    }
  #if USE_FLOAT
    TexturesInstalled = true;
  #endif
}
#endif

int main()
{
#if PRERENDER > 0
    MakeLightMaps();
#endif
#if PRERENDER < 2
  #if USE_FLOAT
    OSMesaContext om = OSMesaCreateContext(OSMESA_RGB, NULL);
    OSMesaMakeCurrent(om, PC::ImageBufferFloat, GL_FLOAT, PC::W, PC::H);
  #else
    OSMesaContext om = OSMesaCreateContext(OSMESA_RGBA, NULL);
    OSMesaMakeCurrent(om, PC::ImageBuffer, GL_UNSIGNED_BYTE, PC::W, PC::H);
  #endif

    ////////
    glEnable(GL_DEPTH_TEST);
    glEnable(GL_CULL_FACE); glCullFace(GL_BACK);

    // Generate a very simple rectangle of a texture.
  #if LIGHTWEIGHT == 0
    const unsigned txW=256, txH=256;
    GLfloat texture[txH*txW];
    for(unsigned y=0; y<txH; ++y)
        for(unsigned x=0; x<txW; ++x)
            texture[y*txW+x] =
             1.0 - std::max(0.40*std::pow((std::rand()%100)/100.0, 2.0),
                            0.20*!(x<8 || y<8 || (x+8)>=txW || (y+8)>=txH))
             ;
    InstallTexture(texture, txW,txH, 1u, GL_LUMINANCE,GL_FLOAT, GL_LINEAR_MIPMAP_LINEAR, GL_REPEAT);
  #else
    const unsigned txW=256, txH=256;
    GLfloat texture[txH*txW];
    for(unsigned y=0; y<txH; ++y)
        for(unsigned x=0; x<txW; ++x)
            texture[y*txW+x] =
             1.0 - 0.70*!(x<8 || y<8 || (x+8)>=txW || (y+8)>=txH);
    InstallTexture(texture, txW,txH, 1u, GL_LUMINANCE,GL_FLOAT, GL_NEAREST, GL_REPEAT);
  #endif

    PC::Init();

    bool need_redraw = true;
    const double below = 1.0, above = 0.3, around = 0.2;

    double look_angle = 170, yaw = 10, move_angle = 0;
    double camera[3] = {  4, 3, 7.25 };
    //double camera[3] = { 21, 3, 2.25 }; yaw = -50;
    double dir[3]    = {  1, 0,    0 };
    double hvel[3]   = {0,0,0}, vvel = 0;
    bool falling = true;
    int moving = 0;
    // Main loop
    for(;;)
    {
        if(need_redraw)
        {
            dir[0] = sin(look_angle * M_PI / 180.0);
            dir[2] = cos(look_angle * M_PI / 180.0);

            glMatrixMode(GL_PROJECTION);
            glLoadIdentity();
            glFrustum(-.1333e-2, .1333e-2, -.1e-2, .1e-2, .001, 30.0);
            glMatrixMode(GL_MODELVIEW);
            glLoadIdentity();
            glRotatef(yaw+vvel*35.f, 1,0,0); // autotilt camera when jumping
            if(moving) // add some camera shake
                glRotatef( drand(.1) * (vlen(hvel) + std::abs(vvel)), 0,0,1);

            gluLookAt( camera[0], camera[1], camera[2],
                       camera[0] + dir[0] + drand(6e-5),
                       camera[1] + dir[1] + drand(6e-5),
                       camera[2] + dir[2] + drand(6e-5),
                       0,-1,0);

//static int flip=0;
//if(++flip % 20 == 0)
//{
        
        #if USE_FLOAT
            for(CurrentBrightnessStep=0; CurrentBrightnessStep<NumBrightnessSteps; ++CurrentBrightnessStep)
            {
        #endif
        
            // Enable depth calculations to work on the new frame.
            glClear(GL_DEPTH_BUFFER_BIT);
        #if LIGHTWEIGHT < 1
            // Create white spheres representing all lightsources.
            glActiveTextureARB(GL_TEXTURE1_ARB); glDisable(GL_TEXTURE_2D); glColor3f(1, 1, 1);
            glActiveTextureARB(GL_TEXTURE0_ARB); glDisable(GL_TEXTURE_2D); glColor3f(1, 1, 1);
            for(const auto& l : lights)
            {
                glTranslatef( l.pos[0], l.pos[1], l.pos[2] );
                GLUquadric* qu = gluNewQuadric();
                 gluSphere(qu, 0.2f/*std::min(1.2f,vlen(l.dif))*2.f*/, 16, 16);
                gluDeleteQuadric(qu);
                glTranslatef( -l.pos[0], -l.pos[1], -l.pos[2] );
            }
        #endif
            ExtractLevelMap();

            IntersectRay( camera, dir );

            glFlush();
        #if USE_FLOAT
            PC::UpdateFloatPixmap();
            }
            PC::ConvertFloatPixmap();
        #endif
//}
            need_redraw = false;
        }

        PC::Render();

        // We do collision tests in two phases. For the purposes of
        // of collision testing, the player is shaped like a letter "E".
        // We first check vertically, whether the player's top or bottom
        // (depending on vertical direction) is colliding with a wall.
        // Then we check horizontally, whether the nearest wall is
        // within the player's cylinder radius.

        /********
        TODO: Use algorithm from the paper
              "Improved Collision detection and Response"
              by Kasper Fauerby.
        Described at:
        http://www.peroxide.dk/papers/collision/collision.pdf

        In principle: Player is represented by an ellipsoid
        with a dimensions in each space. Translate the player
        and the space so that shape of the player an unit sphere.
        I.e. multiply all vectors by
          {1/around, 1/(below+above), 1/around}
       (assuming the ellipse is centered at the player's eye level,
        centering may require more transformations).
        ********/

        // Vertical collision detection
        bool ground = !falling;
        if(falling)
        {
            vvel -= 0.011; // Add gravity
            double down[3] = {0, vvel < 0 ? -1. : 1., 0};
            HitRec vcollision = IntersectRay(camera, down);
            if(!vcollision.set())
            {
                // Absolutely nothing is ahead in our current trajectory.
                // Stop falling, to prevent descending to infinity.
                falling = vvel > 0.f;
                vvel    = 0;
            }
            else if(vvel < 0 && vcollision.distance <= below)
            {
                // The nearest wall is within our feet's distance.
                // Fix to ground.
                // TODO: land on feet sfx if vvel < -0.03
                double fix = below - vcollision.distance;
                camera[1] += fix;
                vvel = 0;
                falling     = false;
                need_redraw = true;
                ground      = true;
            }
            else if(vvel > 0 && vcollision.distance <= above)
            {
                // The nearest wall is within our head's distance.
                // I.e. we hit a ceiling.
                vvel    = 0;
                falling = true;
            }
            if(falling)
            {
                camera[1] += vvel;
                need_redraw = true;
            }
        }
        // Horizontal collision detection
        if(moving)
        {
            // Try to push into the desired direction
            const double maxvel = 0.1*(moving>0), acceleration = (moving>0 ? 0.2 : 0.1);

            double move_vec[3] = {
                dir[0] * std::cos(move_angle) - dir[2]*std::sin(move_angle),
                dir[1],
                dir[0] * std::sin(move_angle) + dir[2]*std::cos(move_angle)
            };
            for(unsigned a=0; a<3; ++a) hvel[a] *= (1-acceleration);
            for(unsigned a=0; a<3; ++a) hvel[a] += move_vec[a]*acceleration*maxvel;

            double vel = vlen(hvel);
            if(std::abs(vel) < 1e-9)
            {
                // Stop moving if our velocity is low enough.
                moving = 0;
            }
            else
            {
                // TODO: walk sfx
                // Do two collision checks.
                //  One for our feet level and one for our head level.
                double heights[2] = { -below+1e-4, above-1e-4 };
                double movedir[3] = { hvel[0]/vel, hvel[1]/vel, hvel[2]/vel };
                double wall_distance = 1e5f, oldvel = vel;
                unsigned nearest_wall = ~0u;
                for(unsigned h=0; h<2; ++h)
                {
                    // Vantage point = camera + (feet or head level)
                    double vantage[3] = { camera[0], camera[1] + heights[h], camera[2] };
                    HitRec hcollision = IntersectRay(vantage, movedir);
                    if(hcollision.set() && hcollision.distance < wall_distance)
                    {
                        // There's a wall ahead us. Keep the nearest wall.
                        wall_distance = hcollision.distance;
                        nearest_wall = hcollision.wallno;
                    }
                }
                // Figure out where the wall is in relation to our body part
                if(wall_distance < around) // already inside the wall, push out
                {
                    const auto& normal = map[nearest_wall].normal;
                    vaddm(camera, normal, around-wall_distance);
                    vel    = 0;
                    moving = 0;
                }
                else if(wall_distance < vel+around) // would go inside wall if we moved now
                {
                    // TODO: Move tangentially to the wall
                    vel = wall_distance+1e-3f;
                    moving = 0;
                }
                // Update the camera position according to our velocity vector.
                for(unsigned a=0; a<3; ++a) camera[a] += (hvel[a] *= vel / oldvel);
                need_redraw = true;
                falling = true;
            }
        }

        moving = -1;
        while(kbhit())
            switch(getch())
            {
                case 'q': case 27: case 'Q': goto done;
                case 'w': moving = 1; move_angle = 0;        break; // forward
                case 's': moving = 1; move_angle = M_PI;     break; // backward
                case 'a': moving = 1; move_angle = M_PI* .5; break; // strafe left
                case 'd': moving = 1; move_angle = M_PI*-.5; break; // strafe right
                case ' ': // jump
                    if(ground)
                    {
                        vvel += 0.18;
                        falling = true;
                        // TODO: jump sfx
                    }
                    break;

            }
      #ifdef __DJGPP__
        // Get mouse relative position (since last poll) and update view
        __dpmi_regs regs = { };
        regs.x.ax = 0xB;
        __dpmi_int(0x33, &regs);
        if(regs.x.cx || regs.x.dx) need_redraw = true;
        look_angle += (short) regs.x.cx;
        yaw        += (short) regs.x.dx;
      #else
        int x,y;
        static int count=2;
        SDL_GetRelativeMouseState(&x,&y);
        if(x || y) need_redraw = true;
        if(count > 0) --count; else {
        look_angle += x;
        yaw    += y; }
      #endif
        if(yaw < -90) yaw = -90; // defer too aggressive tilts
        if(yaw >  90) yaw =  90;
/*
    #if !REPLAY
        static FILE* fp = fopen("movement.log", "wt");
        fprintf(fp, "%d %10.4f %d %10.4f %10.4f %10.4f %d\n",
            moving, move_angle,
            (int)falling, vvel,
            look_angle, yaw,
            (int) need_redraw);
        fflush(fp);
    #else
        static FILE* fp = fopen("movement.log", "rt");
        int f,n;
        fscanf(fp, "%d %lf %d %lf %lf %lf %d\n",
            &moving, &move_angle,
            &f, &vvel,
            &look_angle, &yaw,
            &n);
        falling = f;
        need_redraw = n;
    #endif
*/
        /*if(falling && std::abs(vvel) > 0.01)
            usleep(100000);*/

        static time_t lasttime=0;
        static unsigned nframes=0;
        time_t t = time(NULL);
        if(t != lasttime)
        {
            fprintf(stderr, "fps=%u\n", nframes);
            lasttime=t;
            /*if(nframes > 30)
            {
                double excess = (nframes-30) * (1e6/30);
                usleep(excess);
            }*/
            nframes=0;
        } else ++nframes;
    }
done:;
    PC::Close();

    OSMesaDestroyContext(om);
#endif
}