Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/Changelog
===================================================================
--- trunk/MagicSoft/Simulation/Detector/ReflectorII/Changelog	(revision 5419)
+++ trunk/MagicSoft/Simulation/Detector/ReflectorII/Changelog	(revision 5438)
@@ -1,3 +1,9 @@
 ** Add changes at the beginning! **
+
+21/11/2004 A. Moralejo
+
+Added subdirectory writemagicdef under ReflectorII/. It contains the 
+program to generate the positions and orientations of the mirrors. 
+Directory contains Makefile and writemagicdef.c
 
 03/11/2004 A. Moralejo
Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/writemagicdef/Makefile
===================================================================
--- trunk/MagicSoft/Simulation/Detector/ReflectorII/writemagicdef/Makefile	(revision 5438)
+++ trunk/MagicSoft/Simulation/Detector/ReflectorII/writemagicdef/Makefile	(revision 5438)
@@ -0,0 +1,4 @@
+writemagicdef: writemagicdef.o
+	cc -o writemagicdef writemagicdef.o -lm
+.c.o:
+	cc -c writemagicdef.c
Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/writemagicdef/writemagicdef.c
===================================================================
--- trunk/MagicSoft/Simulation/Detector/ReflectorII/writemagicdef/writemagicdef.c	(revision 5438)
+++ trunk/MagicSoft/Simulation/Detector/ReflectorII/writemagicdef/writemagicdef.c	(revision 5438)
@@ -0,0 +1,127 @@
+#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;
+}
+
+
