/* OpenCL version of the andy galaxy code - based on the OpenCL Nvidia simpleGL example
and originally on the CUDA simpleGL example */

///////////////////////////////////////////////////////////////////////////////
//! version of the galaxy collision code using OpenCL
///////////////////////////////////////////////////////////////////////////////
__kernel void galaxy(__global float4* pos, __global float* pdata, unsigned int width, unsigned int height, unsigned int howManyStars, unsigned int time)
{
// pos contains vec4 positions for numStars followed by vec4 colours for numStars
// pdata contains 7 floats (x,y,z. vx,vy,vz, mass)


#define mass 0.00006
#define stepper (unsigned int) 10

    unsigned int x = get_global_id(0);
    unsigned int y = get_global_id(1);

	unsigned int arrayLoc = y*width*7 + x*7;
	unsigned int posLoc   = y*width+x;
	unsigned int correspondingLoc = width*height + posLoc;
	
	float3 affect = (float3)(0.0, 0.0, 0.0);
	float3 posit;
	float3 dif;

	float diag, invDistCube;
	int velx, vely, velz;
	unsigned int i,j;


if (posLoc < howManyStars)
{	
	j = time % stepper;

	// if first time through then initialize the position and colour data
	if (time < 2)
		{
		pos[posLoc].x = pdata[arrayLoc+0];
		pos[posLoc].y = pdata[arrayLoc+1];
		pos[posLoc].z = pdata[arrayLoc+2];
		pos[posLoc].w = 1.0; // used to be mass = pdata[arrayLoc+6];

		pos[correspondingLoc] = (float4)(1.0, 1.0, 1.0, 1.0);
		}

	posit.xyz = pos[posLoc].xyz;
	
	for (i=j; i <width * height; i+=stepper)
		{
		dif.xyz = pos[i].xyz - posit.xyz;

		diag = (dif.x * dif.x) + (dif.y * dif.y) + (dif.z * dif.z);
		if (diag < 1.0)
			diag = 1.0;

		invDistCube = 1.0f / sqrt(diag * diag * diag); // was sqrtf in CUDA

		affect += invDistCube * dif; // this is twice as fast as 
		
		//affect.x += invDistCube * dif.x; // this
		//affect.y += invDistCube * dif.y;
		//affect.z += invDistCube * dif.z;
		}

	velx =  arrayLoc+3;
	vely =  arrayLoc+4;
	velz =  arrayLoc+5;
	
	
	pdata[velx] += mass * affect.x;
	pdata[vely] += mass * affect.y;
	pdata[velz] += mass * affect.z;

	pos[posLoc].x += 0.001 * pdata[velx]; 
	pos[posLoc].y += 0.001 * pdata[vely];
	pos[posLoc].z += 0.001 * pdata[velz];

		
	float4 c = (float4)(1.0, 1.0, 1.0, 1.0);
	if (pdata[velx] < 0)
		c.x = 0.4;
	if (pdata[vely] < 0)
		c.y = 0.4;
	if (pdata[velz] < 0)
		c.z = 0.4;
	pos[correspondingLoc] = c;
		
	}
}