#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,k,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];
  FILE *fps;
  int p[200],p1,p2;
  char buf[64];
  float xs;

  if(argc!=3 && argc!=4){
    fprintf(stderr,"Usage : %s input reference [pair_tbl]\n",argv[0]);
    exit(0);
  }

  for(i=0;i<200;i++) p[i]=-1;
  if(argc==4){
    if((fps = fopen(argv[3],"r")) == NULL){
      printf("Cannot open %s\n",argv[3]);
      argc=3;
    }
  }
  if(argc==4){
    while(!feof(fps)){
      if(fgets(buf,sizeof(buf),fps)==NULL) break;
      if(strncmp(buf,"#",1)==0) continue;
      sscanf(buf,"%d %d",&p1,&p2);
      p[p1-1]=p2-1;
      p[p2-1]=p1-1;
    }
    fclose(fps);
  }

  status=0;
  if(fits_open_file(&fp,argv[1], READONLY, &status)){
    fprintf(stderr,"Make %s\n",argv[1]);
    for(i=0;i<LINEPIY;i++){
      for(j=0;j<LINEPIX;j++) z[i][j]=1.;
    }
  }else{
    fprintf(stderr,"Read %s  ... ",argv[1]);
    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");
  }

  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)||(naxes[1]>LINEPIY*9)){
    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++){
      if(buffer[j]<-99990){
	xs=0.;
	for(k=0;k<naxes[0];k++){
	  if(buffer[k]>-99990) xs+=fabs(buffer[k]);
	}
	//fprintf(stderr,"%d %d %g %d %d\n",i,j,xs,i/9,p[i/9]);
	if(xs<1e-5){
	  for(k=0;k<naxes[0];k++) z[i/9][k]=1.;
	  if(p[i/9]>=0){
	    for(k=0;k<naxes[0];k++) z[p[i/9]][k]=1.;
	  }
	}else{
	  for(k=0;k<naxes[0];k++) z[i/9][k]=0.;
	  if(p[i/9]>=0){
	    for(k=0;k<naxes[0];k++) z[p[i/9]][k]=0.;
	  }
	}
      }
    }
    fpixel+=naxes[0];
  }
  if(fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK\n");


  remove(argv[1]);
  naxes[1]/=9;
  fprintf(stderr,"Write %s  ... ",argv[1]);
  status=0;
  if (fits_create_file(&fp, argv[1], &status))
    printerror(status);
  if (fits_create_img(fp, bitpix, naxis, naxes, &status))
    printerror(status);
  fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[1]);
  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");
  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 */
}