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

#define LINEPIX 2560
#define LINEPIY 1920

void printerror( int status );

int main(int argc, char *argv[])
{
  fitsfile *fp;
  int i,j,k,status,nfound,anynull;
  long fpixel;
  long naxes[2]={LINEPIX, LINEPIY},nbuffer;
  float nullval;
  static long z0[LINEPIY][LINEPIX];
  static int buffer[LINEPIX],z[LINEPIY][LINEPIX];
  char buf[128],filename[32];
  FILE *fpd;
  unsigned long zamp,zampmin,zsum;
  int dj,di,djmin,dimin,ii,jj,dis,djs;

  if(argc!=2){
    fprintf(stderr,"Usage : %s input_list\n",argv[0]);
    exit(0);
  }

  if((fpd=fopen(argv[1],"r"))==NULL){
    fprintf(stderr,"Can not open file: %s\n",argv[1]);
    exit(1);
  }

  for(i=0;i<LINEPIY;i++){
    for(j=0;j<LINEPIX;j++) z0[i][j]=0;
  }
  k=0;
  while(!feof(fpd)){
    if(fgets(buf,sizeof(buf),fpd)==NULL) break;
    if(sscanf(buf,"%s",filename)!=1) break;
    //fprintf(stderr,"Read %s  ... ",filename);
    status=0;
    if(fits_open_file(&fp,filename, READONLY, &status)){
      fprintf(stderr,"Failed! Exit this process...\n");
      exit(0);
    }
    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)||(naxes[1]>=LINEPIY)){
      fprintf(stderr,"\nIllegal image format\n");
      exit(0);
    }
    fpixel = 1;
    nullval = 0;
    nbuffer = naxes[0];
    for(i=0;i<naxes[1];i++){
      if (fits_read_img(fp, TINT, fpixel, nbuffer, &nullval,
			buffer, &anynull, &status))
	printerror(status);
      for(j=0;j<naxes[0];j++) z0[i][j]+=buffer[j];
      fpixel+=naxes[0];
    }
    if(fits_close_file(fp, &status))
      printerror(status);
    k++;
    //fprintf(stderr,"OK(%d)\n",k);
  }
  fseek(fpd, 0, SEEK_SET);

  zsum=0;
  for(i=0;i<naxes[1];i++){
    for(j=0;j<naxes[0];j++){
      z0[i][j]/=k;
      zsum+=z0[i][j];
    }
  }

  /*
  fprintf(stderr,"Write sum.fits  ... ");
  status=0;
  if (fits_create_file(&fp, "sum.fits", &status))
    printerror(status);
  if (fits_create_img(fp, 16, 2, naxes, &status))
    printerror(status);
  fpixel=1;
  for(i=0;i<naxes[1];i++){
    for(j=0;j<naxes[0];j++) buffer[j]=z0[i][j];
    if (fits_write_img(fp, TINT, fpixel, naxes[0], buffer, &status))
      printerror(status);
    fpixel+=naxes[0];
  }
  if (fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK\n");
  exit(0);
  */

  while(!feof(fpd)){
    if(fgets(buf,sizeof(buf),fpd)==NULL) break;
    if(sscanf(buf,"%s",filename)!=1) break;
    fprintf(stderr,"%s ",filename);
    status=0;
    if(fits_open_file(&fp,filename, READONLY, &status)){
      fprintf(stderr,"Failed! Exit this process...\n");
      exit(0);
    }
    if(fits_read_keys_lng(fp, "NAXIS", 1, 2, naxes, &nfound, &status))
      printerror(status);
    if((naxes[0]>=LINEPIX)||(naxes[1]>=LINEPIY)){
      fprintf(stderr,"\nIllegal image format\n");
      exit(0);
    }
    fpixel = 1;
    nullval = 0;
    nbuffer = naxes[0];
    for(i=0;i<naxes[1];i++){
      if (fits_read_img(fp, TINT, fpixel, nbuffer, &nullval,
			buffer, &anynull, &status))
	printerror(status);
      for(j=0;j<naxes[0];j++) z[i][j]=buffer[j];
      fpixel+=naxes[0];
    }
    if(fits_close_file(fp, &status))
      printerror(status);

    dis=0;djs=0;
    while(1){
      zampmin=1e10;dimin=0;djmin=0;
      for(di=-1;di<=1;di++){
      for(dj=-1;dj<=1;dj++){
	zamp=0;
	for(i=0;i<naxes[1];i++){
	for(j=0;j<naxes[0];j++){
	  ii=i+di+dis;jj=j+dj+djs;
	  if(ii<0) ii=0; else if(ii>=naxes[1]) ii=naxes[1]-1;
	  if(jj<0) jj=0; else if(jj>=naxes[0]) jj=naxes[0]-1;
	  zamp+=abs(z[ii][jj]-z0[i][j]);
	}}
	if(zampmin>zamp){
	  zampmin=zamp;
	  dimin=di;djmin=dj;
	}
	//fprintf(stderr,"%2d %2d %2d %2d zamp=%ld zampmin=%ld\n",dj+djs,di+dis,dj,di,zamp,zampmin);
      }}
      if(dimin==0 && djmin==0) break;

      zamp=zampmin;
      zampmin=1e10;
      while(zampmin>zamp){
	dis+=dimin;djs+=djmin;
	zampmin=zamp;
	zamp=0;
	for(i=0;i<naxes[1];i++){
	for(j=0;j<naxes[0];j++){
	  ii=i+dimin+dis;jj=j+djmin+djs;
	  if(ii<0) ii=0; else if(ii>=naxes[1]) ii=naxes[1]-1;
	  if(jj<0) jj=0; else if(jj>=naxes[0]) jj=naxes[0]-1;
	  zamp+=abs(z[ii][jj]-z0[i][j]);
	}}
	//fprintf(stderr,"%2d %2d %2d %2d zamp=%ld zampmin=%ld\n",djmin+djs,dimin+dis,djmin,dimin,zamp,zampmin);
      }
    }
    fprintf(stderr,"%2d %2d %g\n",-djs,-dis,(double)zampmin/zsum);
    printf("%2d %2d %g\n",-djs,-dis,(double)zampmin/zsum);
  }
  fclose(fpd);
  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 */
}