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

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