#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 ztot,stalksradius,stalksradiuslow,stalksradiushigh,stalkcolor;
bool trapped;

void main(void)
{
	stalksradius=0.05;

	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);

			//stalks orbit trap
			stalksradiushigh=sqrt(c.x*c.x+c.y*c.y)+stalksradius;
   			stalksradiuslow=sqrt(c.x*c.x+c.y*c.y)-stalksradius;
         ztot = 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;

				//stalks orbit traps
            if ((sqrt(magnitude)<=stalksradiushigh) && (sqrt(magnitude)>=stalksradiuslow) && (i>1)) {
					ztot=sqrt(ztot)+(1-abs(sqrt(magnitude)-sqrt(c.x*c.x+c.y*c.y))/stalksradius);
               trapped=true;
				}
				if (trapped==true) break;

			}

			//trapped?
			if (trapped==true) {
				stalkcolor=min(1,sqrt(ztot));
				col=vec4(stalkcolor,stalkcolor,stalkcolor,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));
}