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

#define PI 3.14159265

int main()
{
  double ct_f, k, diameter;
  double f, x, y, z, sx, sy, thetan, phin, xn, yn, zn;
  double x_min, y_min, x_max, y_max, a;
  double norm;
  double step; /* (cm) separation of mirror centers = length of mirror side */

  long x_index, y_index;
  long i_mirror;
  char dummy[256];

  step = 50.;

  ct_f = 16.97;  // In meters-> camera will then be at 17 m to focus at 10 km

  // Convert to cm:
  ct_f *= 100;

  k = 1./(4.*ct_f);


  diameter = 17;  // meters

  x_max = y_max = 8.25; // floor(diameter/2+0.5)-0.25;
  x_min = y_min = -x_max;


  // Convert to cm:
  diameter *= 100;
  x_max *= 100;
  y_max *= 100;
  x_min *= 100;
  y_min *= 100;

  i_mirror = 0;

  for (x = x_min ; x < x_max+1; x += step)  
    for (y = y_min ; y < y_max+1; y += step)    
      {

	if ( (fabs(x)+fabs(y)) > diameter*0.72)
	  continue;

	if ((x >-150 && x < 50) && fabs(y) < 50)
	  continue;

	if ( (x == 825. || x == -825.) && (y == 75. || y == -75.) )
	  continue;

	x_index = x > 0? (int)(x+50)/100 : (int)(x-50)/100 ;
	y_index = y > 0? (int)(y+50)/100 : (int)(y-50)/100 ;

//	printf("%d %d  %.0f %.0f\n", x_index, y_index, x, y);

	i_mirror++;

	z = k*(x*x+y*y);

	// Curvilinear coordinates:

	sx = (2*k*x*sqrt(1+4*k*k*x*x)+asinh(2*k*x))/4/k;
	sy = (2*k*y*sqrt(1+4*k*k*y*y)+asinh(2*k*y))/4/k;

	phin = atan2(y,x);
	if (phin < 0)
	  phin += 2*PI;

//
// OLD: (before chessboarding) Valid only if mirror is ON the parabola
//	xn = - x / sqrt(x*x+y*y+4*ct_f*ct_f);
//	yn = - y / sqrt(x*x+y*y+4*ct_f*ct_f);
//	zn = 2*ct_f / sqrt(x*x+y*y+4*ct_f*ct_f);
//	printf("%.6f %.6f %.6f\n", xn, yn, zn);

//
//	Shift z for chessboarding:
//
	if (!((x_index+y_index)%2 != 0  || 
	      (x_index+y_index) ==  12  || 
	      (x_index+y_index) == -12  || 
	      (x_index-y_index) ==  12  || 
	      (x_index-y_index) == -12 ) )
	  z += 8.;

//
// NEW: Valid also (correct focusing) if we shift mirrors from the parabola,
//      since we use explicitely the z coordinate:
//
	xn = - x / sqrt(x*x+y*y+(ct_f-z)*(ct_f-z));
	yn = - y / sqrt(x*x+y*y+(ct_f-z)*(ct_f-z));
	zn = 1 + (ct_f - z) / sqrt(x*x+y*y+(ct_f-z)*(ct_f-z));

	norm = sqrt(xn*xn+yn*yn+zn*zn);
	xn /= norm;
	yn /= norm;
	zn /= norm;

//	printf("%.6f %.6f %.6f\n\n", xn, yn, zn);


	thetan = acos(zn);

	a = 2*ct_f;

//
// Focal distance of the mirror:
//
	f= 0.5 * (a*sqrt((1.+(x*x+y*y)/(a*a))*
			 (1.+(x*x+y*y)/(a*a))*
			 (1.+(x*x+y*y)/(a*a))) + sqrt(x*x+y*y)/
		  (sin(atan(sqrt(x*x+y*y)/a))))/2.;


	printf("%3d %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.8f %9.8f %9.8f %9.8f %9.8f\n", i_mirror, f, sx, sy, x, y, z, thetan, phin, xn, yn, zn);
      }

  return 0;
}


