#include <stdio.h>
#include <math.h>
#include <glib/gtypes.h> // srsly guys?
#include <gts.h>

void mandelbulb(gdouble **a, GtsCartesianGrid g, guint i, gpointer data)
{
  guint j, k, n;
  gdouble x, y, z = g.z;
  gdouble xi, yi, zi;
  gdouble r, theta, phi, power = 8;
  gdouble *row;

  for(k = 0, y = g.y; k < g.ny; k++, y += g.dy)
  {
    row = a[k];
    for(j = 0, x = g.x; j < g.nx; j++, x += g.dx)
    {
      xi = 0;
      yi = 0;
      zi = 0;
      for(n = 0; n < 6; n++)
      {
        r = sqrt(xi * xi + yi * yi + zi * zi);
        theta = atan2(sqrt(xi * xi + yi * yi), zi);
        phi = atan2(yi, xi);
        xi = pow(r, power) * sin(theta * power) * cos(phi * power) + x;
        yi = pow(r, power) * sin(theta * power) * sin(phi * power) + y;
        zi = pow(r, power) * cos(theta * power) + z;
      }
      row[j] = sqrt(xi * xi + yi * yi + zi * zi);
      if(!isfinite(row[j]))
        row[j] = G_MAXDOUBLE;
    }
  }
}

int main(int argc, char **argv)
{
  GtsCartesianGrid grid;
  GtsSurface *s;

  grid.nx = 1000;
  grid.ny = 1000;
  grid.nz = 1000;
  grid.x = -1.2;
  grid.y = -1.2;
  grid.z = -1.2;
  grid.dx = 2.4 / grid.nx;
  grid.dy = 2.4 / grid.ny;
  grid.dz = 2.4 / grid.nz;
  s = gts_surface_new(gts_surface_class(), gts_face_class(),
    gts_edge_class(), gts_vertex_class());
  gts_isosurface_tetra(s, grid, (GtsIsoCartesianFunc)mandelbulb, NULL, 1e308);
  fprintf(stderr, "Origional surface:\n");
  gts_surface_print_stats(s, stderr);
  gts_surface_write(s, stdout);
  return 0;
}