// constants
int ANTIVOXEL = 1;

// user settings
int VIEWDIST = 240;
float RAYSIZE = 5;
float RELIGHT = 0.025; // amount of light regained by an empty pocket.
int RENDERTYPE = 0;//ANTIVOXEL;

class storm {
	int x, y, z; // temporary variables
	int motes;
	mote m[];

	class mote {
		float x, y, z;
		float xvel, yvel, zvel;
		float opacity = 0.08;
		// smoke : 0-0.2

		mote(float xIn, float yIn, float zIn) {
			x = xIn;
			y = yIn;
			z = zIn;
		}

		void update() {
			x += xvel;
			y += yvel;
			z += zvel;
		}
	}

	storm(int motesIn) {
		motes = motesIn;
		m = new mote[motes];
		createMotes();
	}

	void createMotes() {
		float ah, av, vd;
		float mx, my, mz;

		for(int i = 0; i < motes; i++) {
			ah = random(TWO_PI);
			av = random(TWO_PI);
			vd = random(15);

			// sphere
			mx = (50+vd*cos(ah)*cos(av));
			my = (80+vd*sin(av));//random(40,100);//
			mz = (50+vd*sin(ah)*cos(av));

			m[i] = new mote(mx, my, mz);
		}
	}

	void explode(float ex, float ey, float ez, float force) {
		float ah; // horizontal angle
		float ad; // distance
		float av; // vert angle

		for(int i = 0; i < motes; i++) {
			ah = atan2( (m[i].z-ez), (m[i].x-ex) );
			ad = dist(m[i].x, m[i].y, ex, ey);
			av = atan2( ad, (m[i].y-ey) );

			m[i].xvel += (cos(ah)*cos(av))*(force/ad);
			m[i].yvel += (sin(av))*(force/ad);
			m[i].zvel += (sin(ah)*cos(av))*(force/ad);
		}
	}

	void update() {
		for(int i = 0; i < motes; i++) {
			m[i].update();
		}
	}

	void fillCloud(cloud c) {
		// calculate opacity for all ze pockets.
		for(int i = 0; i < motes; i++) {
			x = (int)m[i].x;
			y = (int)m[i].y;
			z = (int)m[i].z;

			if(x >=0 && y >= 0 && z >= 0 && x < c.xmax && y < c.ymax && z < c.zmax) {
				c.p[x][y][z].motes++;
				c.p[x][y][z].opacity += m[i].opacity / c.p[x][y][z].motes;

				if(c.p[x][y][z].opacity > 1) {
					c.p[x][y][z].opacity = 1;
				}
			}
		}
	}
}

class cloud {
	int xmax, ymax, zmax, area;
	pocket p[][][];

	int x, y, z; // temporary variables

	class pocket {
		// 1-pixel cube area in the cloud
		// opacity is the sum of all the contained motes
		// shade is the sum of all pockets' opacities between it and the light source
		float opacity, shade;
		int motes;
	}

	cloud(int xIn, int yIn, int zIn) {
		xmax = xIn;
		ymax = yIn;
		zmax = zIn;
		area = xmax*ymax*zmax;

		p = new pocket[xmax][ymax][zmax];
		pocketsCreate();
		// pocket coordinate system is (x-horizontal +right),(y-vertical +down), (z-pokey +far);
	}

	void pocketsCreate() {
		for(int i = 0; i < area; i++) {
			p[i%xmax][(i/ymax)%ymax][i/(xmax*ymax)] = new pocket();
		}
	}

	void pocketsClear() {
		for(int i = 0; i < area; i++) {
			x = i%xmax;
			y = (i/ymax)%ymax;
			z = i/(xmax*ymax);
			p[x][y][z].opacity = 0;
			p[x][y][z].shade = 0;
			p[x][y][z].motes = 0;
		}
	}

	void pocketsShade() {
		for(int i = 0; i < area; i++) {
			x = i%xmax;
			y = (i/ymax)%ymax;
			z = i/(xmax*ymax);

			// cycle down through pockets, calculating shading.
			if(p[x][y][z].shade < 1 && y > 0 && x > 0 && z > 0) {
				p[x][y][z].shade = p[x-1][y-1][z].shade + p[x-1][y-1][z].opacity;

				if(p[x][y][z].opacity == 0 && p[x][y][z].shade >= RELIGHT) {
					p[x][y][z].shade -= RELIGHT;
				}

				if(p[x][y][z].shade > 1) {
					p[x][y][z].shade = 1;
				}
				if(y == ymax-1 && x > 1 && z > 1) {
					p[x][y][z].opacity = 1;
				}
			}
		}
	}

	void cheapRender() {
		for(int i = 0; i < area; i++) {
			x = i%xmax;
			y = (i/ymax)%ymax;
			z = (zmax-1)-(i/(xmax*ymax)); // To draw from back > front.
			if(p[x][y][z].opacity > 0) {
				int col = 255-(int)(128 * p[x][y][z].shade);
				stroke(col, col, col, (int)(p[x][y][z].opacity * 255));
				point(x,y,-z);
			}
		}
	}

	boolean validPos(int px, int py, int pz) {
		return (px >= 0 && px < xmax && py >= 0 && py < ymax && pz >= 0 && pz < zmax);
	}

	void fullRender(float cX, float cY, float cZ, float tX, float tY, float tZ, float fov) {
		fov = radians(fov);

		// set up viewing portal
		float angleH = atan2(tZ-cZ, tX-cX);
		float angleV = atan2(tY-cY, dist(cX, cZ, tX, tZ));

		float baseIncX = cos(angleH)*cos(angleV);
		float baseIncY = sin(angleV);
		float baseIncZ = sin(angleH)*cos(angleV);


    float angleHInc = fov/width;
		float angleVInc = angleHInc;//fov/height;

    float angleVMin = angleV - (angleVInc * height/2);//fov/2;
    float angleHMin = angleH - fov/2 * cos(angleV);



		float angleX, angleY, incX, incY, incZ;
		float incRay;
		float posX, posY, posZ;
		int iposX, iposY, iposZ;
		float rayOpacity, rayShade;
		float rayAddOpacity, rayAddShade;

		// components for ray-intersection method
		// boolean crossesX, crossesY, crossesZ;
		// float crsX, crsY, crsZ;

		// components for antialiasing
		float xC, yC, zC, portion;

		for(int i = 0; i < width*height; i++) {
			x = i%width;
			y = i/width;

			angleX = angleHMin + (x * angleHInc);
			angleY = angleVMin + (y * angleVInc);

			incX = cos(angleX)*cos(angleY);
			incY = sin(angleY);
			incZ = sin(angleX)*cos(angleY);

			rayOpacity = 0.0;
			rayShade = 0.0;

			for(int u = 0; u < VIEWDIST*RAYSIZE; u++) {
				incRay = (float)u/RAYSIZE;
				posX = (cX+(incRay*incX));
				posY = (cY+(incRay*incY));
				posZ = (cZ+(incRay*incZ));

				if(posX >= -1 && posX <= xmax && posY >= -1 && posY <= ymax && posZ >= -1 && posZ <= zmax) {
					iposX = (int)posX;
					iposY = (int)posY;
					iposZ = (int)posZ;

					rayAddOpacity = 0.0;
					rayAddShade = 0.0;

					/*
					// Ray Intersection method
					// Ensure RAYSIZE is set to 1 or less.
					// Look at the 1-pixel long ray, figure out what pockets it passes through, and
					// how much of its length lies in those pockets. Then add the opacity/shade of
					// those pockets using the ratio of (amount of ray in pocket):(total ray length).

					crossesX = false;
					crossesY = false;
					crossesZ = false;

					if(iposX != (int)(posX+incX)) {
						crossesX = true;
						crsX = math.round(posX+incX/2);
					}
					if(iposY != (int)(posY+incY)) {
						crossesY = true;
						crsY = math.round(posY+incY/2);
					}
					if(iposZ != (int)(posZ+incZ)) {
						crossesZ = true;
						crsZ = math.round(posZ+incZ/2);
					}
					*/


					// 3d voxel antialiasing
					// Most noticable when close, for long-range shots, it's probably best to use Nearest Pocket
					// with RAYSIZE = 5 or greater
					if(RENDERTYPE == ANTIVOXEL) {
						xC = posX - iposX;
						yC = posY - iposY;
						zC = posZ - iposZ;

						if(validPos(iposX, iposY, iposZ)) {
							portion = abs((1-xC)*(1-yC)*(1-zC));
							rayAddOpacity += p[iposX][iposY][iposZ].opacity * portion;
							rayAddShade += p[iposX][iposY][iposZ].opacity * p[iposX][iposY][iposZ].shade * portion;
						}

						if(validPos(iposX+1, iposY, iposZ)) {
							portion = abs(xC*(1-yC)*(1-zC));
							rayAddOpacity += p[iposX+1][iposY][iposZ].opacity * portion;
							rayAddShade += p[iposX+1][iposY][iposZ].opacity * p[iposX+1][iposY][iposZ].shade * portion;
						}

						if(validPos(iposX, iposY+1, iposZ)) {
							portion = abs((1-xC)*yC*(1-zC));
							rayAddOpacity += p[iposX][iposY+1][iposZ].opacity * portion;
							rayAddShade += p[iposX][iposY+1][iposZ].opacity * p[iposX][iposY+1][iposZ].shade * portion;
						}

						if(validPos(iposX+1, iposY+1, iposZ)) {
							portion = abs(xC*yC*(1-zC));
							rayAddOpacity += p[iposX+1][iposY+1][iposZ].opacity * portion;
							rayAddShade += p[iposX+1][iposY+1][iposZ].opacity * p[iposX+1][iposY+1][iposZ].shade * portion;
						}

						// Lower quad
						if(validPos(iposX, iposY, iposZ+1)) {
							portion = abs((1-xC)*(1-yC)*zC);
							rayAddOpacity += p[iposX][iposY][iposZ+1].opacity * portion;
							rayAddShade += p[iposX][iposY][iposZ+1].opacity * p[iposX][iposY][iposZ+1].shade * portion;
						}

						if(validPos(iposX+1, iposY, iposZ+1)) {
							portion = abs(xC*(1-yC)*zC);
							rayAddOpacity += p[iposX+1][iposY][iposZ+1].opacity * portion;
							rayAddShade += p[iposX+1][iposY][iposZ+1].opacity * p[iposX+1][iposY][iposZ+1].shade * portion;
						}

						if(validPos(iposX, iposY+1, iposZ+1)) {
							portion = abs((1-xC)*yC*zC);
							rayAddOpacity += p[iposX][iposY+1][iposZ+1].opacity * portion;
							rayAddShade += p[iposX][iposY+1][iposZ+1].opacity * p[iposX][iposY+1][iposZ+1].shade * portion;
						}

						if(validPos(iposX+1, iposY+1, iposZ+1)) {
							portion = abs(xC*yC*zC);
							rayAddOpacity += p[iposX+1][iposY+1][iposZ+1].opacity * portion;
							rayAddShade += p[iposX+1][iposY+1][iposZ+1].opacity * p[iposX+1][iposY+1][iposZ+1].shade * portion;
						}
					}
					else {
					// Standard Nearest Pocket
						if(validPos(iposX, iposY, iposZ)) {
							rayAddOpacity = p[iposX][iposY][iposZ].opacity;
							rayAddShade = p[iposX][iposY][iposZ].opacity * p[iposX][iposY][iposZ].shade;
						}
					}

					rayOpacity += rayAddOpacity/RAYSIZE;
					rayShade += rayAddShade/RAYSIZE;

					if(rayOpacity > 1) {
						// 'normalise' ray
						rayShade -= (rayOpacity-1)*rayAddShade/RAYSIZE;
						rayOpacity = 1;
						break;
					}
				}
				else {
					// skip through space;
					if(RAYSIZE > 2) u += (RAYSIZE-1);
				}
			}
			//rayOpacity = constrain(rayOpacity,0,2);
			int col = (int)(255-(rayShade*128));
			stroke(col,col,col);
			modpixel(i, stroke(), (int)((1-rayOpacity)*255));
		}
	}
}

void modpixel(int pos, int col, int alp) {
	int r1=(pixels[pos])&0x0000ff;
	int g1=(pixels[pos])&0x00ff00;
	int b1=(pixels[pos])&0xff0000;
	int r2=(col)&0x0000ff;
	int g2=(col)&0x00ff00;
	int b2=(col)&0xff0000;
	int r3=(((alp*(r1-r2))>>8)+r2)&0x000000ff;
	int g3=(((alp*(g1-g2))>>8)+g2)&0x0000ff00;
	int b3=(((alp*(b1-b2))>>8)+b2)&0x00ff0000;
	pixels[pos]=(r3)|(g3)|(b3);
	}

cloud cld;
storm stm;

void setup() {
	size(320,200);
	smooth();
	background(0xff8D7D5D);
	cld = new cloud(100,100,100);
	stm = new storm(50000);

	stm.fillCloud(cld);
	cld.pocketsShade();
	//c.cheapRender();
}

float a = PI/3;
void loop() {
	stm.explode(60,75,60, 4);
	//a+=PI/120;
	if(a < PI/3+TWO_PI) {
		stm.update();
		cld.pocketsClear();
		stm.fillCloud(cld);
		cld.pocketsShade();
		cld.fullRender(100,-50,100, 50,75,50, 90);
		//cld.fullRender(50+(150*cos(a)),25,50+(150*sin(a)), 50,75,50, 45);
		//screenGrab();
	}

}