/* THIS PROGRAM DOES NOT YET WORK!!! */

#include <stdio.h> /* for input/output */
#define FunctionDataFileName "bessel.dat"
#define X 0
#define Y 1
#define maxData 200
#define fullmaxData 400

/* This program gets slower for higher numbers */

double data[fullmaxData]; /* poor man's "dynamic" array :-/ */
int numData = 0;

void store(double x, double y) {
  /* this is where we would store the value to a dynamic array */
  if (numData >= maxData) {
    printf("EEK! More data than possible!\nI can only store up to %d data pairs!\n", maxData);
    /*+*/
  } else 
    if ((numData > 0) && 
        (x <= data[(numData-1)*2])) {
      printf("EEK! Data file unordered at line %d!\n Data file must be sorted on the x axis.", numData);
      /*+*/
    }
  printf("data[%d]: x=%f; y=%f\n", numData, x, y);
  data[numData*2] = x;
  data[numData*2+1] = y;
  numData++;
}

int lookup(double x, double *x1, double *y1, double *x2, double *y2) {
  /* find values of x and y which are just lower and just above x. */
  /* return 1 if value is found exactly (x2, y2 will be undefined) */
  /*+* assumes data is initialized with at least two data points! *+*/
  int Point;
  if ((x < data[0*2]) || (x > data[(numData-1)*2+1])) {
    printf("EEK! Data point %f is out of range!\nMin: %f; Max: %f.", x, data[0], data[numData-1]);
    /*+*/
  }  
  Point = 0;
  while (data[Point*2] <= x) { Point++; }
  *x1 = x;
  *y1 = data[(Point-1)*2+1];
  if (data[(Point-1)*2] == x) {
    return 1;
  } else {
    *x2 = x;
    *y2 = data[(Point-1)*2+1];
    return 0;
  }
}

void init(char datafilename[255]) {
  /*
     initialisation:

        open file
        while not feof :
           read x
           read y
           store y in array position x
        close file

   */
  FILE *datafile;
  double x, y;
  datafile = fopen(datafilename, "r");
  while (fscanf(datafile, "%f", &x) != EOF) {
    fscanf(datafile, "%f", &y);
    store(x,y);
  }
  close(datafile);
}

main() {
  int Exact;
  double x, x1, y1, x2, y2;
  printf("Lab Session 6, Question 4.4\n"); /* Display program header */
  init(FunctionDataFileName);
  x = 1.0;
  if (lookup(x, &x1, &y1, &x2, &y2)) {
    printf("%f", y1);
  } else {
    printf("%f", y1 + (y2-y1)/(x2-x1) * (x-x1));
  }  
}


/*
Pseudocode For Looking Up A Function
------------------------------------

interpolation:

   select value of x
   if we know y(x)
      set y = y(x)
   otherwise
   if x is lower than lowest x in array
      set s = complaint
      set y = nan
   otherwise
   if x is higher than highest x in array
      set s = complaint
      set y = nan
   otherwise
      get highest value of x in array which is lower than x (x1)
      get lowest value of x in array which is higher than x (x2)
      evaluate y = y(x1) + (y(x2)-y(x1))/(x2-x1) * (x-x1)
   if not y = nan,
      print y
   otherwise
      print complain
*/