density_1d.c

Go to the documentation of this file.

 /*
 *1D density routines

    Copyright (C) 2010 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/>.

    *
    */

#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>
#include <libpostproc/density_1d.h>

//#define DEBUG_DENSITY_1D_C

#define DENSITY_1D_GERROR_DOMAIN 3328

void init_density_1d(struct density_1d* d, double initial_value) {
  for(int i=0; i<d->N; i++) {
    d->storage[i] = initial_value;
  }
}

void density_1d_mul_by_const(struct density_1d* d, double c) {
  for(int i=0; i<d->N; i++) {
    d->storage[i] *= c;
  }
}

struct density_1d* new_density_1d(int N, double dx, double initial_value) {
  struct density_1d *d = (struct density_1d*)malloc(sizeof(struct density_1d));
  if(d == NULL) return NULL;
  d->N = N;
  d->dx = dx;
  d->storage = (double*)malloc(sizeof(double)*N);
  if(d->storage == NULL) {
    free(d);
    return NULL;
  }
  init_density_1d(d, initial_value);
  return d;
}

struct density_1d* new_density_1d_copy(const struct density_1d* orig) {
  struct density_1d *d = (struct density_1d*)malloc(sizeof(struct density_1d));
  if(d == NULL) return NULL;
  d->N = orig->N;
  d->dx = orig->dx;
  d->storage = (double*)malloc(sizeof(double)*orig->N);
  if(d->storage == NULL) {
    free(d);
    return NULL;
  }
  for(int i=0; i<orig->N; i++) d->storage[i] = orig->storage[i];
  return d;
}

void free_density_1d(struct density_1d* d) {
  free(d->storage);
  free(d);
}

int op_density_1d_to_density_1d(struct density_1d* to, enum density_op op, const struct density_1d* from) {
  if(to->N != from->N) return 1;
  #ifdef OPENMP
  #pragma omp parallel for
  #endif
  for(int i=0; i<to->N; i++) {
    switch(op) {
    case ASSIGN:
      to->storage[i]
  += from->storage[i]
  ;
    break;
    case ADD_EQUALS:
      to->storage[i]
  += from->storage[i]
  ;
      break;
    case MINUS_EQUALS:
      to->storage[i]
  -= from->storage[i]
  ;
      break;
    case MUL_EQUALS:
      to->storage[i]
  *= from->storage[i]
  ;
      break;
    case DIV_EQUALS:
      to->storage[i]
  /= from->storage[i]
  ;
      break;
    }
  }
  return 0;
}
 
int op_density_1d_to_density_1d_cb(struct density_1d *restrict to, const struct density_1d *restrict from, void(*callback)(int X, unsigned N, double dx, double *restrict to, double from, void *restrict privdat), void *restrict privdat) {
  if(to->N != from->N) return 1;
  #ifdef OPENMP
  #pragma omp parallel for
  #endif
  for(int i=0; i<to->N; i++) {
    callback(i, to->N, to->dx, &(to->storage[i]), from->storage[i], privdat);
  }
  return 0;
}

double density_1d_sum(const struct density_1d* d, enum density_1d_distribution_hint hint) {
  double total = 0;
  unsigned i;
  switch(hint) {
  case CONSTANT:
    for(i=0; i<d->N; i++) total += d->storage[i];
    break;
  case CENTERED:
    for(i=0; i<(d->N - (d->N%2))/2; i++) {
      total += d->storage[i] + d->storage[d->N-1-i];
    }
    if((d->N%2) != 0) total += d->storage[(d->N - (d->N%2))/2];
    break;
  case ANTICENTERED:
    if((d->N%2) != 0) total = d->storage[(d->N - (d->N%2))/2];
    /*The -1 is because (N-N%2)/2 is the *center* point (if odd; if even, it is the last point on the opposite side!)*/
    for(i=((d->N - (d->N%2))/2-1); i>=0; i--) {
      total += d->storage[i] + d->storage[d->N-1-i];
    }
    break;
  default:
    return nan("NaN");
  }
  return total;
}

struct density_1d_params {
  gsl_spline *d;
  gsl_interp_accel *d_acc;
#ifdef DEBUG_DENSITY_1D_C
  FILE* file;
  long unsigned count;
#endif
};

static double integrand(double x, void *params) {
  struct density_1d_params* p = (struct density_1d_params*)params;
  double val=gsl_spline_eval(p->d, x, p->d_acc);
  #ifdef DEBUG_DENSITY_1D_C
  fprintf(p->file, "%lu\t%g\t%g\n", p->count++, x, val);
  #endif
  return val;
}

double density_1d_integrate_density(const struct density_1d* d, double abserr, double relerr, double* ierr, GError **e, int filename_index) {
  gsl_interp_accel *density_acc = gsl_interp_accel_alloc();
  /*Akima, so that it can cope with dropoffs.*/
  double* x = (double*)malloc(d->N*sizeof(double));
  if(x == NULL) {
    *e = g_error_new(DENSITY_1D_GERROR_DOMAIN, 100, "funkalicious::libpostproc::density_1d.c::density_1d_integrate_density: Failed to allocate array of x values");
    return nan("NaN");
  }
  for(unsigned i=0; i<d->N; i++) x[i] = i*(d->dx);
  gsl_spline *density = gsl_spline_alloc(gsl_interp_akima, d->N);
  gsl_spline_init(density, x, d->storage, d->N);
  free(x);
  struct density_1d_params ig = {
    .d = density, 
    .d_acc = density_acc
    #ifdef DEBUG_DENSITY_1D_C
    ,
    .file = NULL,
    .count = 0
    #endif
  };
  gsl_function gsl_func = {
    .function = &integrand,
    .params = &ig
  };
  double total;
  int status=0;
  gsl_integration_workspace *w = gsl_integration_workspace_alloc(4*d->N);
  double start=0, end=d->N-1;
  for(unsigned i=0; i<(d->N+1); i++) {
    /*fprintf(stderr, "start: %u->%g\n", i, d->storage[i]);*/
    /*If we hit the end, there's no charge in the structure!*/
    if(i == d->N) return 0;
    if(d->storage[i] != 0) {
      if(i == 0) {
  start = 0;
      }else{
  start = (i-1)*d->dx;
      }
      break;
    }
  }
  for(unsigned i=0; i<(d->N+1); i++) {
    /*fprintf(stderr, "end: %u->%g\n", i, d->storage[i]);*/
    if(d->storage[d->N-1-i] != 0) {
      if(i==0) {
  end = (d->N-1)*d->dx;
      }else{  
  end = (d->N-i)*d->dx;
      }
      break;
    }
  }
  /*fprintf(stderr, "start=%g end=%g\n", start, end);*/
  #ifdef DEBUG_DENSITY_1D_C
  char fn[2048];
  sprintf(fn, "density_1d_integrand_log_%d.pdata", filename_index);
  ig.file = fopen(fn, "w");
  ig.count = 0;
  #endif
  gsl_integration_qag(&gsl_func, start, end, abserr, relerr, 4*d->N, GSL_INTEG_GAUSS61, w, &total, ierr);
  #ifdef DEBUG_DENSITY_1D_C
  fclose(ig.file);
  #endif
  if(status != 0) {
    *e = g_error_new(DENSITY_1D_GERROR_DOMAIN, status, "funkalicious::libpostproc::density_1d.c::density_1d_integrate_density: Got error (status=%d) while integrating density to find total density: %s", status, gsl_strerror(status));
  }
  gsl_integration_workspace_free(w);
  gsl_spline_free(density);
  gsl_interp_accel_free(density_acc);
  return total;
}


void density_1d_add_const(struct density_1d* d, double c) {
  #ifdef OPENMP
  #pragma omp parallel for
  #endif
  for(long unsigned i=0; i<d->N; i++) d->storage[i] += c;
}