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

#define LINEPIX 2100
#define LINEPIY 1800

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],p[9][LINEPIX];
  char buf[256],fi[32],fo[32],dmy[32];
  FILE *fpd;
  int xmin,xmax,n,obj[4],sobj,k;
  double w,wtot;

  if(argc!=4){
    fprintf(stderr,"Usage : %s filelist xmin xmax\n",argv[0]);
    fprintf(stderr,"\tFormat of the filelist:\n");
    fprintf(stderr,"\tFits_file1 obj#1 [obj#2 [obj#3 [obj#4]]]\n");
    fprintf(stderr,"\tFits_file2 obj#1 [obj#2 [obj#3 [obj#4]]]\n");
    fprintf(stderr,"\t...\n");
    fprintf(stderr,"Backup file will be created as Fits_file1.bak.\n");
    exit(0);
  }
  sscanf(argv[2],"%d",&xmin);
  sscanf(argv[3],"%d",&xmax);


  if((fpd=fopen(argv[1],"r"))==NULL){
    fprintf(stderr,"Can not open file: %s\n",argv[1]);
    exit(1);
  }
  for(i=0;i<9;i++){
    for(j=xmin-1;j<xmax;j++) p[i][j]=0;
  }

  while(!feof(fpd)){
    if(fgets(buf,sizeof(buf),fpd)==NULL) break;
    n=sscanf(buf,"%s %d %d %d %d %s",fi,&obj[0],&obj[1],&obj[2],&obj[3],dmy);
    if(n<2) continue;
    else if(n>=6){fprintf(stderr,"##### Maximum 4 objects are available. ####\n");n=5;}
    //fprintf(stderr,"%d ",n);
    //for(k=0;k<n-1;k++) fprintf(stderr,"%d ",obj[k]);

    fprintf(stderr,"Read %s  ... ",fi);
    status=0;
    if(fits_open_file(&fp,fi, 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)||(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, TFLOAT, 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);
    fprintf(stderr,"OK\n");
    
    for(k=0;k<n-1;k++){
      if(obj[k]<0) sobj=-1;else sobj=1;
      for(i=0;i<9;i++){
	for(j=xmin-1;j<xmax;j++){
	  p[i][j]+=sobj*z[(obj[k]*sobj-1)*9+i][j];
	  //fprintf(stderr,"%d %d %d %d %d %g %g\n",i,j,k,obj[k],sobj,z[(obj[k]*sobj-1)*9+i][j],p[i][j]);getchar();
	}
      }
    }
  }
  wtot=0;
  for(i=0;i<9;i++){
    for(j=xmin-1;j<xmax;j++) wtot+=p[i][j];
  }
  fprintf(stderr,"Total weight=%g\n",wtot);
  fseek(fpd, 0, SEEK_SET);

  while(!feof(fpd)){
    if(fgets(buf,sizeof(buf),fpd)==NULL) break;
    n=sscanf(buf,"%s %d %d %d %d %s",fi,&obj[0],&obj[1],&obj[2],&obj[3],dmy);
    if(n<2) continue;
    else if(n>=6) n=5;

    fprintf(stderr,"Read %s  ... ",fi);
    status=0;
    if(fits_open_file(&fp,fi, 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)||(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, TFLOAT, 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);
    fprintf(stderr,"OK\n");
    
    for(k=0;k<n-1;k++){
      if(obj[k]<0) sobj=-1;else sobj=1;
      w=0.;
      for(i=0;i<9;i++){
	for(j=xmin-1;j<xmax;j++){
	  w+=z[(obj[k]*sobj-1)*9+i][j];
	}
      }
      fprintf(stderr,"%s Obj#%d: weight=%g\n",fi,obj[k]*sobj,w/wtot);
      for(i=0;i<9;i++){
	for(j=xmin-1;j<xmax;j++){
	  z[(obj[k]*sobj-1)*9+i][j]-=p[i][j]*w/wtot;
	}
      }
    }
    sprintf(fo,"%s.bak",fi);
    remove(fo);
    rename(fi,fo);
    fprintf(stderr,"Write %s  ... ",fi);
    status=0;
    if (fits_create_file(&fp, fi, &status))
      printerror(status);
    if (fits_create_img(fp, bitpix, naxis, naxes, &status))
      printerror(status);
    fpixel=1;
    nelements=naxes[0];
    for(i=0;i<naxes[1];i++){
      if (fits_write_img(fp, TFLOAT, fpixel, nelements, z[i], &status))
	printerror(status);
      fpixel+=naxes[0];
    }
    if (fits_close_file(fp, &status))
      printerror(status);
    fprintf(stderr,"OK\n");
  }
  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 */
}