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

#define NDAT 500
#define RAD   50

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,naxes[2]={NDAT, NDAT};
  float nullval;
  char array[NDAT],*buffer=array,fname[128],buf[256],c[3][64],*ret;
  int s[NDAT],st,n,nmax,np,imax,ii,nn,irs,sp[200],p[2][100],old=0;
  double x,y,r;
  static float t[16][NDAT],z[NDAT][16];
  FILE *fps;

  if(argc!=2 && argc!=3 && argc!=4){
    printf("Usage: %s exposure_list [pair_tbl]\n",argv[0]);
    exit(1);
  }
  if(argc!=2 && (strcmp(argv[2],"1")==0 || strcmp(argv[2],"2"))==0){
    old=1;
    sscanf(argv[2],"%d",&irs);
    irs--;
  }else{
    if(strncmp(argv[1],"irs1",4)==0) irs=0;
    else if(strncmp(argv[1],"irs2",4)==0) irs=1;
    else{
      fprintf(stderr,"The name of exposure_list should be start from 'irs1' or 'irs2' (%s).\n",argv[1]);
      exit(1);
    }
  }
  np=0;
  if(argc==3+old){
    if((fps = fopen(argv[2+old],"r")) == NULL){
      printf("Cannot open %s\n",argv[2+old]);
      exit(1);
    }
    while(!feof(fps)){
      if(fgets(buf,sizeof(buf),fps)==NULL) break;
      if(strncmp(buf,"#",1)==0) continue;
      sscanf(buf,"%d %d",&p[0][np],&p[1][np]);
      np++;
    }
    fclose(fps);
  }

  sprintf(fname,"../irs%d_slit_spine.tbl",irs+1);
  if((fps = fopen(fname,"r")) == NULL){
    sprintf(fname,"irs%d_slit_spine.tbl",irs+1);
    if((fps = fopen(fname,"r")) == NULL){
      printf("Cannot open %s\n",fname);
      exit(1);
    }
  }
  for(i=0;i<200;i++) fscanf(fps,"%d %d",&sp[i],&ii);
  fclose(fps);

  sscanf(argv[1],"%s",fname);
  if((fps = fopen(fname,"r")) == NULL){
    printf("Cannot open %s\n",fname);
    exit(1);
  }
  n=0;imax=0;
  while(!feof(fps)){
    if(fgets(buf,sizeof(buf),fps)==NULL) break;
    if(strncmp(buf,"#",1)==0) continue;
    if(old) sscanf(buf,"%s",fname);
    else {
      if(sscanf(buf,"%s %s %s",c[0],c[1],c[2])!=3){
	fprintf(stderr,"Wrong format: %s",buf);
	exit(1);
      }
      if(strcmp(c[1],"NONE")==0) strcpy(fname,c[2]);
      else strcpy(fname,c[1]);
      strcat(fname,".fits");
    }
    fprintf(stderr,"Read %s  ... ",fname);
    status=0;
    if(fits_open_table(&fp, fname, READONLY, &status))
      printerror(status);
    if(fits_read_keys_lng(fp, "NAXIS", 1, 2, naxes, &nfound, &status)) 
      printerror(status);
    fprintf(stderr,"%ld x %ld format \n",naxes[0],naxes[1]);
    if(fits_read_col(fp, TSTRING, 16, 1, 1, 1LL, &nullval, &buffer, &anynull, &status)){
      fprintf(stderr,"### Header does not contain spineerror information (old format). ###\n");
      exit(0);
    }
    if(imax<naxes[1]) imax=naxes[1];
    nullval = 0;i=0;
    for(j=0;j<naxes[1];j++) {
      if (fits_read_col(fp, TSTRING,  9, j+1, 1, 1LL, &nullval, &buffer, &anynull, &status))
	printerror(status);
      sscanf(buffer,"%lf",&x);
      if (fits_read_col(fp, TSTRING, 10, j+1, 1, 1LL, &nullval, &buffer, &anynull, &status))
	printerror(status);
      sscanf(buffer,"%lf",&y);
      r=sqrt(x*x+y*y)/RAD;
      if (fits_read_col(fp, TSTRING, 12, j+1, 1, 1LL, &nullval, &buffer, &anynull, &status))
	printerror(status);
      sscanf(buffer,"%d",&st);
      if(n>0){
	while(st>s[i]){
	  fprintf(stderr,"Spine No.%d is not allocated.\n",st);
	  t[n][i]=0.;
	  i++;
	}
	if(st<s[i]){
	  fprintf(stderr,"Spine No.%d is newly allocated.\n",st);
	  for(ii=imax-1;ii>i;ii--) {
	    s[ii]=s[ii-1];
	    for(nn=0;nn<n;nn++) t[nn][ii]=t[nn][ii-1];
	  }
	  for(nn=0;nn<n;nn++) t[nn][i]=0.;
	}
      }
      s[i]=st;
      if(r>1.) t[n][i]=0.;
      else{
	t[n][i]=(2*acos(r)-r*sqrt(1-r*r))/M_PI;
      }
      i++;
    }
    if(fits_close_file(fp, &status))
      printerror(status);
    
    for(i=0;i<200;i++){
      z[i][n]=0.;
      for(ii=0;ii<imax;ii++){
	if(sp[i]==s[ii]){
	  z[i][n]=t[n][ii];
	  break;
	}
      }
    }
    for(nn=0;nn<np;nn++){
      z[p[1][nn]-1][n]+=z[p[0][nn]-1][n];
      z[p[1][nn]-1][n]/=2;
      if(p[0][nn]!=p[1][nn]) z[p[0][nn]-1][n]=0.;
    }
    n++;
  }
  fclose(fps);
  
  nmax=n;
  naxes[0]=nmax;naxes[1]=200;
  if(old) sprintf(fname,"irs%d_spineerror_%s.fits",irs+1,argv[1]);
  else{
      sscanf(argv[1],"%s",buf);
      ret = strrchr(buf, '.');
      *(ret)='\0';
      sprintf(fname,"%s_spineerror.fits",buf);
  }
  fprintf(stderr,"\nWrite %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;
  fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[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 */
}