text:
/**
 */
#include <math.h>
#include "fitsio.h"
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_interp.h>
#include "spc_cfg.h"
#include "aXe_utils.h"
#include "aXe_errors.h"
#include "fringe_conf.h"

#define AXE_CONFIG_PATH "AXE_CONFIG_PATH"

/**
 * Function: load_CCD_layer
 * The function creates and returns a structure for an optical
 * CCD layer. The input parameters are interpreted, and
 * the data which defines the CCD layer is loaded.
 *
 * Parameters:
 * @param refr_table - the name of the refractive index table
 * @param thickness  - the value for the thicknes keyword
                       (number or image name)
 *
 * Returns:
 * @return opt_layer  - the optical layer created
 */
ccd_layer *
load_CCD_layer(const char refr_table[], const char thickness[])
{
  ccd_layer *opt_layer;

  float tvalue;

  char **t_err=NULL;

  char refr_table_path[MAXCHAR];
  char thickness_path[MAXCHAR];

  // build the full pathname to the refraction table
  build_path (AXE_CONFIG_PATH, refr_table, refr_table_path);

  // allocate an error array
  t_err = (char **) malloc(sizeof(char *)*1);

  // allocate space for the return structure;
  // complain if this fails
  opt_layer = (ccd_layer *)malloc(sizeof(ccd_layer));

  // initialize both, the single value
  // as well as the matrix to NULL;
  opt_layer->thickness   = 0.0;
  opt_layer->thickness2D = NULL;

  // convert the keyvalue to a float
  // tvalue = atof(thickness);
  tvalue = strtod(thickness,t_err);

  // in case the float value is NULL = 0.0
  // the conversion failed, and it must be a string
  //  if (tvalue)
  if (strcmp(thickness, t_err[0]))
    {
#ifdef DEBUGFCONF
      fprintf(stderr, "Setting the layer thickness to value: %f mum\n",tvalue);
#endif
      // set the fixed thickness value
      opt_layer->thickness = tvalue;
    }
  else
    {
#ifdef DEBUGFCONF
      fprintf(stderr, "Loading layer thickness image: %s\n", thickness);
#endif

      // build the full pathname to the thickness image
      build_path (AXE_CONFIG_PATH, thickness, thickness_path);

      // load the 2D image for the thickness
      opt_layer->thickness2D = FITSimage_to_gsl(thickness_path, 1, 1);
    }

#ifdef DEBUGFCONF
      fprintf(stderr, "Loading refractive index table: %s\n", refr_table);
#endif
  // load the real part of the refraction index
  opt_layer->re_refraction =
    create_interp_ftable(refr_table_path, 2, "WAVELENGTH", "N",
			 REFRAC_INTERP_TYPE);

  // load the real part of the refraction index
  opt_layer->im_refraction =
    create_interp_ftable(refr_table_path, 2, "WAVELENGTH", "K",
			 REFRAC_INTERP_TYPE);

  // return the optical layer
  return opt_layer;
}

/**
 * Function: free_CCD_layer
 * The function releases all the memory allocated
 * in a structure for a CCD layer.
 *
 * Parameters:
 * @param opt_layer - structure for CCD layer
 *
 * Returns:
 * @return -
 */
void
free_CCD_layer(ccd_layer *opt_layer)
{
  // free the thickness matrix
  if (opt_layer->thickness2D)
    gsl_matrix_free(opt_layer->thickness2D);

  // free both interpolators
  free_interp(opt_layer->re_refraction);
  free_interp(opt_layer->im_refraction);

  // free the rest
  free(opt_layer);

  // set the structure to NULL
  opt_layer = NULL;
}

/**
 * Function: load_CCD_layers
 * The function extracts from a fringe configuration file
 * all information on individual CCD layers. It creates
 * the appropriate structure for every CCD layer,
 * and creates and returns a structure which completely
 * describes all layers in a CCD.
 *
 * Parameters:
 * @param fring_conf_path - the name of a fringe configuration file
 *
 * Returns:
 * @return opt_layers - the structure for the optical layers created
 */
ccd_layers *
load_CCD_layers(char fring_conf_path[])
{
  ccd_layers *opt_layers;

  char layer[MAXCHAR];
  char refr_table[MAXCHAR];

  //int nlayers;
  int i=0;

  struct CfgStrings LayerConfig[] = {
    {NULL, NULL},
    {NULL, NULL}
  };

  LayerConfig[0].name = layer;

  // allocate space for the return structure;
  // complain if this fails
  opt_layers = malloc (sizeof (ccd_layers));

  // initialize the substrate
  opt_layers->substrate = NULL;

  // detemrine the number of layers
  // which are described in the
  //configuration file
  opt_layers->num_layers = 0;
  for (i = 0; i < MAX_BEAMS; i++)
    {
      sprintf (layer, "REFR_INDEX_%c", BEAM (i));
      CfgRead (fring_conf_path, LayerConfig);
      if (LayerConfig[0].data != NULL)
	{
	  opt_layers->num_layers += 1;
	  LayerConfig[0].data = NULL;
	}
    }

  //  if (opt_layers->num_layers < NLAYERS_MIN)
  //    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
  //		 "The number of layers in the configuration\n file is %i. "
  //		 "There must be at least: %i!\n",
  //		 opt_layers->num_layers, NLAYERS_MIN);

  // allocate space for the array to
  // the layers described in the
  // configuration file
  opt_layers->opt_layer =
    (ccd_layer **)malloc (opt_layers->num_layers*sizeof(ccd_layers));

  // successively find the
  // description and finally
  // load all layers
  for (i = 0; i < MAX_BEAMS; i++)
    {
      // chekc for the first relevant keyword
      sprintf (layer, "REFR_INDEX_%c", BEAM (i));
      CfgRead (fring_conf_path, LayerConfig);
      if (LayerConfig[0].data != NULL)
	{
	  sprintf (refr_table, "%s", LayerConfig[0].data);
	  LayerConfig[0].data = NULL;

	  // check for the second relevant keyword
	  sprintf (layer, "THICKNESS_%c", BEAM (i));
	  CfgRead (fring_conf_path, LayerConfig);
	  if (LayerConfig[0].data != NULL)
	    {
	      // load the layer
	      opt_layers->opt_layer[i] =
		load_CCD_layer(refr_table, LayerConfig[0].data);

	      LayerConfig[0].data = NULL;
	    }
	  else
	    {
	      // complain that a keyword is missing
	      // for a layer
	      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
			   "Could load index table %s,\n"
			   "but not thickness information "
			   "for beam %c\n", refr_table, BEAM (i));
	    }
	}
    }

  // return the whole structure
  return opt_layers;
}

/**
 * Function: free_CCD_layers
 * The function releases all memory allocated in a
 * CCD layers structure.
 *
 * Parameters:
 * @param opt_layers - the CCD layers structure
 *
 * Returns:
 * @return -
 */
void
free_CCD_layers(ccd_layers *opt_layers)
{
  int index=0;

  // free each layer individually
  for (index=0; index < opt_layers->num_layers; index++)
      free_CCD_layer(opt_layers->opt_layer[index]);

  // free the interpolator for the substrate
  //  free_linint(opt_layers->substrate);
  free_interp(opt_layers->substrate);

  // free everything
  free(opt_layers);

  // set the structure to NULL
  opt_layers = NULL;
}

/**
 * Function: load_fringe_conf
 * The function loads the fringe configuration file given as parameter.
 * All relevant data marked with keywords is extracted. A fringe
 * configuration structure is built up and returned from the
 * data given in the keyvalues.
 *
 * Parameters:
 * @param fring_conf_path - the complete path-name to a fringe configuration file
 *
 * Returns:
 * @return fconf - the fringe configuration structure created
 */
fringe_conf *
load_fringe_conf(char fring_conf_path[])
{
  char beam[MAXCHAR];

  char ffile_name[MAXCHAR];
  char ffile_name_path[MAXCHAR];

  int index=0;

  fringe_conf *fconf=NULL;

  gsl_vector *v=NULL;

  struct CfgStrings FringeConfig[] =
    {
      {"FRINGE_AMPLITUDE", NULL},
      {"FRINGE_PHASE", NULL},
      {"FRINGE_RANGE", NULL},
      {"MAX_DISPERSION", NULL},
      {"FILTER_NAME", NULL},
      {"FRINGE_STEP",NULL},
      {"NUM_STEPS",NULL},
      {"SUBSTR_TRANS",NULL},

      {NULL, NULL},
      {NULL, NULL}		/* array terminator. REQUIRED !!! */
    };

  FringeConfig[8].name = beam;

  // allocate space for the return structure;
  // complain if this fails
  fconf = malloc (sizeof (fringe_conf));
  if (fconf == NULL)
    {
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "Could not allocate memory for fringe configuration");
    }

  // load the optical layers directly from the file
  fconf->opt_layers = load_CCD_layers(fring_conf_path);

  // read in the file
  CfgRead (fring_conf_path, FringeConfig);

  // make some intitializations
  fconf->fringe_amp     = 0.0;
  fconf->fringe_phase   = 1.0e+32;
  fconf->fringe_step    = 0.0;
  fconf->num_steps      = 0;
  fconf->fringe_range   = NULL;
  fconf->max_dispersion = 0.0;
  fconf->filter_through = NULL;

  for (index = 0; index < 9; index++)
    {

      // read in the fringe amplitude
      if (!strcmp (FringeConfig[index].name, "FRINGE_AMPLITUDE"))
	{
	  if (FringeConfig[index].data != NULL)
	    {
	      fconf->fringe_amp = atof(FringeConfig[index].data);
	    }
	}

      // read in the fringe phase
      if (!strcmp (FringeConfig[index].name, "FRINGE_PHASE"))
	{
	  if (FringeConfig[index].data != NULL)
	    {
	      fconf->fringe_phase = atof(FringeConfig[index].data);
	    }
	}

      // read in the fringe step
      if (!strcmp (FringeConfig[index].name, "FRINGE_STEP"))
	{
	  if (FringeConfig[index].data != NULL)
	    {
	      fconf->fringe_step = atof(FringeConfig[index].data);
	    }
	}

      // read in the number of steps
      if (!strcmp (FringeConfig[index].name, "NUM_STEPS"))
	{
	  if (FringeConfig[index].data != NULL)
	    {
	      fconf->num_steps = atoi(FringeConfig[index].data);
	    }
	}

      // read in the fringe range
      if (!strcmp (FringeConfig[index].name, "FRINGE_RANGE"))
	{
	  if (FringeConfig[index].data != NULL)
	    {
	      v = string_to_gsl_array (FringeConfig[index].data);
	      fconf->fringe_range = v;
	    }
	}

      // read in the fringe amplitude
      if (!strcmp (FringeConfig[index].name, "MAX_DISPERSION"))
	{
	  if (FringeConfig[index].data != NULL)
	    {
	      fconf->max_dispersion = atof(FringeConfig[index].data);
	    }
	}

      // read in the filter name
      if (!strcmp (FringeConfig[index].name, "FILTER_NAME"))
	{
	  if (FringeConfig[index].data != NULL)
	    {
	      sprintf (ffile_name, "%s", FringeConfig[index].data);

	      // build the full pathname to the refraction table
	      build_path (AXE_CONFIG_PATH, ffile_name, ffile_name_path);

	      fconf->filter_through =
		create_interp_ftable(ffile_name_path, 2, "WAVELENGTH",
				     "THROUGHPUT", FILTER_INTERP_TYPE);
	    }
	}

      // read in the substrate transmission
      if (!strcmp (FringeConfig[index].name, "SUBSTR_TRANS"))
	{
	  if (FringeConfig[index].data != NULL)
	    {
	      sprintf (ffile_name, "%s", FringeConfig[index].data);

	      // build the full pathname to the refraction table
	      build_path (AXE_CONFIG_PATH, ffile_name, ffile_name_path);

	      fconf->opt_layers->substrate =
		create_interp_ftable(ffile_name_path, 2, "WAVELENGTH",
				     "TRANSMISSION", REFRAC_INTERP_TYPE);
	    }
	}
    }

  return fconf;
}

/**
 * Function: check_fringe_conf
 *
 * Parameters:
 * @param fconf  - the fringe configuration structure
 *
 * Returns:
 * @return
 */
void
check_fringe_conf(fringe_conf *fconf)
{
  // make  sure the fringe amplitude is set
  if (!fconf->fringe_amp)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		 "fringe_conf: The fringing amplitude must\n"
		 "be set to be able determining pixel fringing!\n");

  // make  sure the fringe phase is set
  if (fconf->fringe_phase > 1.0e+31)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		 "fringe_conf: The fringing phase must\n"
		 "be set to be able determining pixel fringing!\n");

  // make  sure the transmission of the substrate
  // layer is given
  if (!fconf->opt_layers->substrate)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		 "fringe_conf: The substrate transmission\n"
		 "must be given to be able determining pixel fringing!\n");

  // make sure that there is a minimum
  // number of layers.
  if (fconf->opt_layers->num_layers < NLAYERS_MIN)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		 "The number of layers in the configuration\n file is %i. "
		 "There must be at least: %i!\n",
		 fconf->opt_layers->num_layers, NLAYERS_MIN);

  // if not defined, give a meaningfull
  // default for the fringing range
   if (!fconf->fringe_range)
    {
      fconf->fringe_range = gsl_vector_alloc(2);
      gsl_vector_set(fconf->fringe_range, 0, DEFAULT_MIN_LAMBDA);
      gsl_vector_set(fconf->fringe_range, 1, DEFAULT_MAX_LAMBDA);
      fprintf(stdout, "Setting the fringing range to defaults: [ %f, %f].\n",
	      DEFAULT_MIN_LAMBDA, DEFAULT_MAX_LAMBDA);
    }

   // give reasonable default for
   // the fringe step
   if (!fconf->fringe_step)
     {
       fconf->fringe_step = DEFAULT_FRINGE_STEP;
       fprintf(stdout, "Setting the fringing step to default: %f.\n",
	       DEFAULT_FRINGE_STEP);
     }

   // give a reasonable default for the
   // number of steps
   if (!fconf->num_steps)
     {
       fconf->num_steps = DEFAULT_NUM_STEPS;
       fprintf(stdout, "Setting the number of steps to default: %i.\n",
	       DEFAULT_NUM_STEPS);
     }

   // give a reasonable default for the
   // maximum dispersion to correct for
   if (!fconf->max_dispersion)
     {
       fconf->max_dispersion = DEFAULT_MAX_DISPERSION;
       fprintf(stdout, "Setting the maximal dispersion to default: %f.\n",
	       DEFAULT_MAX_DISPERSION);
     }
}

/**
 * Function: free_fringe_conf
 * The function deallocates all memory in a
 * fringe configuration structure.
 *
 * Parameters:
 * @param fconf  - the fringe configuration structure
 *
 * Returns:
 * @return -
 */
void
free_fringe_conf(fringe_conf *fconf)
{

  // free the filter throughput
  if (fconf->filter_through)
    free_interp(fconf->filter_through);

  // free the optical layers
  free_CCD_layers(fconf->opt_layers);

  // free the fringe range
  gsl_vector_free(fconf->fringe_range);
  fconf->fringe_range = NULL;

  // free everything
  free(fconf);
}

/**
 * Function: create_interp_ftable
 * This function cretaes an interpolator from data values stored
 * in a fits table. An interpolator is a structure to hold all
 * data to compute the interpolated data values for different
 * interpolation methods. Core of this structure is the interpolator
 * types offered in the gsl-library.
 * This function is a tool to extract the relevant data points from
 * a fits table and to set up and return the interpolator.
 *
 *
 * Parameters:
 * @param table_name  - the name of the fits table
 * @param hdunum      - the extension number with the data
 * @param xcol        - the column name with the indepenent values
 * @param ycol        - the column name with the dependent values
 * @param interp_type - the interpolation type
 *
 * Returns:
 * @return interp - the interpolation structure created
 */
interpolator *
create_interp_ftable(const char table_name[], const int hdunum,
		     char xcol[], char ycol[],
		     const gsl_interp_type *interp_type)
{

  interpolator *interp;

  double *x;
  double *y;
  int    f_status = 0;
  int hdutype, anynul;
  long nrows=0;
  int colnum=0;

  fitsfile *input;

  // allocate space for the return structure;
  // complain if this fails
  //    interp = (interpolator *)malloc (sizeof (interpolator));
  //   if (interp == NULL)
  //     {
  //        aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
  //  		   "Could not allocate memory for interpolator");
  //      }

#ifdef DEBUGFCONF
  fprintf(stderr, "Loading columns %s and %s of fitstable: %s\n", xcol, ycol, table_name);
#endif

  //  Open the file for reading
  fits_open_file (&input, table_name, READONLY, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "get_response_function_fromFITS: "
		   "Could not open" " file: %s",
		   table_name);
    }

  /* Move to the correct hdu */
  fits_movabs_hdu (input, hdunum, &hdutype, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "create_interp_ftable:"
		   "Could not read extention %d from file: %s",
		   hdunum, table_name);
    }

  /* Get number of rows */
  fits_get_num_rows (input, &nrows, &f_status);
  if (f_status) {
    ffrprt (stderr, f_status);
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		 "create_interp_ftable: "
		 "Could not determine the number of rows in"
		 " table %s",table_name);
  }

  /* Allocate temporary memory space */
  x = (double *) malloc(nrows*sizeof(double));
  if (!x) {
    aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
		 "Memory allocation failed");
  }
  y = (double *) malloc(nrows*sizeof(double));
  if (!y) {
    aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
		 "Memory allocation failed");
  }

  /**************************/
  /* Read the X-column      */
  /**************************/
  /* Get column number */
  fits_get_colnum (input, CASEINSEN, xcol, &colnum, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "create_interp_ftable: "
		   "Could not determine column %s in "
		   " table %s", xcol, table_name);
    }

  /* Read the data */
  fits_read_col (input, TDOUBLE, colnum, 1, 1, nrows, NULL, x,
                    &anynul, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "create_interp_ftable: "
		   "Could not read content of WAVELENGTH column "
		   " from BINARY table %s",table_name);
    }

  /**************************/
  /* Read the y-column      */
  /**************************/
  /* Get column number */
  fits_get_colnum (input, CASEINSEN, ycol, &colnum, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "create_interp_ftable: "
		   "Could not determine column %s in "
		   " table %s", ycol, table_name);
    }

  /* Read the data */
  fits_read_col (input, TDOUBLE, colnum, 1, 1, nrows, NULL, y,
                    &anynul, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "create_interp_ftable: "
		   "Could not read column %s"
		   " from BINARY table %s", ycol, table_name);
    }

  fits_close_file(input,&f_status);
  if (f_status) {
      aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
		   "Could not close %s", table_name);
  }

  // create the interpolator
  interp = create_interp(nrows, interp_type, x, y);

  // return the interpolator
  return interp;
}

/**
 * Function: create_interp
 * The function creates, and intializes and returns an interpolator
 * from the basic input, which is the data, the number of data
 * items and the interpolator type requested.
 *
 * Parameters:
 * @param nvals       - the number of data values in the arrays
 * @param interp_type - the interpolation type
 * @param xvals       - the array with the independent data values
 * @param yvals       - the array with the dependent data values
 *
 * Returns:
 * @return interp - the interpolation structure created
 */
interpolator *
create_interp(const int nvals, const gsl_interp_type *interp_type,
	      double *xvals, double *yvals)
{
  interpolator *interp;

  // allocate space for the return structure;
  // complain if this fails
  interp = (interpolator *)malloc (sizeof (interpolator));
  if (interp == NULL)
    {
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "Could not allocate memory for interpolator");
    }

  // store the min and max values
  interp->xmin    = xvals[0];
  interp->xmax    = xvals[nvals-1];

  // check for ascending order in the
  // independent values
  if (interp->xmax < interp->xmin)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		 "The independent data values to be stored\n"
		 " in an interplator must be in INCREASING order!\n");

  // store the number of data values
  interp->nvals   = nvals;

  // store the data arrays
  interp->xvals   = xvals;
  interp->yvals   = yvals;

  // create and intitialize the gsl inteprolator
  interp->acc     = gsl_interp_accel_alloc();
  interp->interp  = gsl_interp_alloc(interp_type, interp->nvals);
  gsl_interp_init(interp->interp, interp->xvals, interp->yvals, interp->nvals);

  // return the new structure
  return interp;
}

/**
 * Function: print_interp
 *
 * Parameters:
 * @param interp - the interpolator structure
 *
 * Returns:
 * @return -
 */
void
print_interp(interpolator *interp)
{
  int index;
  double x_typical;

  x_typical = interp->xmin+(interp->xmax-interp->xmin)/2.0;

  fprintf(stdout, "xmin: %e, xmax: %e\n",interp->xmin, interp->xmax);
  fprintf(stdout, "number of data values: %i\n",interp->nvals);
  fprintf(stdout, "Interpolation type: %s\n", gsl_interp_name(interp->interp));
  fprintf(stdout, "Characteristic value pair: (x,y) = (%e, %e)\n",
	  x_typical, gsl_interp_eval(interp->interp, interp->xvals,
				     interp->yvals, x_typical, interp->acc));

  fprintf(stdout, "Alternative value pair: (x,y) = (%e, %e)\n",
	  x_typical, eval_interp(interp, x_typical));
  for (index=0; index < interp->nvals; index++)
    fprintf(stdout, "xvalue: %e, yvalue: %e\n",interp->xvals[index], interp->yvals[index]);
  fprintf(stdout, "\n");
}

/**
 * Function: eval_interp
 * The function computes and returns the interpolated value
 * at a given position for an interpolater.
 *
 * Parameters:
 * @param interp - the interpolator
 * @param xval   - the position to evaluate the interpolator
 *
 * Returns:
 * @return (value) - the interpolated data value
 */
double
eval_interp(interpolator *interp, const double xval)
{
  // check whether the x-value is within
  // the range spanned by the data;
  // complain if the x-value is outside
  if (xval < interp->xmin || xval > interp->xmax)
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "independent interpolation value %f "
		   "is outside interval (%f, %f)\n", xval,
		   interp->xmin, interp->xmax);

  // evaluate and return the interpolated value
  // in on the spot
  return  gsl_interp_eval(interp->interp, interp->xvals, interp->yvals,
			  xval, interp->acc);
}

/**
 * Function: free_interp
 * the function frees all memory allocated in
 * an interpolator structure.
 *
 * Parameters:
 * @param interp - the interpolator structure
 *
 * Returns:
 * @return -
 */
void
free_interp(interpolator *interp)
{
  // free the data vectors
  free(interp->xvals);
  free(interp->yvals);

  // free the two gsl structures
  gsl_interp_accel_free(interp->acc);
  gsl_interp_free (interp->interp);

  // free the rest
  free(interp);

  // set it to NULL
  interp = NULL;
}

/**
 * Function: create_linint_ftable
 * This function creates a linear interpolator from data values stored
 * in a fits table. An linear interpolator is a structure to hold all
 * data to compute the linear interpolated data values at any
 * point bracketed by the data.
 * This function is a tool to extract the relevant data points from
 * a fits table and to set up and return the linear interpolator.
 *
 * Parameters:
 * @param table_name  - the name of the fits table
 * @param hdunum      - the extension number with the data
 * @param xcol        - the column name with the indepenent values
 * @param ycol        - the column name with the dependent values
 * @param interp_type - the interpolation type
 *
 * Returns:
 * @return lin_int - the linear interplator structure created
 */
linear_interp *
create_linint_ftable(const char table_name[], const int hdunum,
		     char xcol[], char ycol[])
{

  linear_interp *lin_int;

  double *x;
  double *y;
  int    f_status = 0;
  int hdutype, anynul;
  long nrows=0;
  int colnum;
  int index=0;

  gsl_vector *xvals;
  gsl_vector *yvals;

  fitsfile *input;

  // allocate space for the return structure;
  // complain if this fails
  lin_int = (linear_interp *)malloc (sizeof (linear_interp));
  if (lin_int == NULL)
    {
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "Could not allocate memory for linear interpolator");
    }

#ifdef DEBUGFCONF
  fprintf(stderr, "Loading columns %s and %s of fitstable: %s\n", xcol, ycol, table_name);
#endif

  //  Open the file for reading
  fits_open_file (&input, table_name, READONLY, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "get_response_function_fromFITS: "
		   "Could not open" " file: %s",
		   table_name);
    }

  /* Move to the correct hdu */
  fits_movabs_hdu (input, hdunum, &hdutype, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "create_interp_ftable:"
		   "Could not read extention %d from file: %s",
		   hdunum, table_name);
    }

  /* Get number of rows */
  fits_get_num_rows (input, &nrows, &f_status);
  if (f_status) {
    ffrprt (stderr, f_status);
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		 "create_interp_ftable: "
		 "Could not determine the number of rows in"
		 " table %s",table_name);
  }

  /* Allocate temporary memory space */
  x = (double *) malloc(nrows*sizeof(double));
  if (!x) {
    aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
		 "Memory allocation failed");
  }
  y = (double *) malloc(nrows*sizeof(double));
  if (!y) {
    aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
		 "Memory allocation failed");
  }

  /**************************/
  /* Read the X-column      */
  /**************************/
  /* Get column number */
  fits_get_colnum (input, CASEINSEN, xcol, &colnum, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "create_interp_ftable: "
		   "Could not determine column %s in "
		   " table %s", xcol, table_name);
    }

  /* Read the data */
  fits_read_col (input, TDOUBLE, colnum, 1, 1, nrows, NULL, x,
                    &anynul, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "create_interp_ftable: "
		   "Could not read content of WAVELENGTH column "
		   " from BINARY table %s",table_name);
    }

  /**************************/
  /* Read the y-column      */
  /**************************/
  /* Get column number */
  fits_get_colnum (input, CASEINSEN, ycol, &colnum, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "create_interp_ftable: "
		   "Could not determine column %s in "
		   " table %s", ycol, table_name);
    }

  /* Read the data */
  fits_read_col (input, TDOUBLE, colnum, 1, 1, nrows, NULL, y,
                    &anynul, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "create_interp_ftable: "
		   "Could not read column %s"
		   " from BINARY table %s", ycol, table_name);
    }

  fits_close_file(input,&f_status);
  if (f_status) {
      aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
		   "Could not close %s", table_name);
  }

#ifdef DEBUGFCONF
  fprintf(stderr, "Number of rows: %i\n", nrows);
#endif

  // allocate space for the two vectors
  xvals = gsl_vector_alloc(nrows);
  yvals = gsl_vector_alloc(nrows);

  // transport the values from the
  // arrays to the gsl-arrays
  for (index=0; index < nrows; index++)
    {
      gsl_vector_set(xvals, index, x[index]);
      gsl_vector_set(yvals, index, y[index]);
    }

  // reverse both vectors if necessary
  if (gsl_vector_get(xvals, 0) > gsl_vector_get(xvals, nrows-1))
    {
      gsl_vector_reverse(xvals);
      gsl_vector_reverse(yvals);
    }

  // set the two helper elements
  lin_int->act_index = 0;
  lin_int->num_elem  = nrows;

  // set the minimum and maximum
  lin_int->xmin = gsl_vector_get(xvals, 0);
  lin_int->xmax = gsl_vector_get(xvals, nrows-1);

  // set the two vectors
  lin_int->xvals = xvals;
  lin_int->yvals = yvals;

  // release memory in the temp variables
  free(x);
  free(y);

  // return the new structure
  return lin_int;
}

/**
 * Function: free_linint
 * The function frees all memory allocated
 * in a linear interpolator structure.
 *
 * Parameters:
 * @param lin_int - the linear inteprolator
 *
 * Returns:
 * @return -
 */
void
free_linint(linear_interp *lin_int)
{
  // free the two gsl-vectors
  gsl_vector_free(lin_int->xvals);
  gsl_vector_free(lin_int->yvals);

  // free the whole structure
  free(lin_int);

  // set the structure to NULL
  lin_int = NULL;
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
//Gets deltas (1st order) and delta-deltas (2nd order) differences of X,
//which is assumed to be pre-allocated of triple size,
//Thus, for dim==0, C%3 must equal 0, and C/3 is the original ncols of X.
//And,  for dim==1, R%3 must equal 0, and R/3 is the original nrows of X.

//I implement this just like FIR for speed and shorter code,
//except non-causal and mid-sample of B is 0,
//so I don't explicitly make B (e.g., B[n] just equals sc*n).

//Note that this may treat edge samples differently than other code.
//But this could be changed with a few lines of additional code here,
//without changing the super-efficient FIR implementation.
//That is, since out-of-range samps were assumed to be 0, nothing was
//added to Y for them, so can just add something later.

#include <stdio.h>
#include <cblas.h>

#ifdef __cplusplus
namespace ov {
extern "C" {
#endif

int add_delta_deltas_s (float *X, const int iscolmajor, const int R, const int C, const int dim, const int N);
int add_delta_deltas_d (double *X, const int iscolmajor, const int R, const int C, const int dim, const int N);

int add_delta_deltas_s (float *X, const int iscolmajor, const int R, const int C, const int dim, const int N)
{
    const float z = 0.0f;
    const int No = R*C/3;
    int r, c, n;
    float sc = 1.0f;

    //Checks
    if (R<1) { fprintf(stderr,"error in add_delta_deltas_s: R (nrows X) must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in add_delta_deltas_s: C (ncols X) must be positive\n"); return 1; }
    if (N<1) { fprintf(stderr,"error in add_delta_deltas_s: N (delta winlength) must be positive\n"); return 1; }
    if (dim==0 && C%3!=0) { fprintf(stderr,"error in add_delta_deltas_s: for dim==0, C (ncols X) must be divisible by 3\n"); return 1; }
    if (dim==1 && R%3!=0) { fprintf(stderr,"error in add_delta_deltas_s: for dim==1, R (nrows X) must be divisible by 3\n"); return 1; }

    //Get sc (normalizer)
    for (n=2; n<=N; n++) { sc += n*n; }
    sc = 0.5f/sc;

    if (dim==0)
    {
        if (iscolmajor)
        {
            cblas_scopy(2*No,&z,0,&X[No],1);
            for (c=0; c<C/3; c++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[c*R],0,&X[No+c*R],1);             //beg edge samps
                    cblas_saxpy(R-n,-sc*n,&X[c*R],1,&X[No+c*R+n],1);         //past samps
                    cblas_saxpy(R-n,sc*n,&X[c*R+n],1,&X[No+c*R],1);          //future samps
                    cblas_saxpy(n,sc*n,&X[c*R+R-1],0,&X[No+c*R+R-n],1);      //end edge samps
                }
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[No+c*R],0,&X[2*No+c*R],1);        //beg edge samps
                    cblas_saxpy(R-n,-sc*n,&X[No+c*R],1,&X[2*No+c*R+n],1);    //past samps
                    cblas_saxpy(R-n,sc*n,&X[No+c*R+n],1,&X[2*No+c*R],1);     //future samps
                    cblas_saxpy(n,sc*n,&X[No+c*R+R-1],0,&X[2*No+c*R+R-n],1); //end edge samps
                }
            }
        }
        else
        {
            for (c=0; c<C/3; c++)
            {
                cblas_scopy(R,&z,0,&X[C/3+c],C);
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[c],0,&X[C/3+c],C);                      //beg edge samps
                    cblas_saxpy(R-n,-sc*n,&X[c],C,&X[C/3+c+n*C],C);                //past samps
                    cblas_saxpy(R-n,sc*n,&X[c+n*C],C,&X[C/3+c],C);                 //future samps
                    cblas_saxpy(n,sc*n,&X[c+C*(R-1)],0,&X[C/3+c+C*(R-n)],C);       //end edge samps
                }
                cblas_scopy(R,&z,0,&X[2*C/3+c],C);
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[C/3+c],0,&X[2*C/3+c],C);                //beg edge samps
                    cblas_saxpy(R-n,-sc*n,&X[C/3+c],C,&X[2*C/3+c+n*C],C);          //past samps
                    cblas_saxpy(R-n,sc*n,&X[C/3+c+n*C],C,&X[2*C/3+c],C);           //future samps
                    cblas_saxpy(n,sc*n,&X[C/3+c+C*(R-1)],0,&X[2*C/3+c+C*(R-n)],C); //end edge samps
                }
            }
        }
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            for (r=0; r<R/3; r++)
            {
                cblas_scopy(C,&z,0,&X[R/3+r],R);
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[r],0,&X[R/3+r],R);                      //beg edge samps
                    cblas_saxpy(C-n,-sc*n,&X[r],R,&X[R/3+r+n*R],R);                //past samps
                    cblas_saxpy(C-n,sc*n,&X[r+n*R],R,&X[R/3+r],R);                 //future samps
                    cblas_saxpy(n,sc*n,&X[r+R*(C-1)],0,&X[R/3+r+R*(C-n)],R);       //end edge samps
                }
                cblas_scopy(C,&z,0,&X[2*R/3+r],R);
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[R/3+r],0,&X[2*R/3+r],R);                //beg edge samps
                    cblas_saxpy(C-n,-sc*n,&X[R/3+r],R,&X[2*R/3+r+n*R],R);          //past samps
                    cblas_saxpy(C-n,sc*n,&X[R/3+r+n*R],R,&X[2*R/3+r],R);           //future samps
                    cblas_saxpy(n,sc*n,&X[R/3+r+R*(C-1)],0,&X[2*R/3+r+R*(C-n)],R); //end edge samps
                }
            }
        }
        else
        {
            cblas_scopy(2*No,&z,0,&X[No],1);
            for (r=0; r<R/3; r++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[r*C],0,&X[No+r*C],1);             //beg edge samps
                    cblas_saxpy(C-n,-sc*n,&X[r*C],1,&X[No+r*C+n],1);         //past samps
                    cblas_saxpy(C-n,sc*n,&X[r*C+n],1,&X[No+r*C],1);          //future samps
                    cblas_saxpy(n,sc*n,&X[r*C+C-1],0,&X[No+r*C+C-n],1);      //end edge samps
                }
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[No+r*C],0,&X[2*No+r*C],1);        //beg edge samps
                    cblas_saxpy(C-n,-sc*n,&X[No+r*C],1,&X[2*No+r*C+n],1);    //past samps
                    cblas_saxpy(C-n,sc*n,&X[No+r*C+n],1,&X[2*No+r*C],1);     //future samps
                    cblas_saxpy(n,sc*n,&X[No+r*C+C-1],0,&X[2*No+r*C+C-n],1); //end edge samps
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in add_delta_deltas_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int add_delta_deltas_d (double *X, const int iscolmajor, const int R, const int C, const int dim, const int N)
{
    const double z = 0.0;
    const int No = R*C/3;
    int r, c, n;
    double sc = 1.0;

    //Checks
    if (R<1) { fprintf(stderr,"error in add_delta_deltas_d: R (nrows X) must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in add_delta_deltas_d: C (ncols X) must be positive\n"); return 1; }
    if (N<1) { fprintf(stderr,"error in add_delta_deltas_d: N (delta winlength) must be positive\n"); return 1; }
    if (dim==0 && C%3!=0) { fprintf(stderr,"error in add_delta_deltas_d: C (ncols X) must be divisible by 3\n"); return 1; }
    if (dim==1 && R%3!=0) { fprintf(stderr,"error in add_delta_deltas_d: R (nrows X) must be divisible by 3\n"); return 1; }

    //Get sc (normalizer)
    for (n=2; n<=N; n++) { sc += n*n; }
    sc = 0.5/sc;

    if (dim==0)
    {
        if (iscolmajor)
        {
            cblas_dcopy(2*No,&z,0,&X[No],1);
            for (c=0; c<C/3; c++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[c*R],0,&X[No+c*R],1);             //beg edge samps
                    cblas_daxpy(R-n,-sc*n,&X[c*R],1,&X[No+c*R+n],1);         //past samps
                    cblas_daxpy(R-n,sc*n,&X[c*R+n],1,&X[No+c*R],1);          //future samps
                    cblas_daxpy(n,sc*n,&X[c*R+R-1],0,&X[No+c*R+R-n],1);      //end edge samps
                }
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[No+c*R],0,&X[2*No+c*R],1);        //beg edge samps
                    cblas_daxpy(R-n,-sc*n,&X[No+c*R],1,&X[2*No+c*R+n],1);    //past samps
                    cblas_daxpy(R-n,sc*n,&X[No+c*R+n],1,&X[2*No+c*R],1);     //future samps
                    cblas_daxpy(n,sc*n,&X[No+c*R+R-1],0,&X[2*No+c*R+R-n],1); //end edge samps
                }
            }
        }
        else
        {
            for (c=0; c<C/3; c++)
            {
                cblas_dcopy(R,&z,0,&X[C/3+c],C);
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[c],0,&X[C/3+c],C);                      //beg edge samps
                    cblas_daxpy(R-n,-sc*n,&X[c],C,&X[C/3+c+n*C],C);                //past samps
                    cblas_daxpy(R-n,sc*n,&X[c+n*C],C,&X[C/3+c],C);                 //future samps
                    cblas_daxpy(n,sc*n,&X[c+C*(R-1)],0,&X[C/3+c+C*(R-n)],C);       //end edge samps
                }
                cblas_dcopy(R,&z,0,&X[2*C/3+c],C);
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[C/3+c],0,&X[2*C/3+c],C);                //beg edge samps
                    cblas_daxpy(R-n,-sc*n,&X[C/3+c],C,&X[2*C/3+c+n*C],C);          //past samps
                    cblas_daxpy(R-n,sc*n,&X[C/3+c+n*C],C,&X[2*C/3+c],C);           //future samps
                    cblas_daxpy(n,sc*n,&X[C/3+c+C*(R-1)],0,&X[2*C/3+c+C*(R-n)],C); //end edge samps
                }
            }
        }
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            for (r=0; r<R/3; r++)
            {
                cblas_dcopy(C,&z,0,&X[R/3+r],R);
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[r],0,&X[R/3+r],R);                      //beg edge samps
                    cblas_daxpy(C-n,-sc*n,&X[r],R,&X[R/3+r+n*R],R);                //past samps
                    cblas_daxpy(C-n,sc*n,&X[r+n*R],R,&X[R/3+r],R);                 //future samps
                    cblas_daxpy(n,sc*n,&X[r+R*(C-1)],0,&X[R/3+r+R*(C-n)],R);       //end edge samps
                }
                cblas_dcopy(C,&z,0,&X[2*R/3+r],R);
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[R/3+r],0,&X[2*R/3+r],R);                //beg edge samps
                    cblas_daxpy(C-n,-sc*n,&X[R/3+r],R,&X[2*R/3+r+n*R],R);          //past samps
                    cblas_daxpy(C-n,sc*n,&X[R/3+r+n*R],R,&X[2*R/3+r],R);           //future samps
                    cblas_daxpy(n,sc*n,&X[R/3+r+R*(C-1)],0,&X[2*R/3+r+R*(C-n)],R); //end edge samps
                }
            }
        }
        else
        {
            cblas_dcopy(2*No,&z,0,&X[No],1);
            for (r=0; r<R/3; r++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[r*C],0,&X[No+r*C],1);             //beg edge samps
                    cblas_daxpy(C-n,-sc*n,&X[r*C],1,&X[No+r*C+n],1);         //past samps
                    cblas_daxpy(C-n,sc*n,&X[r*C+n],1,&X[No+r*C],1);          //future samps
                    cblas_daxpy(n,sc*n,&X[r*C+C-1],0,&X[No+r*C+C-n],1);      //end edge samps
                }
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[No+r*C],0,&X[2*No+r*C],1);        //beg edge samps
                    cblas_daxpy(C-n,-sc*n,&X[No+r*C],1,&X[2*No+r*C+n],1);    //past samps
                    cblas_daxpy(C-n,sc*n,&X[No+r*C+n],1,&X[2*No+r*C],1);     //future samps
                    cblas_daxpy(n,sc*n,&X[No+r*C+C-1],0,&X[2*No+r*C+C-n],1); //end edge samps
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in add_delta_deltas_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif

ext: c


text:
#pragma once

#include "Periodic.h"

#include <gsl/gsl>

#include <algorithm>
#include <vector>

namespace ad {

class Rythm
{
public:
    explicit Rythm(duration_t aPeriod, int aDuration) :
        mMetronom(aPeriod),
        mDuration(aDuration)
    {}

    Rythm & note(int aTime);

    template <class T_functor, class... VT_args>
    void forEachEvent(duration_t aDelta, const T_functor & aFunctor, VT_args &&... aArgs);

private:
    Periodic mMetronom;
    int mDuration;
    std::vector<int> mScore{};
    int mCurrentNote{0};
};

inline Rythm & Rythm::note(int aTime)
{
    Expects( (aTime >= 0) && (aTime < mDuration) );
    mScore.push_back(aTime);
    return *this;
}

template <class T_functor, class... VT_args>
void Rythm::forEachEvent(duration_t aDelta, const T_functor & aFunctor, VT_args &&... aArgs)
{
    mMetronom.forEachEvent(aDelta, [&aFunctor, this](duration_t aRemainingTime, VT_args &&... aArgs)
    {
        if (std::find(mScore.begin(), mScore.end(), mCurrentNote) != mScore.end())
        {
            aFunctor(aRemainingTime, std::forward<VT_args>(aArgs)...);
        }
        if (++mCurrentNote == mDuration)
        {
            mCurrentNote = 0;
        };
    }, std::forward<VT_args>(aArgs)...);
}

} // namespace ad

ext: h

max_stars_repo_name: FranzPoize/shmurp

max_issues_repo_name: FranzPoize/shmurp

max_forks_repo_name: FranzPoize/shmurp


text:
/**
* Drizzle utilities
*/

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_statistics.h>
#include <math.h>

#include "aXe_grism.h"
#include "aXe_utils.h"
#include "spc_trace_functions.h"
#include "crossdisp_utils.h"
#include "aper_conf.h"
#include "trace_conf.h"
#include "disp_conf.h"
#include "spc_wl_calib.h"
#include "drizzle_utils.h"

d_point
get_refwave_position(dispstruct * disp, trace_func * trace, d_point refpix,
                     aperture_conf *conf)
{
  double a=0.0;
  double b, c, c_star;
  double *cf;
  double rot;
  double dx, dy, dtr;
  d_point res;

  // check for a dispersion solution
  // higher that quadratic order
  if (disp->pol->size > 3)
    // this can not be solved; give an error
    aXe_message(aXe_M_FATAL, __FILE__, __LINE__,
                "Order of dispersion solution: %i!\n"
                "At most qudratic solutions are allowed!\n",
                disp->pol->size-1);
  else if (disp->pol->size > 2)
    // store the quadratic term, if it exists
    a = gsl_vector_get(disp->pol, 2);  // ddlambda-term
  else
    // set the qudratic term to zero
    a=0.0;

  // store the constant and linear terms
  b = gsl_vector_get(disp->pol, 1);  // dlambda-term
  c = gsl_vector_get(disp->pol, 0);  // lambda0-term

  // store the dydx value
  cf = trace->data;

  // dydx -->  rotation angle
  rot = get_rotation_angle(cf[2]);

  // determine the wavelength difference
  c_star = c - conf->drz_lamb0;

  // transform the wavelength difference
  // into a path length difference
  if (a)
    dtr = (-1.0*b+sqrt(b*b - 4.0*a*c_star))/(2.0*a);
  else
    dtr = c_star/b;

  // path length difference --> dx, dy
  dx = dtr*cos(rot);
  dy = dtr*sin(rot);

  // compute the absolute values
  // for the reference position
  // in the stamp images
  res.x = refpix.x + dx;
  res.y = refpix.y + dy;

  // return the result
  return res;
}

double
get_drizzle_width(object *ob, int beamnum,trace_func * trace)
{
  double drizzle_width, drizzle_orient;
  double orig_width, orig_orient;
  double rotation, angle, oangle, factor;
  double *cf;
  beam *b = ob->beams+beamnum;

  cf = trace->data;
  rotation = get_rotation_angle(cf[2]);

  orig_orient = b->orient;
  orig_width  = b->width;

  drizzle_orient = orig_orient-rotation;

  angle  = drizzle_orient / M_PI * 180.;
  oangle = orig_orient / M_PI * 180.;

  drizzle_width = fabs(orig_width * sin (drizzle_orient));
  factor = drizzle_width/orig_width;

  if (angle > 180.0)
    angle = angle-180.0;
  if (angle < 30.0 || angle > 150.0)
     fprintf (stdout,
              "Angle: %4.0f --> %4.0f, Width: %5.1f --> %5.1f, %5.3f\n",
              oangle, angle, orig_width, drizzle_width, factor);

  return drizzle_width;
}

gsl_matrix  *
get_drizzle_coeffs(dispstruct * disp, trace_func * trace,
                   int boxwidth, int boxheight, int trlength,
                   double relx, double rely, aperture_conf *conf,
                   double orient, dispstruct * outdisp, double cdcorr,
                   double sprefreso, double spmeanreso)
{
  double a, b, c;
  double ao, bo, co, bref, cref;
  double a11, a12, a21, a22;
  double xr, yr;
  //double lambda0;
  //double dlambda;
  double tmp, rotation;
  double *cf;
  double shear_term;

  gsl_matrix * ret = gsl_matrix_alloc(2,11);
  gsl_matrix * rotcoeffs;

  gsl_matrix_set_all (ret, 0.0);

  /* get the dispersion at the objects point */
  if (disp->pol->size > 2)
    a = gsl_vector_get(disp->pol, 2);
    //    fprintf (stdout, "quadratic solution:  %f\n", a);
  else
    a = 0.0;

  b = gsl_vector_get(disp->pol, 1);  // dlambda-term
  c = gsl_vector_get(disp->pol, 0);  // lambda0-term

  /* get the dispersion at the mean reference point */
  if (outdisp->pol->size > 2){
    ao = gsl_vector_get(outdisp->pol, 2);
    //    fprintf (stdout, "quadratic solution:  %f, old: %f\n", a, ao);
  }
  else{
    ao = 0.0;
  }
  bo = gsl_vector_get(outdisp->pol, 1);  // dlambda-term
  co = gsl_vector_get(outdisp->pol, 0);  // lambda0-term
  if (conf->drz_resol == 0.0) {
    conf->drz_resol = bo + ao*((double)trlength)/2.0;
  }
  if (conf->drz_lamb0 == 0.0) {
    conf->drz_lamb0 = co;
  }

  //  fprintf (stdout, "lambda terms:  4785.0 <-> %f, 24.0 <-> %f\n", cref, bref);

  cref = conf->drz_lamb0;
  bref = conf->drz_resol;
  bref = sprefreso;

  // Bugfix on Sept. 15th 2010:
  // Check whether the ratio of the average dispersion on the grism images
  // and the dispersion on the drizzled images is outsidethe range 0.95 < ratio < 1.05.
  // If yes, adjust "trlength" which, from now on, is the length of the
  // drizzled image. This fix is necessary to allow a different sampling
  // in axedrizzle. The range was introduced such that the bug-fixed
  // version delivers identical results when using the default sampling.
  if ((spmeanreso / sprefreso > 1.05) || (spmeanreso / sprefreso < 0.95))
    trlength = (int)(spmeanreso / sprefreso * (double)trlength) + 0.5;

  /* put the below lines to go BACK to
   * the old computation of the length of
   * the aXedrizzled images:
  fprintf(stdout, "New spectral length: %i", trlength);
  trlength = (int)ceil(bo/bref*sqrt(pow((double)boxwidth,2.0)
                                    +pow((double)boxheight,2.0)));
  trlength = (int)ceil(spmeanreso/sprefreso*sqrt(pow((double)boxwidth,2.0)
                                                 +pow((double)boxheight,2.0)));
  fprintf(stdout, "<--> old spectral length: %i\n", trlength);
  */

  //fprintf (stdout,"Drizzled to resolution: %e, mean resolution for object: %e\n", sprefreso, spmeanreso);

  gsl_matrix_set(ret, 0,10,trlength);
  //  gsl_matrix_set(ret, 1,10,bref);

  cf = trace->data;
  rotation = get_rotation_angle(cf[2]);
  rotcoeffs = get_coeffs_back(trace);
  a11 = gsl_matrix_get(rotcoeffs,0,0);
  a12 = gsl_matrix_get(rotcoeffs,0,1);
  a21 = gsl_matrix_get(rotcoeffs,1,0);
  a22 = gsl_matrix_get(rotcoeffs,1,1);

  // the following lines have to be changed
  // in order to go from a integer
  // center definition to floating point center
  // definition.
  xr = (relx+1.0)-((double)(boxwidth/2)+1);  // +1.0 to compensate for "sp_sex.c" around line 537
  yr = (rely+1.0)-((double)(boxheight/2)+1); // +1.0 to compensate for "sp_sex.c" around line 537

  if (b < 0.0)
    shear_term = -tan(orient-rotation);
  else
    shear_term = tan(orient-rotation);

  // Transformations for X, without the lambda-terms, only rotation:
  //  tmp = xoffs-(((double)trlength/2.0)+1)-a11*xr-a12*yr;
  //  gsl_matrix_set(ret, 0,0,tmp);  // constant term

  //  tmp = a11;
  //  gsl_matrix_set(ret, 0,1,tmp);  // x-term

  //  tmp = a12;
  //  gsl_matrix_set(ret, 0,2,tmp);  // y-term
  //----------------------------------------------------------------------------------

  //------------------------------------------------------------------
  // Mathematica based run
  tmp = 1.0/bref*(c - cref - a11*b*xr + a*a11*a11*xr*xr - a12*b*yr + 2*a*a11*a12*xr*yr + a*a12*a12*yr*yr) +
    (a21*xr + a22*yr)/shear_term - ((double)(trlength/2)+1.0) + conf->drz_xstart;
  gsl_matrix_set(ret, 0,0,tmp);  // constant term

  tmp = 1.0/bref*(a11*b - 2*a*a11*a11*xr - 2*a*a11*a12*yr) - a21/shear_term;
  gsl_matrix_set(ret, 0,1,tmp);  // x-term

  tmp = 1.0/bref*(a12*b - 2*a*a12*a12*yr - 2*a*a11*a12*xr) - a22/shear_term;
  gsl_matrix_set(ret, 0,2,tmp);  // y-term

  tmp = a*a11*a11/bref;
  gsl_matrix_set(ret, 0,3,tmp);  // x^2-term

  tmp = 2*a*a11*a12/bref;
  gsl_matrix_set(ret, 0,4,tmp); ; // xy-term

  tmp = a*a12*a12/bref;
  gsl_matrix_set(ret, 0,5,tmp);  // y^2-term

  //------------------------------------------------------------------

  // Transformations for Y:
  tmp = -1.0*(a21*xr+a22*yr) * cdcorr; // -1.5  is some kind of a fudge factor, no idea where it comes from
  gsl_matrix_set(ret, 1,0,tmp);  // constant term

  tmp = a21 * cdcorr;
  gsl_matrix_set(ret, 1,1,tmp);  // x-term

  tmp = a22 * cdcorr;
  gsl_matrix_set(ret, 1,2,tmp);  // y-term

  gsl_matrix_free(rotcoeffs);

  return ret;
}

gsl_matrix *
get_coeffs_back(trace_func * trace)
{
  double *cf;
  double rotation;
  gsl_matrix * ret = gsl_matrix_alloc(2,2);

  //  ret = gsl_matrix_alloc (2,2);
  gsl_matrix_set_all (ret, 0.0);

  cf = trace->data;

  rotation = get_rotation_angle(cf[2]);

  gsl_matrix_set(ret, 0, 0,      cos(rotation));
  gsl_matrix_set(ret, 0, 1,      sin(rotation));
  gsl_matrix_set(ret, 1, 0, -1.0*sin(rotation));
  gsl_matrix_set(ret, 1, 1,      cos(rotation));

  return ret;
}

double
get_rotation_angle(double dxdy){
  double rotation=0.0;

  rotation = atan2(dxdy, 1.0);

  return rotation;
}

// Start functions for an alternative approach to store the object information
objectobs **
malloc_objectobs(){
  objectobs **allobjects;

  allobjects = (objectobs **) malloc (NMAXOBJ*sizeof(objectobs *));

  return allobjects;
}

void
free_objectobs(objectobs **allobjects){
  //  int i, nobjs=0;
  //  while (allobjects[nobjs])
  //    nobjs++;
  //  for (i=0; i<nobjs; i++){
  //      free(allobjects[i]);
  //      //    allobjects[i] = NULL;
  //  }
  free(allobjects);
}

int
add_observation(char * filename, char *  conf_file, objectobs ** allobjects, int nobjects,
                object ** oblist, int list_size, px_point pixmax, int sci_numext)
{

  int i, j, dec;
  int beamID = 0;
  int get_scale=0;
  int tlength;
  gsl_matrix * drzcoeffs;
  double cdscale;
  double drzscale;
  dispstruct * disp;
  calib_function *wl_calib;
  double l1, l2, spreso;

  drzcoeffs = get_crossdisp_matrix(filename, sci_numext);

  if (drzcoeffs->size1 > 1 && drzcoeffs->size2)
    get_scale = 1;

  if (get_scale){
    drzscale =  (double)get_float_from_keyword(filename, 1, "DRZSCALE");
  }
  for (i=0; i < list_size; i++){
    if (oblist[i]->beams[beamID].ignore == 0){
      dec = 0;

      if (get_scale)
        cdscale = drzscale * get_crossdisp_scale(oblist[i]->beams[beamID].spec_trace,
                                                 oblist[i]->beams[beamID].refpoint, drzcoeffs, pixmax);
      else
        cdscale=1.0;

      disp = get_dispstruct_at_pos(conf_file, 1, beamID,oblist[i]->beams[beamID].refpoint);
      wl_calib = create_calib_from_gsl_vector(1, disp->pol);
      l1 =  wl_calib->func (0.0,  wl_calib->order,  wl_calib->coeffs);
      l2 =  wl_calib->func (1.0,  wl_calib->order,  wl_calib->coeffs);
      spreso = fabs(l2-l1);

      // compute the tracelength
      tlength = get_beam_trace_length(oblist[i]->beams[beamID]);

      for (j=0; j < nobjects; j++){
        if (allobjects[j]->OBJID == oblist[i]->ID){
          add_obs_to_allobj(allobjects[j], oblist[i], pixmax, cdscale, spreso, tlength);
          dec = 1;
        }
      }
      if (dec == 0) {
        nobjects = add_obj_to_allobj(allobjects, nobjects, oblist[i], pixmax, cdscale, spreso, tlength);
      }
      free_dispstruct(disp);
      free_calib(wl_calib);
    }
  }
  gsl_matrix_free(drzcoeffs);
  return nobjects;
}

void
add_obs_to_allobj(objectobs *actobject, object * actobs, px_point pixmax, double cdscale, double spreso, int tlength)
{

  int xmin, xmax, ymin, ymax;
  //double m, b;
  gsl_vector_int * xvec;
  gsl_vector_int * yvec;
  px_point bbox;
  px_point mins;

  double *gaga;

  xvec = gsl_vector_int_alloc (4);
  yvec = gsl_vector_int_alloc (4);

  /* Store the refpoint */
  actobject->refpoint[actobject->nobs].x = actobs->beams[0].refpoint.x;
  actobject->refpoint[actobject->nobs].y = actobs->beams[0].refpoint.y;

  //*************************************************
  // patch to correct the reference point in case
  // that the the trace descritpion
  // does have a non negligeable first order term!
  gaga = actobs->beams[0].spec_trace->data;
  actobject->refpoint[actobject->nobs].y = actobject->refpoint[actobject->nobs].y + gaga[1];
  //*************************************************

  gsl_vector_int_set(xvec, 0, actobs->beams[0].corners[0].x);
  gsl_vector_int_set(xvec, 1, actobs->beams[0].corners[1].x);
  gsl_vector_int_set(xvec, 2, actobs->beams[0].corners[2].x);
  gsl_vector_int_set(xvec, 3, actobs->beams[0].corners[3].x);
  xmin = gsl_stats_int_min(xvec->data, 1, 4);
  xmax = gsl_stats_int_max(xvec->data, 1, 4);
  gsl_vector_int_set(yvec, 0, actobs->beams[0].corners[0].y);
  gsl_vector_int_set(yvec, 1, actobs->beams[0].corners[1].y);
  gsl_vector_int_set(yvec, 2, actobs->beams[0].corners[2].y);
  gsl_vector_int_set(yvec, 3, actobs->beams[0].corners[3].y);
  ymin = gsl_stats_int_min(yvec->data, 1, 4);
  ymax = gsl_stats_int_max(yvec->data, 1, 4);

  bbox = recalc_bbox(xmin, xmax, ymin, ymax, pixmax);
  mins = recalc_mins(xmin, xmax, ymin, ymax, pixmax);

  actobject->width[actobject->nobs]    = bbox.x;
  actobject->height[actobject->nobs]   = bbox.y;
  actobject->tlength[actobject->nobs]  = tlength;
  actobject->objwidth[actobject->nobs] = actobs->beams[0].width;
  actobject->orient[actobject->nobs]   = actobs->beams[0].orient;
  actobject->cdscale[actobject->nobs]   = cdscale;
  actobject->spreso[actobject->nobs]   = spreso;

  actobject->relrefpt[actobject->nobs].x = actobs->beams[0].refpoint.x - (double)mins.x;
  actobject->relrefpt[actobject->nobs].y = actobs->beams[0].refpoint.y - (double)mins.y;

  actobject->nobs = actobject->nobs + 1;

  gsl_vector_int_free(xvec);
  gsl_vector_int_free(yvec);
}

int
add_obj_to_allobj(objectobs **allobjects, int nobjects, object * actobs,
                   px_point pixmax, double cdscale, double spreso, int tlength)
{
  int xmin, xmax, ymin, ymax;
  //double m, b;

  gsl_vector_int * xvec;
  gsl_vector_int * yvec;

  px_point bbox;
  px_point mins;

  double *gaga;

  //  object *ob = malloc (sizeof (object));
  objectobs *objobs = malloc (sizeof (objectobs));

  xvec = gsl_vector_int_alloc (4);
  yvec = gsl_vector_int_alloc (4);

  allobjects[nobjects] = objobs;
  /* Store the refpoint */
  allobjects[nobjects]->refpoint[0].x = actobs->beams[0].refpoint.x;
  allobjects[nobjects]->refpoint[0].y = actobs->beams[0].refpoint.y;

  //*************************************************
  // patch to correct the reference point in case
  // that the the trace descritpion
  // does have a non negligeable first order term!
  gaga = actobs->beams[0].spec_trace->data;
  allobjects[nobjects]->refpoint[0].y = allobjects[nobjects]->refpoint[0].y + gaga[1];
  //*************************************************

  allobjects[nobjects]->OBJID = actobs->ID;
  allobjects[nobjects]->nobs = 1;
  allobjects[nobjects]->pointer = 0;

  gsl_vector_int_set(xvec, 0, actobs->beams[0].corners[0].x);
  gsl_vector_int_set(xvec, 1, actobs->beams[0].corners[1].x);
  gsl_vector_int_set(xvec, 2, actobs->beams[0].corners[2].x);
  gsl_vector_int_set(xvec, 3, actobs->beams[0].corners[3].x);
  xmin = gsl_stats_int_min(xvec->data, 1, 4);
  xmax = gsl_stats_int_max(xvec->data, 1, 4);

  gsl_vector_int_set(yvec, 0, actobs->beams[0].corners[0].y);
  gsl_vector_int_set(yvec, 1, actobs->beams[0].corners[1].y);
  gsl_vector_int_set(yvec, 2, actobs->beams[0].corners[2].y);
  gsl_vector_int_set(yvec, 3, actobs->beams[0].corners[3].y);
  ymin = gsl_stats_int_min(yvec->data, 1, 4);
  ymax = gsl_stats_int_max(yvec->data, 1, 4);

  bbox = recalc_bbox(xmin, xmax, ymin, ymax, pixmax);
  mins = recalc_mins(xmin, xmax, ymin, ymax, pixmax);

  allobjects[nobjects]->width[0]   = bbox.x;
  allobjects[nobjects]->height[0]  = bbox.y;
  allobjects[nobjects]->tlength[0] = tlength;
  allobjects[nobjects]->objwidth[0] = actobs->beams[0].width;
  allobjects[nobjects]->orient[0]   = actobs->beams[0].orient;
  allobjects[nobjects]->cdscale[0]  = cdscale;
  allobjects[nobjects]->spreso[0]   = spreso;

  allobjects[nobjects]->relrefpt[0].x = actobs->beams[0].refpoint.x - (double)mins.x;
  allobjects[nobjects]->relrefpt[0].y = actobs->beams[0].refpoint.y - (double)mins.y;

  ++nobjects;

  gsl_vector_int_free(xvec);
  gsl_vector_int_free(yvec);

  return nobjects;
}

void print_objectobs(objectobs **allobjects, int nobjects){
  int i, j;

  for (i = 0; i < nobjects; i++){
   fprintf(stdout, "OBJID: %i, no of observ.: %i\n", allobjects[i]->OBJID, allobjects[i]->nobs);
    for (j = 0; j < allobjects[i]->nobs; j++){
      fprintf(stdout, "No %i: width %f, height: %f, \n", j, allobjects[i]->relrefpt[j].x, allobjects[i]->relrefpt[j].y);
    }
  }
}
void print_objectobs2(objectobs allobjects[], int nobjects){
  int i, j;

  for (i = 0; i < nobjects; i++){
   fprintf(stdout, "OBJID: %i, no of observ.: %i\n", allobjects[i].OBJID, allobjects[i].nobs);
    for (j = 0; j < allobjects[i].nobs; j++){
      fprintf(stdout, "No %i: width %f, height: %f, \n", j, allobjects[i].relrefpt[j].x, allobjects[i].relrefpt[j].y);
    }
  }
}
int make_refpoints( char * conf_file, char * filename, px_point pixmax, objectobs **allobjects, int nobjects)
{
  int i, j;
  int nobs;
  int max_width, max_height, omax;
  double xmean, ymean;
  //double xdata[NMAXOBS];
  //double ydata[NMAXOBS];
  //  objectobs  oneobject;
  objectobs *oneobject; // = malloc (sizeof (objectobs));

  gsl_vector * xvec;
  gsl_vector * yvec;
  aperture_conf *conf;
  gsl_matrix * drzcoeffs;
  //double cdscale;
  double drzscale;
  trace_func  *trace;
  int beamID=0;
  dispstruct * disp;
  calib_function *wl_calib;
  double l1, l2;

  conf = get_aperture_descriptor(conf_file);
  get_extension_numbers(filename, conf,conf->optkey1,conf->optval1);
  drzcoeffs = get_crossdisp_matrix(filename,conf->science_numext);
  drzscale =  (double)get_float_from_keyword(filename, 1, "DRZSCALE");

  for (i=0; i < nobjects; i++){
    oneobject = allobjects[i];
    nobs = oneobject->nobs;

    xvec = gsl_vector_alloc (nobs);
    yvec = gsl_vector_alloc (nobs);
    for (j=0; j < nobs; j++){
      gsl_vector_set(xvec, j, oneobject->refpoint[j].x);
      gsl_vector_set(yvec, j, oneobject->refpoint[j].y);
    }
    xmean = gsl_stats_mean(xvec->data, 1, nobs);
    ymean = gsl_stats_mean(yvec->data, 1, nobs);

    max_width  = gsl_stats_int_max(oneobject->width,  1, nobs);
    max_height = gsl_stats_int_max(oneobject->height, 1, nobs);
    omax       = gsl_stats_max(oneobject->objwidth, 1, nobs);

    allobjects[i]->mean_refpoint.x = xmean;
    allobjects[i]->mean_refpoint.y = ymean;
    allobjects[i]->max_width  = max_width;
    allobjects[i]->max_height = max_height;
    allobjects[i]->owidthmax  = omax;

    allobjects[i]->max_tlength = gsl_stats_int_max(oneobject->tlength, 1, nobs);

    /*
     * Look whether the cross dispersion scale is given in the configuration.
     * If not determine the cross dispersion scale at the mean reference point an store it
     * store also the mean correction in for to correct the width.
     *
     */
    if (conf->drz_scale < 1.0e-16){
      trace = get_tracefunc_at(conf_file, allobjects[i]->mean_refpoint);
      allobjects[i]->cdrefscale = drzscale * get_crossdisp_scale(trace, allobjects[i]->mean_refpoint, drzcoeffs, pixmax);
    }
    else{
      allobjects[i]->cdrefscale = conf->drz_scale;
    }
    for (j=0; j < nobs; j++){
      gsl_vector_set(xvec, j, oneobject->cdscale[j]);
    }
    xmean = gsl_stats_mean(xvec->data, 1, nobs);
    allobjects[i]->cdmeanscale  = xmean;

    /*
     * look whether the wavelength dispersion is given in the conig file
     * if not determine the wavelength dispersion at the mean reference point amd store
     * it. Store also the mean wavelength dispersion to correct the length
     */
    if (conf->drz_resol < 1.0e-16){
      disp = get_dispstruct_at_pos(conf_file, 1, beamID, allobjects[i]->mean_refpoint);
      wl_calib = create_calib_from_gsl_vector(1, disp->pol);
      l1 =  wl_calib->func (0.0,  wl_calib->order,  wl_calib->coeffs);
      l2 =  wl_calib->func (1.0,  wl_calib->order,  wl_calib->coeffs);
      allobjects[i]->sprefreso = fabs(l2-l1);
    }
    else{
       allobjects[i]->sprefreso = conf->drz_resol;
    }
    for (j=0; j < nobs; j++){
      gsl_vector_set(xvec, j, oneobject->spreso[j]);
    }
    xmean = gsl_stats_mean(xvec->data, 1, nobs);
    allobjects[i]->spmeanreso  = xmean;

    gsl_vector_free(xvec);
    gsl_vector_free(yvec);
  }

  gsl_matrix_free(drzcoeffs);
  free(conf);

  return 0;
}

d_point
get_mean_refpoint(objectobs **allobjects, int nobjects, int ID,int * boxwidth,
                  int * boxheight, double * relx, double * rely,
                  double * objwidth, double *orient, double * cdref,
                  double *cdscale, double *cdmeanscale, double *sprefreso, double *spreso,
                  double *spmeanreso, int *tlength)
{
  int i, dec=0;
  d_point mpoint;

  for (i=0; i < nobjects; i++){
    if (allobjects[i]->OBJID == ID){
      mpoint.x = allobjects[i]->mean_refpoint.x;
      mpoint.y = allobjects[i]->mean_refpoint.y;
      *boxwidth  = allobjects[i]->max_width;
      *boxheight = allobjects[i]->max_height;
      *tlength   = allobjects[i]->max_tlength;
      *objwidth  = allobjects[i]->owidthmax;
      *cdref     = allobjects[i]->cdrefscale;
      *cdmeanscale     = allobjects[i]->cdmeanscale;
      *spmeanreso= allobjects[i]->spmeanreso;
      *sprefreso = allobjects[i]->sprefreso;
      *orient    = allobjects[i]->orient[allobjects[i]->pointer];
      *cdscale   = allobjects[i]->cdscale[allobjects[i]->pointer];
      *spreso    = allobjects[i]->spreso[allobjects[i]->pointer];
      *relx      = allobjects[i]->relrefpt[allobjects[i]->pointer].x;
      *rely      = allobjects[i]->relrefpt[allobjects[i]->pointer].y;
      ++allobjects[i]->pointer;
      dec=1;
    }
  }
  return mpoint;
}
px_point
recalc_bbox(int xmin, int xmax, int ymin, int ymax, px_point pixmax)
{
  px_point ret;

  if (xmax > pixmax.x)
    xmax = pixmax.x;
  if (ymax > pixmax.y)
    ymax = pixmax.y;
  if (xmin < 1)
    xmin = 1;
  if (ymin < 1)
    ymin = 1;

  ret.x = xmax - xmin + 2;
  ret.y = ymax - ymin + 2;

  return ret;
}
/*
I am not sure what this routine is doing,
but one should make it better.
*/
px_point
recalc_mins(int xmin, int xmax, int ymin, int ymax, px_point pixmax)
{
  double m, b;
  px_point ret;

  m = ((double)ymin - (double)ymax)/((double)xmax - (double)xmin);
  b = (double)ymin - m * (double)xmax;
  if (xmin < 1){
    ymax = (int)ceil(1.0*m + b);
    xmin=1;
  }
  if (xmax > pixmax.x){
    ymin = (int)ceil((float)pixmax.x*m + b);
    xmax = pixmax.x;
  }
  if (ymin < 1){
    xmax = (int)ceil((1-b)/m);
    ymin=1;
  }
  if (ymax > pixmax.y){
    xmin = (int)ceil(((float)pixmax.y-b)/m);
    ymax=pixmax.y;
  }
  ret.x = xmin;
  ret.y = ymin;

  return ret;
}

int
get_beam_trace_length(const beam actbeam)
{
  int trlength=0;
  double diagonal_1 = 0.0;
  double diagonal_2 = 0.0;
  double xdiff, ydiff;

  // compute the x- and y- differences for one corner pair
  xdiff = (float)actbeam.corners[0].x - (float)actbeam.corners[2].x;
  ydiff = (float)actbeam.corners[0].y - (float)actbeam.corners[2].y;

  // compute the diagonal distance
  diagonal_1 = xdiff*xdiff + ydiff*ydiff;

  // compute the x- and y- differences for the other corner pair
  xdiff = (float)actbeam.corners[1].x - (float)actbeam.corners[3].x;
  ydiff = (float)actbeam.corners[1].y - (float)actbeam.corners[3].y;

  // compute the diagonal distance
  diagonal_2 = xdiff*xdiff + ydiff*ydiff;

  // compute the tracelength
  // from the larger diagonal
  if (diagonal_1 > diagonal_2)
    trlength = (int)ceil(sqrt(diagonal_1)+3.0);
  else
    trlength = (int)ceil(sqrt(diagonal_2)+3.0);

  // return the result
  return trlength;
}


text:
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
int rosenbrock_f (const gsl_vector * x, void *params,
		  gsl_vector * f);
int print_state (size_t iter, gsl_multiroot_fsolver * s);

struct rparams
  {
    double a;
    double b;
  };

int rosenbrock_f (const gsl_vector * x, void *params,
              gsl_vector * f)
{
  double a = ((struct rparams *) params)->a;
  double b = ((struct rparams *) params)->b;

  const double x0 = gsl_vector_get (x, 0);
  const double x1 = gsl_vector_get (x, 1);

  const double y0 = a * (1 - x0);
  const double y1 = b * (x1 - x0 * x0);

  gsl_vector_set (f, 0, y0);
  gsl_vector_set (f, 1, y1);

  return GSL_SUCCESS;
}

int print_state (size_t iter, gsl_multiroot_fsolver * s)
{
  printf ("iter = %3u x = % .3f % .3f "
          "f(x) = % .3e % .3e\n",
          iter,
          gsl_vector_get (s->x, 0),
          gsl_vector_get (s->x, 1),
          gsl_vector_get (s->f, 0),
          gsl_vector_get (s->f, 1));
  return 0;
}

int main (void)
{
  const gsl_multiroot_fsolver_type *T;
  gsl_multiroot_fsolver *s;

  int status;
  size_t iter = 0;

  const size_t n = 2;
  struct rparams p = {1.0, 10.0};
  gsl_multiroot_function f = {&rosenbrock_f, n, &p};

  double x_init[2] = {-10.0, -5.0};
  gsl_vector *x = gsl_vector_alloc (n);

  gsl_vector_set (x, 0, x_init[0]);
  gsl_vector_set (x, 1, x_init[1]);

  T = gsl_multiroot_fsolver_hybrids;
  s = gsl_multiroot_fsolver_alloc (T, 2);
  gsl_multiroot_fsolver_set (s, &f, x);

  print_state (iter, s);

  do
    {
      iter++;
      status = gsl_multiroot_fsolver_iterate (s);

      print_state (iter, s);

      if (status)   /* check if solver is stuck */
        break;

      status =
        gsl_multiroot_test_residual (s->f, 1e-7);
    }
  while (status == GSL_CONTINUE && iter < 1000);

  printf ("status = %s\n", gsl_strerror (status));

  gsl_multiroot_fsolver_free (s);
  gsl_vector_free (x);
  return 0;
}

max_stars_repo_name: turtlewangbin/Meanfield

max_issues_repo_name: turtlewangbin/Meanfield

max_forks_repo_name: turtlewangbin/Meanfield


text:
/**
* @file likelihood.c
* @brief base functions to calculate the likelihood for general multitype OU processes on a phylogenetic tree
* @author Carl Boettiger <cboettig@gmail.com>
* @version 0.1
* @date 2011-04-22
*
* @section LICENSE
*
* Copyright (C)
* 2011 - Carl Boettiger
* 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.
*
* 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, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_errno.h>
#include "optimizers.h"
#include "mvn.h"

typedef struct {
	int n_nodes;
	size_t n_regimes;
	int * ancestor;
	int * regimes;
	double * branch_length;
	double * traits;
	double * alpha;
	double * theta;
	double * sigma;
	double * Xo;
	int * lca_matrix;
} tree;

int get_lca (int i, int j, int n_nodes, const int * ancestor,
             const double * branch_length, double * sep);

void simulate (const gsl_rng * rng, tree * mytree);

void calc_lik (const double *Xo, const double alpha[],
               const double theta[], const double sigma[],
               const int regimes[], const int ancestor[],
               const double branch_length[], const double traits[],
               int *n_nodes, int lca_matrix[], double *llik);

max_stars_repo_name: cboettig/wrightscape

max_issues_repo_name: cboettig/wrightscape

max_forks_repo_name: cboettig/wrightscape


text:
#pragma once

#include <gsl/gsl_vector.h>
#include <gsl/gsl_statistics_double.h>
#include <cmath>
#include <memory>
#include <cassert>

namespace halo {

// A simple wrapper for gsl_vector to manage its memory, and also some user-friendly additions.
class GSLVector {
public:
  GSLVector(size_t Sz, bool Zeroed=true)
    : Vec(gsl_vector_alloc(Sz), gsl_vector_free) {
      if (Zeroed)
        gsl_vector_set_all(Vec.get(), 0);
    }

  // number of elements the vector can handle
  size_t capacity() const { return Vec->size; }

  // The stride is the step-size from one element to the next in physical memory,
  // measured in units of the appropriate datatype.
  size_t stride() const { return Vec->stride; }

  // clears the vector, setting everything to zero
  void clear() {
    gsl_vector_set_all(Vec.get(), 0);
  }

  // access underlying vector, without taking ownership.
  gsl_vector* vec() { return Vec.get(); }
  gsl_vector const* vec() const { return Vec.get(); }

  // access the underlying vector's data
  double* data() { return Vec->data; }
  double const* data() const {return Vec->data; }

  // read an element from the vector
  double get(size_t idx) const { return gsl_vector_get(Vec.get(), idx); }

  // write an element of the vector
  void set(size_t idx, double val) { gsl_vector_set(Vec.get(), idx, val); }

private:
  // a managed gsl_vector, which must be constructed with an explicit deleter.
  std::unique_ptr<gsl_vector, decltype(&gsl_vector_free)> Vec;
};

/// A value whose quantity depends on some random phenomenon.
/// Not all observations are used to compute properties of the quantity,
/// with priority given to the most recent N observations, for some N.
class RandomQuantity {
public:
  RandomQuantity(size_t N=10) : Obs(N) {}

  /// Records an observation of the random quantity.
  /// If the history becomes full, the oldest observation is dropped to make room
  /// for this one.
  void observe(double NewSample) {
    Count += 1;
    Obs.set(NextFree, NewSample);

    const size_t Len = Obs.capacity();
    NextFree = (NextFree + 1) % Len; // bump with wrap-around

    // If this is the first observation, fill the vector with this value
    // if (Count == 1)
    //   for (size_t i = NextFree; i < Len; i++)
    //     Obs.set(i, NewSample);
  }

  // adds observations from another random quantity into this one.
  void merge(RandomQuantity const& RQ) {
    for (size_t i = 0; i < RQ.size(); i++)
      observe(RQ.Obs.get(i));
  }

  // resets the random quantity, forgetting all previous observations.
  void clear() {
    Obs.clear();
    NextFree = 0;
    Count = 0;
  }

  // returns the most recent observation
  double last() const {
    assert(Count > 0);
    size_t LastIdx = NextFree == 0 ? Obs.capacity()-1 : NextFree-1;
    return Obs.get(LastIdx);
  }

  /// returns the maximum capacity of this random quantity's ability to
  /// remember observations.
  size_t capacity() const {
    return Obs.capacity();
  }

  /// returns the number of observations currently remembered
  /// by this random quantity.
  size_t size() const { return std::min(Count, Obs.capacity()); }

  /// returns the mean of the RQ's observations.
  double mean() const {
    return gsl_stats_mean(Obs.data(), Obs.stride(), size());
  }

  /// Given the mean of this RQ's observations, returns the variance.
  double variance(double mean) const {
    return gsl_stats_variance_m(Obs.data(), Obs.stride(), size(), mean);
  }

  private:
    GSLVector Obs;
    size_t NextFree{0};
    size_t Count{0};
};

} // end namespace halo

max_stars_repo_name: halo-project/halo

max_issues_repo_name: halo-project/halo

max_forks_repo_name: halo-project/halo


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef rose_fbed1fbf_2ba3_46bc_97da_b1df11752358_h
#define rose_fbed1fbf_2ba3_46bc_97da_b1df11752358_h

#include <gslib/utility.h>
#include <ariel/config.h>
#include <ariel/rendersys.h>
#include <ariel/painter.h>
#include <ariel/batch.h>
#include <ariel/texbatch.h>

__ariel_begin__

class rose;

/* this block was used for synchronize with the constant buffer */
struct rose_configs
{
    vec4                screen;
    float               mapscreen[3][4];
    /* more.. */
};

struct rose_bind_info_cr
{
    vec4                color;      /* RGBA */
};

struct rose_bind_info_tex
{
    texture2d*          img;
    vec2                tex;
};

struct vertex_info_cr
{
    vec2                pos;
    vec4                cr;
};

struct vertex_info_klm_cr
{
    vec2                pos;
    vec3                klm;
    vec4                cr;
};

struct vertex_info_coef_cr
{
    vec2                pos;
    vec4                coef;       /* float3 coef & float tune */
    vec4                cr;
};

struct vertex_info_klm_tex
{
    vec2                pos;
    vec3                klm;
    vec2                tex;
};

struct vertex_info_coef_tex
{
    vec2                pos;
    vec4                coef;
    vec2                tex;
};

class rose_batch;
typedef list<rose_bind_info_cr> rose_bind_list_cr;
typedef list<rose_bind_info_tex> rose_bind_list_tex;
typedef vector<rose_batch*> rose_batch_list;
typedef vector<vertex_info_cr> vertex_stream_cr;
typedef vector<vertex_info_klm_cr> vertex_stream_klm_cr;
typedef vector<vertex_info_coef_cr> vertex_stream_cf_cr;
typedef vector<vertex_info_klm_tex> vertex_stream_klm_tex;
typedef vector<vertex_info_coef_tex> vertex_stream_cf_tex;
typedef vector<int> index_stream;
typedef bat_type rose_batch_tag;

class __gs_novtable rose_batch abstract
{
public:
    typedef rendersys::vertex_buffer vertex_buffer;
    typedef rendersys::index_buffer index_buffer;

public:
    rose_batch(int index);
    virtual ~rose_batch();
    virtual rose_batch_tag get_tag() const = 0;
    virtual void create(bat_batch* bat) = 0;
    virtual int buffering(rendersys* rsys) = 0;
    virtual void draw(rendersys* rsys) = 0;
    virtual void tracing() const = 0;

protected:
    int                 _bat_index;
    vertex_shader*      _vertex_shader;
    pixel_shader*       _pixel_shader;
    vertex_format*      _vertex_format;
    vertex_buffer*      _vertex_buffer;

public:
    int get_batch_index() const { return _bat_index; }
    void set_vertex_shader(vertex_shader* p) { _vertex_shader = p; }
    void set_pixel_shader(pixel_shader* p) { _pixel_shader = p; }
    void set_vertex_format(vertex_format* p) { _vertex_format = p; }
    void setup_vs_and_ps(rendersys* rsys);
    void setup_vf_and_topology(rendersys* rsys, uint topo);

protected:
    template<class stream_type>
    int template_buffering(stream_type& stm, rendersys* rsys);
};

class rose_fill_batch_cr:
    public rose_batch
{
public:
    rose_fill_batch_cr(int index): rose_batch(index) {}
    rose_batch_tag get_tag() const override { return bf_cr; }
    void create(bat_batch* bat) override;
    int buffering(rendersys* rsys) override;
    void draw(rendersys* rsys) override;
    void tracing() const override;

protected:
    vertex_stream_cr    _vertices;

private:
    void create_from_fill(bat_fill_batch* bat);
    void create_from_stroke(bat_stroke_batch* bat);
};

class rose_fill_batch_klm_cr:
    public rose_batch
{
public:
    rose_fill_batch_klm_cr(int index) : rose_batch(index) {}
    rose_batch_tag get_tag() const override { return bf_klm_cr; }
    void create(bat_batch* bat) override;
    int buffering(rendersys* rsys) override;
    void draw(rendersys* rsys) override;
    void tracing() const override;

protected:
    vertex_stream_klm_cr _vertices;
};

class rose_fill_batch_klm_tex:
    public rose_batch
{
    friend class rose_stroke_batch_assoc_with_klm_tex;

public:
    rose_fill_batch_klm_tex(int index, render_sampler_state* ss);
    ~rose_fill_batch_klm_tex() { destroy(); }
    rose_batch_tag get_tag() const override { return bf_klm_tex; }
    void create(bat_batch* bat) override;
    int buffering(rendersys* rsys) override;
    void draw(rendersys* rsys) override;
    void tracing() const override;
    void destroy();

protected:
    vertex_stream_klm_tex _vertices;
    tex_batcher         _texbatch;
    render_sampler_state* _sstate;
    shader_resource_view* _srv;
    render_texture2d*   _tex;

private:
    void create_from_fill(bat_fill_batch* bat);
    void create_from_stroke(bat_stroke_batch* bat);
};

class rose_stroke_batch_coef_cr:
    public rose_batch
{
public:
    rose_stroke_batch_coef_cr(int index) : rose_batch(index) {}
    rose_batch_tag get_tag() const override { return bs_coef_cr; }
    void create(bat_batch* bat) override;
    int buffering(rendersys* rsys) override;
    void draw(rendersys* rsys) override;
    void tracing() const override;

protected:
    vertex_stream_cf_cr _vertices;
};

class rose_stroke_batch_coef_tex:
    public rose_batch
{
public:
    rose_stroke_batch_coef_tex(int index, render_sampler_state* ss);
    ~rose_stroke_batch_coef_tex() { destroy(); }
    rose_batch_tag get_tag() const override { return bs_coef_tex; }
    void create(bat_batch* bat) override;
    int buffering(rendersys* rsys) override;
    void draw(rendersys* rsys) override;
    void tracing() const override;
    void destroy();

protected:
    vertex_stream_cf_tex _vertices;
    tex_batcher         _texbatch;
    render_sampler_state* _sstate;
    shader_resource_view* _srv;
    render_texture2d*   _tex;

protected:
    void create_vertices(bat_batch* bat);
    void create_vertices(bat_lines& lines, const tex_batcher& bat);
};

class rose_stroke_batch_assoc_with_klm_tex:
    public rose_stroke_batch_coef_tex
{
public:
    rose_stroke_batch_assoc_with_klm_tex(int index, rose_fill_batch_klm_tex* assoc);
    ~rose_stroke_batch_assoc_with_klm_tex();
    void create(bat_batch* bat) override;
    int buffering(rendersys* rsys) override;

protected:
    rose_fill_batch_klm_tex* _assoc;

protected:
    void create_vertices(bat_batch* bat);
};

class rose_bindings
{
public:
    rose_bind_list_cr& get_cr_bindings() { return _cr_bindings; }
    rose_bind_list_tex& get_tex_bindings() { return _tex_bindings; }
    void clear_binding_cache();

protected:
    rose_bind_list_cr   _cr_bindings;
    rose_bind_list_tex  _tex_bindings;
};

class graphics_obj_entity:
    public loop_blinn_processor
{
public:
    graphics_obj_entity(float w, float h) : loop_blinn_processor(w, h) {}
    void proceed_fill(const painter_path& path) { __super::proceed(path); }
    void proceed_stroke(const painter_path& path);

protected:
    int create_from_path(const painter_path& path, int start);

    struct path_seg
    {
        lb_line*        first;
        lb_line*        last;
        path_seg() { first = last = nullptr; }
    };
    void add_line_seg(path_seg& seg, const painter_path::line_to_node* node);
    void add_quad_seg(path_seg& seg, const painter_node* node1, const painter_path::quad_to_node* node2);
    void add_cubic_seg(path_seg& seg, const painter_node* node1, const painter_path::cubic_to_node* node2);
};

class graphics_obj:
    public std::shared_ptr<graphics_obj_entity>
{
public:
    graphics_obj(float w, float h);
};

extern void rose_paint_non_picture_brush(graphics_obj& gfx, rose_bindings& bindings, const painter_brush& brush);
extern void rose_paint_solid_brush(graphics_obj& gfx, rose_bind_list_cr& bind_cache, const painter_brush& brush);
extern void rose_paint_picture_brush(graphics_obj& gfx, const rectf& bound, rose_bind_list_tex& bind_cache, const painter_brush& brush);
extern void rose_paint_pen(graphics_obj& gfx, const painter_path& path, rose_bindings& bindings, const painter_pen& pen);
extern void rose_paint_solid_pen(graphics_obj& gfx, rose_bind_list_cr& bind_cache, const painter_pen& pen);
extern void rose_paint_picture_pen(graphics_obj& gfx, const rectf& bound, rose_bind_list_tex& bind_cache, const painter_pen& pen);
/* more to come. */

class rose:
    public painter
{
public:
    typedef render_constant_buffer constant_buffer;
    typedef render_sampler_state sampler_state;
    typedef list<graphics_obj> graphics_obj_cache;

public:
    rose();
    virtual ~rose();
    virtual void resize(int w, int h) override {}
    virtual void draw_path(const painter_path& path) override;
    virtual void on_draw_begin() override;
    virtual void on_draw_end() override;

public:
    void setup(rendersys* rsys);
    rendersys* get_rendersys() const { return _rsys; }
    void fill_non_picture_graphics_obj(graphics_obj& gfx, uint brush_tag);
    bat_batch* fill_picture_graphics_obj(graphics_obj& gfx);
    void stroke_graphics_obj(graphics_obj& gfx, uint pen_tag);

protected:
    template<class _addline>
    void meta_stroke_graphics_obj(graphics_obj& gfx, _addline fn_add);

protected:
    rendersys*          _rsys;
    rose_configs        _cfgs;
    batch_processor     _bp;
    rose_batch_list     _batches;
    rose_bindings       _bindings;
    graphics_obj_cache  _gocache;
    float               _nextz;

protected:
    void setup_configs();
    void prepare_fill(const painter_path& path, const painter_brush& brush);
    void prepare_picture_fill(const painter_path& path, const painter_brush& brush);
    void prepare_stroke(const painter_path& path, const painter_pen& pen);
    rose_batch* create_fill_batch_cr(int index);
    rose_batch* create_fill_batch_klm_cr(int index);
    rose_batch* create_fill_batch_klm_tex(int index);
    rose_batch* create_stroke_batch_cr(int index);
    rose_batch* create_stroke_batch_tex(int index);
    rose_batch* create_stroke_batch_assoc(int index, rose_fill_batch_klm_tex* assoc);
    void clear_batches();
    void prepare_batches();
    void draw_batches();

#if use_rendersys_d3d_11
protected:
    void initialize();
    void destroy_miscs();
    sampler_state* acquire_default_sampler_state();

protected:
    vertex_shader*      _vsf_cr;
    vertex_shader*      _vsf_klm_cr;
    vertex_shader*      _vsf_klm_tex;
    pixel_shader*       _psf_cr;
    pixel_shader*       _psf_klm_cr;
    pixel_shader*       _psf_klm_tex;
    vertex_format*      _vf_cr;
    vertex_format*      _vf_klm_cr;
    vertex_format*      _vf_klm_tex;
    vertex_shader*      _vss_coef_cr;
    vertex_shader*      _vss_coef_tex;
    pixel_shader*       _pss_coef_cr;
    pixel_shader*       _pss_coef_tex;
    vertex_format*      _vf_coef_cr;
    vertex_format*      _vf_coef_tex;
    sampler_state*      _sampler_state;
    constant_buffer*    _cb_configs;
    uint                _cb_config_slot;
#endif
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib