#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <strings.h>

#include "3d.h"
#include "trig.h"

#define SQRT_3_4 .8660254037844386467637231707529361834714026269051903140279034897
#define SQRT_1_2 .707106781186547524400844362104849039284835937688474036588339869

static PT3 apts[2][2][2];
static PT3 opts[2][2][2];
static XFORM xf_o2;
static PT3 xo2[2][2][2];
static XFORM xf_a;
static PT3 xa[2][2][2];

static void gen_pts(double a, PT3 pts[2][2][2])
{
 double sa;
 double ca;
 double r;
 double d;
 double cg;
 double sg;
 PT3 cc;
 PT3 u;

 sa = sindeg(a);
 ca = cosdeg(a);
 pts[0][0][0] = (PT3){0,0,0};
 pts[1][0][0] = (PT3){1,0,0};
 pts[0][1][0] = (PT3){ca,sa,0};
 pts[1][1][0] = add3(pts[1][0][0],pts[0][1][0]);
 cc = scale3(pts[1][1][0],.5);
 d = dist3(pts[0][0][0],cc);
 r = SQRT_3_4 * dist3(pts[1][0][0],pts[0][1][0]);
 // r² = 1² + d² - 2 × 1 × d × cos gamma
 // cos gamma = (1 + d² - r²) / (2 d)
 // we actually don't need gamma, just its sine
 cg = (1 + (d * d) - (r * r)) / (2 * d);
 if ((cg < -1) || (cg > 1)) abort();
 sg = sqrt(1-(cg*cg));
 u = unit3(pts[1][1][0]);
 pts[0][0][1] = (PT3){.x=cg*u.x,.y=cg*u.y,.z=sg};
 pts[1][0][1] = add3(pts[0][0][1],pts[1][0][0]);
 pts[0][1][1] = add3(pts[0][0][1],pts[0][1][0]);
 pts[1][1][1] = add3(pts[0][0][1],pts[1][1][0]);
}

static void gen(double a)
{
 gen_pts(a,&apts[0]);
 gen_pts(180-a,&opts[0]);
}

#if 0
static void dump_locs(PT3 p[2][2][2])
{
 int i;

 for (i=0;i<8;i++)
  { printf("%g %g %g\n",
	p[i&1][(i>>1)&1][(i>>2)&1].x,
	p[i&1][(i>>1)&1][(i>>2)&1].y,
	p[i&1][(i>>1)&1][(i>>2)&1].z);
  }
}
#endif

static double tryfit(void)
{
 int i;
 double d1;
 double d2;

 xf_o2 = find_transform(opts[1][1][1],opts[1][0][1],opts[1][1][0]);
 xf_a = find_transform(opts[1][0][1],opts[1][1][1],opts[0][0][1]);
 for (i=0;i<8;i++)
  { xo2[i&1][(i>>1)&1][(i>>2)&1] = transform(opts[i&1][(i>>1)&1][(i>>2)&1],xf_o2);
    xa[i&1][(i>>1)&1][(i>>2)&1] = transform(apts[i&1][(i>>1)&1][(i>>2)&1],xf_a);
  }
 d1 = dist3(xa[0][1][0],xa[1][0][1]);
 d2 = dist3(xa[0][1][0],xo2[0][0][1]);
 printf("# %g %g %g\n",d1,d2,d1-d2);
 return(d1-d2);
}

int main(void);
int main(void)
{
 double l;
 double m;
 double h;
 double d;

 l = 75;
 h = 61;
 while (1)
  { m = (l + h) / 2;
    printf("# %30.26f  %02x %02x %02x %02x %02x %02x %02x %02x\n",m,
	((unsigned char *)&m)[0],
	((unsigned char *)&m)[1],
	((unsigned char *)&m)[2],
	((unsigned char *)&m)[3],
	((unsigned char *)&m)[4],
	((unsigned char *)&m)[5],
	((unsigned char *)&m)[6],
	((unsigned char *)&m)[7]);
    if ((m == l) || (m == h)) break;
    gen(m);
    d = tryfit();
    if (d < 0) l = m; else h = m;
  }
 return(0);
}