#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "comphys.h"
#include "comphys.c"

int main(int argc, char *argv[])
{
  float **array;
  int err,i,j,N;
  double sum1,sum12,av,var,sum,totalS,counts,tsum;
  char infile[80],buf[10000],buffer[10000];
  char *tmp;
  FILE *in;

  array = matrix(1,1024,1,256);

  strcpy(infile,argv[1]);

  sprintf(buf,"imlist %s > temp.dat",infile);
  err = system(buf);

  in = fopen("temp.dat","r");
  if(in == NULL) {
    printf("Cannot open temp.dat\n");
    exit(1);
  }

  fgets(buf,10000,in);
  fgets(buf,10000,in);

  for(j=256;j>=1;j--) {
    fgets(buffer,10000,in);
    tmp = strtok(buffer," ");
    for(i=1;i<=1024;i++) {
      tmp = strtok(NULL," ");
      array[i][j] = atof(tmp);
    }
  }
  fclose(in);

  /* Dark area:  X 220 - 240, Y  20 -  70
     Spectrum:   X 220 - 240, Y 110 - 160 */

  // Get mean and std deviation of dark area

  sum1 = 0.0;
  sum12 = 0.0;
  N = 0;
  for(i=220;i<=240;i++) {
    for(j=20;j<=70;j++) {
      sum1 += array[i][j];
      sum12 += array[i][j]*array[i][j];
      N++;
    }
  }

  av = sum1/(double)N;
  var = sum12/(double)N - av*av;

  //printf("mean = %f sig = %f\n",av,sqrt(var));

  // Now find total counts in spectrum area, and subtract off dark

  sum = 0.0;
  N = 0;
  for(i=220;i<=240;i++) {
    for(j=110;j<=160;j++) {
      sum += array[i][j];
      N++;
    }
  }

  totalS = sum - N*av;  // total signal in spectrum area
  counts = totalS/21.0;  // total counts per dispersion pixel

  //printf("sum = %f\n",sum);
  //printf("Total signal = %f\n",totalS);
  printf("%f\n",counts);

  /*  tsum = 0.0;
  for(j=110;j<=160;j++) {
    //printf("%d %f\n",j,array[220][j]-av);
    tsum += array[220][j]-av;
    } */


  free_matrix(array,1,1024,1,256);

  return(0);
}