/*                                                                              *
 *   This file is part of the ESO UVES Pipeline                                 *
 *   Copyright (C) 2004,2005 European Southern Observatory                      *
 *                                                                              *
 *   This library 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.                                        *
 *                                                                              *
 *   This program is distributed in the hope that it will be useful,            *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of             *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
 *   GNU General Public License for more details.                               *
 *                                                                              *
 *   You should have received a copy of the GNU General Public License          *
 *   along with this program; if not, write to the Free Software                *
 *   Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA       *
 *                                                                              */
/* ===========================================================================
 * $Id: uves_physmod_utils.c,v 1.19 2010-09-24 09:32:07 amodigli Exp $
 * $Name: not supported by cvs2svn $
 * $Revision: 1.19 $
 *
 * ===========================================================================
 */
/*----------------------------------------------------------------------------*/
/**
 * @defgroup uves_physmod_utils   Substep: UVES physical model utilities
 */
/*----------------------------------------------------------------------------*/

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

/*+++ uves_physmod_utils.c +++++++++++++++++++++++++++++++++++++++++++++++++++
 
 .COPYRIGHT    (C) 1998  European Southern Observatory
 
 .IDENT        
 
 .KEYWORDS    uves physical model, spectroscopy, echelle,
 
 .USAGE        .
 
 .INPUT        .
        
 .OUTPUT    .
 
 .RETURN    Q1:  0: successful return
             -1: a successful return failed
 
 .PURPOSE    Library of UVES Physical Model Functions. 
 
 .ALGORITHM    UVES Physical Model.

 .ENVIRON    MIDAS, UVES context
    
 .LANGUAGE    C
 
 .AUTHOR    Pascal Ballester, Olivier BOITQUIN, ESO-DMD
 
 .VERSION    1.0 1999/07/01      Creation
                    1999/09/30-SW Binning implemented.
            1999/10/27-OB Bugs with n index fixed.

 .COMMENT       
----------------------------------------------------------------------------*/


/*-----------------------------------------------------------------------------
                            Includes
 -----------------------------------------------------------------------------*/
#include <uves_physmod_utils.h>

#include <uves_utils.h>
#include <uves_utils_wrappers.h>
#include <uves_error.h>
#include <uves_msg.h>

#include <cpl.h>

#include <stdio.h>
#include <math.h>

/*-----------------------------------------------------------------------------
                            Functions prototypes
 -----------------------------------------------------------------------------*/
static void 
beta2lamb(double uves_beta_ech, double uves_beta_cd, double* plambda, int m);

static double 
cameraFocal(double lm);

static void 
uves_physmod_find_alpha_beta(double lm, 
                             int m, 
                             double k, 
                             double theta, 
                             double *alpha, 
                             double *beta);

static double 
uves_physmod_find_lambda(double k, double alpha, double beta);

static double 
uves_physmod_find_order_lambda(double k, double alpha, double beta);
/* not used:
static void 
uves_physmod_lambda_model(double* plambda, 
                          int m, 
                          double fc, 
                          double x, 
                          double y);
*/
static void 
uves_physmod_lambda_order_beta(double* plambda, 
                               int* pm, 
                               double x, 
                               double y, 
                               double* puves_beta_ech, 
                               double* puves_beta_cd);

static void 
uves_physmod_lambda_order_focus_model(double* plambda, 
                                      double* pdm, 
                                      double fc, 
                                      double x, 
                                      double y);

static double uves_physmod_wave_bin(double l, int m);

static void uves_physmod_xy2beta(double* puves_beta_ech, 
                                 double* puves_beta_cd, 
                                 double fc, 
                                 double x, 
                                 double y);

static void uves_physmod_set_binning(float binx, float biny);

static double dsqrarg;

/*-----------------------------------------------------------------------------
                                Defines
 -----------------------------------------------------------------------------*/

#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)

#define     PROID "physmod.c" 




enum uves_arm_ident {UVES_ARM_UNDEF,UVES_ARM_BLUE,UVES_ARM_RED};
enum uves_arm_ident uves_arm_ident = UVES_ARM_UNDEF;


/* for messout function: no "static" (clear) -
   no "extern" (otherwise redeclaration in all files neccessary) */

/* globals declared as extern in physmod.h */

float    uves_bin[2]      = {1, 1}; /* binning of exposure in x and y */
static double   uves_physmod_pix_size[2]  = {15e-6, 15e-6}; /* meters */
static double   delta[2]  = {75.04, 76.0};

static double   cdgroov[4]  = {1000.0e-6, 660.0e-6, 600.0e-6, 312.5e-6};  /* grooves/nm */
static double   uves_ech_groov[2] = {31.6e-6,  41.59e-6};  /* grooves/nm */
/*double   uves_ech_blaze[2]  = {75.04, 75.9};  degrees */
/*double   uves_ech_blaze[2]  = {74.74,  76.1};*/
static double   uves_ech_blaze[2]  = {74.57, 75.9}; /* degrees */
/*double   cdbeam[2]   = {45.336,  45.0};  degrees */
/*double   cdbeam[2]   = {45.336,  45.9};*/
/*double   cdbeam[2]   = {45.336,  45.0};  degrees */
static double   cdbeam[2]   = {45.336,  46.0}; /* degrees */
/*double   uves_ccd_rot[2]   = {0.7, -0.45};  degrees */
/*double   uves_ccd_rot[2]   = {-1.0,  -1.50}; */
/*double   uves_ccd_rot[2]   = {0.2, -0.45};   degrees */
double   uves_ccd_rot[2]   = {0.3, -0.55};   /* degrees old */
/* blue detector upgrade 2004-10-13 */
/* double   uves_ccd_rot[2]   = {0.3, -0.10}; degrees  new */
static int      imsize[2]   = {4096, 3000}; /* pixels */
static int      uves_physmod_row_size[2]   = {2048, 2048}; /* pixels */

/* Configurations 
   1:Blue CD1,     2:Blue CD2, 
   3:Red CD 3 EEV, 4:Red CD 4 EEV, 
   5:Red CD3 MIT,  6: Red CD4 MIT 
*/

/* old */
double uves_physmod_offsetx[6]={1391.0,1385.0,1852.0,1835.0,1852.0,1834.0};
double uves_physmod_offsety[6]={1030.0,1025.0,2098.0,2104.0,-52.0,-49.0}; 

/* blue detector upgrade 2004-10-13 */
/* double uves_physmod_offsetx[6]={1355.0,1350.0,1852.0,1835.0,1852.0,1834.0}; new */
/* double uves_physmod_offsety[6]={1030.0,1023.0,2098.0,2104.0,-52.0,-49.0}; new */

static double   flipx       = 1.; 
static double   flipy       = -1.;
double   uves_airt        = 25;  /* Celsius */
double   uves_airp        = 760; /* mmHg */
double   uves_airw        = 3;   /* mmHg */

int      uves_cfg_indx     = 0; 
int      uves_x_disp_id     = 0;
static int      uves_ech_id       = 0;
double   uves_alpha0_cd, uves_beta0_cd;
char     uves_arm_id       = 'x';

static double   uves_deg2rad     = M_PI/180.;

/* For history, please keep track HERE of the previous offset used: 
Garching 1: 
double uves_physmod_offsetx[6]={1470.0,1450.0,2130.0,2140.0,2130.0,2140.0}; 
double uves_physmod_offsetx[6]={1515.0,1514.0,2010.0,2000.0,2010.0,2000.0};
 
Comissioning 1: 
double uves_physmod_offsetx[6]={1474.0,1471.0,1960.0,1948.0,1961.0,1949.0}; 

Comissioning 2, Dec 2000:
double uves_physmod_offsetx[6]={1390.0,1386.0,1849.0,1840.0,1854.0,1840.0};
 
Paranal 2, Feb 2000: 
double uves_physmod_offsetx[6]={1390.0,1384.0,1851.0,1840.0,1851.0,1839.0}; 

Garching 1: 
double uves_physmod_offsety[6]={1030.0,1030.0,2020.0,2020.0,-125.0,-115.0}; 
double uves_physmod_offsety[6]={1029.0,1025.0,2072.0,2080.0,-74.0,-65.0}; 

Comissioning 1: 
double uves_physmod_offsety[6]={1027.0,1024.0,2069.0,2077.0,-74.0,-65.0};

Comissioning 2, Dec 2000: 
double uves_physmod_offsety[6]={1027.0,1025.0,2084.0,2094.0,-65.0,-54.0}; 

Paranal 2, Feb 2000:
double uves_physmod_offsety[6]={1030.0,1025.0,2088.0,2094.0,-57.0,-54.0}; 

The gap between the EEV and MIT chips is estimated to 95 pixels */
/*
Mean pixel-scales for the 6 configurations: 
   pixscale : 0.252396, 0.246, 0.182, 0.175266, 0.182, 0.175266 */
 
/*---------------------------------------------------------------------*/


/**@{*/

/**
@brief set blaze incidence angles
@param echred  echelle offset angle (red arm)
@param echblue echelle offset angle (blue arm)
@param xred    Cross Disperser offset (red arm)
@param xblue   Cross Disperser offset (blue arm)


*/
void 
uves_physmod_set_incidence(double echred, 
                           double echblue, 
                           double xred, 
                           double xblue)
{
     uves_ech_blaze[0] += echred;
     uves_ech_blaze[1] += echblue;
     cdbeam[0]  += xred;
     cdbeam[1]  += xblue;
     uves_msg_debug("uves_ech_blaze=%f,%f ccdbeam=%f,%f",
                   uves_ech_blaze[0],uves_ech_blaze[1],cdbeam[0],cdbeam[1]);
}

/**
@brief set CCD rotation angles
@param ccdrot   CCD rotation angle
@param uves_ccd_rot_off_red  CCD rotation angle offset (red arm)
@param uves_ccd_rot_off_blue CCD rotation angle offset (blue arm)
 */
void 
uves_set_ccd_rot(double* ccdrot, 
                 double uves_ccd_rot_off_red, 
                 double uves_ccd_rot_off_blue)
{
  uves_ccd_rot[0] =ccdrot[0];
  uves_ccd_rot[1] =ccdrot[1];

  uves_ccd_rot[0] += uves_ccd_rot_off_red;
  uves_ccd_rot[1] += uves_ccd_rot_off_blue;

  uves_msg_debug("uves_ccd_rot[0,1]=%f,%f uves_ccd_rot_off: red,blue=%f,%f",
  uves_ccd_rot[0],uves_ccd_rot[1],uves_ccd_rot_off_red,uves_ccd_rot_off_blue);

}

/**
@brief set x,y offset to physical model
@param offset_x x offset (hardcoded, instrument setting dependent)
@param offset_y y offset (hardcoded, instrument setting dependent)
@param uves_physmod_x_off additional x offset
@param yoff               additional y offset


 */
void 
uves_physmod_set_offset(double offset_x, 
                        double offset_y,
                        double uves_physmod_x_off, 
                        double yoff)
{
  uves_physmod_offsetx[uves_cfg_indx-1]=offset_x;
  uves_physmod_offsety[uves_cfg_indx-1]=offset_y;
/*
  uves_msg("offset_x=%f offset_y=%f",offset_x,offset_y);
     uves_msg("uves_physmod_offsetx=%f 
               uves_physmod_offsety=%f 
               uves_physmod_x_off=%f 
               yoff=%f",
           uves_physmod_offsetx[uves_cfg_indx-1],
               uves_physmod_offsety[uves_cfg_indx-1],
               uves_physmod_x_off,
               yoff);
*/
     uves_physmod_offsetx[uves_cfg_indx-1] += uves_physmod_x_off;
     uves_physmod_offsety[uves_cfg_indx-1] += yoff;
     /*
     uves_msg("uves_physmod_offsetx=%f 
               uves_physmod_offsety=%f 
               uves_physmod_x_off=%f 
               yoff=%f",
           uves_physmod_offsetx[uves_cfg_indx-1],
               uves_physmod_offsety[uves_cfg_indx-1],
               uves_physmod_x_off,
               yoff);
     */
}


/**
@brief set binning
@param x binning
@param y binning

 */
void uves_physmod_set_binning(float binx, float biny)
{
  uves_bin[0] = binx;
  uves_bin[1] = biny;
}


/**
@brief set parameters describing atmospheric conditions
@param p pressure
@param t temperature
@param w humidity
@note Initialises atmospheric conditions with  p in mmHg, t in Celsius, w in 
      mmHg. Typical values, in Garching: t=25, p=760, w=3; 
      in Paranal: t=15, p=560, w=3

 */
void uves_air_config(double p, double t, double w)
{
    uves_airt = t;
    uves_airp = p;
    uves_airw = w;
    uves_msg_debug("uves_airt=%f uves_airp=%f uves_airw=%f",
                    uves_airt,uves_airp,uves_airw);
}

/**
@brief computes atmospheric index of refraction
@param lambda wavelength
 */
double uves_air_index(double lambda)
{
   /* wavelength is expected in nanometers. */
   double t1, t2, t3, airdx;
 
   t1 = 1.0e-6/lambda/lambda;  /* 1e-6 for nm to um conversion, squared */

   t2 = (64.328+29498.1/(146.0-t1)+255.4/(41.0-t1))*uves_airp*
               (1.0+1.0e-6*(1.049-0.0157*uves_airt)*uves_airp)/
               (720.883*(1.0+0.003661*uves_airt));

   t3 = t2 - uves_airw*(0.0624 - 0.000680*t1)/(1+0.003661*uves_airt);

   airdx = 1.0+ t3*1.0e-6;
   /*   airdx = 1.0; */
   /*
   uves_msg_debug("uves_airt=%f uves_airp=%f uves_airw=%f",
                   uves_airt,uves_airp,uves_airw);
   uves_msg_debug("lambda=%f t1=%g t2=%g t3=%g airdx=%g",
                  lambda,t1,t2,t3,airdx);
   */
   return(airdx);
}


/**
@brief find alpha and beta angles in echelle dispersion relation
@param lm    wavelength
@param m     order
@param k     kappa
@param theta theta angle
@param alpha alpha angle
@param beta  beta angle
 */
void 
uves_physmod_find_alpha_beta(double lm, 
                             int m, 
                             double k, 
                             double theta, 
                             double *alpha, 
                             double *beta)
{

 /* Solves sin(alpha) + sin(beta) = m*k*lm, given beta - alpha = theta */
 uves_msg_debug("lm, m, k, theta : %f %d %f %f",lm,m,k,theta);

 lm /= uves_air_index(lm);

 *alpha = 0.5* ( 2*asin( m*k*lm/2/cos(theta*uves_deg2rad/2) ) + 
                 theta*uves_deg2rad );
 *beta  = 0.5* ( 2*asin( m*k*lm/2/cos(theta*uves_deg2rad/2) ) - 
                 theta*uves_deg2rad );
 *alpha /= uves_deg2rad;
 *beta  /= uves_deg2rad;
}

/**
@brief set UVES configuration
@param uves_arm     UVES arm setting
@param uves_ccd_id  UVES CCD ID
@param disp         dispersion
@param waveCent     central wavelength
@param binx         x binning
@param biny         y binning

ArmId = 'r'   for the red uves_arm, ArmId = 'b' for the blue uves_arm of Uves 
uves_ccd_id = 'e' for the EEV detector, uves_ccd_id = 'm' for the MIT detector
uves_x_disp_id = 1 to 4 for cross-dispersers CD#1 to CD#4 

*/

int  
uves_config(char uves_arm, 
            char uves_ccd_id, 
            int disp, 
            double waveCent, 
            float binx, 
            float biny)
{

  int cfg = 0;
  uves_ech_id   = 2;
  uves_x_disp_id = disp;
  uves_arm_id   = uves_arm;
  
/*  uves_msg("Configuring: 
              Arm %c 
              CCD %c  
              Xdisp %d 
              Wave %f",
              uves_arm,
              uves_ccd_id,
              disp,
              waveCent); 
*/
  uves_msg("Cfg: Arm %c CCD %c  Xdisp %d Wave %f",
                  uves_arm,uves_ccd_id,disp,waveCent);

  if (uves_arm == 'b' && disp == 1)  cfg = 1;
  if (uves_arm == 'b' && disp == 2)  cfg = 2;

  if (uves_arm == 'r' && disp == 3) {
          uves_ech_id  = 1;
          if (uves_ccd_id == 'e') cfg = 3;
          if (uves_ccd_id == 'm') cfg = 5;
  }
  if (uves_arm == 'r' && disp == 4) {
          uves_ech_id  = 1;
          if (uves_ccd_id == 'e') cfg = 4;
          if (uves_ccd_id == 'm') cfg = 6;
  }

  if (cfg == 0) {
    cpl_msg_error(__func__,"Wrong configuration!");
    return -1;
  }

  uves_cfg_indx = cfg;

  (void) uves_physmod_set_binning(binx, biny);
  (void) uves_physmod_find_alpha_beta(waveCent, 1, cdgroov[uves_x_disp_id-1], 
                     cdbeam[uves_ech_id-1], &uves_alpha0_cd, &uves_beta0_cd); 


  uves_msg("alpha, beta for Xdisp: %f %f\nin config %d and CCD-ID %c",
                 uves_alpha0_cd, uves_beta0_cd,cfg, uves_ccd_id);
  
  return(cfg);

}
/**
@brief new way to set UVES instrument configuration
@param uves_arm     UVES arm setting
@param upper        UVES CCD ID switch
@param disp         dispersion
@param waveCent     central wavelength
@param binx         x binning
@param biny         y binning
*/
int  
uves_config_cpl_new(int uves_arm,
                    int upper, 
                    int disp, 
                    double waveCent, 
                    float binx, 
                    float biny)
{
  int cfg = 0;
  uves_ech_id   = 2;
  uves_x_disp_id = disp;
  
/*  uves_msg("Configuring: 
              Arm %c 
              CCD %c  
              Xdisp %d 
              Wave %f",
              uves_arm,
              uves_ccd_id,
              disp,
              waveCent); 
*/
  uves_msg("New Cfg: Arm [b/r] %d CCD eev/mit %d  Xdisp %d Wave %f",
                  uves_arm,upper,disp,waveCent);

  if (uves_arm == ARM_BLUE && disp == 1)  cfg = 1;
  if (uves_arm == ARM_BLUE && disp == 2)  cfg = 2;

  if (uves_arm == ARM_RED && disp == 3) {
          uves_ech_id  = 1;
          if (upper == 0) cfg = 3;
          if (upper == 1) cfg = 5;
  }
  if (uves_arm == ARM_RED && disp == 4) {
          uves_ech_id  = 1;
          if (upper == 0) cfg = 4;
          if (upper == 1) cfg = 6;
  }

  if (cfg == 0) {
    cpl_msg_error(__func__,"Wrong configuration!");
    return -1;
  }

  uves_cfg_indx = cfg;

  (void) uves_physmod_set_binning(binx, biny);
  (void) uves_physmod_find_alpha_beta(waveCent, 1, cdgroov[uves_x_disp_id-1], 
                    cdbeam[uves_ech_id-1], &uves_alpha0_cd, &uves_beta0_cd); 

 
  uves_msg("alpha, beta for Xdisp: %f %f\nin config %d and CCD-ID %c",
	   uves_alpha0_cd, uves_beta0_cd,cfg, upper == 0 ? 'e' : 'm');
 
  return(cfg);

}

/**
@brief set UVES instrument configuration: 
          # 1 346
          # 2 437
          # 3 520
          # 4 564
          # 5 580
          # 6 860

@param blue     blue/red arm switch
@param upper    upper/lower chip ID switch
@param disp     X dispersion
@param waveCent central wavelength
@param binx     x binning
@param biny     y binning


*/
int  
uves_config_cpl(int blue, 
                int upper, 
                int disp, 
                double waveCent, 
                float binx, 
                float biny)
{


  int cfg = 0;
  uves_ech_id   = 2;
  uves_x_disp_id = disp;
  
/*  uves_msg("Configuring: 
              Arm %c 
              CCD %c  
              Xdisp %d 
              Wave %f",
              uves_arm,
              uves_ccd_id,
              disp,
              waveCent); 
*/
  uves_msg("Cfg cpl: Arm [b/r] %d CCD eev/mit %d  Xdisp %d Wave %f",
                  blue,upper,disp,waveCent);

  if (blue == 1 && disp == 1)  cfg = 1;
  if (blue == 1 && disp == 2)  cfg = 2;

  if (blue == 0 && disp == 3) {
          uves_ech_id  = 1;
          if (upper == 0) cfg = 3;
          if (upper == 1) cfg = 5;
  }
  if (blue == 0 && disp == 4) {
          uves_ech_id  = 1;
          if (upper == 0) cfg = 4;
          if (upper == 1) cfg = 6;
  }

  if (cfg == 0) {
    uves_msg_error("Wrong configuration!");
    return -1;
  }

  uves_cfg_indx = cfg;

  (void) uves_physmod_set_binning(binx, biny);
  (void) uves_physmod_find_alpha_beta(waveCent, 1, cdgroov[uves_x_disp_id-1], 
                    cdbeam[uves_ech_id-1], &uves_alpha0_cd, &uves_beta0_cd); 

  
  uves_msg("alpha, beta for Xdisp: %f %f\nin config %d and CCD-ID %c",
		 uves_alpha0_cd, uves_beta0_cd,cfg, upper == 0 ? 'e': 'u');
 
  return(cfg);

}

/**
@brief computes the focal length at a given wavelength
@param lm wavelength
@return focal length

 Blue Arm  

 computes the focal length of the blue camera (in mm), lm is in nm.

 nfk=5.3961886e-7*lm*lm*lm-0.00079597882e0*lm*lm+0.41122805e0*lm+287.89644e0;

 Red Arm

 nfk=6.0172051e-13*lm*lm*lm*lm*lm-2.5623231e-9*lm*lm*lm*lm+
     4.3585543e-6*lm*lm*lm-0.0037286381e0*lm*lm+
     1.6289971e0*lm+210.06767e0;


 */
double cameraFocal(double lm)

{
   double nfk=0.;
   /*uves_msg("lm = %f ", lm);*/
   /* uves_msg_debug("lm=%g uves_air_index(lm)=%g",lm,uves_air_index(lm)); */

   lm /= uves_air_index(lm);
   /* uves_msg_debug("uves_arm=%d",uves_arm_ident); */
   if (uves_arm_id == 'b' ) { 
     /* uves_msg_debug("uves_arm blue"); */
    nfk=5.3961886e-7*lm*lm*lm-0.00079597882*lm*lm+0.41122805*lm+287.89644;
   }

   if (uves_arm_id == 'r' ) { 
     /* uves_msg_debug("uves_arm red"); */
    nfk=6.0172051e-13*lm*lm*lm*lm*lm-2.5623231e-9*lm*lm*lm*lm+
        4.3585543e-6*lm*lm*lm -0.0037286381*lm*lm+
        1.6289971*lm + 210.06767;

   }
   /*  uves_msg_debug("lm=%g nfk=%g",lm,nfk); */

 /*uves_msg("Nfk = %f, lm = %f",nfk/1000.,lm); */
return(nfk/1000.);

}


/**
@brief Solves sin(alpha) + sin(beta) = m*k*lm to find m 
@param lm wavelength
@return order
 */
int uves_physmod_find_order(double lm)
{
  int order;
  double k, alpha, beta;

  lm /= uves_air_index(lm);

  k     = uves_ech_groov[uves_ech_id-1];
  alpha = uves_ech_blaze[uves_ech_id-1];
  beta  = uves_ech_blaze[uves_ech_id-1];
  order = (int)((sin(alpha*uves_deg2rad) + sin(beta*uves_deg2rad))/k/lm + 0.5);
  /*
  uves_msg_debug("uves_ech_id=%d lm %g airindex %g k=%g alpha=%g beta=%g 
           order=%d",uves_ech_id,lm,uves_air_index(lm),k,alpha,beta,order);
  */
  return order;
}

/**
@brief Solves 
       sin(alpha*uves_deg2rad) + sin(beta*uves_deg2rad) = m*k*lm to find m*lm 
@param k     grating constant
@param alpha grating alpha
@param beta  grating beta
 */
double uves_physmod_find_order_lambda(double k, double alpha, double beta)
{
  double ML;
  ML =  ( ( sin(alpha*uves_deg2rad) + sin(beta*uves_deg2rad) )/k);
  ML *= uves_air_index(ML);
  return(ML);
}


/**
@brief Solves sin(alpha) + sin(beta) = m*k*lm to find m*lm 
@param k     grating constant
@param alpha grating alpha
@param beta  grating beta
 */
double uves_physmod_find_lambda(double k, double alpha, double beta)
{
  double L;
  L = (sin(alpha*uves_deg2rad) + sin(beta*uves_deg2rad) )/k ;
  L *= uves_air_index(L);
  return(L);
}

/**
@brief Solves sin(alpha) + sin(beta) = m*k*lm to find m*lm 
@param m     grating order
@param k     grating constant
@param l     grating wavelength
@param alpha grating alpha
 */
double uves_physmod_find_beta(int m, double k, double l, double alpha)
{
   /* uves_msg ("m, k, l, alpha: %d %f %f %f",m, k, l, alpha); */
 
  l /= uves_air_index(l);
  
  /*  I check that the argument of asin is <=0 otherwhise return
      a dummy angle (~89 deg) which will produce (x,y) pos out
      of detector so the line will be discarted
   */
  /*
  uves_msg_debug("l=%g m*k*l=%g alpha=%g",l,m*k*l,alpha);
  */
  if ( (m*k*l - sin(alpha*uves_deg2rad)) <=1.0 )
  {
     return( (asin(m*k*l - sin(alpha*uves_deg2rad)))/uves_deg2rad );
  }
  else
  {
     return( (asin(0.999))/uves_deg2rad );
  }
}

/**
@brief determines beta ech, CD angle, focal length
@param lambda          wavelength
@param m               order
@param puves_beta_ech  echelle beta angle
@param puves_beta_cd   CD angle
@param pfc             focal length
 */
void 
uves_physmod_lambda_order2beta(double lambda, 
                               int m, 
                               double* puves_beta_ech, 
                               double* puves_beta_cd, 
                               double* pfc)
{
    /* uves_msg ("Disp: %d Ech: %d",uves_x_disp_id,uves_ech_id-1); */

/*   lambda /= uves_air_index(lambda);  bog fixed ! */
   
   *pfc       = cameraFocal(lambda);

    /* uves_msg ("New Camera focal (m) : %f",fc); */

   *puves_beta_ech = uves_physmod_find_beta(m, uves_ech_groov[uves_ech_id-1], 
                     lambda, uves_ech_blaze[uves_ech_id-1]);

   *puves_beta_cd  = uves_physmod_find_beta(1, cdgroov[uves_x_disp_id-1], 
                     lambda, uves_alpha0_cd); 
   /* uves_msg_debug("fc=%g uves_beta_ech=%g uves_beta_cd=%g"
      ,*pfc,*puves_beta_ech,*puves_beta_cd); */

}

/**
@brief finds imaging points after flip and rotation and offset
@param uves_beta_cd   Cross Disperser beta angle
@param uves_beta_ech  Echelle beta angle
@param fc             focal length
@param px             new x position
@param py             new y position
 */
void 
uves_beta2xy(double uves_beta_cd, 
             double uves_beta_ech, 
             double fc, 
             double* px, 
             double* py)
{
   double  xd, yd, xr, yr, angle;   
    
   xd       = fc*tan( (uves_beta_ech - uves_ech_blaze[uves_ech_id-1])*
                       uves_deg2rad )/
                       uves_physmod_pix_size[0]/uves_bin[0];    
   yd       = fc*tan( (uves_alpha0_cd - uves_beta_cd - 
                       cdbeam[uves_ech_id-1])*uves_deg2rad )/
                       uves_physmod_pix_size[1]/uves_bin[1];

   //CHECK
   uves_msg_debug("beta(CD), yorg: %f %f", uves_beta_cd, yd); 

   xd = xd*flipx;
   yd = yd*flipy;   
   uves_msg_debug ("Positions after flip: %f %f",xd,yd); 
    
   angle    = uves_ccd_rot[uves_ech_id-1]*uves_deg2rad;
   xr       =  xd*cos(angle) + yd*sin(angle);
   yr       = -xd*sin(angle) + yd*cos(angle);

   uves_msg_debug ("Rotated positions %f %f",xr,yr);

   *px = uves_physmod_offsetx[uves_cfg_indx-1] / uves_bin[0] + xr;
   *py = uves_physmod_offsety[uves_cfg_indx-1] / uves_bin[1] + yr;
    

}

/**
@brief computes the blaze function
@param lambda               wavelength
@param uves_beta_ech        echelle beta angle
@param uves_beta_cd         Cross Disperser beta angle
@param puves_physmod_rech   echelle beta angle in radians
@param puves_physmod_rcd    Cross Disperser beta angle in radians
@param pblz                 blaze

 */
void 
uves_physmod_photo_beta(double lambda, 
                        double uves_beta_ech, 
                        double uves_beta_cd, 
                        double* puves_physmod_rech, 
                        double* puves_physmod_rcd, 
                        double* pblz)
{
   double gam;
   /*
   uves_msg("uves_ech_id = %d, uves_ech_blaze = %f",
             uves_ech_id,uves_ech_blaze[uves_ech_id-1]); 
   uves_msg("uves_deg2rad=%f uves_beta_ech=%f uves_alpha0_cd=%f",
             uves_deg2rad,uves_beta_ech,uves_alpha0_cd); 
   */
   *puves_physmod_rech = cos(uves_ech_blaze[uves_ech_id-1]*uves_deg2rad)/
                         cos(uves_beta_ech*uves_deg2rad);
   
   *puves_physmod_rcd = cos(uves_alpha0_cd*uves_deg2rad)/
                        cos(uves_beta_cd*uves_deg2rad);
   /*
   uves_msg("puves_physmod_rech=%f *puves_physmod_rcd=%f",
             *puves_physmod_rech,*puves_physmod_rcd);
   */
   /* Computes the blaze function */
   gam = M_PI / uves_ech_groov[uves_ech_id-1] * 
                cos(delta[uves_ech_id-1]*uves_deg2rad) * 
         (sin((uves_beta_ech-delta[uves_ech_id-1])*uves_deg2rad)-
          sin((uves_ech_blaze[uves_ech_id-1]-delta[uves_ech_id-1])*
              uves_deg2rad))
     / lambda/ uves_air_index(lambda);
   *pblz = DSQR(sin(gam)/gam);
   /*
   uves_msg("gamma = %f, Blaze function = %g ",gam,blz); 
   */  
}

/**
@brief Computation of the true pixel-scale along echelle dispersion
@param wave  wavelength
@param order order
@param uves_physmod_rech echelle angle in radians 
@param uves_physmod_rcd  Cross Disperser angle in radians
@param binx              x binning
@param biny              y binning
@param fc                focal length
@param slitwidth         slit width
@param slitlength        slit length
@param pbinsize          output pixel size
@param ppixscale         output pixel scale
@param ppixscalCD        output pixel scale on Cross Disperser direction
@param plinewidpx        output line width in pixels
@param plinelenpx        output line length in pixels
@param plinewidth        output line width
@param presol            output resolution

     Computation of the true pixel-scale along echelle dispersion 
     scale = (K*FocalRatioNasmyth*PixelSize*ApertureCamera*1e-3)/FocalNasmyth
     K = 206265 : conversion radian to arcsec (1 radian = 206265 arcsec)
 */
void uves_physmod_pixscl(
                         double wave, 
                         int order, 
                         double uves_physmod_rech, 
                         double uves_physmod_rcd, 
                         float binx, 
                         float biny, 
                         double fc, 
                         double slitwidth, 
                         double slitlength,
                         double* pbinsize, 
                         double* ppixscale, 
                         double* ppixscalCD, 
                     double* plinewidpx, 
                         double* plinelenpx, 
                         double* plinewidth, 
                         double* presol)
{
   double binsize, pixscale, pixscalCD, linewidpx, linelenpx, linewidth, resol;
   static double scale;

   if (!(scale)) scale = (206265.0*15.0*0.015*200*1e-3*binx)/120000;

   /*
   uves_msg("scale=%f 
             wave=%f 
             order=%d 
             fc=%f 
             uves_physmod_rech=%f 
             binx=%f 
             biny=%f 
             uves_physmod_rcd=%f",
         scale,    
             wave,   
             order,   
             fc,   
             uves_physmod_rech,   
             binx,   
             biny,   
             uves_physmod_rcd);
   */
/* Computes the width (in pixel and A) and resolution lines */
   binsize = uves_physmod_wave_bin(wave, order) * 1e4;  /* in mA */
   pixscale = scale/(fc*uves_physmod_rech); /* in arcsec/pixel */
   pixscalCD = pixscale *(biny/binx) *uves_physmod_rech/uves_physmod_rcd; /* in arcsec/pixel */ 
   linewidpx = slitwidth / pixscale;  /* pixel */
   linelenpx = slitlength /pixscalCD;  /* pixel */    
   linewidth = binsize * linewidpx * 1e-3; /* in A (* pixel) */
   resol = wave * 10.0 / linewidth;  /* without unit */
                              /* (10.0: conversion factor from nm to A)*/

   /*
   uves_msg("slitwidth=%f slitlength=%f binsize=%f pixscale=%f pixscaleCD=%f",
      slitwidth,   slitlength,   binsize,   pixscale,   pixscalCD);

   uves_msg("linewidpx=%f linewidth=%f resol=%f",
      linewidpx,    linewidth,   resol);
   */
   *pbinsize = binsize;
   *ppixscale = pixscale;
   *ppixscalCD = pixscalCD;
   *plinewidpx = linewidpx;
   *plinelenpx = linelenpx;
   *plinewidth = linewidth;
   *presol = resol;

}

/**
@brief computes x,y positions predicted by the model
@param lambda wavelength
@param m      order
@param px     predicted x pos
@param py     predicted y pos
 */
void 
uves_physmod_xy_model(double lambda, int m, double* px, double* py)
{
   double fc, uves_beta_ech, uves_beta_cd;
   //CHECK
   uves_physmod_lambda_order2beta(lambda, m, &uves_beta_ech, &uves_beta_cd, &fc);
   uves_msg_debug("lambda=%f m=%d uves_beta_ech=%f,uves_beta_cd=%f,fc=%f",
		    lambda,m,uves_beta_ech,uves_beta_cd,fc);

   uves_beta2xy(uves_beta_cd, uves_beta_ech, fc, px, py);
   uves_msg_debug("px=%f py=%f",*px,*py);

/* exemple : 
   uves_physmod_photo_beta(lambda, uves_beta_ech, uves_beta_cd, 
                           puves_physmod_rech, puves_physmod_rcd, pblz); 
*/

}

/**
@brief computes predicted echelle beta angle and cross disperser angle
@param puves_beta_ech predicted echelle beta angle
@param puves_beta_cd  predicted cross disperser angle
@param fc             focal length
@param x              x pos
@param y              y pos
 */
void 
uves_physmod_xy2beta(double* puves_beta_ech, 
                     double* puves_beta_cd, 
                     double fc, 
                     double x, 
                     double y)
{
   double xr, yr, xd, yd, angle;

   angle    = uves_ccd_rot[uves_ech_id-1]*uves_deg2rad;

   xr = (x - uves_physmod_offsetx[uves_cfg_indx-1]/uves_bin[0]);
   yr = (y - uves_physmod_offsety[uves_cfg_indx-1]/uves_bin[1]);

   xd = xr*cos(angle) - yr*sin(angle);
   yd = xr*sin(angle) + yr*cos(angle);

   xd  /= flipx;
   yd  /= flipy;

   *puves_beta_ech = atan(xd*uves_physmod_pix_size[0]*uves_bin[0]/fc)/
                         uves_deg2rad + uves_ech_blaze[uves_ech_id-1];
   *puves_beta_cd  =  uves_alpha0_cd - cdbeam[uves_ech_id-1] - 
                      atan(yd*uves_physmod_pix_size[1]*uves_bin[1]/fc)/
                      uves_deg2rad; 

}

/**
@brief computes the wavelength corresponding to a given echelle beta and cross
disperser beta and a given order
@param uves_beta_ech  input beta echelle angle
@param uves_beta_cd   input beta cross disperser angle
@param plambda        output wavelength
@param m              input order
 */
static void 
beta2lamb(double uves_beta_ech, 
          double uves_beta_cd, 
          double* plambda, 
          int m)
{
  uves_beta_cd=uves_beta_cd;
   *plambda = uves_physmod_find_lambda(m*uves_ech_groov[uves_ech_id-1], 
              uves_ech_blaze[uves_ech_id-1], uves_beta_ech);

}

/**
@brief
     Finds the couple lambda,m corresponding to a given position (x,y), 
     assuming the focal of the camera for this wavelength is fc. 

@param plambda   output wavelength
@param pdm       output order
@param fc        input focal length
@param x         input x
@param y         input y

 */
void 
uves_physmod_lambda_order_focus_model(double* plambda, 
                                      double* pdm, 
                                      double fc, 
                                      double x, 
                                      double y)
{

   double uves_beta_ech, uves_beta_cd;
   
   uves_physmod_xy2beta(&uves_beta_ech, &uves_beta_cd, fc, x, y);
   
   *plambda = uves_physmod_find_order_lambda(cdgroov[uves_x_disp_id-1], 
               uves_alpha0_cd, uves_beta_cd);
   *pdm      = uves_physmod_find_order_lambda(uves_ech_groov[uves_ech_id-1], 
                uves_ech_blaze[uves_ech_id-1], uves_beta_ech)/(*plambda);
/*
   uves_msg("  m= %f, lambda= %f, position (x,y)= (%f , %f)",
            *pdm,*plambda, x,y);
*/
}

/************************************************************************* 
 *  Finds the couple lambda,m corresponding to a given position (x,y),
 *  assuming the focal of the camera for this wavelength is fc. 
 *************************************************************************
 */
/* Not used:
void 
uves_physmod_lambda_model(double* plambda, 
                          int m, 
                          double fc, 
                          double x, 
                          double y)
{

   double uves_beta_ech, uves_beta_cd;

   uves_physmod_xy2beta(&uves_beta_ech, &uves_beta_cd, fc, x, y);


   beta2lamb(uves_beta_ech, uves_beta_cd, plambda, m);

}
*/

/**
@brief
  Iteratively finds the couple lambda,m corresponding to a given position 
  (x,y), by starting with an initial guess and iteratively refining it. 
  assuming the focal of the camera for this wavelength is fc. 

@param plambda          output wavelength
@param pm               output order
@param x                input x
@param y                input y
@param puves_beta_ech   output echelle beta angle
@param puves_beta_cd    output ecross disperser beta angle

 */
void 
uves_physmod_lambda_order_beta(double* plambda, 
                               int* pm, 
                               double x, 
                               double y, 
                               double* puves_beta_ech, 
                               double* puves_beta_cd)
{

double fcguess=0., wave=0., mdbl=0., xe=0., ye=0., xd=0., yd=0.;
int  i=0;

if (uves_arm_id == 'b') fcguess = 0.360;/*0.35722;*/
if (uves_arm_id == 'r') fcguess = 0.500;

uves_physmod_lambda_order_focus_model(&wave,&mdbl,fcguess,x,y);
/*
if (x == 1500.) {
  uves_msg("m= %f, lambda= %f, position (x,y)= (%f , %f)",mdbl,wave, x,y);
    uves_msg("focal guess= %f",fcguess);
}
*/
do {
   
   fcguess = cameraFocal(wave);
   /*
     if (x == 1500.) { 
       uves_msg("i= %d, focal = %f, m= %f, lambda= %f", i,fcguess,mdbl,wave);
     }
   */
   if (*pm <= 0) *pm = (int)(mdbl+0.5);
/*   uves_physmod_lambda_model(&wave,*pm,fcguess,x,y); */
   uves_physmod_xy2beta(puves_beta_ech, puves_beta_cd, fcguess, x, y);
   beta2lamb(*puves_beta_ech, *puves_beta_cd, &wave, *pm);

   /*
     if (x == 1500.) { 
       uves_msg("i= %d, focal = %f, m= %d, lambda= %f", i,fcguess,*pm,wave);
     }
   */
   
   uves_physmod_xy_model(wave,(int)(mdbl+0.5), &xe, &ye);
   /*
      if (x == 1500.) { 
        uves_msg("m= %f, lambda= %f, position (xe,ye)= (%f , %f)",
                  mdbl,wave, xe,ye);
        uves_msg("focal = %f",fcguess);
      }
   */
   i++;
   
   xd = fabs(x-xe);
   yd = fabs(y-ye);
   
   } while (!((xd < 1.) && (yd < 1.)) && (i <= 4));

*plambda = wave;

/*uves_msg("i= %d", i);*/
/* uves_physmod_photo_beta(wave, *puves_beta_ech, *puves_beta_cd, 
                 puves_physmod_rech, puves_physmod_rcd, pblz); 
uves_msg("uves_physmod_lambda_order_beta(%f, %d, %f, %f, %f, %f)",
          wave, *pm, x, y, *puves_beta_ech, *puves_beta_cd);
*/

}

/**
@brief
   Iteratively finds the couple lambda,m corresponding to a given position 
   (x,y), by starting with an initial guess and iteratively refining it. 
   assuming the focal of the camera for this wavelength is fc. 


@param plambda          output wavelength
@param pm               output order
@param x                input x
@param y                input y

 */
void 
uves_physmod_lambda_order_model(double* plambda, int* pm, double x, double y)
{

  double uves_beta_ech, uves_beta_cd;
  
  uves_physmod_lambda_order_beta(plambda, pm, x, y, 
                                 &uves_beta_ech, &uves_beta_cd);
  
/* exemple : 
  uves_physmod_photo_beta(*plambda, &uves_beta_ech, &uves_beta_cd, 
                   puves_physmod_rech, puves_physmod_rcd, pblz); 

uves_msg("uves_physmod_lambda_order_beta(%f, %d, %f, %f, %f, %f)",
         *plambda, *pm, x, y, uves_beta_ech, uves_beta_cd);
*/
}

/**
@brief Finds start and end of the free-spectral range 
@param m            order
@param lambdaC      central wavelength
@param fsrStart     start wavelength free spectral range
@param fsrEnd       end   wavelength free spectral range

 */
void 
uves_physmod_find_FSR(int m, double* lambdaC, double* fsrStart, double* fsrEnd)
{

  double tmp_delta;

  *lambdaC = 2*sin(uves_ech_blaze[uves_ech_id-1]*uves_deg2rad)/
                   m/uves_ech_groov[uves_ech_id-1];

  tmp_delta = *lambdaC/m;

  *fsrStart = *lambdaC - tmp_delta/2.;
  *fsrEnd   = *lambdaC + tmp_delta/2.;

}

/**
@brief Finds start and end of the free-spectral range 
@param l wavelength
@param m order
 */
double uves_physmod_wave_bin(double l, int m)
{

  double dl, x0,y_0,x1,y_1;

  dl = 20e-4; /* Wavelength increment in nm (nearly one pixel)*/

  uves_physmod_xy_model( l,    m,&x0,&y_0);
  uves_physmod_xy_model((l+dl),m,&x1,&y_1);

  return( dl/(x1-x0) );
}

/**
@brief Finds x and y dimension of the selected CCD 
@param nx x dimension
@param ny y dimension
 */
void uves_ccd_size(int* nx, int* ny)
{
  /*
  uves_msg("imsize[uves_ech_id-1]=%d 
            uves_physmod_row_size[uves_ech_id-1]=%d 
            uves_bin[0]=%f uves_bin[1]=%f",
        imsize[uves_ech_id-1],   
            uves_physmod_row_size[uves_ech_id-1],   
            uves_bin[0],   
            uves_bin[1]);
  */
  *nx = imsize[uves_ech_id-1] / uves_bin[0];
  *ny = uves_physmod_row_size[uves_ech_id-1] / uves_bin[1];
}

/**
@brief Computes corrected x,y positions
@param x  input x
@param y  input y
@param px predicted x
@param py predicted y



 */
void uves_physmod_xy_regres(double x,double y,double* px,double* py)
{
  double xdiff=0;
  double ydiff=0;
  int    xnpix=0;
  int    ynpix=0;
  
  goto simplified;
  /*  We comment the following to remove a compilation warning
      Anyway the code would not be executed due to goto statement      
  xdiff = xcoef[uves_cfg_indx-1][8]*DSQR(x*y) + 
          xcoef[uves_cfg_indx-1][7]*x*DSQR(y) + 
          xcoef[uves_cfg_indx-1][6]*DSQR(y)   + 
          xcoef[uves_cfg_indx-1][5]*DSQR(x)*y + 
          xcoef[uves_cfg_indx-1][4]*x*y       + 
          xcoef[uves_cfg_indx-1][3]*y         + 
          xcoef[uves_cfg_indx-1][2]*DSQR(x)   + 
          xcoef[uves_cfg_indx-1][1]*x         + 
          xcoef[uves_cfg_indx-1][0];

  ydiff = ycoef[uves_cfg_indx-1][8]*DSQR(x*y) + 
          ycoef[uves_cfg_indx-1][7]*x*DSQR(y) + 
          ycoef[uves_cfg_indx-1][6]*DSQR(y)   + 
          ycoef[uves_cfg_indx-1][5]*DSQR(x)*y + 
          ycoef[uves_cfg_indx-1][4]*x*y       + 
          ycoef[uves_cfg_indx-1][3]*y         + 
          ycoef[uves_cfg_indx-1][2]*DSQR(x)   + 
          ycoef[uves_cfg_indx-1][1]*x         + 
          ycoef[uves_cfg_indx-1][0];
  */
  
  /* New, simplified correction */ 
  simplified: {
  uves_ccd_size(&xnpix, &ynpix);
  /* uves_msg("xnpix=%d ynpix=%d",xnpix,ynpix); */
  xdiff = (-7.)*(x-(double)xnpix/2.)/((double)xnpix/2.);
  /* ydiff = (5.)*pow((x-(double)xnpix/2.)/((double)xnpix/2.),2.); */
  ydiff = (5.)*DSQR((x-(double)xnpix/2.)/((double)xnpix/2.));
  /* uves_msg("xdiff=%f ydiff=%f",xdiff,ydiff); */
  *px = x + xdiff;
  *py = y + ydiff;
  }
}

/**@}*/