#include <vector>

#include <assert.h>
#include <math.h>
#include <stdio.h>

struct coord { double x, y, z; };

static const double C = 1.0;
static const double D = 0.5;

double compute_f(double phi) {
  return M_PI * ((exp(phi/M_PI/C + 0.5/C) - 1) / (exp(1/C) - 1) - 0.5);
}

double compute_r(double theta, double phi) {
  double delta_r;

  if (phi >= M_PI/2) {
    delta_r = 1;
  } else if (phi <= -M_PI/2) {
    delta_r = -1;
  } else {
    const double f1 = sin(compute_f(phi));
    const double f2 = -sin(compute_f(-phi));
    const double tt = tan(theta);
    const double a = 1 / (1 + tt*tt*(1 - f2*f2)/(1 - f1*f1));
    delta_r = a*f1 + (1 - a)*f2;
  }

  return 1.0 + D * delta_r;
}

void facet(coord v1, coord v2, coord v3) {
  coord n;
  n.x = (v1.y-v3.y) * (v2.z-v3.z) - (v1.z-v3.z) * (v2.y-v3.y);
  n.y = (v1.z-v3.z) * (v2.x-v3.x) - (v1.x-v3.x) * (v2.z-v3.z);
  n.z = (v1.x-v3.x) * (v2.y-v3.y) - (v1.y-v3.y) * (v2.x-v3.x);
  double norm = sqrt(n.x * n.x + n.y * n.y + n.z * n.z);
  if (norm > 1e-20) {
    printf("facet normal %g %g %g\n", n.x / norm, n.y / norm, n.z / norm);
    printf("  outer loop\n");
    printf("    vertex %g %g %g\n", v1.x, v1.y, v1.z);
    printf("    vertex %g %g %g\n", v2.x, v2.y, v2.z);
    printf("    vertex %g %g %g\n", v3.x, v3.y, v3.z);
    printf("  endloop\n");
    printf("endfacet\n");
  }
}

int main() {
  printf("solid gomboc\n");

  std::vector<coord> last, next;
  for (double phi = -M_PI/2; phi <= M_PI/2; phi += M_PI/36) {
    for (double theta = 0; theta < 2*M_PI; theta += M_PI/36) {
      double r = compute_r(theta, phi);
      coord xyz;
      xyz.x = cos(theta) * cos(phi) * r;      
      xyz.y = sin(theta) * cos(phi) * r;
      xyz.z = sin(phi) * r;
      next.push_back(xyz);
    }

    if (!last.empty()) {
      assert(last.size() == next.size());
      for (size_t i = 0; i < last.size(); ++i) {
        size_t j = (i + 1) % last.size();
        facet(last[i], next[j], next[i]);
        facet(next[j], last[i], last[j]);
      }
    }

    last.swap(next);
    next.clear();
  }

  printf("endsolid gomboc\n");
  return 0;
}