// clang a.c -O3
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
int const m=20; int const n=20;
typedef double R;
typedef struct{R x; R y; R xd; R yd;} pt;
pt a[m][n];
R t = 0;
void prn() {printf("t = %20.16f\n", t);
    for(int i=0; i<m; ++i) for(int j=0; j<n; ++j) {pt p = a[i][j];
      printf("%2d %2d %20.16f %20.16f %20.16f %20.16f\n",
      i, j, p.x, p.y, p.xd, p.yd);}};
int main(){
{int i=m; while(i--){int j=n; while(j--) 
   a[i][j] = (pt){i + 0.5*j, (sqrt(3)/2)*j, 0, 0};}}
R dt = 0.005, tp = -dt/2, Xm = a[m-1][n-1].x, Ym = a[m-1][n-1].y;
    while(t < 10.55) {if(t > tp) {tp += 1; prn();}
// Symmetries??
if(1) {R c6 = sqrt(3)/2; R eps=0000001;
if(m!=n) exit(printf("Sym!!\n"));
for(int i=0; i<m; ++i) for(int j=0; j<m; ++j) {
  pt p=a[i][j], q=a[j][i], r=a[m-1-i][n-1-j];
  // Verify symmetry about long diagonal:
  if(fabs(p.x - ( 0.5*q.x + c6*q.y))> eps) exit(printf("Fox %d %d\n", i, j));
  if(fabs(p.y - (-0.5*q.y + c6*q.x))> eps) exit(printf("Foy %d %d\n", i, j));
  // Verify central symmetry.
  if(fabs(p.x + r.x - Xm) > eps) exit(printf("Fux %d %d\n", i, j));
  if(fabs(p.y + r.y - Ym) > eps) exit(printf("Fuy %d %d\n", i, j));
  }}
   R dx[n-1], dy[n-1];
   for(int j=0; j<n-1; ++j) {
      dx[j] = a[0][j+1].x - a[0][j].x;
      dy[j] = a[0][j+1].y - a[0][j].y;}
   for(int i=0; i<m-1; ++i) {R
      dex = a[i+1][0].x - a[i][0].x,
      dey = a[i+1][0].y - a[i][0].y;
   for(int j=0; j<n-1; ++j) {R
      p0 = dt*(1/(dy[j]*dex - dx[j]*dey)-2),
      ddx = a[i][j+1].x - a[i+1][j].x,
      ddy = a[i][j+1].y - a[i+1][j].y;
      a[i][j].xd -= p0*ddy; a[i][j].yd += p0*ddx;
      a[i][j+1].xd -= p0*dey; a[i][j+1].yd += p0*dex;
      a[i+1][j].xd += p0*dy[j]; a[i+1][j].yd -= p0*dx[j];
      dx[j] = a[i+1][j+1].x - a[i+1][j].x;
      dy[j] = a[i+1][j+1].y - a[i+1][j].y;
      dex = a[i+1][j+1].x - a[i][j+1].x;
      dey = a[i+1][j+1].y - a[i][j+1].y;
    R p1 = dt*(1/(dy[j]*dex - dx[j]*dey)-2);
      a[i][j+1].xd -= p1*dy[j]; a[i][j+1].yd += p1*dx[j];
      a[i+1][j].xd += p1*dey;  a[i+1][j].yd -= p1*dex;
      a[i+1][j+1].xd += p1*ddy; a[i+1][j+1].yd -= p1*ddx;
      a[i][j].x += dt*a[i][j].xd; a[i][j].y += dt*a[i][j].yd;}
    a[i][n-1].x += dt*a[i][n-1].xd; a[i][n-1].y += dt*a[i][n-1].yd;}
    for(int j=0; j<n; ++j) {
      a[m-1][j].x += dt*a[m-1][j].xd; a[m-1][j].y += dt*a[m-1][j].yd;}
t += dt;}
prn(); return 0;}