#version 400

//the following uniform values are set by Visions of Chaos prior to shader execution
uniform vec2 resolution;
uniform vec3 palette[256];
uniform double xmin;
uniform double xmax;
uniform double ymin;
uniform double ymax;
uniform double bailout;
uniform int maxiters;
uniform int samplepixels;

double sqrsamplepixels=double(samplepixels*samplepixels);
double bailout_squared=double(bailout*bailout);
double magnitude=0.0,r1,r2,g1,g2,b1,b2,tweenval;
float realiters;
vec4 finalcol,col;
int superx,supery;
double stepx=(xmax-xmin)/resolution.x/double(samplepixels);
double stepy=(ymax-ymin)/resolution.y/double(samplepixels);
int index,colval,colval2;
dvec2 z,c;

double trapsize,sxmin,sxmax,symin,symax,trapdist,trapx,trapy,distpercentage;
bool trapped;

void main(void)
{
	//center of trap
   trapx=0;
   trapy=0.5;
	//square trap sizes
   trapsize=0.4;
   sxmin=trapx-trapsize;
   sxmax=trapx+trapsize;
   symin=trapy-trapsize;
   symax=trapy+trapsize;

	finalcol=vec4(0,0,0,0);
	for (supery=0;supery<samplepixels;supery++)
	{
		for (superx=0;superx<samplepixels;superx++)
		{
			c.x = xmin+gl_FragCoord.x/resolution.x*(xmax-xmin)+(stepx*double(superx));
			c.y = ymin+gl_FragCoord.y/resolution.y*(ymax-ymin)+(stepy*double(supery));
			int i;
			z = dvec2(0.0,0.0);

         trapped=false;

			for(i=0; i<maxiters; i++) 
			{
				//START OF FRACTAL FORMULA
				double x = (z.x * z.x - z.y * z.y) + c.x;
				double y = (z.y * z.x + z.x * z.y) + c.y;
				//END OF FRACTAL FORMULA
				
				magnitude=(x * x + y * y);
				if(magnitude>bailout_squared) break;
				z.x = x;
				z.y = y;

				//squares orbit traps
            if ((z.x>sxmin) && (z.x<sxmax) && (z.y>symin) && (z.y<symax)) {
               trapdist=(abs(z.x-trapx)+abs(z.y-trapy))/2;
					trapped=true;
				}
				if (trapped==true) break;

			}

			//trapped?
			if (trapped==true) {
				distpercentage=1-trapdist/trapsize;
				col=vec4(distpercentage,distpercentage,distpercentage,1);
			} else 
			if (i==maxiters) {
				col=vec4(0.0,0.0,0.0,1.0);
			} else {
				//CPM smooth colors
				//note that double precision does not support log so it needs to be cast as float
				realiters=float(i+1-((log(log(sqrt(float(magnitude))))/log(2.0))));
				colval=int(mod(realiters,255));
				colval2=int(mod(colval+1,255));
				tweenval=realiters-int(realiters);
				r1=palette[colval].r;
				g1=palette[colval].g;
				b1=palette[colval].b;
				r2=palette[colval2].r;
				g2=palette[colval2].g;
				b2=palette[colval2].b;
				col=vec4(r1+((r2-r1)*tweenval),g1+((g2-g1)*tweenval),b1+((b2-b1)*tweenval),1.0);
			}

			finalcol+=col;
		}
	}
	gl_FragColor = vec4(finalcol/double(sqrsamplepixels));
}