I changed cartesian coordinate to spherical coordinate because i want to improve this code by General relativity, but when i tried to check this kernel when compile in Mathematica it has been error
this is my kernel
__kernel
void
nbody_sim(
__global float4* pos ,
__global float4* vel,
mint numBodies,
float deltaTime,
float epsSqr,
__local float4* localPos,
__global float4* newPosition,
__global float4* newVelocity)
{
unsigned int tid = get_local_id(0);
unsigned int gid = get_global_id(0);
unsigned int localSize = get_local_size(0);
// Number of tiles we need to iterate
unsigned int numTiles = numBodies / localSize;
// position of this work-item
float4 myPos = pos[gid];
float myPosR = sqrt(myPos.x * myPos.x + myPos.y * myPos.y + myPos.z * myPos.z);
float myPosTheta = acos(myPos.z / sqrt(myPos.x * myPos.x + myPos.y * myPos.y + myPos.z * myPos.z));
float myPosPhi = atan(myPos.y / myPos.z);
float4 myPosSphere = (float4)(myPosR, myPosTheta, myPosPhi, 1.0f);
float4 acc = (float4)(0.0f, 0.0f, 0.0f, 0.0f);
for(int i = 0; i < numTiles; ++i)
{
// load one tile into local memory
int idx = i * localSize + tid;
localPos[tid] = pos[idx];
float localPosR = sqrt(localPos[tid].x * localPos[tid].x + localPos[tid].y * localPos[tid].y + localPos[tid].z * localPos[tid].z);
float localPosTheta = acos(localPos[tid].z / sqrt(localPos[tid].x * localPos[tid].x + localPos[tid].y * localPos[tid].y + localPos[tid].z * localPos[tid].z));
float localPosPhi = atan(localPos[tid].y / localPos[tid].z);
float4 localPosSphere[tid] = (float4)(localPosR, localPosTheta, localPosPhi, 1.0f);
// Synchronize to make sure data is available for processing
barrier(CLK_LOCAL_MEM_FENCE);
// calculate acceleration effect due to each body
// a[i->j] = m[j] * r[i->j] / (r^2 + epsSqr)^(3/2)
for(int j = 0; j < localSize; ++j)
{
// Calculate acceleration caused by particle j on particle i
float4 r = localPosSphere[j] - myPosSphere;
float sineTheta = sin(localPosSphere[j].y);
float sineThetaPrime = sin(myPosSphere[j].y);
float cosineTheta = cos(localPosSphere[j].y);
float cosineThetaPrime = cos(myPosSphere.y);
float cosineDeltaPhi = cos(r.z);
float distSqr = localPosSphere.x * localPosSphere.x + myPosSphere.x * myPosSphere.x - (2.0f * localPosSphere.x * myPosSphere.x)*(sineTheta * sineThetaPrime * cosineDeltaPhi + (cosineTheta * cosineThetaPrime));
float invDist = 1.0f / sqrt(distSqr + epsSqr);
float invDistCube = invDist * invDist * invDist;
localPos[j].w = 1.0f;
float s = localPos[j].w * invDistCube;
// accumulate effect of all particles
acc += s * r;
}
// Synchronize so that next tile can be loaded
barrier(CLK_LOCAL_MEM_FENCE);
}
float4 oldVel = vel[gid];
float oldVelR = sqrt(oldVel.x * oldVel.x + oldVel.y * oldVel.y + oldVel.z * oldVel.z);
float oldVelTheta = acos(oldVel.z / sqrt(oldVel.x * oldVel.x + oldVel.y * oldVel.y + oldVel.z * oldVel.z);
float oldVelPhi = atan(oldVel.y / oldVel.z);
float4 oldVelSphere = (float4)(oldVelR, oldVelTheta, oldVelPhi, 0.0f);
// updated position and velocity
float4 newPosSphere = myPosSphere + oldVelSphere * deltaTime + acc * 0.5f * deltaTime * deltaTime;
newPosSphere.w = 1.0f;
float4 newVelSphere = oldVelSphere + acc * deltaTime;
// write to global memory
newPosition[gid] = newPosSphere;
newVelocity[gid] = newVelSphere;
}