#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)    
      {
	// Skip non existent mirrors:

	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 ;

	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;

	//
	// Shift z if chessboarding is desired:
	//
	//	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.; // cm

	//
	// Calculate orientation of vector normal to the mirror:
	//
	// Add the versor pointing from the mirror element center (x,y,z)
	// to the camera center (0,0,ct_f), and versor (0,0,1). The 
        // direction of the resulting vector is that of the bisector of 
	// the two versors. If the normal of the mirror element is 
	// oriented along that direction, light rays paralel to the
	// telescope axis and impinging on the center of the mirror
	// element will be reflected to the center of the camera.
	// 
	//
	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));

	//
	// Normalize the vector to obtain the versor:
	//
	norm = sqrt(xn*xn+yn*yn+zn*zn);
	xn /= norm;
	yn /= norm;
	zn /= norm;

	//
	// Get angle between the normal to the mirror and the 
	// telescope axis:
	//
	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;
}