// clang sph.c test.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;}
R const pi = 3.141592653589793;
R rnd(){static int f = 0; R static const scale = 0x1.p-30;
  static R m;
  if(f) {f = 0; return m;}
  while (1){R x = scale*random()-1, y = scale*random()-1, s = x*x +y*y; 
    if(s<1) {R r = sqrt(- 2*log(s)/s); f = 1; m = x*r; return y*r;}}}
static void chk(char * m){printf("chk %s\n", m);
  {int x=0; for(int c = 0; c<3*SN; ++c) {assert(L[c] < 3*SN && L[c] > -2);
                  if(L[c]==-1) ++x;} if(x!=VN) exit(printf("Box\n"));}
  for(int c = 0; c<  VN; ++c) assert(h[c] < 3*SN && h[c] > -2);}
void P(Pt t){R v = 0, p = 0;
     for(int j = 0; j<SN; ++j) {T u = Q[j]; R r = len(sub(t, centr(j, q)));
       R d = det(sub(q[u.a], t), sub(q[u.b], t), sub(q[u.c], t));
       v += d; p += d/r;}
     {R r = len(t); printf("v=%14.10f, r=%10.6f, p=%14.10f, tp=%12.8f\n",
            v/6, r, p/4, (pi*4/3)*(r<1? (3-r*r)/2: 1/r));}}
int main(){sph(); printf("n = %d vertices\n", n); chk("a");
  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(ip(c, c) > E) {
             prP("D", c); exit(printf("damn %e 😠\n", len(c)));}}
  if(1) {int k = 102%VN; prP("vec", q[k]); prP("wind", puff(k));}
  if(0) for(int y=0; y<SN; ++y) printf("%d %d %d\n", Q[y].a, Q[y].b, Q[y].c);
  for(int j=0; j<30; ++j) P((Pt){0.7*rnd(), 0.7*rnd(), 0.7*rnd()});
  for(int j=0; j<21; ++j) P((Pt){0, j*0.1, 0});
  P((Pt){0,1,0});  P((Pt){0,0,0.999});  P((Pt){0.6,0.8001,0});
  P((Pt){0,0,0}); P((Pt){0,0.01,0});
  return 0;}

/*
qx = 4002, VN = 4002
Qx = 8000, SN = 8000
total mass =  4.18298294339, 4/3 π =  4.18879020479
n = 20 vertices
chk a
 0.740783169626  0.584370049471  0.331288304774 vert[k]
 0.001256676132  0.000939563067  0.000598086754 nvv
 0.001235663009  0.001000582014  0.000620212974 nvv
 0.001241486659  0.000906632556  0.000531060873 nvv
 0.001203641584  0.001025374769  0.000565309600 nvv
 0.001210558587  0.000931489370  0.000479945892 nvv
 0.001190117640  0.000989749374  0.000499313073 nvv
 0.022014430832  0.017380173452  0.009881787497 sum
 0.368117876835  0.699466179729  0.612565337063 va
 0.346860766779  0.664646410445  0.661764880870 vb
 0.404391307602  0.654318880455  0.639010542180 vc
 0.373123317072  0.672810490209  0.637780253371 sliv
 0.000650205294  0.001173392441  0.001111366800 nv
-0.000000000000  0.000000000000 -0.000000000000 vec sum area
 0.641055147461  0.613066296619  0.461733704492 vec
 0.007055220685  0.006748364413  0.005089226787 wind
v=  4.1829829434, r=  0.930575, p=  4.4638378474, tp=  4.46950214
v=  4.1829829434, r=  1.043076, p=  4.0102115840, tp=  4.01580619
v=  4.1829829434, r=  2.020515, p=  2.0702551706, tp=  2.07312956
v=  4.1829829434, r=  1.622347, p=  2.5783529991, tp=  2.58193287
v=  4.1829829434, r=  1.164962, p=  3.5907001241, tp=  3.59564685
v=  4.1829829434, r=  0.793423, p=  4.9592723378, tp=  4.96472302
v=  4.1829829434, r=  1.190370, p=  3.5140301235, tp=  3.51889741
v=  4.1829829434, r=  1.223792, p=  3.4180911212, tp=  3.42279623
v=  4.1829829434, r=  1.615251, p=  2.5896795909, tp=  2.59327552
v=  4.1829829434, r=  0.864158, p=  4.7136098909, tp=  4.71915503
v=  4.1829829434, r=  1.261335, p=  3.3162970015, tp=  3.32091897
v=  4.1829829434, r=  1.136241, p=  3.6814087052, tp=  3.68653196
v=  4.1829829434, r=  0.722344, p=  5.1850259217, tp=  5.19036955
v=  4.1829829434, r=  0.949516, p=  4.3891884494, tp=  4.39491863
v=  4.1829829434, r=  0.544618, p=  5.6568430187, tp=  5.66196906
v=  4.1829829434, r=  1.581574, p=  2.6448247086, tp=  2.64849378
v=  4.1829829434, r=  1.513427, p=  2.7639053472, tp=  2.76775124
v=  4.1829829434, r=  0.517463, p=  5.7172730958, tp=  5.72237298
v=  4.1829829434, r=  0.734096, p=  5.1491648810, tp=  5.15452314
v=  4.1829829434, r=  1.470827, p=  2.8439584489, tp=  2.84791432
v=  4.1829829434, r=  1.011105, p=  4.1370683084, tp=  4.14278599
v=  4.1829829434, r=  0.775085, p=  5.0195422667, tp=  5.02496219
v=  4.1829829434, r=  0.984362, p=  4.2481767504, tp=  4.25378413
v=  4.1829829434, r=  0.841111, p=  4.7959437261, tp=  4.80146801
v=  4.1829829434, r=  0.881032, p=  4.6518771184, tp=  4.65747890
v=  4.1829829434, r=  0.903342, p=  4.5684825260, tp=  4.57410115
v=  4.1829829434, r=  1.142119, p=  3.6625032852, tp=  3.66755954
v=  4.1829829434, r=  1.509684, p=  2.7707583189, tp=  2.77461474
v=  4.1829829434, r=  1.260913, p=  3.3174250830, tp=  3.32202954
v=  4.1829829434, r=  1.282586, p=  3.2613604792, tp=  3.26589439
v=  4.1829829434, r=  0.000000, p=  6.2783447092, tp=  6.28318531
v=  4.1829829434, r=  0.100000, p=  6.2573910674, tp=  6.26224136
v=  4.1829829434, r=  0.200000, p=  6.1945301461, tp=  6.19940950
v=  4.1829829434, r=  0.300000, p=  6.0897619778, tp=  6.09468975
v=  4.1829829434, r=  0.400000, p=  5.9430866528, tp=  5.94808209
v=  4.1829829434, r=  0.500000, p=  5.7545042508, tp=  5.75958653
v=  4.1829829434, r=  0.600000, p=  5.5240144782, tp=  5.52920307
v=  4.1829829434, r=  0.700000, p=  5.2516158643, tp=  5.25693171
v=  4.1829829434, r=  0.800000, p=  4.9373050226, tp=  4.94277244
v=  4.1829829434, r=  0.900000, p=  4.5810789196, tp=  4.58672527
v=  4.1829829434, r=  1.000000, p=  4.1829901067, tp=  4.18879020
v=  4.1829829434, r=  1.100000, p=  3.8027062537, tp=  3.80799110
v=  4.1829829434, r=  1.200000, p=  3.4858140442, tp=  3.49065850
v=  4.1829829434, r=  1.300000, p=  3.2176745934, tp=  3.22214631
v=  4.1829829434, r=  1.400000, p=  2.9878410427, tp=  2.99199300
v=  4.1829829434, r=  1.500000, p=  2.7886520813, tp=  2.79252680
v=  4.1829829434, r=  1.600000, p=  2.6143617549, tp=  2.61799388
v=  4.1829829434, r=  1.700000, p=  2.4605761422, tp=  2.46399424
v=  4.1829829434, r=  1.800000, p=  2.3238777753, tp=  2.32710567
v=  4.1829829434, r=  1.900000, p=  2.2015686659, tp=  2.20462642
v=  4.1829829434, r=  2.000000, p=  2.0914904289, tp=  2.09439510
v=  4.1829829434, r=  1.000000, p=  4.1829901067, tp=  4.18879020
v=  4.1829829434, r=  0.999000, p=  4.1871373213, tp=  4.19297690
v=  4.1829829434, r=  1.000080, p=  4.1826515723, tp=  4.18845512
v=  4.1829829434, r=  0.000000, p=  6.2783447092, tp=  6.28318531
v=  4.1829829434, r=  0.010000, p=  6.2781351728, tp=  6.28297587
*/