#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.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],s[LINEPIY][LINEPIX];
  FILE *fpd;
  char buf[256];
  int i1,i2,nl=9,nb,nt,ok;
	float dmy;

  if(argc!=5 && argc!=6){
    fprintf(stderr,"Usage : %s input_file pair_tbl output_file align_param_file\n",argv[0]);
		fprintf(stderr,"        or\n");
    fprintf(stderr,"        %s input_file pair_tbl output_file cut_bottom cut_top\n",argv[0]);
    exit(0);
  }
	if(argc==6){
		sscanf(argv[4],"%d",&nb);
		sscanf(argv[5],"%d",&nt);
	}else{
		if((fpd=fopen(argv[4],"r"))==NULL){
			fprintf(stderr,"Can not open file: %s\n",argv[4]);
			exit(1);
		}
		nb=0;nt=-1;
		while(!feof(fpd)){
			if(fgets(buf,sizeof(buf),fpd)==NULL) break;
			if(sscanf(buf,"%f %f %d",&dmy,&dmy,&ok)!=3){
				printf("Wrong format\n");
				exit(1);
			}
			if(ok==0){
				if(nt<0)nb++;
				else nt++;
			}
			else nt=0;
		}
		fclose(fpd);	
	}
	fprintf(stderr,"cut_bottom=%d cut_top=%d\n",nb,nt);

  fprintf(stderr,"Read %s  ... ",argv[1]);
  status=0;
  if(fits_open_file(&fp,argv[1], 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 && 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++) z[i][j]=buffer[j];
    fpixel+=naxes[0];
  }
  if(fits_close_file(fp, &status)) printerror(status);
  fprintf(stderr,"OK\n");
  nl=naxes[1]*9/LINEPIY;

  for(i=0;i<naxes[1];i++){
    for(j=0;j<naxes[0];j++) s[i][j]=-99999.;
  }

  if((fpd=fopen(argv[2],"r"))==NULL){
    fprintf(stderr,"Can not open file: %s\n",argv[2]);
    exit(1);
  }
  while(!feof(fpd)){
    if(fgets(buf,sizeof(buf),fpd)==NULL) break;
    if(strncmp(buf,"#",1)==0) continue;
    if(sscanf(buf,"%d %d",&i1,&i2)!=2) continue;
    if(i1<=nb || i1>200-nt) continue;
    for(i=0;i<nl;i++){
      for(j=0;j<naxes[0];j++){
	if(i2<=nb || i2>200-nt) z[(i2-1)*nl+i][j]=-99999.;
	if(fabs(z[(i1-1)*nl+i][j]+99999) > 0.01) s[(i2-1)*nl+i][j]=(5-nl)/4*z[(i1-1)*nl+i][j];
	if(i1==i2) z[(i2-1)*nl+i][j]=-99999.;
      }
    }
  }
  fclose(fpd);

  remove(argv[1]);
  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);
  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");

  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];
  for(i=0;i<naxes[1];i++){
    if (fits_write_img(fp, TFLOAT, fpixel, nelements, s[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 */
}