#include <glib-2.0/glib.h>
#include <gio/gio.h>
#include <libpostproc/density_common.h>

Include dependency graph for density_1d.h:

This graph shows which files directly or indirectly include this file:

Classes
struct	density_1d
Enumerations
enum	density_1d_distribution_hint { CONSTANT, CENTERED, ANTICENTERED }
Functions
void	init_density_1d (struct density_1d *d, double initial_value)
void	density_1d_mul_by_const (struct density_1d *d, double c)
void	density_1d_add_const (struct density_1d *d, double c)
struct density_1d *	new_density_1d (int N, double dx, double initial_value)
struct density_1d *	new_density_1d_copy (const struct density_1d *orig)
void	free_density_1d (struct density_1d *d)
int	op_density_1d_to_density_1d (struct density_1d to, enum density_op op, const struct density_1d from)
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)
double	density_1d_sum (const struct density_1d *d, enum density_1d_distribution_hint hint)
double	density_1d_integrate_density (const struct density_1d d, double abserr, double relerr, double ierr, GError **e, int filename_index)

Enumeration Type Documentation

enum density_1d_distribution_hint

Density distribution hints. These help ensure numerical accuracy. constant: there is roughly an even distribution across the structure centered: there is very very little at the edges and much more at the center. anticentered: there is very little at the center and much more at the edges.

Enumerator:

CONSTANT
CENTERED
ANTICENTERED

Definition at line 45 of file density_1d.h.

                                    {
    CONSTANT,
    CENTERED,
    ANTICENTERED
  };

Function Documentation

void density_1d_add_const	(	struct density_1d *	d,
		double	c
	)

Definition at line 245 of file density_1d.c.

References density_1d::N, and density_1d::storage.

Referenced by potential_1d_get_potential_integrate().

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

Here is the caller graph for this function:

double density_1d_integrate_density	(	const struct density_1d *	d,
		double	abserr,
		double	relerr,
		double *	ierr,
		GError **	e,
		int	filename_index
	)

Definition at line 171 of file density_1d.c.

References density_1d_params::d, DENSITY_1D_GERROR_DOMAIN, density_1d::dx, integrand(), density_1d::N, and density_1d::storage.

                                                                                                                                            {
  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;
}

Here is the call graph for this function:

void density_1d_mul_by_const	(	struct density_1d *	d,
		double	c
	)

Definition at line 40 of file density_1d.c.

References c, density_1d::N, and density_1d::storage.

Referenced by get_potential_for_density_1d().

                                                             {
  for(int i=0; i<d->N; i++) {
    d->storage[i] *= c;
  }
}

Here is the caller graph for this function:

double density_1d_sum	(	const struct density_1d *	d,
		enum density_1d_distribution_hint	hint
	)

Gets the sum of the density, using distribution hint.

Definition at line 127 of file density_1d.c.

References ANTICENTERED, CENTERED, CONSTANT, density_1d::N, and density_1d::storage.

                                                                                          {
  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;
}

void free_density_1d ( struct density_1d * d )

Definition at line 74 of file density_1d.c.

References density_1d::storage.

Referenced by get_potential_for_density_1d(), potential_1d_get_dfield_integrate(), potential_1d_get_potential_from_dfield(), and potential_1d_get_potential_integrate().

                                           {
  free(d->storage);
  free(d);
}

Here is the caller graph for this function:

void init_density_1d	(	struct density_1d *	d,
		double	initial_value
	)

Definition at line 34 of file density_1d.c.

References density_1d::N, and density_1d::storage.

Referenced by new_density_1d().

                                                                 {
  for(int i=0; i<d->N; i++) {
    d->storage[i] = initial_value;
  }
}

Here is the caller graph for this function:

struct density_1d* new_density_1d	(	int	N,
		double	dx,
		double	initial_value
	)		`[read]`

Definition at line 46 of file density_1d.c.

References density_1d::dx, init_density_1d(), density_1d::N, and density_1d::storage.

Referenced by get_potential_for_density_1d(), potential_1d_get_dfield_integrate(), and potential_1d_get_potential_from_dfield().

                                                                          {
  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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

struct density_1d* new_density_1d_copy ( const struct density_1d * orig ) [read]

Definition at line 60 of file density_1d.c.

References density_1d::dx, density_1d::N, and density_1d::storage.

Referenced by get_potential_for_density_1d().

                                                                      {
  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;
}

Here is the caller graph for this function:

int op_density_1d_to_density_1d	(	struct density_1d *	to,
		enum density_op	op,
		const struct density_1d *	from
	)

Definition at line 79 of file density_1d.c.

References ADD_EQUALS, ASSIGN, DIV_EQUALS, MINUS_EQUALS, MUL_EQUALS, density_1d::N, and density_1d::storage.

Referenced by get_potential_for_density_1d().

                                                                                                          {
  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;
}

Here is the caller graph for this function:

int op_density_1d_to_density_1d_cb	(	struct density_1d *restrict	to,
		const struct density_1d *restrict	from,
		void()(int X, unsigned N, double dx, double restrict to, double from, void *restrict privdat)	callback,
		void *restrict	privdat
	)

Allows the caller to specify a callback to perform a more advanced manipulation.

Parameters:

to	density to which values are saved and from which values are drawn
from	density from which values are drawn
callback	function to call at each point
privdat	private data to send to the callback

Note:: CALLBACK MUST NOT SAVE STATE; callback may be called in numerous threads!; that to and from CANNOT be the same or overlap! (restrict keyword)

Definition at line 116 of file density_1d.c.

                                                                                                                                                                                                                                            {
  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;
}

density_1d.h File Reference

Classes

Enumerations

Functions

Enumeration Type Documentation

Function Documentation