/* 
 * Copyright (C) 2006 Robert Geist
 * All rights reserved.
 */
#define GL_GLEXT_PROTOTYPES 1
#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glut.h>
#include <GL/glx.h>
#include <GL/glext.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <fcntl.h>
#include "xrgba.h"

#define G (0.6)
#define SIGMA_A (0.008)
#define SIGMA_S (0.792)
#define SIGMA_T (SIGMA_A + SIGMA_S)
#define SUN_AMBIENT 0.65
#define SUN_DIFFUSE 4.0

#define DIRECTIONS 19
#define FINAL_TIME (2*(EDGE))
#define EPSILON 0.00001

struct ivector {
	int x;	/* col */
	int y;	/* row */
	int z;  /* dep */
	}; 

struct ivector ci[DIRECTIONS] = {
	 /* index 0 is rest density == absorbed/transmitted */
	 0, 0, 0,   
	 /* axis directions are indices 1 - 6    */
	 1, 0, 0,   0, 1, 0,   0, 0, 1,
	-1, 0, 0,   0,-1, 0,   0, 0, -1,
	 /* non-corner cube lattice points for 7 - 18 */
	 1, 1, 0,   1,-1, 0,  -1, 1, 0,  -1,-1, 0,
	 1, 0, 1,   1, 0,-1,  -1, 0, 1,  -1, 0,-1,
	 0, 1, 1,   0, 1,-1,   0,-1, 1,   0,-1,-1
	};

struct rvector {
	float x;
	float y;
	float z;
	};

float dot(struct rvector *uptr, struct rvector *vptr)
{
return(uptr->x*vptr->x + uptr->y*vptr->y + uptr->z*vptr->z);
}

void normalize(struct rvector *pv)
{
float length;

length = sqrt(dot(pv,pv));
if(length<=0.0){
	fprintf(stderr,"eek\n");
	exit(1);
	}
pv->x /= length;
pv->y /= length;
pv->z /= length;
}

int vec_eq(struct rvector *uptr, struct rvector *vptr)
{
float sum;
sum = fabs(uptr->x - vptr->x)+fabs(uptr->y - vptr->y)+fabs(uptr->z - vptr->z);
if(sum<EPSILON) return(1);
return(0);
}

struct rvector vec_smult(struct rvector *uptr, float s)
{
struct rvector p;
p.x = uptr->x*s;
p.y = uptr->y*s;
p.z = uptr->z*s;
return(p);
}

struct rvector vec_sub(struct rvector *uptr, struct rvector *vptr)
{
struct rvector p;
p.x = uptr->x - vptr->x;
p.y = uptr->y - vptr->y;
p.z = uptr->z - vptr->z;
return(p);
}

/*
 * This is needed for both init_lattice and collision matrix setup.
 */
struct rvector univ[DIRECTIONS];

fix_univ()
{
/* 
 * compute unit-length direction vectors; these are
 * used for the collision matrix too, hence the global
 */ 
int i;
univ[0].x = univ[0].y = univ[0].z = 0.0;
for(i=1;i<DIRECTIONS;i++){
	univ[i].x = (float)(ci[i].x);
	univ[i].y = (float)(ci[i].y);
	univ[i].z = (float)(ci[i].z);
	normalize(&univ[i]);
	}
}

/* 
 * Photon density and cloud density textures are too big for stack passing.
 * We don't need them once we hand them off to the graphics card, but we need 
 * the space again when we're done to readback and re-assemble the lighting 
 * values.  So, we use global statics and just avoid the re-allocation time. 
 */
float photon_d[(EDGE/4)*EDGE][EDGE*DIRECTIONS][4];
float cloud_d[(EDGE/4)*EDGE][EDGE][4];
float omega[DIRECTIONS][DIRECTIONS];
float directions[3*DIRECTIONS];

void init_lattice(char *filename)
{
/* 
 * Load cloud densities from file and init photon densities on cube BORDER. 
 * These will remain fixed. For each face, we need to find dot products with 
 * the sun direction. For those that are positive, allocate source density 
 * accordingly. Ambient density is added to each BORDER node in each direction.
 * Init interior photon densities to zero. Overall brightness is controlled by 
 * sun_ambient and sun_diffuse parameters.
 *
 * The flattening map for photon density is:
 * f(i,j,k,m) -> pingpong[mod(k,EDGE/4)*EDGE+i][j*DIRECTIONS+m][floor(k,EDGE/4)]
 * so cube coordinates can be recovered from texture coordinates (s,t)
 * within the fragment program by 
 *
 * i = mod(t,EDGE)
 * j = floor(s/DIRECTIONS)
 * mod(k,(EDGE/4)) = floor(t/EDGE), 
 *	+0 (red) = 0*(EDGE/4) + mod(k,EDGE/4)
 *	+1 (green) = 1*(EDGE/4) + mod(k,EDGE/4)
 *	+2 (blue) = 2*(EDGE/4) + mod(k,EDGE/4)
 *	+3 (alpha) = 3*(EDGE/4) + mod(k,EDGE/4)
 * m = mod(s,DIRECTIONS)
 *
 * and that for cloud density is:
 * cd(i,j,k) -> cloud_d[mod(k,EDGE/4)*EDGE+i][j][floor(k,EDGE/4)]
 *
 * where dimensions are pingpong[(EDGE/4)*EDGE][EDGE*DIRECTIONS][4] and
 * cloud_d[(EDGE/4)*EDGE][EDGE][4].  If EDGE==128, EDGE*EDGE would exceed 
 * the 4096 limit on texture dimensions, so we fold the arrays into RGBA 
 * components to reduce the row dimension from EDGE*EDGE to (EDGE/4)*EDGE.
 *
 */
int i,j,k,m;
struct rvector sun_dir, dyn_sun_dir, this_much;
float sun_ambient, sun_diffuse;
FILE *fptr;
int dim[3], size;
float *densities;
float max_component;
float sun_component[DIRECTIONS];
int max_comp_index;

sun_ambient = SUN_AMBIENT;
sun_diffuse = SUN_DIFFUSE;
/* sun goes somewhere like (-1000,1000,-1000) and points back at us */
sun_dir.x = 1.0; 
sun_dir.y = -1.0; 
sun_dir.z = 1.0;

normalize(&sun_dir);

/* get cloud densities here ... assume Karl's format */
fptr=fopen(filename,"r");
fread(dim,sizeof(int),3,fptr);
if(dim[0]!=EDGE || dim[1]!=EDGE || dim[2]!=EDGE){
	fprintf(stderr,"can't handle dimensions\n");
	fclose(fptr);
	exit(1);
	}
size=dim[0]*dim[1]*dim[2];
densities = (float *)calloc(sizeof(float),size);
fread(densities,sizeof(float),size,fptr);
fclose(fptr);

/* load cloud densities and zip out photon densities */
for(i=0;i<EDGE;i++){
	for(j=0;j<EDGE;j++){
		for(k=0;k<EDGE;k++){
			cloud_d[(k%(EDGE/4))*EDGE+i][j][k/(EDGE/4)] = 
				densities[i*EDGE*EDGE+j*EDGE+k];
			for(m=0;m<DIRECTIONS;m++){
				photon_d[(k%(EDGE/4))*EDGE+i][j*DIRECTIONS+m]
					[k/(EDGE/4)] = 0.0;
				}
			}
		}
	}
free(densities);

/*
 * resolve the sun direction into lattice direction components for
 * each face; all BORDER nodes in that face are allocated source
 * photon density accordingly
 */
for(m=0;m<DIRECTIONS;m++) sun_component[m] = 0.0;
dyn_sun_dir = sun_dir;

while(!vec_eq(&dyn_sun_dir,&univ[0])){
	max_component = -1.0;
	for(m=1;m<DIRECTIONS;m++){
		if(dot(&univ[m],&dyn_sun_dir)>=(max_component)){
			max_component = dot(&univ[m],&dyn_sun_dir);
			max_comp_index = m;
			}
		}
	if(max_component<=0.0){
		fprintf(stderr,"eek\n");
		exit(1);
		}
	sun_component[max_comp_index] += max_component;
	this_much = vec_smult(&univ[max_comp_index],max_component);
	dyn_sun_dir = vec_sub(&dyn_sun_dir,&this_much);
	}


/* now load up the cube boundary values */
for(m=1;m<DIRECTIONS;m++){
	/* left and right */
	for(j=0;j<EDGE;j++){
		for(k=0;k<EDGE;k++){
			photon_d[(k%(EDGE/4))*EDGE+0][j*DIRECTIONS+m]
				[k/(EDGE/4)] = sun_ambient + 
				sun_diffuse*sun_component[m]; 
			photon_d[(k%(EDGE/4))*EDGE+(EDGE-1)][j*DIRECTIONS+m]
				[k/(EDGE/4)] = sun_ambient + 
				sun_diffuse*sun_component[m]; 
			}
		}
	/* top and bottom */
	for(i=0;i<EDGE;i++){
		for(k=0;k<EDGE;k++){
			photon_d[(k%(EDGE/4))*EDGE+i][0*DIRECTIONS+m]
				[k/(EDGE/4)] = sun_ambient + 
				sun_diffuse*sun_component[m]; 
			photon_d[(k%(EDGE/4))*EDGE+i][(EDGE-1)*DIRECTIONS+m]
				[k/(EDGE/4)] = sun_ambient +
				sun_diffuse*sun_component[m]; 
			}
		}
	/* front and back */
	for(i=0;i<EDGE;i++){
		for(j=0;j<EDGE;j++){
			photon_d[0*EDGE+i][j*DIRECTIONS+m][0] = sun_ambient + 
				sun_diffuse*sun_component[m]; 
			photon_d[(EDGE/4-1)*EDGE+i][j*DIRECTIONS+m][3] = 
				sun_ambient + sun_diffuse*sun_component[m]; 
			}
		}
	}

// have to pass directions to the shader too
for(j=0;j<3*DIRECTIONS;j += 3){
	directions[j+0] = (float)(ci[j/3].x);
	directions[j+1] = (float)(ci[j/3].y);
	directions[j+2] = (float)(ci[j/3].z);
	}

}

unsigned int t_unit[] = {GL_TEXTURE0_ARB, GL_TEXTURE1_ARB, GL_TEXTURE2_ARB, 
			GL_TEXTURE3_ARB, GL_TEXTURE4_ARB};
unsigned int t_id[] = {1, 2, 3, 4, 5};
unsigned int t_attach[] = {GL_COLOR_ATTACHMENT0_EXT, GL_COLOR_ATTACHMENT1_EXT};

void load_textures()
{
int tex;
int tcount;

/* set a current framebuffer id; we only need 1 */
glBindFramebufferEXT(GL_FRAMEBUFFER_EXT,1);

/* init both pingpong buffers */
for(tex=0;tex<2;tex++){
	glActiveTextureARB(t_unit[tex]);
	glBindTexture(GL_TEXTURE_RECTANGLE_NV,t_id[tex]);

	glTexImage2D(GL_PROXY_TEXTURE_RECTANGLE_NV,0,GL_FLOAT_RGBA32_NV,
		EDGE*DIRECTIONS,EDGE*EDGE/4,0,GL_RGBA,GL_FLOAT,NULL);

	glGetTexLevelParameteriv(GL_PROXY_TEXTURE_RECTANGLE_NV,0,
		GL_TEXTURE_WIDTH,&tcount);
	fprintf(stderr,"check width of allocation for texture: %d\n",tcount);
	glGetTexLevelParameteriv(GL_PROXY_TEXTURE_RECTANGLE_NV,0,
		GL_TEXTURE_HEIGHT,&tcount);
	fprintf(stderr,"check height of allocation for texture: %d\n",tcount);
	glGetTexLevelParameteriv(GL_PROXY_TEXTURE_RECTANGLE_NV,0,
		GL_TEXTURE_RED_SIZE,&tcount);
	fprintf(stderr,"check red allocation for texture: %d\n",tcount);
	glGetTexLevelParameteriv(GL_PROXY_TEXTURE_RECTANGLE_NV,0,
		GL_TEXTURE_GREEN_SIZE,&tcount);
	fprintf(stderr,"check green allocation for texture: %d\n",tcount);
	glGetTexLevelParameteriv(GL_PROXY_TEXTURE_RECTANGLE_NV,0,
		GL_TEXTURE_BLUE_SIZE,&tcount);
	fprintf(stderr,"check blue allocation for texture: %d\n",tcount);
	glGetTexLevelParameteriv(GL_PROXY_TEXTURE_RECTANGLE_NV,0,
		GL_TEXTURE_ALPHA_SIZE,&tcount);
	fprintf(stderr,"check alpha allocation for texture: %d\n",tcount);

	glTexImage2D(GL_TEXTURE_RECTANGLE_NV,0,GL_FLOAT_RGBA32_NV,
		EDGE*DIRECTIONS,(EDGE/4)*EDGE,0,GL_RGBA,GL_FLOAT,
		&photon_d[0][0][0]);
	glTexEnvf(GL_TEXTURE_ENV,GL_TEXTURE_ENV_MODE,GL_REPLACE);
	glTexParameterf(GL_TEXTURE_RECTANGLE_NV,GL_TEXTURE_MIN_FILTER,
		GL_NEAREST);
	glTexParameterf(GL_TEXTURE_RECTANGLE_NV,GL_TEXTURE_MAG_FILTER,
		GL_NEAREST);

	/* 
	 * Bind texture to current fbo (GL_FRAMEBUFFER_EXT)
	 */
	glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT,t_attach[tex],
		GL_TEXTURE_RECTANGLE_NV,t_id[tex],0);

	tcount = glCheckFramebufferStatusEXT(GL_FRAMEBUFFER_EXT);
	if(tcount==GL_FRAMEBUFFER_COMPLETE_EXT) fprintf(stderr,"fb ok\n");
	else fprintf(stderr,"fb oops %x\n",tcount);
	}

/* cloud density texture is static, not bound */
glActiveTextureARB(t_unit[2]);
glBindTexture(GL_TEXTURE_RECTANGLE_NV,t_id[2]);
glTexImage2D(GL_TEXTURE_RECTANGLE_NV,0,GL_FLOAT_RGBA32_NV,EDGE,(EDGE/4)*EDGE,0,
	GL_RGBA, GL_FLOAT,&cloud_d[0][0][0]);
glTexEnvf(GL_TEXTURE_ENV,GL_TEXTURE_ENV_MODE,GL_REPLACE);
glTexParameterf(GL_TEXTURE_RECTANGLE_NV,GL_TEXTURE_MIN_FILTER,GL_NEAREST);
glTexParameterf(GL_TEXTURE_RECTANGLE_NV,GL_TEXTURE_MAG_FILTER,GL_NEAREST);

/* omega can't be built into fragment shader due to dynamic index */
glActiveTextureARB(t_unit[3]);
glBindTexture(GL_TEXTURE_RECTANGLE_NV,t_id[3]);
glTexImage2D(GL_TEXTURE_RECTANGLE_NV,0,GL_FLOAT_R32_NV,DIRECTIONS,DIRECTIONS,
		0,GL_RED, GL_FLOAT,&omega[0][0]);
glTexEnvf(GL_TEXTURE_ENV,GL_TEXTURE_ENV_MODE,GL_REPLACE);
glTexParameterf(GL_TEXTURE_RECTANGLE_NV,GL_TEXTURE_MIN_FILTER,GL_NEAREST);
glTexParameterf(GL_TEXTURE_RECTANGLE_NV,GL_TEXTURE_MAG_FILTER,GL_NEAREST);

/* 3-component textures here make it easier for the fragment shader */
glActiveTextureARB(t_unit[4]);
glBindTexture(GL_TEXTURE_RECTANGLE_NV,t_id[4]);
glTexImage2D(GL_TEXTURE_RECTANGLE_NV,0,GL_FLOAT_RGB32_NV,DIRECTIONS,1,
		0,GL_RGB, GL_FLOAT,&directions[0]);
glTexEnvf(GL_TEXTURE_ENV,GL_TEXTURE_ENV_MODE,GL_REPLACE);
glTexParameterf(GL_TEXTURE_RECTANGLE_NV,GL_TEXTURE_MIN_FILTER,GL_NEAREST);
glTexParameterf(GL_TEXTURE_RECTANGLE_NV,GL_TEXTURE_MAG_FILTER,GL_NEAREST);
glEnable(GL_TEXTURE_RECTANGLE_NV);
}

/* 
 * the Henyey-Greenstein scattering function; 
 * see Astrophysics Journal 93,70; g>0 is forward,
 * g<0 is backward, g=0 is isotropic 
 */
float hg(float g, float cos_theta)
{
float val, dval;

dval = 1.0 + g*g - 2.0*g*cos_theta;
val = (1.0 - g*g)/pow(dval,1.5);
return(val);
}


void load_collision_matrix()
{
/* omega[i,j] controls scattering from dir j into dir i */
float hg_row[DIRECTIONS];
float hg_sum;
float row_wt[DIRECTIONS];
int i,j;

row_wt[0] = 1.0;
for(i=1;i<=6;i++) row_wt[i] = 1.0/12.0;
for(i=7;i<DIRECTIONS;i++) row_wt[i] = 1.0/24.0;

/* row 0 controls absorption; column 0 controls retransmission */
omega[0][0] = -1.0;
for(i=1;i<DIRECTIONS;i++) omega[i][0] = row_wt[i];

/* make sure columns sum to zero */
for(j=1;j<DIRECTIONS;j++){ 
	omega[0][j] = SIGMA_A;
	hg_sum = 0.0;
	for(i=1;i<=6;i++){
		hg_row[i] = hg(G,dot(&univ[i],&univ[j]));
		hg_sum += hg_row[i];
		}
	for(i=1;i<=6;i++){
		omega[i][j] = 6.0*row_wt[i]*SIGMA_S*hg_row[i]/hg_sum;
		}
	hg_sum = 0.0;
	for(i=7;i<DIRECTIONS;i++){
		hg_row[i] = hg(G,dot(&univ[i],&univ[j]));
		hg_sum += hg_row[i];
		}
	for(i=7;i<DIRECTIONS;i++){
		omega[i][j] = 12.0*row_wt[i]*SIGMA_S*hg_row[i]/hg_sum;
		}
	omega[j][j] -= SIGMA_T;
	}
}

char *read_shader_program(char *filename) 
{
FILE *fp;
char *content = NULL;
int fd, count;
fd = open(filename,O_RDONLY);
count = lseek(fd,0,SEEK_END);
close(fd);
content = (char *)calloc(1,(count+1));
fp = fopen(filename,"r");
count = fread(content,sizeof(char),count,fp);
content[count] = '\0';
fclose(fp);
return content;
}

unsigned int shader_std, shader_rs, shader_rn;

void load_shaders()
{
GLint vertCompiled, fragCompiled;
char *vs, *fs;
GLuint v, f, p;
int logsize, written;
char *logptr;

v = glCreateShader(GL_VERTEX_SHADER);
f = glCreateShader(GL_FRAGMENT_SHADER);
vs = read_shader_program("xrgba.vert");
fs = read_shader_program("xrgba.frag");
glShaderSource(v,1,(const char **)&vs,NULL);
glShaderSource(f,1,(const char **)&fs,NULL);
free(vs);
free(fs); 
glCompileShader(v);
glCompileShader(f);
glGetShaderiv(f,GL_COMPILE_STATUS,&fragCompiled);
if(!fragCompiled){
	fprintf(stderr,"xrgba.frag compile status is %d\n",fragCompiled);
        glGetShaderiv(f,GL_INFO_LOG_LENGTH,&logsize);
        logptr = (char *)(calloc(1,logsize));
        glGetShaderInfoLog(f,logsize,&written,logptr);
        fprintf(stderr,"log is:%s\n",logptr);
        free(logptr);
        exit(1);
        }

p = glCreateProgram();
glAttachShader(p,f);
glAttachShader(p,v);
glLinkProgram(p);
shader_std = p;

v = glCreateShader(GL_VERTEX_SHADER);
f = glCreateShader(GL_FRAGMENT_SHADER);
vs = read_shader_program("xrgba.vert");
fs = read_shader_program("xrgba_rs.frag");
glShaderSource(v,1,(const char **)&vs,NULL);
glShaderSource(f,1,(const char **)&fs,NULL);
free(vs);
free(fs); 
glCompileShader(v);
glCompileShader(f);
glGetShaderiv(f,GL_COMPILE_STATUS,&fragCompiled);
if(!fragCompiled){
	fprintf(stderr,"xrgba_rs.frag compile status is %d\n",fragCompiled);
        glGetShaderiv(f,GL_INFO_LOG_LENGTH,&logsize);
        logptr = (char *)(calloc(1,logsize));
        glGetShaderInfoLog(f,logsize,&written,logptr);
        fprintf(stderr,"log is:%s\n",logptr);
        free(logptr);
        exit(1);
        }

p = glCreateProgram();
glAttachShader(p,f);
glAttachShader(p,v);
glLinkProgram(p);
shader_rs = p;

v = glCreateShader(GL_VERTEX_SHADER);
f = glCreateShader(GL_FRAGMENT_SHADER);
vs = read_shader_program("xrgba.vert");
fs = read_shader_program("xrgba_rn.frag");
glShaderSource(v,1,(const char **)&vs,NULL);
glShaderSource(f,1,(const char **)&fs,NULL);
free(vs);
free(fs); 
glCompileShader(v);
glCompileShader(f);
glGetShaderiv(f,GL_COMPILE_STATUS,&fragCompiled);
if(!fragCompiled){
	fprintf(stderr,"xrgba_rn.frag compile status is %d\n",fragCompiled);
        glGetShaderiv(f,GL_INFO_LOG_LENGTH,&logsize);
        logptr = (char *)(calloc(1,logsize));
        glGetShaderInfoLog(f,logsize,&written,logptr);
        fprintf(stderr,"log is:%s\n",logptr);
        free(logptr);
        exit(1);
        }

p = glCreateProgram();
glAttachShader(p,f);
glAttachShader(p,v);
glLinkProgram(p);
shader_rn = p;
}

unsigned int do_stuff()
{
int i,j,k,m;
int dim[3];
int size;
int cloud_sampler_addr, photon_sampler_addr, omega_sampler_addr, 
	directions_sampler_addr;
int tcount;

glGetProgramiv(shader_std,GL_VALIDATE_STATUS,&tcount);
fprintf(stderr,"status is %d\n",tcount);

glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glOrtho(0.0,(double)(EDGE*DIRECTIONS),0.0,(double)((EDGE/4)*EDGE),-1.0,1.0);
glViewport(0,0,EDGE*DIRECTIONS,(EDGE/4)*EDGE);
glReadBuffer(t_attach[0]);

for(i=0;i<FINAL_TIME;i++){
	j = i%2;
	glUseProgram(shader_std);
	cloud_sampler_addr = glGetUniformLocation(shader_std,"cloud");
	glUniform1i(cloud_sampler_addr,2);
	omega_sampler_addr = glGetUniformLocation(shader_std,"omega");
	glUniform1i(omega_sampler_addr,3);
	directions_sampler_addr = glGetUniformLocation(shader_std,"directions");
	glUniform1i(directions_sampler_addr,4);
	photon_sampler_addr = glGetUniformLocation(shader_std,"photon");
	glUniform1i(photon_sampler_addr,j);
	glDrawBuffer(t_attach[1-j]);

	for(k=1;k<(EDGE/4)-1;k++){
		// write strip of width DIRECTIONS*(EDGE-2), height (EDGE-2)
		glBegin(GL_QUADS);
		glTexCoord2i(1*DIRECTIONS,k*EDGE+(EDGE-1));
		glVertex2i(1*DIRECTIONS,k*EDGE+(EDGE-1));

		glTexCoord2i(1*DIRECTIONS,k*EDGE+1);
		glVertex2i(1*DIRECTIONS,k*EDGE+1);

		glTexCoord2i((EDGE-1)*DIRECTIONS,k*EDGE+1);
		glVertex2i((EDGE-1)*DIRECTIONS,k*EDGE+1);

		glTexCoord2i((EDGE-1)*DIRECTIONS,k*EDGE+(EDGE-1));
		glVertex2i((EDGE-1)*DIRECTIONS,k*EDGE+(EDGE-1));
		glEnd();
		glFlush();
		glFinish();
		}

	// adjust lower boundary in layers r, g, b; 
	// read from upper boundary in layers g, b, a below
	glUseProgram(shader_rs);
	cloud_sampler_addr = glGetUniformLocation(shader_rs,"cloud");
	glUniform1i(cloud_sampler_addr,2);
	omega_sampler_addr = glGetUniformLocation(shader_rs,"omega");
	glUniform1i(omega_sampler_addr,3);
	directions_sampler_addr = glGetUniformLocation(shader_rs,"directions");
	glUniform1i(directions_sampler_addr,4);
	photon_sampler_addr = glGetUniformLocation(shader_rs,"photon");
	glUniform1i(photon_sampler_addr,j);
	glDrawBuffer(t_attach[1-j]);
	glBegin(GL_QUADS);
	k = (EDGE/4)-1; 
	glTexCoord2i(1*DIRECTIONS,k*EDGE+(EDGE-1));
	glVertex2i(1*DIRECTIONS,k*EDGE+(EDGE-1));
	glTexCoord2i(1*DIRECTIONS,k*EDGE+1);
	glVertex2i(1*DIRECTIONS,k*EDGE+1);
	glTexCoord2i((EDGE-1)*DIRECTIONS,k*EDGE+1);
	glVertex2i((EDGE-1)*DIRECTIONS,k*EDGE+1);
	glTexCoord2i((EDGE-1)*DIRECTIONS,k*EDGE+(EDGE-1));
	glVertex2i((EDGE-1)*DIRECTIONS,k*EDGE+(EDGE-1));
	glEnd();
	glFlush();
	glFinish();

	// adjust upper boundary in layers b, g, a; 
	// read from lower boundary in layers r, g, b above
	glUseProgram(shader_rn);
	cloud_sampler_addr = glGetUniformLocation(shader_rn,"cloud");
	glUniform1i(cloud_sampler_addr,2);
	omega_sampler_addr = glGetUniformLocation(shader_rn,"omega");
	glUniform1i(omega_sampler_addr,3);
	directions_sampler_addr = glGetUniformLocation(shader_rn,"directions");
	glUniform1i(directions_sampler_addr,4);
	photon_sampler_addr = glGetUniformLocation(shader_rn,"photon");
	glUniform1i(photon_sampler_addr,j);
	glDrawBuffer(t_attach[1-j]);
	glBegin(GL_QUADS);
	k =  0;
	glTexCoord2i(1*DIRECTIONS,k*EDGE+(EDGE-1));
	glVertex2i(1*DIRECTIONS,k*EDGE+(EDGE-1));
	glTexCoord2i(1*DIRECTIONS,k*EDGE+1);
	glVertex2i(1*DIRECTIONS,k*EDGE+1);
	glTexCoord2i((EDGE-1)*DIRECTIONS,k*EDGE+1);
	glVertex2i((EDGE-1)*DIRECTIONS,k*EDGE+1);
	glTexCoord2i((EDGE-1)*DIRECTIONS,k*EDGE+(EDGE-1));
	glVertex2i((EDGE-1)*DIRECTIONS,k*EDGE+(EDGE-1));
	glEnd();
	glFlush();
	glFinish();
	}
// return the correct readbuffer attachment for readback and writeout
return(1-j);
}

void dump_lights(int attach_index)
{
int i, j, k, m, dim[3], size;
float *lights, rho;

glBindFramebufferEXT(GL_FRAMEBUFFER_EXT,1);
glReadBuffer(t_attach[attach_index]);

/* with any luck, this puts them back in the correct positions */
glReadPixels(0,0,EDGE*DIRECTIONS,(EDGE/4)*EDGE,GL_RGBA,GL_FLOAT,
	&photon_d[0][0][0]);

dim[0] = dim[1] = dim[2] = EDGE;
size = EDGE*EDGE*EDGE;
lights = (float *)(calloc(size,sizeof(float)));
for(i=0;i<EDGE;i++){
	for(j=0;j<EDGE;j++){
		for(k=0;k<EDGE;k++){
			rho = 0.0;
			for(m=0;m<DIRECTIONS;m++){
				rho += photon_d[(k%(EDGE/4))*EDGE+i][j*DIRECTIONS+m][k/(EDGE/4)];
				}
			lights[i*EDGE*EDGE+j*EDGE+k]= rho/((float)(DIRECTIONS));
			}
		}
	}
fwrite(dim,sizeof(int),3,stdout);
fwrite(lights,sizeof(float),size,stdout);
free(lights);
}

int main(int argc, char **argv)
{
int shader_program;
int photon_d_addr;
unsigned int readbuffer_attach;

glutCreateWindow("zoom lighting");
fix_univ();
load_collision_matrix();
init_lattice(argv[1]);
load_textures();
load_shaders();
readbuffer_attach = do_stuff();
dump_lights(readbuffer_attach);
return 0;
}