#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,naxes0,naxes1;
  long naxes[2]={LINEPIX, LINEPIY},nbuffer;
  float nullval;
  static float buffer[LINEPIX],z[LINEPIY][LINEPIX];
  int iw,ndat,n,nframe,i1,i2,nb,nt,ok,y='y';
  double z2[LINEPIY],zav1,zsg1,zav0,zsg0,zs1,zs2,convfac;
  FILE *fpd;
  char buf[256],fname[128],dmy[16];

  if(argc<2 || argc>5){
    fprintf(stderr,"Usage : %s input_file [maskedge_file [pair_tbl [align_param_file]]]\n",argv[0]);
    fprintf(stderr,"        The name of the output file is 'input_file_prefix'_nosq.fits\n");
    exit(0);
  }
  
  if(strncmp(argv[1],"irs2",4)==0) convfac=4.;
  else convfac=2.;

  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");
  iw=naxes[1]*9/LINEPIY;
  naxes0=naxes[0];
  naxes1=naxes[1];

  sscanf(argv[1],"%s",fname);
  ok=strrchr(fname,y)-fname;
  fname[ok]='\0';
  sprintf(fname,"%sxpb.tbl",fname);
  if((fpd=fopen(fname,"r"))==NULL){
    fprintf(stderr,"##### Can not open file: %s #####\n",fname);
    fprintf(stderr,"##### The Number of the combined frames is assumed ONE. #####\n");
    nframe=1;
  }else{
    fprintf(stderr,"Read %s ...\n",fname);
    nframe=0;
    while(!feof(fpd)){
      if(fgets(buf,sizeof(buf),fpd)==NULL) break;
      if(strncmp(buf,"#",1)==0) continue;
      if(strlen(buf)<=1) continue;
      nframe++;
    }
    fclose(fpd);
    if(nframe==0){
      fprintf(stderr,"##### No lines are contained in the file: %s #####\n",fname);
      fprintf(stderr,"##### The Number of the combined frames is assumed ONE. #####\n");
      nframe=1;
    }
  }

  for(j=0;j<naxes[0];j++){
    if(j%100==0) fprintf(stderr,"*");
    zav0=0.;zsg0=10000.;
    for(n=0;n<10;n++){
      zs1=0.;zs2=0.;ndat=0;
      for(i=50;i<naxes[1]-50;i++){
	if(fabs(z[i][j]-zav0)<3*zsg0 && fabs(z[i][j])>1e-7){
	  zs1+=z[i][j];
	  zs2+=z[i][j]*z[i][j];
	  ndat++;
	}
      }
      if(ndat==0){zav0=0.;zsg0=0.;}
      else{
	zav0=zs1/ndat;
	zsg0=sqrt(zs2/ndat-zav0*zav0);
      }
    }
    //if(j==1228){fprintf(stderr,"%d(%d,%g,%g) ",j,ndat,zav0,zsg0);getchar();}
    if(ndat<200){
      for(i=0;i<naxes[1];i++){
	if(fabs(z[i][j]-zav0)<3*zsg0) z[i][j]=zsg0*zsg0;
	else z[i][j]=(fabs(z[i][j]-zav0)-3*zsg0)/convfac/nframe+zsg0*zsg0;
      }
    }else{
      for(i=0;i<naxes[1];i++){
	zs1=0.;zs2=0.;ndat=0;n=0;
	while(ndat<200 && n<naxes[1]-100){
	  if(i-n>=50){
	    //if(j==1228 && i==1621)fprintf(stderr,"%g %g %g\n",z[i-n][j],zav0,3*zsg0);
	    if(fabs(z[i-n][j]-zav0)<3*zsg0){
	      zs1+=z[i-n][j];
	      zs2+=z[i-n][j]*z[i-n][j];
	      ndat++;
	    }
	  }
	  if(n!=0 && i+n<naxes[1]-50){
	    //if(j==1228 && i==1621)fprintf(stderr,"%g %g %g\n",z[i+n][j],zav0,3*zsg0);
	    if(fabs(z[i+n][j]-zav0)<3*zsg0){
	      zs1+=z[i+n][j];
	      zs2+=z[i+n][j]*z[i+n][j];
	      ndat++;
	    }
	  }
	  //if(j==1228 && i==1621){fprintf(stderr,"%d %d %d %d %g %g\n",j,i,n,ndat,zs1,zs2);getchar();}
	  n++;
	}
	zav1=zs1/ndat;
	zsg1=sqrt(zs2/ndat-zav1*zav1);
	if(fabs(z[i][j]-zav1)<3*zsg1) z2[i]=zsg1*zsg1;
	else z2[i]=(fabs(z[i][j]-zav1)-3*zsg1)/convfac/nframe+zsg1*zsg1;
      }
      for(i=0;i<naxes[1];i++) z[i][j]=z2[i];
    }
  }
  fprintf(stderr,"\n");

  if(argc>2){
    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]!=naxes0 || naxes[1]!=naxes1){
      fprintf(stderr,"\nImage size of %s is different from that of %s.\n",argv[2],argv[1]);
      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]*buffer[j]);
      fpixel+=naxes[0];
    }
    if(fits_close_file(fp, &status)) printerror(status);
    fprintf(stderr,"OK\n");
  }

  if(argc==5){
    if((fpd=fopen(argv[4],"r"))==NULL){
      fprintf(stderr,"Can not open file: %s\n",argv[4]);
      exit(1);
    }
    fprintf(stderr,"Read %s ...\n",argv[4]);
    nb=0;nt=-1;
    while(!feof(fpd)){
      if(fgets(buf,sizeof(buf),fpd)==NULL) break;
      if(sscanf(buf,"%s %s %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);
  }else{nb=0;nt=0;}

  if(argc>3){
    if((fpd=fopen(argv[3],"r"))==NULL){
      fprintf(stderr,"Can not open file: %s\n",argv[3]);
      exit(1);
    }
    fprintf(stderr,"Read %s ...\n",argv[3]);
    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<iw;i++){
	for(j=0;j<naxes[0];j++){
	  if(i2>nb && i2<=200-nt && i1!=i2) z[(i2-1)*iw+i][j]=(z[(i1-1)*iw+i][j]+z[(i2-1)*iw+i][j])/4;
	}
      }
    }
    fclose(fpd);
  }

  sscanf(argv[1],"%s",fname);
  fname[strlen(argv[1])-5]='\0';
  sprintf(fname,"%s_nosq.fits",fname);
  remove(fname);
  fprintf(stderr,"Write %s  ... ",fname);
  status=0;
  if (fits_create_file(&fp, fname, &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");
  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 */
}