#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "fitsio.h"

#define LINEPIX 2100
#define LINEPIY 200

void printerror( int status );

int main(int argc, char *argv[])
{
  fitsfile *fp;
  int i,j,status,nfound,anynull,bitpix=-32;
  long fpixel, nelements,naxis=2;
  long naxes[2]={LINEPIX, LINEPIY},nbuffer;
  float nullval;
  static float buffer[LINEPIX],z[LINEPIY][LINEPIX],f[LINEPIY][LINEPIX],zm[LINEPIY];
  int iw,nj,nn,jw,bo[LINEPIY];
  float ww,ws,we,wwmin;
  double zmav,zmsg,zamp,zampmin;
  FILE *fpd;
  char fname[128];
  
  if(argc!=4 && argc!=5){
    fprintf(stderr,"Usage : %s input_file reduced_flat output_file [-]\n",argv[0]);
    exit(0);
  }
  
  fprintf(stderr,"Read %s  ... ",argv[2]);
  status=0;
  if(fits_open_file(&fp,argv[2], READONLY, &status))
    printerror(status);
  if(fits_read_keys_lng(fp, "NAXIS", 1, 2, naxes, &nfound, &status))
    printerror(status);
  fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[1]);
  if(naxes[0]>LINEPIX){
    fprintf(stderr,"\nNumber of columns must be < %d\n",LINEPIX);
    exit(0);
  }
  if(naxes[1]!=LINEPIY){
    fprintf(stderr,"\nNumber of lines must be %d\n",LINEPIY);
    exit(0);
  }
  fpixel = 1;
  nullval = 0;
  nbuffer = naxes[0];
  for(i=0;i<naxes[1];i++){
    if (fits_read_img(fp, TFLOAT, fpixel, nbuffer, &nullval,
                      buffer, &anynull, &status))
      printerror(status);
    zm[i]=0;
    for(j=0;j<naxes[0];j++){
      if(buffer[j]<0) buffer[j]=0;
      f[i][j]=buffer[j];
      zm[i]+=buffer[j]/naxes[0];
    }
    fpixel+=naxes[0];
  }
  if(fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK\n");
  nj=naxes[0];
  for(i=0;i<naxes[1];i++){
    for(j=0;j<naxes[0];j++) f[i][j]/=zm[i];
  }

  sscanf(argv[1],"%s",fname);
  fprintf(stderr,"Read %s  ... ",fname);
  status=0;
  if(fits_open_file(&fp,fname, READONLY, &status))
    printerror(status);
  if(fits_read_keys_lng(fp, "NAXIS", 1, 2, naxes, &nfound, &status))
    printerror(status);
  fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[1]);
  if(naxes[0]!=nj){
    fprintf(stderr,"\nNumber of columns is different from that of %s{%d)\n",argv[2],nj);
    exit(0);
  }
  if(naxes[1]!=LINEPIY&&naxes[1]!=LINEPIY*9){
    fprintf(stderr,"\nNumber of lines must be %d or %d\n",LINEPIX,LINEPIY*9);
    exit(0);
  }
  iw=naxes[1]/LINEPIY;
  for(i=0;i<LINEPIY;i++){
    for(j=0;j<LINEPIX;j++) z[i][j]=0.;
  }
  fpixel = 1;
  nullval = 0;
  nbuffer = naxes[0];
  for(i=0;i<naxes[1];i++){
    if (fits_read_img(fp, TFLOAT, fpixel, nbuffer, &nullval,
                      buffer, &anynull, &status))
      printerror(status);
    for(j=0;j<naxes[0];j++) z[i/iw][j]+=buffer[j];
    fpixel+=naxes[0];
  }
  if(fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK\n");

  zmav=0.;zmsg=0.;
  for(i=0;i<LINEPIY;i++){
    zm[i]=0.;
    for(j=0;j<naxes[0];j++) zm[i]+=z[i][j]/naxes[0];
    zmav+=zm[i]/LINEPIY;
    zmsg+=zm[i]*zm[i]/LINEPIY;
  }

  zmsg=sqrt(zmsg-zmav*zmav);
  for(i=0;i<LINEPIY;i++){
    bo[i]=0;
    if(fabs(zm[i]-zmav)>zmsg*3){
      fprintf(stderr,"Bright object #%d is found.(%g %g %g %g)\n",i+1,zm[i-1]-zmav,zm[i]-zmav,zm[i+1]-zmav,zmsg*3);
      bo[i]=1;
    }
  }

  ww=-49.;ws=1;we=0.001;zampmin=1e100;nn=0;
  while(1){
    jw=(int)floor(ww);
    zamp=0.;
    for(i=0;i<LINEPIY;i++){
      if(bo[i]){
	for(j=50;j<naxes[0]-50;j++){
	  zmav=f[i][j+jw+1]*(ww-jw)+f[i][j+jw]*(1-ww+jw);
	  if(zmav<0.1) zmav=0.1;
	  buffer[j]=z[i][j]/zmav;
	  if(j>50 && f[i][j+jw]>0.5 && f[i][j+jw+1]>0.5) {
	    zamp+=fabs(buffer[j]-buffer[j-1]);
	    //if(j==903 && i==99){fprintf(stderr,"%d %d %g %g %g %g %g %g %g %g\n",i,j,buffer[j-1],buffer[j],z[i][j],zmav,f[i][j+jw],f[i][j+jw+1],(f[i][j+jw+1]*(ww-jw)+f[i][j+jw]*(1-ww+jw)),zamp);}
	  }
	}
      }
    }
    fprintf(stderr,"dx=%g zamp=%g zampmin=%g\n",ww,zamp,zampmin);
    if(fabs(ws-1)<1e-5){
      if(zampmin>zamp){
	zampmin=zamp;
	wwmin=ww;
      }
      ww+=ws;
      if(ww>49.5){
	ws/=2;
	ww=wwmin;
      }
    }else{
      if(zampmin>zamp){
	zampmin=zamp;
	ww+=ws;
      }
      else{
	ww-=ws;
	if(nn==0){
	  ws*=-1;
	  ww+=ws;
	  nn=1;
	}
	else{
	  if(fabs(ws)<we) break;
	  ws/=2;
	  ww+=ws;
	  nn=0;
	}
      }
    }
  }
  fprintf(stderr,"X-shift=%g\n",ww);
  jw=(int)ww;
  fname[strlen(argv[1])-5]='\0';
  sprintf(fname,"%s.dat",fname);
  if((fpd=fopen(fname,"w"))==NULL){
    fprintf(stderr,"Can not open file: %s\n",fname);
    exit(1);
  }
  fprintf(stderr,"Write %s\n",fname);  
  fprintf(fpd,"%f\n",ww);
  fclose(fpd);

  for(i=0;i<LINEPIY;i++){
    for(j=0;j<naxes[0];j++) buffer[j]=f[i][j];
    for(j=-jw;j<naxes[0]-jw-1;j++){
      if(j>=0 && j<naxes[0]) f[i][j]=buffer[j+jw+1]*(ww-jw)+buffer[j+jw]*(1-ww+jw);
    }
  }

  if(argc!=5){
    for(j=0;j<naxes[0];j++){
      buffer[j]=0;
      for(i=0;i<LINEPIY;i++) buffer[j]+=f[i][j]/LINEPIY;
      for(i=0;i<LINEPIY;i++){
	if(buffer[j]<1e-4) f[i][j]=1;
	else f[i][j]/=buffer[j];
      }
    }
  }

  naxes[1]=LINEPIY;
  fprintf(stderr,"Write %s  ... ",argv[3]);
  status=0;
  if (fits_create_file(&fp, argv[3], &status))
    printerror(status);
  if (fits_create_img(fp, bitpix, naxis, naxes, &status))
    printerror(status);
  fpixel=1;
  nelements=naxes[0];
  fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[1]);
  for(i=0;i<naxes[1];i++){
    if (fits_write_img(fp, TFLOAT, fpixel, nelements, f[i], &status))
      printerror(status);
    fpixel+=naxes[0];
  }
  if (fits_close_file(fp, &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 */
}