#include <stdio.h>
#include <math.h>
#include "blob.h"
#include <stdlib.h>
Pt centr(int j){T t = Q[j];// return centroid of three verts that define sliver.
  return sm(1./3, add(q[t.a], add(q[t.b], q[t.c])));}
R ip(Pt a, Pt b){return a.x*b.x + a.y*b.y + a.z*b.z;}
R det(Pt a, Pt b, Pt c){return ip(a, cp(b, c));}
Pt sm(R a, Pt b){return (Pt){a*b.x, a*b.y, a*b.z};}
R rf(){return random()*0x1.p-31;}
R len(Pt x){return sqrt(ip(x, x));}
Pt cp(Pt a, Pt b){return (Pt){ // Cross product
  a.y*b.z - a.z*b.y,
  a.z*b.x - a.x*b.z,
  a.x*b.y - a.y*b.x};}
Pt nv(int j){ // compute the normal vector for the triangle for this sliver.
  if(0) {prP("d1", sub(q[Q[j].b], q[Q[j].a])); prP("d2", sub(q[Q[j].c], q[Q[j].b]));}
  return sm(0.5, cp(sub(q[Q[j].b], q[Q[j].a]), sub(q[Q[j].c], q[Q[j].b])));}
Pt puff(int n){Pt s = {0, 0, 0}; for(int c = h[n]; c<(1<<16)-1; c = L[c]) {
  unsigned int I = c/3, i = c - 3*I; T x = Q[I];
  unsigned int jk = i==0?(x.b<<16|x.c):i==2?(x.a<<16|x.b):(x.c<<16|x.a);
  s = add(s, cp(q[jk>>16], q[jk&0xffff]));} return sm(0.5, s);}
Pt add(Pt a, Pt b){return (Pt){a.x + b.x, a.y + b.y, a.z + b.z};}
Pt sub(Pt a, Pt b){return (Pt){a.x - b.x, a.y - b.y, a.z - b.z};}
void prP(char * s, Pt sl) {
  printf("%15.12f %15.12f %15.12f %s\n", sl.x,  sl.y, sl.z, s);}
void pR(R x){printf("%14.11f ", x);}
void rm(char * c, R m0, R m1[3], R m2[3][3]){  // print moments
  printf("m0 = %13.10f, %s\n", m0, c);
  for(int j=0; j<3; ++j) pR(m1[j]); printf("m1\n");
  for(int j=0; j<3; ++j) {for(int k=0; k<3; ++k) pR(m2[j][k]); printf("m2\n");}}
R moments(R m1[3], R m2[3][3]){// Moments:
  R m0 = 0; for(int j=0; j<3; ++j) {m1[j] = 0; for(int k=0; k<3; ++k) m2[j][k] = 0;}
  for(int j=0; j<SN; ++j) {Pt c = centr(j); R * C = (R *)&c; R m = ip(nv(j), c)/3;
    m0 += m; for(int i=0; i<3; ++i) {m1[i] += m * C[i];
      for(int z=0; z<3; ++z) m2[z][i] += m*C[z]*C[i];}}
    for(int j=0; j<3; ++j) {m1[j] *= 3./4;
      for(int k=0; k<3; ++k) m2[j][k] *= (3./5);}
  return m0;}