// clang sph.c fac.c gen.c -Wmost -O2
// (define π (* 4 (atan 1)))
#include <stdio.h>
#include <math.h>
#include "blob.h"
// #include <stdlib.h>
static R pt[VN]; // remember potential at vertices.
static R P(Pt w){R tp = 0;
  for(int j = 0; j<SN; ++j) {Pt lv = sub(centr(j, q), w);
     tp += ip(lv, nv(j))/len(lv);}
  return (1.5/3)*tp;}
R * mv(Pt * x){return &x->x;} // ugly!
static void tst(Pt w){prP("at:", w); printf("pot is %10.7f\n", P(w));}
int main(){sph();
  if(0) {tst((Pt){0, 0, 0}); tst((Pt){0, 0.5, 0});
    tst((Pt){0, sqrt(0.5), 0}); tst((Pt){1, 0, 0});
    tst((Pt){0, 0, 1}); tst((Pt){0, 0, 2});
    tst((Pt){0, 4, 3}); tst((Pt){5, 0, 0});}
  if(1) {R tu = 0, td = 1000, ts = 0; for(int j=0; j<VN; ++j) {
      R t = P(q[j]); pt[j] = t; tu = tu<t?t:tu; td = td<t?td:t; ts += t;}
    printf("minMax = %10.7f %10.7f, average = %10.7f\n", td, tu, ts/VN);}
  {// Moments:
    R m0 = 0; Pt m1 = {0, 0, 0}, m2[3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
    for(int j=0; j<SN; ++j) {Pt c = centr(j, q); R m = ip(nv(j), c);
      Pt r = sm(m, c); m1 = add(m1, r);
      for(int z=0; z<3; ++z) m2[z] = add(m2[z], sm(mv(&c)[z], r)); m0 += m;}
    m0 *= 1./3; m1 = sm(3./4, m1);
    for(int j=0; j<3; ++j) m2[j] = sm(0.4, m2[j]);
    printf("m0 = %13.10f\n", m0); prP("m1", m1);
    for(int j=0; j<3; ++j) prP("Inertia", m2[j]);}
  return 0;}