/***************************************************************************
                          main.cpp  -  description
                             -------------------
    begin                : Sun Aug  1 01:31:51 UTC 2004
    copyright            : (C) 2004 by Martin Merck
    email                : merck@astro.uni-wuerzburg.de
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <iostream.h>
#include <iomanip.h>
#include <stdlib.h>
#include <math.h>

#define RADDEG	(180. / M_PI)
#define DEGRAD	(M_PI / 180.)

int main(int argc, char *argv[])
{

	for( int m_i = -8; m_i <=8; m_i++)
		for( int m_j = -8; m_j <=8; m_j++)
		{
    	// First we define some constants used in the calculation
    	const double kDistMirror = 17000.0;  		// Distance from camera to mirror center
    																					// Nominal focal length for infinity focusing
    	const double kDistRoque = 980000.0;     // Distance from the MAGIC telescope to the
    																					// Roque lamp.
    	const double kFocPlaneRoque = 300.1;    // change of focal plane for roque lamp focusing.
    	const double kDistWinstonCones = 35.0;	// Distance from the plane of focusing ( white
    																					// paper or styropor surface) to the PMT camera
    																					// plane (front of winston cone surface)

    	/*
       * First we calculate the distance of the center of this panel to the center of the
    	 * mirror.
    	 */
      double dMirrorRadius = sqrt( pow( m_i*1000.0, 2 ) + pow( m_j*1000.0, 2 ) );

    	/*
       * dHeightMirror is the height of the center of the selected mirror panel over the
       * tangential plane to the parabolas vertex.
       * Using the equation for a parabola open to the right:
       *                         y^2 = 4ax
    	 *    a is the distance from the vertex to the focal point), or focal length of the parabola.
       *
    	 * This height is calculated to:
    	 *                        x = y^2 / 4a;
       * or
    	 *                        dHeightMirror = SQR( dMirrorRadius ) / dLatusRectum;
       *
       * where we used the latus rectum of the parabooa: dLatusRectum = 4 * a.
       *
       * See: http://mathworld.wolfram.com/Parabola.html
       */
    	double dLatusRectum = 4.0 * kDistMirror;
    	double dHeightMirror = pow( dMirrorRadius, 2. ) / dLatusRectum;

      /*
       * Alpha is the angle betweeen the optical axis and the connection line
    	 * between the mirror center and the focal point of the parabola.
    	 * (The focal point of the parabola)
       *
    	 *
    	 *                                             ____-\         -
    	 *                                     ____----               ^
    	 *                             ____----                   dMirrorRadius
    	 *										 ____----      dAlpha                   v
    	 *                 x------------------------------------|     -
    	 *
       *                 <------------ kDistMirror ------------>
       *                                                   <--->   dHeightMirror
       */
       double dAlpha = atan( dMirrorRadius / ( kDistMirror - dHeightMirror ) );

    	/*
    	 * For the Roque Lamp focusing the ray from the roque lamp is seen at the panels
       * center under an angle dDelta.
    	 *                   dDelta = atan( dMirrorRadius / (kDistRoque - dHeightMirror) );
    	 *
       *                                                           ________-\         -
    	 *                                           ________--------     __--          ^
    	 *                           ________--------                 __--        dMirrorRadius
    	 *					 ________--------      dDelta                __--     dGamma        v    	
       *   x-------------------------------------------------x--*-------------|       -
    	 *
       *   <--------------------------- kDistRoque --------------------------->
       *                                                     <- dFocusRoque -->
       *                                                     		<-kDistMirror->
       *                                                                    <->   dHeightMirror
       *
       * The reflected ray intersects the optical axis at a focal distance dFocusRoque
       * which is longer then focal length of the parabola.
       * For a spherical mirror (Davis-Cotton design like CT1) this difference would be
       * 300.1 mm which is the distance we shift the camera plane back for the roque lamp
       * focusing.
       * The angle at which this reflected ray intersects the optical axis is:
    	 *      dGamma = dAlpha - dDelta;
       *
       * For the parabola this distance is only correct for the vertex and bigger as farer the
       * mirror is away from the vertex. We calculate now the focal length at roque focusing
       * for this mirror dFocusRoque and the difference dDiffFocusRQ between the roque focusing
       * plane amd this focal length.
       */
    	double dDelta = atan( dMirrorRadius / ( kDistRoque - dHeightMirror) );
    	double dGamma = dAlpha - dDelta;
      double dFocusRoque = (dMirrorRadius / tan( dGamma )) + dHeightMirror;
    	double dDiffFocusRQ = dFocusRoque - kFocPlaneRoque - kDistMirror;

      /*
       * The correction we have to apply results from the error we made by forcing the spots of
       * the individual mirrors onto the focal spot of the vertex for the roque lamp distance.
    	 *
       *                                                        ____-\              -
    	 *                                                ____---- ..                 ^
    	 *                                        ____----     ....                   |
    	 *					                      ____----________________________            |
    	 *                        ____----    |        ....              ^       dMirrorRadius
    	 *                ____----            |    ....              dCorrection      |
    	 *			  ____----        	          |....                      V            V
       *   x--------------------------------*---------------------------------|     -
       *
       *   <--------------------------- dFocusRoque -------------------------->
       *                                    <- kDistMirror + kFocPlaneRoque -->
       *   <---  dDiffFocusRQ ------------->                         <--------> dHeightMirror
       *   <-------------- dFocusRoque - dHeightMirror ------------->
       *
       * Using the rule of three we get
       *
    	 *         dCorrection / dDiffFocusRQ = dMirrorRadius / ( dFocusRoque - dHeightMirror )
       * or:
       *         dCorrection = (dMirrorRadius / (dFocusRoque- dHeightMirror) ) * dDiffFocusRQ
       */
    	double dCorrection1 = (dMirrorRadius / (dFocusRoque - dHeightMirror) ) * dDiffFocusRQ;

      /*
       * The second correction needs to be applied because we can not use the front of the
       * winston cones as reference plane but a plane defined py the styropor panel put on
       * top of the plexiglas window of the camera.
       * This reference plane is 35 mm (kDistWinstonCones) infornt of the real focal plane we
       * should have used.
       * Looking at the scetch below we forced the ray to follow the line with stars, but
       * should have focus the mirror in such a way that it would have followed the dashed
       * line.
    	 *
       *                                                        ____--\              -
    	 *                                                ____---- ***                 ^
    	 *                                        ____----     ****                   |
    	 *					                      ____----         ****         ---   |
    	 *                        ____----    |        ****              ^       dMirrorRadius
    	 *                ____----            |    ****              dCorrection      |
    	 *			  ____----        	          |****                      V            V
       *   x--------------------------------*---------------------------------|     -
       *
       *                                    <- kDistMirror + kFocPlaneRoque -->
       *   <----- kDistWinstonCones ------->                             <----> dHeightMirror
       *   																	<- dFocus - dHeightMirror -->
       *
       * Using the rule of three we get
       *
    	 *         dCorrection / kDistWinstonCones = dMirrorRadius / ( dFocus - dHeightMirror - kDistWinstonCones )
       * or:
       *         dCorrection = (dMirrorRadius / (dFocus - dHeightMirror - kDistWinstonCones) ) * kDistWinstonCones
       */
      double dFocus = kDistMirror + kFocPlaneRoque;
      double dCorrection2 = (dMirrorRadius / (dFocus - dHeightMirror - kDistWinstonCones) ) * kDistWinstonCones;


      cout << "I: " << setw(2) << m_i;
      cout << "  J: " << setw(2) << m_j;
      cout << "  R: " << setiosflags(ios::showpoint | ios::right) << setw(6) << setprecision(5) << dMirrorRadius;
      cout << "  F_RQ: " << setiosflags(ios::showpoint | ios::right) << setw(6) << setprecision(5) << dFocusRoque;
      cout << "  C1: " << setiosflags(ios::showpoint | ios::right) << setw(4) << setprecision(3) << dCorrection1;
      cout << "  C2: " << setiosflags(ios::showpoint | ios::right) << setw(4) << setprecision(3) << dCorrection2 << endl;
			double dCorrection = ( dCorrection1 + dCorrection2 ) / 2.0;

			/*
			 * The angle zeta is the angle of the mirror panel relative to the mirror center.
			 * 0 degrees corresponds to panels (0,x) with x being a positve number.
			 * The other angles are definedclockwise (positive) if looking onto the
			 * mirror plane from the camera.
			 */
			double dZeta = atan2( m_i * 1000.0, m_j * 1000.0 );

			/*	
			 * We assumethat the video camera is horizontally alligned with mirror dish.
			 * If we have a rotation of the video image in respect to the mirror frame we
			 * must add an extra correction here.
			 * As the coordinate system of the video image has a reversed Y axis ( 0rign being
			 * the upper left corner) we have to multiply the Y correction with -1. to get the right
			 * sign.
			 */
			double dCorrSpotX = sin( dZeta ) * dCorrection;
			double dCorrSpotY = -1.0 * cos( dZeta ) * dCorrection;
      cout << "  Alpha: " << setiosflags(ios::showpoint | ios::right) << setw(6) << setprecision(4) << dZeta * RADDEG;
      cout << "  X: " << setiosflags(ios::showpoint | ios::right) << setw(6) << setprecision(4) << dCorrSpotX;
      cout << "  Y: " << setiosflags(ios::showpoint | ios::right) << setw(6) << setprecision(4) << dCorrSpotY << endl;
		}

  return EXIT_SUCCESS;
}
