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

#define LINEPIX 2048
#define SATULATE 60000
#define MARK_SATULATE -2000000000
#define MARK_BAD      -2100000000
#define BAD_PIXEL     0
#define COSMIC_RAY    100
#define ZSHIFT        1000000

void printerror(int status);

int main(int argc,char *argv[])
{
  fitsfile *fp,*fps;
  int i,j,status,nfound,anynull,bitpix=-32;
  long fpixel, nelements,npixels,naxis=2;
  long naxes[2]={LINEPIX, LINEPIX},nbuffer;
  float nullval; 
  static unsigned short buffer[LINEPIX],buffers[LINEPIX];
  static int zs[LINEPIX][LINEPIX],xzs[LINEPIX][LINEPIX],cr[LINEPIX][LINEPIX];
  static double zzs[LINEPIX][LINEPIX];
  int ihdu,ihdumax,ihdumaxs,ihdus,itype,sx,sx2,sigma,zex,zre,irs,star=0,scr=0,ncr;
  char buf[256],fname[3][128],*ret,mod[8],mod2[4];
  double a,b,stddev;
  FILE *fpt;

  if(argc!=4 && argc!=5){
    fprintf(stderr,"Usage: %s input input_sky output [starimage or js/jl/hs/hl/lr/...]\n",argv[0]);
    exit(0);
  }

  if((fpt=fopen(argv[3],"rb"))!=NULL){
    fprintf(stderr,"Output file (%s) exists.\n",argv[3]);
    fclose(fpt);
    exit(0);
  }

  if(argc==5){
    if(strstr(argv[4], ".fits") == NULL){
      strcpy(mod,argv[4]);
      if(strncmp(argv[4],"lr",2)==0) strcpy(mod2,"s");
      else strcpy(mod2,"");
      scr=1;
    }
    else star=1;    
  }

  fprintf(stderr,"Read %s  ... ",argv[1]);
  npixels=naxes[0]*naxes[1];
  status=0;
  if(fits_open_file(&fp, argv[1], READONLY, &status))
    printerror(status);
  if(fits_read_keyword(fp, "B_SPECID", buf, NULL, &status))
    printerror(status);
  //fprintf(stderr,"B_SPECID=%s\n",buf);
  if(strncmp(buf+1,"SPEC1",5)==0) irs=0;
  else if(strncmp(buf+1,"SPEC2",5)==0) irs=1;
  else{
    fprintf(stderr,"B_SPECID is wrong (%s).\n",buf);
    exit(0);
  }

  if(irs){
    if(fits_read_key(fp, TSTRING, "DET-NSMP", buf, NULL, &status))
      printerror(status);
    sscanf(buf,"%d",&ihdumax);
  }else ihdumax=1;

  fprintf(stderr,"Read %s  ... \n",argv[2]);
  npixels=naxes[0]*naxes[1];
  status=0;
  if(fits_open_file(&fps, argv[2], READONLY, &status))
    printerror(status);
  if(irs){
    if(fits_read_key(fp, TSTRING, "DET-NSMP", buf, NULL, &status))
      printerror(status);
    sscanf(buf,"%d",&ihdumaxs);
  }else ihdumaxs=1;

  //fprintf(stderr,"ihdumax=%d,%d\n",ihdumax,ihdumaxs);
  if(ihdumaxs!=ihdumax){
    fprintf(stderr,"NUM_EXPS is not the same (%d,%d).\n",ihdumax,ihdumaxs);
    exit(0);
  }
  if(ihdumax==2) ihdumax=1;
  for(i=0;i<LINEPIX;i++){
    for(j=0;j<LINEPIX;j++){
      zs[i][j]=0;
      xzs[i][j]=0;
      zzs[i][j]=0.;
      cr[i][j]=0;
    }
  }
  sx=0;sx2=0;sigma=0;
  if(ihdumax>1) ihdus=2; else ihdus=1;
  fprintf(stderr,"ihdu=");
  for(ihdu=0;ihdu<ihdumax;ihdu++){
    fprintf(stderr,"%d ",ihdu+ihdus);
    if(fits_movabs_hdu(fp, ihdu+ihdus, &itype, &status) )
      printerror(status);
    if(fits_movabs_hdu(fps, ihdu+ihdus, &itype, &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,"Wrong format in images #1(%d %ld %ld).\n",ihdu,naxes[0],naxes[1]);
      exit(0);
    }
    if(fits_read_keys_lng(fps, "NAXIS", 1, 2, naxes, &nfound, &status)) 
      printerror(status);
    //fprintf(stderr,"%ld x %ld format\n",naxes[0],naxes[1]);
    if(naxes[0]!=LINEPIX || naxes[1]!=LINEPIX){
      fprintf(stderr,"Wrong format in images #2(%d %ld %ld).\n",ihdu,naxes[0],naxes[1]);
      exit(0);
    }
    fpixel = 1;
    nullval = 0;
    nbuffer = naxes[0];
    for(i=0;i<naxes[1];i++){
      if(ihdumax==1){
	if (fits_read_img(fp, TINT, fpixel, nbuffer, &nullval,
			  xzs[0], &anynull, &status))
	  printerror(status);
	if (fits_read_img(fps, TINT, fpixel, nbuffer, &nullval,
			  xzs[1], &anynull, &status))
	  printerror(status);
	for(j=0;j<naxes[0];j++) zs[i][j]=xzs[0][j]-xzs[1][j];
      }else{
	if (fits_read_img(fp, TUSHORT, fpixel, nbuffer, &nullval,
			  buffer, &anynull, &status))
	  printerror(status);
	if (fits_read_img(fps, TUSHORT, fpixel, nbuffer, &nullval,
			  buffers, &anynull, &status))
	  printerror(status);
	for(j=0;j<naxes[0];j++){
	  if(xzs[i][j]>MARK_SATULATE){
	    a=((double)xzs[i][j]*ihdu-(double)zs[i][j]*sx)/sigma;
	    b=(zs[i][j]-a*sx)/ihdu;
	    if(buffer[j]<SATULATE && buffers[j]<SATULATE){
	      zex=(int)(a*ihdu+b);
	      zre=(int)buffer[j]-(int)buffers[j];
	      stddev=sqrt((zzs[i][j]-2*a*xzs[i][j]-2*b*zs[i][j]+a*a*sx2+2*a*b*sx)/ihdu+b*b);
	      //if(j==1100 && i==1100) fprintf(stderr,"%d %d %d %d %d %g %d %d %d %g %g %d %d %d %d %d\n",ihdu,buffer[j],buffers[j],zs[i][j],xzs[i][j],zzs[i][j],sx,sx2,sigma,a,b,zre,zex,cr[i][j],zex-zre,(int)stddev);
	      if(ihdu>3 && cr[i][j]<=ZSHIFT/2){
		if(ihdu>4 && abs(zre-zex)>COSMIC_RAY && abs((zre-zex)-cr[i][j])>COSMIC_RAY && abs(zre-zex)>stddev*10 && abs((zre-zex)-cr[i][j])>stddev*10){
		  //fprintf(stderr,"%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\n",ihdu,j,i,(int)(a*ihdu),(int)b,zre,zex,cr[i][j],zre-zex,(int)stddev);
		  cr[i][j]=zre-zex+ZSHIFT;
		}
		else
		  cr[i][j]=zre-zex;
	      }
	      //if(i==0&&j==0)getchar();
	      if(cr[i][j]>ZSHIFT/2) zre=(int)buffer[j]-(int)buffers[j]-cr[i][j]+ZSHIFT;
	      else zre=(int)buffer[j]-(int)buffers[j];
	      zs[i][j]+=zre;
	      xzs[i][j]+=ihdu*zre;
	      zzs[i][j]+=zre*zre;
	    }
	    else{
	      if(sigma==0){
		zs[i][j]=BAD_PIXEL;
		xzs[i][j]=MARK_BAD;
	      }
	      else{
		zs[i][j]=(int)(a*(ihdumax-1));
		xzs[i][j]=MARK_SATULATE;
	      }
	    }
	  }
	}
      }
      fpixel+=naxes[0];
    }
    if(ihdumax!=1){
      sx+=ihdu;
      sx2+=ihdu*ihdu;
      sigma=(ihdu+1)*sx2-sx*sx;
      //fprintf(stderr,"%d %d %d %d %d %d %f\n",ihdu,zs[0][1],xzs[0][1],sx,sx2,sigma,(float)(((double)xzs[0][1]*(ihdu+1)-zs[0][1]*sx)*(ihdumax-1)/sigma));
    }
    //j=xxx;i=xxx;
    //fprintf(stderr,"%d\t%d\t%d\t%d\t%d\n",sx,sx2,zs[i][j],xzs[i][j],(int)(((double)xzs[i][j]*(ihdu+1)-(double)zs[i][j]*sx)*ihdu/sigma));
  }fprintf(stderr,"\n");
  if(ihdumax!=1){
    if(sigma==0) sigma=1;
    ncr=0;
    for(i=0;i<naxes[1];i++){
      for(j=0;j<naxes[0];j++){
	if(xzs[i][j]!=MARK_BAD && xzs[i][j]!=MARK_SATULATE)
	  zs[i][j]=(int)(((double)xzs[i][j]*ihdumax-(double)zs[i][j]*sx)*(ihdumax-1)/sigma);
	if(cr[i][j]>ZSHIFT/2){
	  cr[i][j]-=ZSHIFT;
	  ncr++;
	}else cr[i][j]=0;
      }
    }
    fprintf(stderr,"CR rejected : %d pixels\n",ncr);
  }
  if(fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK ");
  if(fits_close_file(fps, &status))
    printerror(status);
  fprintf(stderr,"OK\n");

  if(ihdumax!=1){
    remove("cr.fits");
    fprintf(stderr,"Write cr.fits ... ");
    status=0;
    if (fits_create_file(&fp, "cr.fits", &status))
      printerror(status);
    if (fits_create_img(fp, bitpix, naxis, naxes, &status))
      printerror(status);
    fpixel=1;
    fprintf(stderr,"%d x %d format ",LINEPIX,LINEPIX);
    nelements=LINEPIX;
    for(i=0;i<LINEPIX;i++){
      if(fits_write_img(fp, TINT, fpixel, nelements, cr[i], &status))
	printerror(status);
      fpixel+=LINEPIX;
    }
    if (fits_close_file(fp, &status))
      printerror(status);
    fprintf(stderr,"OK\n");
  }

  if(star){
    fprintf(stderr,"Read %s  ... ",argv[4]);
    status=0;
    if(fits_open_file(&fp,argv[4], READONLY, &status))
      printerror(status);
    if(fits_read_keys_lng(fp, "NAXIS", 1, 2, naxes, &nfound, &status))
      printerror(status);
    status=0;
    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, TINT, fpixel, nbuffer, &nullval,
			xzs[i], &anynull, &status))
	printerror(status);
      for(j=0;j<naxes[0];j++) zs[i][j]+=xzs[i][j];
      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;
  fprintf(stderr,"%d x %d format ",LINEPIX,LINEPIX);
  nelements=LINEPIX;
  for(i=0;i<LINEPIX;i++){
    if(fits_write_img(fp, TINT, fpixel, nelements, zs[i], &status))
      printerror(status);
    fpixel+=LINEPIX;
  }
  if (fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK\n");

  if(scr){
    fprintf(stderr,"Write subtract.cl ... ");
    if((fpt = fopen("subtract.cl","w")) == NULL){
      fprintf(stderr,"Cannot open subtract.cl\n");
      exit(1);
    }

    sscanf(argv[1],"%s",fname[0]);
    sscanf(argv[2],"%s",fname[1]);
    sscanf(argv[3],"%s",fname[2]);
    if((ret = strrchr(fname[0], '.')) != NULL) *(ret)='\0';
    if((ret = strrchr(fname[0], '/')) != NULL) strcpy(fname[0],ret+1);
    if((ret = strrchr(fname[1], '.')) != NULL) *(ret)='\0';
    if((ret = strrchr(fname[1], '/')) != NULL) strcpy(fname[1],ret+1);
    if((ret = strrchr(fname[2], '.')) != NULL) *(ret)='\0';
    
    if(irs){
      fprintf(fpt,"!sed -e 's/FMSAfileno/%s/g' ../redscr1%s_ | sed -e 's/BAND/%s/g' | sed -e 's/SHIFT/0/g' > tmpscr\n",fname[2],mod2,mod);
      fprintf(fpt,"!sed -e 's/FMSAfileno/%s/g' ../redscr2%s | sed -e 's/irs1/irs2/g' | sed -e 's/BAND/%s/g' | sed -e 's/SHIFT/0/g' >> tmpscr\n",fname[2],mod2,mod);
    }else{
      fprintf(fpt,"!sed -e 's/FMSAfileno/%s/g' ../redscr1%s | sed -e 's/BAND/%s/g' | sed -e 's/SHIFT/0/g' > tmpscr\n",fname[2],mod2,mod);
      fprintf(fpt,"!sed -e 's/FMSAfileno/%s/g' ../redscr2%s | sed -e 's/BAND/%s/g' | sed -e 's/SHIFT/0/g' >> tmpscr\n",fname[2],mod2,mod);
    }
    fprintf(fpt,"cl < tmpscr\n");
    fprintf(fpt,"# If you want to detect fainter stars, replace the 3rd argument of \"inverse\" command to a smaller absolute value (-4 => -3 etc...)\n");
    fprintf(fpt,"#   and restart from here.\n");
    fprintf(fpt,"imdel tmp0,check,%ssub1,%ssub\n",fname[0],fname[0]);
    fprintf(fpt,"imcopy %sxp tmp0.fits\n",fname[2]);
    if(strcmp(mod2,"s")==0) fprintf(fpt,"!../inverse tmp0.fits %ssub1.fits -4 irs%d_cont_%sa.dat irs%d_cont_%se1.dat irs%d_cont_%se2.dat\n",fname[0],irs+1,mod,irs+1,mod,irs+1,mod);
    else fprintf(fpt,"!../inverse tmp0.fits %ssub1.fits -4 irs%d_cont_%sa.dat irs%d_cont_%se.dat\n",fname[0],irs+1,mod,irs+1,mod);
    fprintf(fpt,"imcomb %ssub?.fits %ssub\n",fname[0],fname[0]);
    fprintf(fpt,"displ %s 1\n",fname[2]);
    if(irs) fprintf(fpt,"imarith %ssub * ../flat_ %ssub\n",fname[0],fname[0]);
    else fprintf(fpt,"imarith %ssub * ../flat %ssub\n",fname[0],fname[0]);
    fprintf(fpt,"displ %ssub 2\n",fname[0]);
    fprintf(fpt,"imarith tmp0 * -1 tmp0\n");
    fprintf(fpt,"imdel check,%ssub2,%ssub\n",fname[1],fname[1]);
    if(strcmp(mod2,"s")==0) fprintf(fpt,"!../inverse tmp0.fits %ssub2.fits -4 irs%d_cont_%sa.dat irs%d_cont_%se1.dat irs%d_cont_%se2.dat\n",fname[1],irs+1,mod,irs+1,mod,irs+1,mod);
    else fprintf(fpt,"!../inverse tmp0.fits %ssub2.fits -4 irs%d_cont_%sa.dat irs%d_cont_%se.dat\n",fname[1],irs+1,mod,irs+1,mod);
    fprintf(fpt,"imcomb %ssub?.fits %ssub\n",fname[1],fname[1]);
    if(irs) fprintf(fpt,"imarith %ssub * ../flat_ %ssub\n",fname[1],fname[1]);
    else fprintf(fpt,"imarith %ssub * ../flat %ssub\n",fname[1],fname[1]);
    fprintf(fpt,"displ %ssub 3\n",fname[1]);
  
    fclose(fpt);
    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 */
}