gimme-density.c

Go to the documentation of this file.

/*
** gimme-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  Fri Feb 27 15:03:31 2009 Johnny Q. Hacker
** Last update Sun May 12 01:17:25 2002 Speed Blue
*/

#include <stdio.h>
#include <glib.h>
#include <gio/gio.h>
#include <stdlib.h>
#include <math.h>
#include <gimmeprogs/process-gimme-density-args.h>
#include <libnextnano/read_nextnano_inputfile.h>
#include <libnextnano/read_nextnano_density.h>
#include <gimmeprogs/argdefaults.h>
#include <libcommon/density.h>

struct integrand_range {
  double xmin, xmax;
};

/*Faster, if you're not gonna worry about roundoff*/
static gint coordinate_density_sorter(gconstpointer a, gconstpointer b) {
  double *c;
  double *d;
  if((c=g_hash_table_lookup((GHashTable*)a, "x")) == NULL) {
    if((d=g_hash_table_lookup((GHashTable*)b, "x")) == NULL) return 0;
    return -1;
  }else if((d=g_hash_table_lookup((GHashTable*)b, "x")) == NULL) {
    return 1;
  }
  if(*c < *d) return -1;
  if(*c > *d) return 1;
  return 0;
}

/*Slower, but gets the answer right if you have roundoff*/
static gint approx_coordinate_density_sorter(gconstpointer a, gconstpointer b) {
  double *c;
  double *d;
  if((c=g_hash_table_lookup((GHashTable*)a, "x")) == NULL) {
    if((d=g_hash_table_lookup((GHashTable*)b, "x")) == NULL) return 0;
    return -1;
  }else if((d=g_hash_table_lookup((GHashTable*)b, "x")) == NULL) {
    return 1;
  }
  //fprintf(stderr, "got x values c=%g(%p) and d=%g(%p) (a=%p,b=%p)\n", *c, c, *d, d, a, b);
  if(fabs(*c - *d) < fabs(CLOSE_ENOUGH*(*c))) return 0;
  if(*c < *d) return -1;
  if(*c > *d) return 1;
  return 0;
}

int main (int argc, char* argv[]) {
  g_type_init();
  GError *err = NULL;
  struct gimme_density_args* arg = process_gimme_density_args(argc, argv);
  if(arg == NULL) exit(2);
  GList *l = NULL;
  struct integrand_range *s;
  if(arg->material != NULL) {
    printf("MATERIAL-SPECIFIED RANGES NOT IMPLEMENTED YET\n");
    exit(1);
  }else{
    s = (struct integrand_range*)malloc(sizeof(struct integrand_range));
    s->xmin = arg->xmin;
    s->xmax = arg->xmax;
    l = g_list_append(l, s);
  }
  /*Read in the density information*/
  //printf("Got xmin=%g xmax=%g\n", ((struct integrand_range*)(l->data))->xmin, ((struct integrand_range*)(l->data))->xmax);

  GFile *keywords = g_file_get_child(arg->base_dir, "keywords.in");
  GFile *database = g_file_get_child(arg->base_dir, "database.in");
  GFileInputStream *db_istream = g_file_read(database, NULL, &err);
  if(err != NULL) {
    fprintf(stderr, "Got error opening database file %s: %s\n", arg->dir, err->message);
    exit(1);
  }
  GHashTable *db = read_nextnano_dbfile(G_INPUT_STREAM(db_istream), &err);
  if(err != NULL) {
    fprintf(stderr, "Got error reading database file %s: %s\n", arg->dir, err->message);
    exit(1);
  }
  GFile *input_file = g_file_get_child(arg->base_dir, "input_file1.in");
  GFileInputStream *input_istream = g_file_read(input_file, NULL, &err);
  if(err != NULL) {
    fprintf(stderr, "Got error opening input file %s: %s\n", arg->dir, err->message);
    exit(1);
  }
  GHashTable *ifile_hash = read_nextnano_inputfile_stringy(G_INPUT_STREAM(input_istream), &err);

  GList *density_func = read_nextnano_density_el(arg->base_dir, ifile_hash);
  if(density_func == NULL) {
    fprintf(stderr, "ERROR: error reading in the density function.\n");
    exit(4);
  }
  /*Ensure that we're ordered by point*/
  density_func = g_list_sort(density_func, approx_coordinate_density_sorter);
  int len = g_list_length(density_func);
  struct density_point *dp;
  if((dp = (struct density_point*)malloc(len*sizeof(struct density_point))) == NULL) {
    fprintf(stderr, "ERROR: failed to allocate enough memory to store the density information\n");
    exit(5);
  }
  int i;
  GList *density_pt;
  double *v;
  for(i=0, density_pt=density_func; i<len; i++, density_pt = g_list_next(density_pt)) {
    if(unlikely((v = g_hash_table_lookup((GHashTable*)(density_pt->data), "x")) == NULL)) {
      fprintf(stderr, "INTERNAL ERROR: a density function point's x coordinate was NULL!\n");
    }
    dp[i].x = *v;
      if(unlikely((v = g_hash_table_lookup((GHashTable*)(density_pt->data), "total")) == NULL)) {
      fprintf(stderr, "INTERNAL ERROR: a density function point's total density was NULL!\n");
    }
    dp[i].rho = *v;
  }
  /*
  printf("#x\trho\n");
  for(i=0; i<len; i++) {
    printf("%g\t%g\n", dp[i].x, dp[i].rho);
  }
  printf("GSL:%g\n", get_density_within_range(
  */
  printf("%ge11 (cm^-2)\n", get_density_within_range(
            ((struct integrand_range*)l->data)->xmin, 
            ((struct integrand_range*)l->data)->xmax, 
            dp, 
            len, 0));
  /*
  printf("const interpolation: %g\n", get_density_within_range(
            ((struct integrand_range*)l->data)->xmin, 
            ((struct integrand_range*)l->data)->xmax, 
            dp, 
            len, 1));
  printf("summation: %g\n", get_density_within_range(
            ((struct integrand_range*)l->data)->xmin, 
            ((struct integrand_range*)l->data)->xmax, 
            dp, 
            len, 2));
  */
  return 0;
}