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

#define LINEPIX 2100

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, LINEPIX},nbuffer;
  float nullval;
  static float buffer[LINEPIX],z[LINEPIX][LINEPIX];
  double zav[LINEPIX][2],ave,std,sum,sum2,bs[LINEPIX/2];
  int k,p,n,n0,loop,m[LINEPIX][LINEPIX];


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

  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]>LINEPIX)){
    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(loop=0;loop<5;loop++){
    for(j=0;j<naxes[0];j++){
      for(i=0;i<naxes[1];i++) m[i][j]=1;
    }
    for(j=0;j<naxes[0];j++){
      for(i=0;i<naxes[1];i++){
	if(i>=6 && i<=naxes[1]-6 && (z[i-1][j]+z[i][j]+z[i+1][j])/3>(z[i-5][j]+z[i-4][j]+z[i-3][j]+z[i+3][j]+z[i+4][j]+z[i+5][j])/6+50.){
	  for(p=-4;p<=4;p++) m[i+p][j]=0;
	}
      }
    }

    for(i=0;i<naxes[1];i++){
      for(k=0;k<2;k++){
	ave=0;std=1e5;n=100;n0=0;
	while(n!=n0 && n>50){
	  n0=n;
	  sum=0;sum2=0;n=0;
	  for(j=naxes[0]/2*k;j<naxes[0]/2*(k+1);j++){
	    if(z[i][j] < ave + 1.5*std && m[i][j]){
	      sum+=z[i][j];
	      sum2+=z[i][j]*z[i][j];
	      n++;
	    }
	  }
	  if(n==0){
	    ave=-999.;std=-999.;
	    fprintf(stderr,"No available data (%d %d)\n",i,k);
getchar();
	    break;
	  }else{
	    ave=sum/n;
	    std=sqrt(sum2/n-ave*ave);
	  }
	}
	zav[i][k]=ave;
	//fprintf(stderr,"%d %d %g %g\n",loop,n,ave,std);
      }
    }

    sum=0;sum2=0;n=0;
    for(i=0;i<naxes[1]/2;i++){
      bs[i]=(zav[naxes[1]/2-i-1][1]-zav[naxes[1]/2-i-1][0]+zav[naxes[1]/2+i][0]-zav[naxes[1]/2+i][1])/2;
      sum+=bs[i];sum2+=bs[i]*bs[i];n++;
      //fprintf(stderr,"%d %g\n",i,bs[i]);
    }

    ave=sum/n;
    std=sqrt(sum2/n-ave*ave);
    fprintf(stderr,"%g %g\n",ave,std);

    for(i=0;i<naxes[1]/2;i++){
      for(j=0;j<naxes[0]/2;j++){
	z[i][j]-=bs[naxes[0]/2-j-1];
	z[i][j+naxes[0]/2]-=bs[naxes[1]/2-i-1];
	z[i+naxes[1]/2][j]-=bs[i];
	z[i+naxes[1]/2][j+naxes[0]/2]-=bs[j];
      }
    }
  }

  //for(i=0;i<naxes[1];i++){
  //  for(j=0;j<naxes[0];j++) z[i][j]*=m[i][j];
  //}

  if(strcmp(argv[1],argv[2])==0) remove(argv[1]);
  fprintf(stderr,"Write %s  ... ",argv[2]);
  status=0;
  if (fits_create_file(&fp, argv[2], &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 */
}