// gcc conic.c p.c -lm -std=c99 -Wall; ./a.out
#include "h.h"
static R pi = 3.1415926535897932384626434L;
d r(d x, d y){return sqrt(x*x + y*y);}
int main(){assert(N==2 && D == 2);
 if(0) {pt z[] = {{{{1,0},{0,1}}, 1.e-20},{{{0,0},{0,0}},1}}; // Circle
   for(R t=0; t < 10.; t += 0.7) {
     printf("t, dx, dy = %e, %e, %e\n", (d)t,
     (d)(z[0].dv[0].x - cosl(t)), (d)(z[0].dv[1].x - sinl(t)));
      spin(z, 0.7, 1, 1000);}}
 if(0) {pt z[] = {{{{1,0},{0, sqrtl(2)}},1.e-20},{{{0,0},{0,0}},1}}; // Parabola
   int j; for(j=0; j < 16; ++j) {
    d x = z[0].dv[0].x,   y = z[0].dv[1].x,
     xd = z[0].dv[0].xd, yd = z[0].dv[1].xd;
     d pe = -1/r(x, y), ke = s(r(xd, yd))/2;
     printf("x, y = %13.10f, %13.10f; xd,yd = %15.12f, %15.12f;\n 1-x-yy/4= %e; "
       "pe, ke = %13.10f, %13.10f\n", x, y, xd, yd, (d)(1-x-y*y/4), pe, ke);
       spin(z, 1., 1, 1000);}}
 if(1) {R a = 2, e = 0.5, v = a/(1-e*e), ts = v*sqrt(v), af = -2*e*v, ax = 1/sqrt(a);
   // same angular momentum as circle.
   pt z[] = {{{{a/(1+e),0},{0, ax*(1+e)}},1.e-20},{{{0,0},{0,0}},1}}; // Ellipse
   R g(R x){return acosl(x) - sqrtl(1-x*x)*e;}
   assert(fabs(pi - g(-1)) < 1.e-13);
   if(1) {R t[] = {g(1), g(e), g(0), g(-e), g(-1), 2*pi-g(-e), 2*pi};
   int n = sizeof(t)/sizeof(R);
   for(int j=0; j < n; ++j) {
    d x = z[0].dv[0].x,   y = z[0].dv[1].x,
     xd = z[0].dv[0].xd, yd = z[0].dv[1].xd;
     d pe = -1/r(x, y), ke = s(r(xd, yd))/2;
     d td = sqrt(x*x+y*y) + sqrt(s(x-af)+y*y);
     printf("x,y = %13.10f, %13.10f; xd,yd = %13.12f, %13.12f;\n"
       "pe+ke = %13.10f, td = %13.10f\n", x, y, xd, yd, pe+ke, td);
     printf("atan=%14.11f, tt=%16.13f\n", atan2(z[0].dv[1].x, z[0].dv[0].x), (d)t[j]);
     if(j<n-1) spin(z, ts*(t[j+1] - t[j]), 1, 1000);}}}
return 0;}