END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

//A function to check if the choosen point is in the mandelbrot set
bool isInMSet(
  float cr,
  float ci,
  const unsigned int maxIter,
  const float escapeOrbit
)
{
    int iter = 0;
    float zr = 0.0;
    float zi = 0.0;
    float zr2 = zr * zr;
    float zi2 = zi * zi;
    float temp = 0.0;

    //Quick rejection check if c is in the 2nd order period bulb.
    if( sqrt( ((cr+1.0) * (cr+1.0)) + (ci * ci) ) < 0.25 )
    {
        return true;
    }

    //Quick rejection check if c is in the main cardiod
    //IF ((ABS( 1.0 - SQRT(1-(4*c)) ))  < 1.0 ) THEN RETURN TRUE (main cardioid)
    float tempi = ci*(-4.0);
    float tempr = 1.0 - cr*4.0;
    float theta = atan2(tempi,tempr)/2.0;
    float r = pow((tempr*tempr + tempi*tempi),0.25);
    tempr = 1.0 - r * cos(theta);
    tempi = -r * sin(theta);
    if( (tempr * tempr + tempi * tempi) < 1.0)
    {
        return true;
    }


    //Bruteforce check if c is escaping escapeOrbit (with a good old iteration up to maxIter)
    while( (iter < maxIter) && ((zr2 + zi2) < escapeOrbit) )
    {
        temp = zr * zi;
        zr2 = zr * zr;
        zi2 = zi * zi;
        zr = zr2 - zi2 + cr;
        zi = temp + temp + ci;
        iter++;
    }

    if ( iter < maxIter )
    {
        return false;
    } else {
        return true;
    }

}   

//Main kernel

__kernel void buddhabrot(
  const float realMax,
  const float imaginaryMax,
  const float realMin,
  const float imaginaryMin,
  const unsigned int maxIter,
  const unsigned int escapeOrbit,
  const unsigned int hRes,
  const float offset,
  __global int* outputi
) {

  const int xId = get_global_id(0);
  const int yId = get_global_id(1);
  const int offsetStep = get_global_id(2);

  const int maxX = get_global_size(0);
  const int maxY = get_global_size(1);

  const float deltaReal = (realMax - realMin) / (maxX - 1);
  const float deltaImaginary = (imaginaryMax - imaginaryMin) / (maxY - 1);

  float cr = realMin + (xId * deltaReal) + (offsetStep * offset);
  float ci = imaginaryMin + (yId * deltaImaginary);

  int iter = 0;
  float zr = 0.0;
  float zi = 0.0;
  float zr2 = zr * zr;
  float zi2 = zi * zi;
  float temp = 0.0;
  int x, y;
   
  if(isInMSet(cr, ci, maxIter, escapeOrbit) == false)
  {
    iter = 0;
    zr = 0.0;
    zi = 0.0;
    zr2 = zr * zr;
    zi2 = zi * zi;
    temp = 0.0;
    while ((iter < maxIter) && ((zr2 + zi2) < escapeOrbit) )
    {
      temp = zr * zi;
      zr2 = zr * zr;
      zi2 = zi * zi;
      zr = zr2 - zi2 + cr;
      zi = temp + temp + ci;
      x = (maxX * (zr - realMin) / (realMax - realMin));
      y = (maxY * (zi - imaginaryMin) / (imaginaryMax - imaginaryMin));
       
      if( (x > 0) && (y > 0) && (x < maxX) && (y < maxY) && (iter > 2))
      {
        outputi[(y * hRes) + x] += 1;
      }
      iter++;
    } //EndWhile
  } //EndIf
} //EndKernel