source: trunk/MagicSoft/Simulation/Detector/ReflectorII/writemagicdef/writemagicdef.c@ 5438

#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;
}