density.c

Go to the documentation of this file.

/*
** density.c
** 
** Made by (Johnny Q. Hacker)
** Login   <solarion@borkborkbork>

    Copyright (C) 2009 Joseph Pingenot

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU Affero General Public License as published by
    the Free Software Foundation, either version 3 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 Affero General Public License for more details.

    You should have received a copy of the GNU Affero General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

** Started on  Tue Mar  3 15:40:58 2009 Johnny Q. Hacker
** Last update Sun May 12 01:17:25 2002 Speed Blue
*/

#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <libcommon/density.h>
#include <libcommon/gcc_hints.h>

double get_density_within_range2(double xmin, double xmax, double x[], double rho[], int len) {
  /*Use GSL to integrate using a spline interpolator*/
  gsl_interp_accel *acc = gsl_interp_accel_alloc();
  gsl_spline *spline = gsl_spline_alloc(gsl_interp_akima_periodic, len);
  gsl_spline_init(spline, x, rho, len);
  double density = gsl_spline_eval_integ(spline, xmin, xmax, acc);
  gsl_spline_free(spline);
  gsl_interp_accel_free(acc);
  return density;
}

double get_density_within_range2_stairstep(double xmin, double xmax, double x[], double rho[], int len) {
  /*Use simple stair step (const interpolation) to calculate the density*/
  double density=0;
  int prev=0;
  for(int i=1; i<len; i++) {
    if((x[i] >= xmin) && (x[i] <= xmax)) {
      if(unlikely((x[i-1] < xmin) && (x[i] != xmin))) {
  /*If we stepped over xmin, we only have a partial distance*/
  density+=rho[i-1]*(xmin-x[i]);
      }else{
  density+=rho[i-1]*(x[i]-x[i-1]);
      }
    }
    if(x[i] > xmax) {
      if(x[i-1] < xmax) {
  density += rho[i-1]*(xmax-x[i-1]);
      }
      /*No point in continuing the loop past this exit point.*/
      break;
    }
  }
  return density;
}

double get_density_within_range2_sum(double xmin, double xmax, double x[], double rho[], int len) {
  /*Use simple stair step (const interpolation) to calculate the density*/
  double density=0;
  int prev=0;
  for(int i=1; i<len; i++) {
    if((x[i] >= xmin) && (x[i] <= xmax)) {
      if(unlikely((x[i-1] < xmin) && (x[i] != xmin))) {
  /*If we stepped over xmin, we only have a partial distance*/
  density+=rho[i-1];
      }else{
  density+=rho[i-1];
      }
    }
    if(x[i] > xmax) {
      if(x[i-1] < xmax) {
  density += rho[i-1];
      }
      /*No point in continuing the loop past this exit point.*/
      break;
    }
  }
  return density;
}

double get_density_within_range(double xmin, double xmax, struct density_point *p, int len, short int type) {
  int i;
  double x[len];
  double rho[len];
  for(i=0; i<len; i++) {
    x[i] = p[i].x;
    rho[i] = p[i].rho;
    //fprintf(stderr, "i=%d: x=%g rho=%g -> x=%g rho=%g\n", i, p[i].x, p[i].rho, x[i], rho[i]);
  }
  switch (type) {
  case 1:
    return get_density_within_range2_stairstep(xmin, xmax, x, rho, len);
    break;
  case 2:
    return get_density_within_range2_sum(xmin, xmax, x, rho, len);
    break;
  default:
    return get_density_within_range2(xmin, xmax, x, rho, len);
  }
}