// clang sph.c blobC.c mcP.c gen.c -Wmost
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include "blob.h"
#include <stdlib.h>
static R sq(R x){return x*x;}
static R P(Pt A, Pt B){R a2 = ip(A, A), ai = 1/a2, x = ai*ip(A, B);
  Pt CP = cp(A, B); R y2 = ai*ai*ip(CP, CP), x2 = x*x;
  R f = sqrt(sq(x-1) + y2), g = sqrt(x2+y2);
  return 0.5*a2*(3*x*(f - g) + f + (2*x2 - y2)*log((f-x+1)/(g-x)));}
static void chk(char * m){printf("chk %s\n", m);
  if(0) for(int c = 0; c<3*SN; ++c) if(L[c]==-1) printf("%d\n", c);
  for(int c = 0; c<3*SN; ++c) assert(L[c] < 3*SN && L[c] > -2);
  for(int c = 0; c<  VN; ++c) assert(h[c] < 3*SN && h[c] > -2);}
void mv(Pt om, Pt tp){R p = 0, M=0; Pt f = {0,0,0}, qn[VN], C = {0,0,0};
    const Pt B = add(tp, om);
    for(int j = 0; j<VN; ++j) qn[j] = add(q[j], om);
    for(int j = 0; j<SN; ++j) if(j==300%SN) {
        R m = det(qn[Q[j].a], qn[Q[j].b], qn[Q[j].c])/6;
        Pt c = centr(j, qn); R dp = 3*m*P(c, B);
        {T t = Q[j]; printf("mP = %e, m*mcP = %e\n",
                            dp/3, mcP(q[t.a], q[t.b], q[t.c], B));}
        M += m; C = add(C, c);
        f = add(nv(j), f); p += dp;}
    printf("M=%12.9f, r = %12.9f, M/r = %12.9f, p=%12.9f\n", M, len(tp), M/len(tp), p);
    prP(" = force ", f); prP(" = center", sm(1./SN, C));}
int main(){sph(); printf("n = %d vertices\n", n); chk("a");
  if(0) { // several test cases
  if(1) printf("zz %20.17f\n", P((Pt){1, 0, 0}, (Pt){0.43, 0.61, 0})); // => 0.472826963
  if(1) {int k = 85%VN; Pt s = {0,0,0}; prP("vert[k]", q[k]); // k selects vertex
        for(int c = h[k]; c>=0; c = L[c]) {Pt w = nv(c/3); prP("nvv", w); s=add(s,w);}
        prP("sum", sm(3, s));}
  if(1) {int k = 318%SN; T t = Q[k]; // k selects sliver
           prP("va", q[t.a]); prP("vb", q[t.b]); prP("vc", q[t.c]);
           prP("sliv", centr(k, q)); prP("nv", nv(k));}
  if(1) {Pt c = {0, 0, 0}; for(int j=0; j<VN; ++j) c = add(c, puff(j));
           prP("vec sum area", c); if(len(c) > 1e-8) printf("damn 😠\n"); }
  if(1) {int k = 102%VN; prP("vec", q[k]); prP("wind", puff(k));}}
  else if(0) { // this test code is not needed under "easy sliver" plans, and is wrong.
        mv((Pt){0,  0,  0}, (Pt){5, 9, 3});
        mv((Pt){0.6, 0, 0}, (Pt){5, 9, 3});
        mv((Pt){0, 0.6, 0}, (Pt){5, 9, 3});
        mv((Pt){0, 0, 0.6}, (Pt){5, 9, 3});
        mv((Pt){0, 0,-0.6}, (Pt){5, 9, 3});}
  else {R min=100, max=0; int cnt = 0;
     for(int j=0; j<1000; ++j) {Pt r = {8.*rf()-4., 8.*rf()-4., 8.*rf()-4.};
     R rd = len(r); if(rd > 1.1) {R p = 0; 
      for(int j = 0; j<SN; ++j) p += 3*m[j]*P(centr(j, q), r);
      {R mr =  p*rd; if(mr<min) min=mr; if(mr>max) max=mr; ++cnt;
       if(0) printf("xs %13.10f %9.6f\n", mr, rd);}}}
    printf("%d %13.10f %13.10f\n", cnt, min, max);}
  return 0;}

/*
qx = 4002, VN = 4002
Qx = 8000, SN = 8000
tm =  4.18298294339, 4/3 π =  4.18879020479
n = 20 vertices
chk a
zz  0.47282696269790586
 0.7407832  0.5843700  0.3312883 vert[k]
 0.0012567  0.0009396  0.0005981 nvv
 0.0012357  0.0010006  0.0006202 nvv
 0.0012415  0.0009066  0.0005311 nvv
 0.0012036  0.0010254  0.0005653 nvv
 0.0012106  0.0009315  0.0004799 nvv
 0.0011901  0.0009897  0.0004993 nvv
 0.0220144  0.0173802  0.0098818 sum
 0.3681179  0.6994662  0.6125653 va
 0.3468608  0.6646464  0.6617649 vb
 0.4043913  0.6543189  0.6390105 vc
 0.3731233  0.6728105  0.6377803 sliv
 0.0006502  0.0011734  0.0011114 nv
 0.6410551  0.6130663  0.4617337 vec
-0.0037245  0.0098449  0.0318008 wind
*/