#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
#include <sys/time.h>
#include <fcntl.h>
#include <assert.h>
#include <errno.h>
#include <sys/types.h>
#include <sys/wait.h>
#include "cement.h"
#include "pnmutils.h"
#include "file.h"


/* Global Variables */

extern int errno;
extern char **environ;

/* Function prototypes needed in executable and in library */
static int my_system (char *command);
static inline double createValue(double a, unsigned char b, double weight);

static int verbose=0;
static double powLookupTable[256];
static double exponent;

static int convert = 0;
static int grayscale = 0;
static int newfile = 0;


#ifndef OBJECT_COMPILE
/* Functions needed only in executable */
void usage(char * string);
void parse_commandline(int argc, char ** argv,
                       char**input_filenameA,
                       char**input_filenameB,
                       double *R,
                       double *G,
                       double *B,
                       double *yR,
                       double *yG,
                       double *yB);

void mystatus(char *action, char *filename, double R,double G, double B, double complete);
#endif

void readlut(char *lutfilename);

/* Functions */
#ifndef CEMENT_REMOVE
int cementi(char *input_filenameA,
            char *input_filenameB,
            double R,
            double G,
            double B,
	    double yR,
	    double yG,
	    double yB,
            char* command,
            void (*status)(char *action, char *filename, double R,double G, double B, double complete)
           )
#else
int cemento(char *input_filenameA,
            char *input_filenameB,
            double R,
            double G,
            double B,
	    double yR,
	    double yG,
	    double yB,
            char* command,
            void (*status)(char *action, char *filename, double R,double G, double B, double complete)
           )
#endif
{
    struct image_params a_params={NULL,0,0,0,0,0.0};
    struct image_params b_params={NULL,0,0,0,0,0.0};
    File * a, * b, * dest, * new;
    double p;  /* exponent for pow */
    double        a_pixel[ARRAY_SIZE];
    unsigned char b_pixel[ARRAY_SIZE];
    double out_pixel[ARRAY_SIZE];
    int i, j, k, exp;
    int channels_per_pixel;
    char * lut_filename;
//    double R, G, B;  /* weights for cement */
    double new_pixel;
    double * newimage;

    /* Timing measurments */
    struct timeval timeValStart;
    struct timeval timeValEnd;
    struct timezone timeZone;
    int minutes;
    double seconds=0.0;
#ifdef GET_TIME_MEASUREMENTS
    struct timeval runTimeStart;
    struct timeval runTimeStop;
    double fileReadTime=0.0;
    double fileWriteTime=0.0;
    double computeTime=0.0;
#endif
    unsigned int pixels=0;
    unsigned int numberOfPixels;

    char string[MAX_STRING_LENGTH];
    unsigned int image_ptr=0;

    a_params.filename=input_filenameA;
    b_params.filename=input_filenameB;
    
    //    parse_commandline(argc, argv, &a_params, &b_params, &R, &G, &B);

    if ((a=fileopen(a_params.filename, "r"))==NULL)
    {
        fprintf(stderr, "Unable to open %s.\n", a_params.filename);
	fprintf(stderr, "Creating a new file.\n");
	newfile = 1;
    }

    if ((b=fileopen(b_params.filename, "r"))==NULL)
    {
        fprintf(stderr, "Unable to open %s.\n", b_params.filename);
        exit(EXIT_FAILURE);
    }

    sprintf(string,"tmp.%s",a_params.filename);
    if ((dest=fileopen(string, "w"))==NULL)
    {
        fprintf(stderr, "Unable to open %s.\n", string);
        exit(EXIT_FAILURE);
    }

    if(newfile == 1) 
      {
	a_params.type='A';
	channels_per_pixel=3;
	sprintf(string, "P%c\n", a_params.type);
	filewrite(string,sizeof(unsigned char),strlen(string),dest);
      }
    else 
      { 
	get_image_type(a, &a_params);
	if(a_params.type=='A') {
	  channels_per_pixel=3;
	  sprintf(string, "P%c\n", a_params.type);
	  filewrite(string,sizeof(unsigned char),strlen(string),dest);
	  get_image_params_preserve_comments(a,&a_params,dest);
	}
	else
	  {
	    fprintf(stderr, "Portable Lightspace Map (plm) must be type A and it is type %c\n",a_params.type);
	    exit(EXIT_FAILURE);
	  }
      }

    if(isJpegFile(b)==0)
    {
        get_image_params(b, &b_params);
	if(b_params.type=='5')
	  grayscale = 1;
    }
    else
    {
      b_params.type='6';
      b_params.max_val=255;
      
      p = exponent = a_params.exponent;
      lut_filename = (char*) malloc(50 * sizeof(char));
      exp = rint(p * 10.0);
      
      sprintf(lut_filename, "powLookup%d.txt", exp);
      printf("p: %6g  exp: %d  filename: %s\n", p, exp, lut_filename);
      
      
      if(lut_filename !=NULL) 
	{
	  readlut(lut_filename);
	}      if(read_JPEG_file (b_params.filename)==0)
	  {
	    printf("ERROR: File %s does not exist\n",b_params.filename);
	    exit(0);
	  }
      //        write_JPEG_file("test.jpg",95);
      //        exit(0);
      
      image_ptr=0;
      b_params.width=image_width;
      b_params.height=image_height;
    }
    
    if(newfile == 1) {

      /* Create a new plm file and initialize it with zeros. */
      /* This is a kludge, but a lot easier than trying to rewrite the
	 code to get rid of all the file reads. */
      new = fileopen(a_params.filename, "w");
      if (new->fp == NULL) {
	fprintf(stderr,"Couldn't open new file either!\n");
	exit(EXIT_FAILURE);
      }
      fprintf(new->fp, "PA\n");
      fprintf(new->fp, "%d %d\n", b_params.width, b_params.height);
      fprintf(new->fp, "255\n");
      fprintf(new->fp, "2\n"); /* don't know what this field is */
      fprintf(new->fp, "%2g\n", EXPONENT);
      newimage = calloc(b_params.width*b_params.height, sizeof(double));
      if (newimage == NULL) {
	fprintf(stderr, "Unable to allocate.\n");
	exit(EXIT_FAILURE);
      }
      fwrite(newimage, sizeof(double), b_params.width*b_params.height, new->fp);
      fclose(new->fp);

      /* Now open it again and pretend it's the original plm file we were given. */
      a = fileopen(a_params.filename, "r");
      if (a->fp==NULL) {
	fprintf(stderr, "Unable to open %s after supposed creation.\n", a_params.filename);
	exit(EXIT_FAILURE);
      }
      get_image_type(a,&a_params); /* to move the pointer along */
      get_image_params_preserve_comments(a,&a_params,dest);
    }

    p = exponent = a_params.exponent;
    lut_filename = (char*) malloc(50 * sizeof(char));
    exp = rint(p * 10.0);

    sprintf(lut_filename, "powLookup%d.txt", exp);

    if(lut_filename !=NULL) 
    {
        readlut(lut_filename);
    }

    if ((a_params.width!=b_params.width) || (a_params.height!=b_params.height)
        || (a_params.max_val!=b_params.max_val))
    {
        printf("%f %f\n",a_params.exponent,EXPONENT);
        printf("%d %d\n",a_params.width,b_params.width);
        printf("%d %d\n",a_params.height,b_params.height);
        printf("%d %d\n",a_params.max_val,b_params.max_val);
        fprintf(stderr, "Lightspace must have the same dimensions, maximum "
                "values, and exponent.\n");
        exit(EXIT_FAILURE);
    }

    sprintf(string, "# %s %s %s %f %f %f\n",command,input_filenameA,input_filenameB,R,
            G,B);
    filewrite(string,sizeof(unsigned char),strlen(string),dest);
    
#ifndef CEMENT_REMOVE
    sprintf(string, "%d %d\n%d\n%d\n%f\n",a_params.width,
            a_params.height,a_params.max_val, ++a_params.numberOfImages,exponent);
#else
    sprintf(string, "%d %d\n%d\n%d\n%f\n",a_params.width,
            a_params.height, a_params.max_val,--a_params.numberOfImages,exponent);
#endif
    filewrite(string,sizeof(unsigned char),strlen(string),dest);

    
#ifdef SILENT_RUN

#ifndef CEMENT_REMOVE
//    printf("Add %s %d %d %d\r",b_params.filename,(int)R,(int)G,(int)B);
#else
//    printf("Del %s %d %d %d\r",b_params.filename,(int)R,(int)G,(int)B);
#endif
#else

#ifndef CEMENT_REMOVE
    printf("Adding file: \"%s\", ",b_params.filename);
#else
    printf("Removing file: \"%s\", ",b_params.filename);
#endif
    printf("P%c, %dx%d, ",b_params.type,b_params.width,b_params.height);

    printf("R=%7.3f ",R);
    printf("G=%7.3f ",G);
    printf("B=%7.3f to ",B);

    printf("lightspace \"%s\", ",a_params.filename);
    printf("P%c\n",a_params.type);
#endif
    
    /* Start Timer */
    if(gettimeofday(&timeValStart,&timeZone)==-1)
    {
        fprintf(stderr,"Error getting the time\n");
        exit(0);
    }

    numberOfPixels=a_params.width*a_params.height;

    if(((float)numberOfPixels/NUMBER_PIXELS_PER_READ)!=(int)((float)numberOfPixels/NUMBER_PIXELS_PER_READ))
    {
        if(verbose)
            printf("\nERROR: The number of pixels %d is not evenly divisible by %d. ",numberOfPixels,NUMBER_PIXELS_PER_READ);

        numberOfPixels=(int)((float)numberOfPixels/NUMBER_PIXELS_PER_READ)*NUMBER_PIXELS_PER_READ;
        if(verbose)
            printf("We will go up to %d pixels.\n",numberOfPixels);
    }

    for(i=0;i<numberOfPixels;i+=NUMBER_PIXELS_PER_READ)
    {
        if(i%(3*NUMBER_PIXELS_PER_READ)==0)
        {
#ifndef CEMENT_REMOVE
            if(status!=NULL)
                status("Add",b_params.filename,R,G,B,(double)100*i/numberOfPixels);
#else
            if(status!=NULL)
                status("Del",b_params.filename,R,G,B,(double)100*i/numberOfPixels);
#endif
        }
        
#ifdef GET_TIME_MEASUREMENTS
        gettimeofday(&runTimeStart,&timeZone);
#endif

        fileread(a_pixel, sizeof(double),        ARRAY_SIZE, a);

        if(isJpegFile(b)==0)
            fileread(b_pixel, sizeof(unsigned char), ARRAY_SIZE, b);
        else
        {
            memcpy(b_pixel,&image_buffer[image_ptr],ARRAY_SIZE);
            image_ptr+=ARRAY_SIZE;
        }
#ifdef GET_TIME_MEASUREMENTS
        gettimeofday(&runTimeStop,&timeZone);

        seconds=(double)runTimeStop.tv_sec+(double)runTimeStop.tv_usec/1000000;
        seconds-=(double)runTimeStart.tv_sec+(double)runTimeStart.tv_usec/1000000;
        fileReadTime+=seconds;

        gettimeofday(&runTimeStart,&timeZone);
#endif

	k=0;
        for(j=0;j<ARRAY_SIZE;j+=3)
        {
	  /* We want to convert the RGB to grayscale first. */
	  if(convert==1){

	    /* cement the three channels together */
	    new_pixel = (yR*powLookupTable[b_pixel[j]])+(yG*powLookupTable[b_pixel[j+1]])+(yB*powLookupTable[b_pixel[j+2]]);

#ifndef CEMENT_REMOVE
	    /* add the cemented channels (now plm) to the new image. */
	    /* this is essentially the same as create value, but cements 2 plms. */
	    out_pixel[j  ]=a_pixel[j  ] + R*new_pixel;
	    out_pixel[j+1]=a_pixel[j+1] + G*new_pixel;
	    out_pixel[j+2]=a_pixel[j+2] + B*new_pixel;
#else
	    /* sub the cemented channels (now plm) from the new image. */
	    out_pixel[j  ]=a_pixel[j  ] - R*new_pixel;
	    out_pixel[j+1]=a_pixel[j+1] - G*new_pixel;
	    out_pixel[j+2]=a_pixel[j+2] - B*new_pixel;
#endif
	  } 
	  /* The image was grayscale to start with. */
	  else if(grayscale == 1) {
	    out_pixel[j  ]=createValue(a_pixel[j  ],b_pixel[k],R);
	    out_pixel[j+1]=createValue(a_pixel[j+1],b_pixel[k],G);
	    out_pixel[j+2]=createValue(a_pixel[j+2],b_pixel[k],B);
	    k++;
	  }
	  else {
	    out_pixel[j  ]=createValue(a_pixel[j  ],b_pixel[j  ],R);
	    out_pixel[j+1]=createValue(a_pixel[j+1],b_pixel[j+1],G);
	    out_pixel[j+2]=createValue(a_pixel[j+2],b_pixel[j+2],B);
	  }
	}

#ifdef GET_TIME_MEASUREMENTS
        gettimeofday(&runTimeStop,&timeZone);

        seconds=(double)runTimeStop.tv_sec+(double)runTimeStop.tv_usec/1000000;
        seconds-=(double)runTimeStart.tv_sec+(double)runTimeStart.tv_usec/1000000;

        computeTime+=seconds;

        gettimeofday(&runTimeStart,&timeZone);
#endif
        
        filewrite(out_pixel, sizeof(double), ARRAY_SIZE, dest);

#ifdef GET_TIME_MEASUREMENTS
        gettimeofday(&runTimeStop,&timeZone);

        seconds=(double)runTimeStop.tv_sec+(double)runTimeStop.tv_usec/1000000;
        seconds-=(double)runTimeStart.tv_sec+(double)runTimeStart.tv_usec/1000000;
        fileWriteTime+=seconds;
#endif
    }

    if(numberOfPixels!=a_params.width*a_params.height)
    {
        if(verbose)
            printf("Cleanup on pixels %d to %d\n",numberOfPixels,a_params.width*a_params.height);

        for(i=numberOfPixels;i<a_params.width*a_params.height;i++)
        {

            fileread(a_pixel, sizeof(double),CHANNELS_PER_PIXEL, a);

            if(isJpegFile(a)==0)
                fileread(b_pixel, sizeof(unsigned char), CHANNELS_PER_PIXEL, b);
            else
            {
                memcpy(b_pixel,&image_buffer[image_ptr],CHANNELS_PER_PIXEL);
                image_ptr+=CHANNELS_PER_PIXEL;
            }

	    
	    /* We want to convert the RGB to grayscale first. */
	    /* Have to cement the three channels together. */
	    if(convert==1){
	      /* cement the three channels together */
	      new_pixel = (yR*powLookupTable[b_pixel[0]])+(yG*powLookupTable[b_pixel[1]])+(yB*powLookupTable[b_pixel[2]]);

#ifndef CEMENT_REMOVE      
	      /* add the cemented channels (now plm) to the new image. */
	      out_pixel[0]=a_pixel[0] + R*new_pixel;
	      out_pixel[1]=a_pixel[1] + G*new_pixel;
	      out_pixel[2]=a_pixel[2] + B*new_pixel;
#else
	      /* sub the cemented channels (now plm) from the new image. */
	      out_pixel[0]=a_pixel[0] - R*new_pixel;
	      out_pixel[1]=a_pixel[1] - G*new_pixel;
	      out_pixel[2]=a_pixel[2] - B*new_pixel;
#endif
	    }
	    /* The image was grayscale to start with. */
	    else if(grayscale == 1) {
	      out_pixel[0]=createValue(a_pixel[0],b_pixel[0],R);
	      out_pixel[1]=createValue(a_pixel[1],b_pixel[0],G);
	      out_pixel[2]=createValue(a_pixel[2],b_pixel[0],B);
	    }
	    else {
	      out_pixel[0]=createValue(a_pixel[0],b_pixel[0],R);
	      out_pixel[1]=createValue(a_pixel[1],b_pixel[1],G);
	      out_pixel[2]=createValue(a_pixel[2],b_pixel[2],B);
	    }
            filewrite(out_pixel, sizeof(double), CHANNELS_PER_PIXEL, dest);
            pixels++;    
        }
    }
    
    /* End timer */
    if(gettimeofday(&timeValEnd,&timeZone)==-1)
    {
        printf("Error getting the time\n");
        exit(0);
    }

    pixels=i;

    while (!fileeof(a))
    {
        fileread(a_pixel, sizeof(double),CHANNELS_PER_PIXEL, a);

        if(isJpegFile(a)==0)
            fileread(b_pixel, sizeof(unsigned char), CHANNELS_PER_PIXEL, b);
        else
        {
            memcpy(b_pixel,&image_buffer[image_ptr],CHANNELS_PER_PIXEL);
            image_ptr+=CHANNELS_PER_PIXEL;
        }

        if(fileeof(a))
            break;
        
	    
	/* We want to convert the RGB to grayscale first. */
	if(convert==1){	      
	  /* cement the three channels together */
	  new_pixel = (yR*powLookupTable[b_pixel[0]])+(yG*powLookupTable[b_pixel[1]])+(yB*powLookupTable[b_pixel[2]]);
	  
#ifndef CEMENT_REMOVE
	  /* add the cemented channels (now plm) to the new image. */
	  out_pixel[0]=a_pixel[0] + R*new_pixel;
	  out_pixel[1]=a_pixel[1] + G*new_pixel;
	  out_pixel[2]=a_pixel[2] + B*new_pixel;
#else
	  /* sub the cemented channels (now plm) from the new image. */
	  out_pixel[0]=a_pixel[0] - R*new_pixel;
	  out_pixel[1]=a_pixel[1] - G*new_pixel;
	  out_pixel[2]=a_pixel[2] - B*new_pixel;
#endif
	}  
	/* The image was grayscale to start with. */
	else if(grayscale == 1) {
	  out_pixel[0]=createValue(a_pixel[0],b_pixel[0],R);
	  out_pixel[1]=createValue(a_pixel[1],b_pixel[0],G);
	  out_pixel[2]=createValue(a_pixel[2],b_pixel[0],B);
	}
	else {
	  out_pixel[0]=createValue(a_pixel[0],b_pixel[0],R);
	  out_pixel[1]=createValue(a_pixel[1],b_pixel[1],G);
	  out_pixel[2]=createValue(a_pixel[2],b_pixel[2],B);
	}
	
        printf("EP A: %f %f %f\t",a_pixel[0],a_pixel[1],a_pixel[2]);
        printf("B: %d %d %d\t",b_pixel[0],b_pixel[1],b_pixel[2]);
        printf("OUT: %f %f %f\n",out_pixel[0],out_pixel[1],out_pixel[2]);
        filewrite(out_pixel, sizeof(double), CHANNELS_PER_PIXEL, dest);
        pixels++;
    }
    
    fileclose(a);
    fileclose(b);
    fileclose(dest);

    seconds=(double)timeValEnd.tv_sec+(double)timeValEnd.tv_usec/1000000;
    seconds-=(double)timeValStart.tv_sec+(double)timeValStart.tv_usec/1000000;

#ifdef GET_TIME_MEASUREMENTS
    printf("%f s doing File Read %f%%\n",fileReadTime,100*fileReadTime/seconds);
    printf("%f s doing File Writes %f%%\n",fileWriteTime,100*fileWriteTime/seconds);
    printf("%f s doing Computation %f%%\n",computeTime,100*computeTime/seconds);
    printf("%f s doing Timing and looping\n",seconds-fileReadTime-fileWriteTime-computeTime);
#endif
    
    if(seconds>60)
    {
        minutes=seconds/60;
        seconds=seconds-minutes*60;
    }
    else
        minutes=0;

    if(verbose || GET_RUN_TIME)
    {
#ifndef SILENT_RUN
        if(minutes==1)
            printf("Run Time is about %d minute %f seconds.\n",minutes,seconds);
        else
            printf("Run Time is about %d minutes %f seconds.\n",minutes,seconds);

#else


#ifndef CEMENT_REMOVE
//        printf("Add %s %d %d %d %4.2f(s)\n",b_params.filename,(int)R,(int)G,(int)B,minutes*60+seconds);
            if(status!=NULL)
                status("Add",b_params.filename,R,G,B,(double)100);
#else
//        printf("Del %s %d %d %d %4.2f(s)\n",b_params.filename,(int)R,(int)G,(int)B,minutes*60+seconds);
            if(status!=NULL)
                status("Del",b_params.filename,R,G,B,(double)100);
#endif

#endif
    }

    if(isJpegFile(b)==1)
        free(image_buffer);

    sprintf(string,"rm %s; mv tmp.%s  %s",a_params.filename,a_params.filename,a_params.filename);
    if(verbose)
        printf("%s\n",string);

    if(my_system (string)==-1)
    {
        printf("Error calling %s\n",string);
        exit(EXIT_FAILURE);
    }

    return(EXIT_SUCCESS);
}


/**************************************************************************/

static inline double createValue(double a, unsigned char b, double weight)
{
#ifndef CEMENT_REMOVE
    return(a+weight*powLookupTable[b]);
#else
    return(a-weight*powLookupTable[b]);
#endif
}

/**************************************************************************/


static int my_system (char *command)
{
    int pid, status;

    if (command == 0)
        return 1;
    pid = fork();
    if (pid == -1)
        return -1;
    if (pid == 0)
    {
        char *argv[4];
        argv[0] = "sh";
        argv[1] = "-c";
        argv[2] = command;
        argv[3] = 0;
        execve("/bin/sh", argv, environ);
        exit(127);
    }
    do {
        if (waitpid(pid, &status, 0) == -1) {
            if (errno != EINTR)
                return -1;
        } else
            return status;
    } while(1);
}

#ifndef OBJECT_COMPILE
int main(int argc, char ** argv)
{
    char * input_filenameA=NULL;
    char * input_filenameB=NULL;
    double R,G,B,yR,yG,yB;

    parse_commandline(argc, argv, &input_filenameA,&input_filenameB,&R,&G,&B,&yR,&yG,&yB);

#ifndef CEMENT_REMOVE
    return(cementi(input_filenameA, input_filenameB,R,G,B,yR,yG,yB,argv[0],mystatus));
#else
    return(cemento(input_filenameA, input_filenameB,R,G,B,yR,yG,yB,argv[0],mystatus));
#endif
}

void mystatus(char *action, char *filename, double R,double G, double B, double complete)
{
    printf("%s %s %d %d %d %d%%\r",action,filename,(int)R,(int)G,(int)B,(int)complete);
    fflush(stdout);
    if((int)complete==100)
      printf("\n");
}



void usage(char * string)
{
    fprintf(stderr, "Use: %s [Lightspace file] [Imagespace file] <R> <G> <B>\n",string);
    fprintf(stderr, " [Lightspace file] is Portable Lightspace Map (plm) file image to be added to\n");
    fprintf(stderr, " [Imagespace file] is ppm file of image to add to the lightspace.\n");
    fprintf(stderr, " <R>=Float for Red weight.\n");
    fprintf(stderr, " <G>=Float for Green weight.\n");
    fprintf(stderr, " <B>=Float for Blue weight.\n");
    fprintf(stderr, " <-y>=Convert ppm input to pgm before cementing.\n");
    exit(EXIT_SUCCESS);
}


void parse_commandline(int argc, char ** argv,
                       char**input_filenameA,
                       char**input_filenameB,
                       double *R,
                       double *G,
                       double *B,
                       double *yR,
                       double *yG,
                       double *yB)
{
  if (argc<6) 
    usage(argv[0]);

  *input_filenameA=(char*)strdup(argv[1]);
  *input_filenameB=(char*)strdup(argv[2]);
  
  sscanf(argv[3],"%lf",R); /* sscanf does not recognize %g, must be %lf */
  sscanf(argv[4],"%lf",G); /* sscanf does not recognize %g, must be %lf */
  sscanf(argv[5],"%lf",B); /* sscanf does not recognize %g, must be %lf */

  if(argc==10) {
    if (strcmp("-y", argv[6]) == 0)
      {
	convert = 1;
	printf("averaging color image to get grayscale\n");
	sscanf(argv[7],"%lf",yR); /* sscanf does not recognize %g, must be %lf */
	sscanf(argv[8],"%lf",yG); /* sscanf does not recognize %g, must be %lf */
	sscanf(argv[9],"%lf",yB); /* sscanf does not recognize %g, must be %lf */
      }
    else
      usage(argv[0]);
  }
}


/* ===================================================================== */
/* ===================================================================== */
/* ===================================================================== */
/* Readlut - reads a given lookup table. */

void readlut(char *lutfilename)
{

    // some variables used in reading a text file lut
    FILE *lutfile;
    int lutsize=256;
    float p;
    int i;

    if(lutfilename !=NULL) 
      {
        if ((lutfile = fopen(lutfilename, "r")) == NULL)
	  {
            fprintf(stderr, "Unable to open %s.\n", lutfilename);
	    fclose(lutfile);
            exit(EXIT_FAILURE);
	  } 
	else 
	  {
	    fscanf(lutfile, "%f", &p);
	    for (i=0; i<lutsize; i++)
	      fscanf(lutfile, "%lf", &(powLookupTable[i]));
		   
	    if (i != lutsize) {
	      fprintf(stderr, "Incomplete read of LUT. Numread: %d\n", i);
	      fclose(lutfile);
	      exit(1);
	    }
	  }
	fclose(lutfile);
      }
}

#endif