# include <stdio.h>
# include <stdlib.h>

int main ( );
__global__ void julia_set ( int w, int h, float xl, float xr, float yb, 
  float yt, int julia_gpu[] );
__device__ int julia_value ( int w, int h, float xl, float xr, float yb, 
  float yt, int i, int j );
void tga_write ( int w, int h, unsigned char rgb[], const char *filename );

/******************************************************************************/

int main ( ) 

/******************************************************************************/
/*
  Purpose:

    MAIN is the main program for JULIA_GPU.

  Discussion:

    Consider points (X,Y) in a rectangular domain R = [XL,XR]x[YB,YT].

    Let Z be the complex number X+Yi, and let C be some complex constant.

    Let Z(0) = Z, Z(k+1) = Z(k)^2 + C

    The Julia set is the set of points Z in R with the property that
    the sequence of points Z(k) remain within R.

    To compute a picture of the Julia set, we choose a discrete array
    of WxH points in R.  We carry out up to 200 steps of the iteration for
    each point Z.  If 1000 < |Z| at any time, we assume Z is not in the
    Julia set.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    07 March 2017

  Parameters:

    Local, int H, W, the height and width of the region in pixels.

    Local, float XL, XR, YB, YT, the left and right X limits, the
    bottom and top Y limits, of the region.

    Local, unsigned char *RGB_CPU, will hold W*H*3 values between 0 and 255,
    specifying the pixel color values.

    Local, unsigned char *RGB_GPU, a pointer to the RGB data on the GPU.
*/
{
  int h = 1000;
  int i;
  int j;
  int *julia_cpu;
  int *julia_gpu;
  int julia_size;
  int k;
  unsigned char *rgb;
  int w = 1000;
  float xl = - 1.5;
  float xr = + 1.5;
  float yb = - 1.5;
  float yt = + 1.5;

  printf ( "\n" );
  printf ( "JULIA_GPU:\n" );
  printf ( "  CUDA version.\n" );
  printf ( "  Plot a version of the Julia set for Z(k+1)=Z(k)^2-0.8+0.156i\n" );
/*
  Allocate GPU memory for the JULIA array.
*/
  julia_size = w * h * sizeof ( int );
  cudaMalloc ( ( void ** ) &julia_gpu, julia_size );
/*
  JULIA_SET runs on the GPU.
  Run on a WxHx1 grid of blocks, and 1 thread.
*/
  dim3 blocks ( w, h, 1 );
  int threads = 1;
  julia_set <<< blocks, threads >>> ( w, h, xl, xr, yb, yt, julia_gpu );
/*
  Copy JULIA data from GPU to CPU.
*/
  julia_cpu = ( int * ) malloc ( julia_size );
  cudaMemcpy ( julia_cpu, julia_gpu, julia_size, cudaMemcpyDeviceToHost );
/*
  Create an RGB image of the data.
*/
  k = 0;
  rgb = ( unsigned char * ) malloc ( w * h * 3 * sizeof ( unsigned char ) );

  for ( j = 0; j < h; j++ )
  {
    for ( i = 0; i < w; i++ )
    {
/*
  TARGA format orders colors B/G/R!
*/
      rgb[k]   = 255 * ( 1 - julia_cpu[i+j*w] );
      rgb[k+1] = 255 * ( 1 - julia_cpu[i+j*w] );
      rgb[k+2] = 255;

      k = k + 3;
    }
  }
/*
  Write RGB data to TARGA file.
*/
  tga_write ( w, h, rgb, "julia_gpu.tga" );
/*
  Free CPU memory.
*/
  free ( julia_cpu );
  free ( rgb );
/*
  Free GPU memory.
*/
  cudaFree ( julia_gpu );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "JULIA_GPU:\n" );
  printf ( "  Normal end of execution.\n" );

  return 0;
}
/******************************************************************************/

__global__ void julia_set ( int w, int h, float xl, float xr, float yb, 
  float yt, int julia_gpu[] )

/******************************************************************************/
/*
  Purpose:

    JULIA_SET checks all the points.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    06 March 2017

  Parameters:

    Input, int W, H, the width and height of the region in pixels.

    Input, float XL, XR, YB, YT, the left, right, bottom and top limits.

    Output, int JULIA_GPU[W*H], contains 1 for points in the Julia set,
    0 otherwise.
*/
{
  int i;
  int j;
  int k;

  i = blockIdx.x;
  j = blockIdx.y;
  k = i + j * w;
/*
  Determine if point Z=(X,Y) with indices (I,J) is in the Julia set.
*/
  julia_gpu[k] = julia_value ( w, h, xl, xr, yb, yt, i, j );

  return;
}
/******************************************************************************/

__device__ int julia_value ( int w, int h, float xl, float xr, float yb, 
  float yt, int i, int j )

/******************************************************************************/
/*
  Purpose:

    JULIA_VALUE checks one point.

  Discussion:

    The iteration Z(k+1) = Z(k) + C is used, with C=-0.8+0.156i.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    06 March 2017

  Parameters:

    Input, int W, H, the width and height of the region in pixels.

    Input, float XL, XR, YB, YT, the left, right, bottom and top limits.

    Input, int I, J, the indices of the point to be checked.

    Ouput, int JULIA, is 1 if the point is in the Julia set.
*/
{
  float ai;
  float ar;
  float ci = 0.156;
  float cr = -0.8;
  int k;
  float t;
  int value;
  float x;
  float y;
/*
  Convert (I,J) indices to (X,Y) coordinates.
*/
  x = ( ( float ) ( w - i - 1 ) * xl
      + ( float ) (     i     ) * xr ) 
      / ( float ) ( w     - 1 );

  y = ( ( float ) ( h - j - 1 ) * yb
      + ( float ) (     j     ) * yt ) 
      / ( float ) ( h     - 1 );
/*
  Think of (X,Y) as real and imaginary components of
  a complex number A = x + y*i.
*/
  ar = x;
  ai = y;
/*
  A -> A * A + C
*/
  value = 1;

  for ( k = 0; k < 200; k++ )
  {
    t  = ar * ar - ai * ai + cr;
    ai = ar * ai + ai * ar + ci;
    ar = t;
/*
  if 1000 < ||A||, reject the point.
*/
    if ( 1000.0 < ar * ar + ai * ai )
    {
      value = 0;
      return value;
    }
  }

  return value;
}
/******************************************************************************/

void tga_write ( int w, int h, unsigned char rgb[], const char *filename )

/******************************************************************************/
/*
  Purpose:

    TGA_WRITE writes a TGA or TARGA graphics file of the data.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    06 March 2017

  Parameters:

    Input, int W, H, the width and height of the image.

    Input, unsigned char RGB[W*H*3], the pixel data.

    Input, const char *FILENAME, the name of the file to contain the screenshot.
*/
{
  FILE *file_unit;
  unsigned char header1[12] = { 0,0,2,0,0,0,0,0,0,0,0,0 };
  unsigned char header2[6] = { w%256, w/256, h%256, h/256, 24, 0 };
/* 
  Create the file.
*/
  file_unit = fopen ( filename, "wb" );
/*
  Write the headers.
*/
  fwrite ( header1, sizeof ( unsigned char ), 12, file_unit );
  fwrite ( header2, sizeof ( unsigned char ), 6, file_unit );
/*
  Write the image data.
*/
  fwrite ( rgb, sizeof ( unsigned char ), 3 * w * h, file_unit );
/*
  Close the file.
*/
  fclose ( file_unit );

  printf ( "\n" );
  printf ( "TGA_WRITE:\n" );
  printf ( "  Graphics data saved as '%s'\n", filename );

  return;
}