#include <math.h>
#include <sys/types.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <sys/wait.h>        //wait
#include <sys/shm.h>         //shm
#include <fitsio.h>
#include <malloc.h>

#define EXP_TIMER_RESOURCE  "EXP_TIMER"
#define SHM_IMAGE_RESOURCE  "SHM_IMAGE"

#define SATULATE 65000 //satulation value
#define MARK_SATULATE -2000000000
#define MARK_BAD      -2100000000
#define BAD_PIXEL     0       //sigma = 0
#define NBYTE (2048 * 2048 * sizeof(int))
#define NPIXEL 2048
#define COSMIC_RAY 200
#define ZSHIFT 1000000

void printerror(int status);

int main(int argc, char *argv[] ){
  int i, j;
  double sec_exp;
  int times;

  char command[256], buf[256], comment[256];
  fitsfile *fp0;
  FILE *in_pipe;
  //FLIE *fp;
  int read_times;

  double sum_xi, sum_xi2, sigma;

  static int imbuf1[NPIXEL][NPIXEL],imbuf2[NPIXEL][NPIXEL],imbuf3[NPIXEL][NPIXEL];
  static double imbufd[NPIXEL][NPIXEL];
  long naxes[2], fpixel, naxis = 2;
  int status, nfound, anynull;

  int nullval; 
  int buffer[NPIXEL];

  ///cosmic ray
  static int cr[NPIXEL][NPIXEL];
  int fitcomp, tmpbuf, ncr;
  double a, b, stddev;

  if( argc != 3 && argc != 4 ){
    fprintf( stderr, "Usage: %s imputimage(FMSA00001235) outimage.fits\n", argv[0]);
    fprintf( stderr, "   or if you need CR rejected ramp images...\n");
    fprintf( stderr, "Usage: %s imputimage(FMSA00001235) outimage.fits suffix\n", argv[0]);
    exit(1);
  }

  ///get sample number
  sprintf(command, "ls %s_???.fits | wc -l", argv[1]);
  if((in_pipe = popen(command, "r")) == NULL){
    fprintf(stderr, "ERROR: %s\n", command);
   exit(1);
  }
  fgets(buf, 256, in_pipe);
  sscanf(buf, "%d", &read_times);
  read_times--;
  fprintf(stderr,"ramp sample set: 000 -- %03d\n", read_times);
  
  pclose(in_pipe);

  for(j = 0; j < NPIXEL; j++){
    for(i = 0; i < NPIXEL; i++){
      imbuf1[i][j] = 0;
      imbuf2[i][j] = 0;
      imbufd[i][j] = 0.;
      cr[i][j] = 0;
    }
  }
  //initiallization
  sum_xi=0.0;
  sum_xi2=0.0;
  sigma=0.0;

  for(times = 0; times<= read_times; times++){
    sprintf(buf, "%s_%03d.fits", argv[1], times);
    fprintf(stderr,"Read %s ... ",buf);
    status = 0;
    if(fits_open_file(&fp0, buf, READONLY, &status))
      printerror(status);
    if(fits_read_keys_lng(fp0, "NAXIS", 1, 2, naxes, &nfound, &status)) 
      printerror(status);
    fprintf(stderr, "%ld x %ld format\n", naxes[0], naxes[1]);
    if((naxes[0] != NPIXEL) || (naxes[1] != NPIXEL)){
      fprintf(stderr,"\nIllegal image format: %ld(%d)  %ld(%d)\n",naxes[0], NPIXEL, naxes[1], NPIXEL);
      exit(1);
    }
    if(fits_read_keyword(fp0, "EXPOSURE", buf, comment, &status))
      printerror(status);
    sscanf(buf, "%lf", &sec_exp);
    fpixel = 1;
    nullval = 0;
    for(j = 0;j < naxes[1];j++){
      if (fits_read_img(fp0, TINT, fpixel, naxes[0], &nullval, buffer, &anynull, &status))
	printerror(status);
      for(i = 0;i < naxes[0];i++){
	imbuf3[i][j] = buffer[i];
      }      
      fpixel += naxes[0];
    }
    if(fits_close_file(fp0, &status))
      printerror(status);
    
    ///fitting
    for(j = 0; j < naxes[1]; j++){
      for(i = 0; i < naxes[0]; i++){
	if(imbuf2[i][j] > MARK_SATULATE){
	  a = ((double)(imbuf2[i][j]) * times - imbuf1[i][j] * sum_xi) / sigma;
	  b = (imbuf1[i][j] - a * sum_xi) / times;
	  if (imbuf3[i][j] < SATULATE){
	    fitcomp = (int)(a * sec_exp + b);
	    //if(i == 1676 && j == 1754){
	    //fprintf(stdout,"%d %d %d %d %d\n", imbuf3[i][j], fitcomp, cr[i][j], imbuf3[i][j] - fitcomp, (int)(a * sec_exp));
	    //}
	    if(times > 3 && cr[i][j] <= ZSHIFT/2){
	      stddev = sqrt((imbufd[i][j]-2*a*imbuf2[i][j]-2*b*imbuf1[i][j]+a*a*sum_xi2+2*a*b*sum_xi)/times+b*b);
	      if(times > 4 && abs(imbuf3[i][j] - fitcomp) > COSMIC_RAY && abs(imbuf3[i][j] - fitcomp - cr[i][j]) > COSMIC_RAY && abs(imbuf3[i][j] - fitcomp) > stddev * 10 && abs(imbuf3[i][j] - fitcomp - cr[i][j]) > stddev * 10){
		cr[i][j] = imbuf3[i][j] - fitcomp + ZSHIFT;
		//fprintf(stdout,"%d %d %d %d %d %d %d %d %d %d\n", times, j, i, (int)(a * sec_exp), (int)b, imbuf3[i][j], fitcomp, cr[i][j], imbuf3[i][j] - fitcomp, (int)stddev);
	      } else{
		cr[i][j] = imbuf3[i][j] - fitcomp;
	      }
	    }
	    //if(i==0 && j==0)getchar();
	    if(cr[i][j] > ZSHIFT/2) tmpbuf = imbuf3[i][j] - cr[i][j] + ZSHIFT;
	    else tmpbuf = imbuf3[i][j];
	    imbuf1[i][j] += tmpbuf;  //sum_yi
	    imbuf2[i][j] += (int)(tmpbuf * sec_exp);  //sum_xi*yi
	    imbufd[i][j] += (double)tmpbuf * tmpbuf;  //sum_yi*yi
	  } else{
	    //Satulate:estimate final value from previous data
	    //To keep resolution, set exptime as 10000sec
	    if(sigma < 1e-7){
	      imbuf2[i][j] = MARK_BAD;
	    }else{
	      imbuf1[i][j] = (int)(a * 10000.0);
	      imbuf2[i][j] = MARK_SATULATE;
	    }
	  }
	}
      }
    }
    sum_xi += sec_exp; //sum_xi
    sum_xi2 += sec_exp * sec_exp; //sum_xi^2
    sigma = (double)(times + 1) * sum_xi2 - sum_xi * sum_xi; ///sigma
    ///    fprintf(stderr,"%d %lf %d %d %lf %lf %lf %d\n", imbuf3[1][1], sec_exp, imbuf1[1][1], imbuf2[1][1], sum_xi, sum_xi2, sigma, cr[1][1]);

    if(argc == 4){
      sprintf(buf, "%s_%03d%s.fits", argv[1], times, argv[3]);
      fprintf(stderr,"Write %s ... ",buf);
      status=0;
      if (fits_create_file(&fp0, buf, &status))
	printerror(status);
      if (fits_create_img(fp0, LONG_IMG, naxis, naxes, &status))
	printerror(status);
      fpixel=1;
      fprintf(stderr,"%d x %d format ", NPIXEL, NPIXEL);
      for(j = 0; j < naxes[1]; j++){
	for(i = 0;i < naxes[0];i++){
	  if(cr[i][j] > ZSHIFT/2) buffer[i] = imbuf3[i][j] - cr[i][j] + ZSHIFT;
	  else buffer[i] = imbuf3[i][j];
	}
	if(fits_write_img(fp0, TINT, fpixel, naxes[0], buffer, &status))
	  printerror(status);
	fpixel += naxes[0];
      }
      if (fits_close_file(fp0, &status))
	printerror(status);
    fprintf(stderr,"OK\n");
    }
  }

    ///final calclation
  if(sigma < 1e-5){
    sigma = 1e-5;
  }
  ncr=0;

  for(j = 0; j < naxes[1]; j++){
    for(i = 0; i < naxes[0]; i++){
      if (imbuf2[i][j] == MARK_BAD){
	imbuf3[i][j] = BAD_PIXEL;
      }else if(imbuf2[i][j] == MARK_SATULATE){
	//convert 10000sec to sec_exp
	imbuf3[i][j] = (int)(imbuf1[i][j] * sec_exp / 10000.0);
      }else{
	//estimate count for sec_exp
	imbuf3[i][j] = (int)(((double)(imbuf2[i][j]) * times - imbuf1[i][j] * sum_xi)* sec_exp / sigma);
      }
      if(cr[i][j] > ZSHIFT/2){
	cr[i][j]-= ZSHIFT;
	ncr++;
      }
      else cr[i][j]=0;
    }
  }
  ///  fprintf(stderr,"%d %lf %lf %lf\n", imbuf3[1][1], (double)(imbuf2[1][1]) * times - (double)(imbuf1[1][1] * sum_xi), (double)(imbuf2[1][1] )* times, imbuf1[1][1] * sum_xi);
  fprintf(stderr,"CR rejected : %d pixels\n",ncr);

  ///write fitted data
  fprintf(stderr,"Write %s ... ",argv[2]);
  status=0;
  if (fits_create_file(&fp0, argv[2], &status))
    printerror(status);
  if (fits_create_img(fp0, LONG_IMG, naxis, naxes, &status))
    printerror(status);
  fpixel=1;
  fprintf(stderr,"%d x %d format ", NPIXEL, NPIXEL);
  for(j = 0; j < naxes[1]; j++){
    for(i = 0;i < naxes[0];i++){
      buffer[i]= imbuf3[i][j];
    }
    if(fits_write_img(fp0, TINT, fpixel, naxes[0], buffer, &status))
      printerror(status);
    fpixel += naxes[0];
  }
  if (fits_close_file(fp0, &status))
    printerror(status);
  fprintf(stderr,"OK\n");

  ///write cr.fits ... number for cr detected
  remove("cr.fits");
  fprintf(stderr,"Write cr.fits ... ");
  status=0;
  if (fits_create_file(&fp0, "cr.fits", &status))
    printerror(status);
  if (fits_create_img(fp0, LONG_IMG, naxis, naxes, &status))
    printerror(status);
  fpixel=1;
  fprintf(stderr,"%d x %d format ", NPIXEL, NPIXEL);
  for(j = 0; j < naxes[1]; j++){
    for(i = 0;i < naxes[0];i++){
      buffer[i]= cr[i][j];
    }
    if(fits_write_img(fp0, TINT, fpixel, naxes[0], buffer, &status))
      printerror(status);
    fpixel += naxes[0];
  }
  if (fits_close_file(fp0, &status))
    printerror(status);
  fprintf(stderr,"OK\n");

  return 0;
 
}


void printerror(int status){
  /*****************************************************/
  /* Print out cfitsio error messages and exit program */
  /*****************************************************/

  char status_str[FLEN_STATUS], errmsg[FLEN_ERRMSG];
  
  if (status)
    fprintf(stderr, "\n*** Error occurred during program execution ***\n");

  fits_get_errstatus(status, status_str);   /* get the error description */
  fprintf(stderr, "\nstatus = %d: %s\n", status, status_str);

  /* get first message; null if stack is empty */
  if ( fits_read_errmsg(errmsg) ){
    fprintf(stderr, "\nError message stack:\n");
    fprintf(stderr, " %s\n", errmsg);
    
    while ( fits_read_errmsg(errmsg) )  /* get remaining messages */
      fprintf(stderr, " %s\n", errmsg);
  }

  exit( status );       /* terminate the program, returning error status */
}