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

#define LINEPIX 2048
#define LINEPIY 2200
#define NOPIX  -1e10

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,naxes1;
  float nullval;
  static float buffer[LINEPIX],zo[LINEPIY][LINEPIX],z[LINEPIY][LINEPIX],zs[200][LINEPIX],yc[200],bg[81][LINEPIX],bpo[LINEPIY][LINEPIX],bp[LINEPIY][LINEPIX],bps[200][LINEPIX];

  float zav[LINEPIY],dyld,dylu,dyrd,dyru,dycd,dycu,dy,dyl,dyr,dyd,dyu,dya,dyb,dyc,dx,ddy,ddyl,ddyr,ddyc,su,sc,sd,acc,accmin,zt,zt1,zt2,ys=0.,bpt,bpt1,bpt2;
  float w[9]={0,1,1,1,1,1,1,1,0};
  //float w[9]={0.5,0.75,1,1,1,1,1,0.75,0.5};
  //float w[9]={0.06,0.21,0.5,0.84,1,0.84,0.5,0.21,0.06};
  double zamp,zampmax,zampmax2,yi,zmav,zmsg,zm[4][200];
  int di,k,kmin,kmax,mode=9,dyi,loop,dk,kk,nd,bpm;
  int inc[10][4]={{0,0,0,0},{-1,-1,0,0},{1,1,0,0},{-1,1,0,0},{1,-1,0,0},{0,0,-1,-1},{0,0,1,1},{0,0,-1,1},{0,0,1,-1},{0,0,0,0}};
  int inc2[4]={0,-1,1,0},bo[200];
  char buf[128],fname[64];
  FILE *fpd;

  if(argc<3 || argc>5){
    fprintf(stderr,"Usage : %s input output\n",argv[0]);
    fprintf(stderr,"  or    %s input output param_file [yshift]\n",argv[0]);
    exit(0);
  }
  
  if(argc>=4){
    if((fpd=fopen(argv[3],"r"))==NULL){
      fprintf(stderr,"Can not open file: %s\n",argv[3]);
      exit(1);
    }
    fgets(buf,sizeof(buf),fpd);
    if(sscanf(buf,"%f %f %f %f",&dyd,&dyu,&dya,&dyb)!=4){
      fprintf(stderr,"Wrong format\n");
      exit(1);
    }
    while(!feof(fpd)){
      if(fgets(buf,sizeof(buf),fpd)==NULL) break;
      if(sscanf(buf,"%d %f",&k,&sc)!=2){
	fprintf(stderr,"Wrong format\n");
	exit(1);
      }
      yc[k]=sc;
    }
    fclose(fpd);
    if(argc==5){
      sscanf(argv[4],"%f",&ys);
    }
  }else{
    dyd=0.;dyu=0.;dya=0.;dyb=0.;mode=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]>LINEPIY)){
    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++) zo[i][j]=buffer[j];
    fpixel+=naxes[0];
  }
  if(fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK\n");
  
  strcpy(fname,argv[1]);
  fname[strlen(argv[1])-5]='\0';
  sprintf(fname,"%s_bpm.fits",fname);
  status=0;
  fprintf(stderr,"Read %s  ... ",fname);
  if(fits_open_file(&fp, fname, READONLY, &status)){
    fprintf(stderr,"Failed! Bad Pixel Map is neglected.\n");
    for(i=0;i<naxes[1];i++){
      for(j=0;j<naxes[0];j++) bpo[i][j]=1;
    }
    bpm=0;
  }else{
    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)){
      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++) bpo[i][j]=buffer[j];
      fpixel+=naxes[0];
    }
    if(fits_close_file(fp, &status))
      printerror(status);
    fprintf(stderr,"OK\n");
    bpm=1;
  }

  acc=10.;accmin=0.001;
  while(1){
    dyd+=acc*inc[mode][0];
    dyu+=acc*inc[mode][1];
    dya+=acc*inc[mode][2];
    dyb+=acc*inc[mode][3];
    dyrd= dyd;dyru= dyu;
    dyld=-dyd;dylu=-dyu;
    dycd= dya;dycu= dyb;
    dy = dyld;
    if(dyrd < dy) dy = dyrd;
    dy=floor(dy);
    if(dycd<0)dy+=dycd;
    dyld-=dy;
    dylu-=dy;
    dyrd-=dy;
    dyru-=dy;
    for(i=0;i<LINEPIY;i++){
      for(j=0;j<LINEPIX;j++){z[i][j]=0.;bp[i][j]=0.;}
    }
 
    for(i=-100;i<naxes[1]+100;i++){
      dx=(float)i/(naxes[1]-1);
      dyl=(1-dx)*dyld+dx*dylu;
      dyr=(1-dx)*dyrd+dx*dyru;
      dyc=(1-dx)*dycd+dx*dycu;
      ddyl=(dylu-dyld)/(naxes[1]-1);
      ddyr=(dyru-dyrd)/(naxes[1]-1);
      ddyc=(dycu-dycd)/(naxes[1]-1);
      for(j=0;j<naxes[0];j++){
	dx=(float)j/(naxes[0]-1);
	dy=(1-dx)*dyl+dx*dyr+4*dyc*(dx-0.5)*(dx-0.5);
	ddy=((1-dx)*ddyl+dx*ddyr+4*ddyc*(dx-0.5)*(dx-0.5))*0.5;
	dyi=(int)dy;
	if(dy<0) dyi--;
	if(i<0 || i>=naxes[1]){
	  if(i+dyi>=0 && i+dyi<LINEPIY){z[i+dyi][j]=NOPIX;bp[i+dyi][j]=0.;}
	  if(i+dyi+1>=0 && i+dyi+1<LINEPIY){z[i+dyi+1][j]=NOPIX;bp[i+dyi+1][j]=0.;}
	}
	else{
	  if(i+dyi>=0 && i+dyi<LINEPIY){
	    if(z[i+dyi][j]>NOPIX*0.99 && zo[i][j]>NOPIX*0.99){
	      if(dy-dyi-ddy<0){
		z[i+dyi][j]+=zo[i][j]/(1+2*ddy);
		bp[i+dyi][j]+=bpo[i][j]/(1+2*ddy);
		if(i+dyi-1>=0){
		  if(z[i+dyi-1][j]>NOPIX*0.99){
		    z[i+dyi-1][j]+=zo[i][j]*(ddy-(dy-dyi))/(1+2*ddy);
		    bp[i+dyi-1][j]+=bpo[i][j]*(ddy-(dy-dyi))/(1+2*ddy);
		    //if((i+dyi-1==1097||i+dyi-1==1098||i+dyi-1==1099) && j==14) fprintf(stderr,"2: %d %d %d %g %g %g %g\n",i+dyi-1,i,dyi,dy,ddy,bpo[i][j],bp[i+dyi-1][j]);
		  }
		}
	      }else if(dy-dyi-ddy<1){
		z[i+dyi][j]+=zo[i][j]*(1-(dy-dyi-ddy))/(1+2*ddy);
		bp[i+dyi][j]+=bpo[i][j]*(1-(dy-dyi-ddy))/(1+2*ddy);
	      }
	      //if((i+dyi==1097||i+dyi==1098||i+dyi==1099) && j==14) fprintf(stderr,"1: %d %d %d %g %g %g %g\n",i+dyi,i,dyi,dy,ddy,bpo[i][j],bp[i+dyi][j]);
	    }
	    else{
	      z[i+dyi][j]=NOPIX;
	      bp[i+dyi][j]=0.;
	    }
	  }
	  if(i+dyi+1>=0 && i+dyi+1<LINEPIY){
	    if(z[i+dyi+1][j]>NOPIX*0.99 && zo[i][j]>NOPIX*0.99){
	      if(dy-dyi+ddy>1){
		z[i+dyi+1][j]+=zo[i][j]/(1+2*ddy);
		bp[i+dyi+1][j]+=bpo[i][j]/(1+2*ddy);
		if(i+dyi+2<LINEPIY){
		  if(z[i+dyi+2][j]>NOPIX*0.99){
		    z[i+dyi+2][j]+=zo[i][j]*(dy-dyi+ddy-1)/(1+2*ddy);
		    bp[i+dyi+2][j]+=bpo[i][j]*(dy-dyi+ddy-1)/(1+2*ddy);
		    //if((i+dyi+2==1097||i+dyi+2==1098||i+dyi+2==1099) && j==14) fprintf(stderr,"2: %d %d %d %g %g %g %g\n",i+dyi+2,i,dyi,dy,ddy,bpo[i][j],bp[i+dyi+2][j]);
		  }
		}
	      }else if(dy-dyi+ddy>0){
		z[i+dyi+1][j]+=zo[i][j]*(dy-dyi+ddy)/(1+2*ddy);
		bp[i+dyi+1][j]+=bpo[i][j]*(dy-dyi+ddy)/(1+2*ddy);
	      }
	      //if((i+dyi+1==1097||i+dyi+1==1098||i+dyi+1==1099) && j==14) fprintf(stderr,"2: %d %d %d %g %g %g %g\n",i+dyi+1,i,dyi,dy,ddy,bpo[i][j],bp[i+dyi+1][j]);
	    }else{
	      z[i+dyi+1][j]=NOPIX;
	      bp[i+dyi+1][j]=0.;
	    }
	  }
	}
      }
    }
    if(dylu>dyru) di=(int)dylu+1;
    else di=(int)dyru+1;
    if(dycu>0)di+=dycu;
    //fprintf(stderr,"%g %g %g %g %g %g %d\n",dyld,dyrd,dycd,dylu,dyru,dycu,di);getchar();
    
    for(i=0;i<naxes[1]+di;i++){
      zav[i]=0.;
      for(j=50;j<naxes[0]-50;j++){
	if(z[i][j]>NOPIX*0.99) zav[i]+=z[i][j]/(naxes[0]-100);
      }
    }
    zamp=0;
    for(i=5;i<naxes[1]+di;i++)
      zamp+=(zav[i]-zav[i-5])*(zav[i]-zav[i-5])/(naxes[1]+di-5);
    zamp=sqrt(zamp);
    
    fprintf(stderr,"%d %g %g %g %g Amplitude=%g\n",mode,dyd,dyu,dya,dyb,zamp);
    
    if(mode==9) break;
    else if(mode==0){
      mode=1;
      zampmax=zamp;
      zampmax2=zampmax;
    }
    else if(zampmax>zamp){
      dyd-=acc*inc[mode][0];
      dyu-=acc*inc[mode][1];
      dya-=acc*inc[mode][2];
      dyb-=acc*inc[mode][3];
      if(mode==8){
	if(zampmax-zampmax2>1e-5){mode=0;zampmax2=zampmax;}
	else if(acc>accmin){acc/=10.;mode=0;}
      }
      mode++;
    }else zampmax=zamp;
  }
  naxes[1]+=di;
  if(argc==3) printf("%g %g %g %g\n",dyd,dyu,dya,dyb);

  if(1){
    fprintf(stderr,"Write check.fits  ... ");
    status=0;
    if (fits_create_file(&fp, "check.fits", &status))
      printerror(status);
    if (fits_create_img(fp, bitpix, naxis, naxes, &status))
      printerror(status);
    fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[1]);
    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");
  }

  k=0;
  if(argc==3){
    for(i=4;i<naxes[1]-4;i++){
      if(zav[i]>zav[i-1] && zav[i]>zav[i+1] && k<200){
	sd=zav[i-4]+zav[i-3]+zav[i-2];
	sc=zav[i-1]+zav[i]+zav[i+1];
	su=zav[i+2]+zav[i+3]+zav[i+4];
	//fprintf(stderr,"%d %d (%g %g %g)\n",k,i,sd,sc,su);getchar();
	if(sc<sd || sc<su) continue;
	if((sc<sd+zamp/4 && sc<su+zamp/4)||(k==0 && sc<sd+zamp/2 && sc<su+zamp/2)) continue;
	dy=(su-sd)/(2*sc-sd-su)*3/2;
	if(dy<-1.5 || dy>1.5) {fprintf(stderr,"Detection error: %d %g %g %g %g %g %g %g %g %g\n",i,zav[i-4],zav[i-3],zav[i-2],zav[i-1],zav[i],zav[i+1],zav[i+2],zav[i+3],zav[i+4]);exit(1);}
	yc[k]=i+dy;
	//fprintf(stderr,"%d %d %g (%g %g %g)\n",k,i,yc[k],sd,sc,su);getchar();
	k++;i+=5;
      }
    }
    kmax=k;
    if(yc[1]-yc[0]>15.){
      for(j=0;j<kmax-1;j++) yc[j]=yc[j+1];
      kmax--;
    }
    
    for(k=0;k<kmax-1;k++){
      //fprintf(stderr,"%d %d %g %g\n",k+1,k,yc[k+1],yc[k]);getchar();
      if(yc[k+1]-yc[k]>14.){
	if(yc[k+1]-yc[k]<22.){
	  for(kk=kmax-1;kk>k;kk--) yc[kk+1]=yc[kk];
	  kmax++;
	  yc[k+1]=(yc[k]+yc[k+2])/2;
	  fprintf(stderr,"Position of the spectrum @Y=%g in check.fits is added.\n",yc[k+1]);
	}else{
	  if(k<20) kmin=19-k;
	  //fprintf(stderr,"%d %d %d\n",kmin,k,kmax);getchar();
	  dk=(k+kmin+1)%20;
	  if(dk!=0){
	    fprintf(stderr,"Number of the fiber in the group[%d] is %d!\n",(k+kmin+1)/20,dk);
	    exit(1);
	  }
	}
      }
    }
    //fprintf(stderr,"%d %d\n",kmin,kmax);
    if((120<kmax && kmax<180 && yc[0]<100 && (kmax+kmin+1)%20==0)||yc[0]>800) kmin+=(200-kmax)/20*20;
    for(k=kmax-1;k>=0;k--) yc[k+kmin]=yc[k];
    kmax+=kmin;
    if(kmin!=0){
      for(k=kmin-1;k>=0;k--) yc[k]=2*yc[k+20]-yc[k+40];
    }
    //fprintf(stderr,"%d %d\n",kmin,kmax);
    if(kmax!=200){
      for(k=kmax;k<200;k++) yc[k]=2*yc[k-20]-yc[k-40];
    }
  }

  if(argc==5) mode=0; else mode=3;
  acc=0.1;accmin=0.001;yi=yc[100];
  for(k=0;k<200;k++){
    yc[k]+=ys;
    bo[k]=0;
  }

  if((fpd=fopen("extract.dat","r"))==NULL){
    fprintf(stderr,"Can not open file: extract.dat\n");
    exit(1);
  }

  loop=0;
  naxes1=naxes[1];
  naxes[1]=200;
  while(!feof(fpd)){
    if(fgets(buf,sizeof(buf),fpd)==NULL) break;
    if(sscanf(buf,"%f %f %f %f %f %f %f %f %f",&w[0],&w[1],&w[2],&w[3],&w[4],&w[5],&w[6],&w[7],&w[8])!=9){
      fprintf(stderr,"Wrong format in extract.dat (%s)\n",buf);
      continue;
    }
    
    while(1){
      for(k=0;k<200;k++) yc[k]+=acc*inc2[mode];
      zamp=0;
      for(k=0;k<200;k++){
	zm[0][k]=0;zm[1][k]=0;zm[2][k]=0;
	for(j=0;j<naxes[0];j++){zs[k][j]=0.;bps[k][j]=0.;}
	if(argc==3) printf("%d %f\n",k,yc[k]);
	ddy=yc[k]-floor(yc[k]);
	i=(int)(yc[k]-ddy+0.5);
	for(j=0;j<naxes[0];j++){
	  for(di=-4;di<=4;di++){
	    if(i+di>=0 && i+di+1<naxes1){zt1=z[i+di][j];zt2=z[i+di+1][j];bpt1=bp[i+di][j];bpt2=bp[i+di+1][j];}
	    else{zt1=NOPIX;zt2=NOPIX;bpt1=0.;bpt2=0.;}
	    if(zt1>NOPIX*0.99 && zt2>NOPIX*0.99){zt=zt1*(1-ddy)+zt2*ddy;bpt=bpt1*(1-ddy)+bpt2*ddy;}
	    else{zt=NOPIX;bpt=0.;}
	    if(zt>NOPIX*0.99 && zs[k][j]>NOPIX*0.99){zs[k][j]+=zt*w[di+4];bps[k][j]+=bpt*w[di+4];}
	    else{zs[k][j]=NOPIX;bps[k][j]=0.;}
	    if(j>=200 && j<naxes[0]-200){
	      if(di<=-2 && zt>NOPIX*0.99) zm[1][k]+=zt/(naxes[0]-400)/3;
	      else if(di>=2 && zt>NOPIX*0.99) zm[2][k]+=zt/(naxes[0]-400)/3;
	      else if(di>=-1 && di<=1 && zt>NOPIX*0.99) zm[0][k]+=zt/(naxes[0]-400)/3;
	    }
	  }
	}
	zm[3][k]=zm[0][k]-(zm[1][k]+zm[2][k])/2;
      }
      for(k=10;k<190;k++){
	if(bo[k]){
	  zamp+=fabs(zm[3][k]-(zm[3][k-1]+zm[3][k+1])/2);
	  //fprintf(stderr,"%d %g %g %g %g %g\n",k+1,yc[k],zm[3][k-1],zm[3][k],zm[3][k+1],fabs(zm[3][k]-(zm[3][k-1]+zm[3][k+1])/2));
	}
	//fprintf(stderr,"%d %d %g %g\n",k,j,zs[k][j],zamp);getchar();
      }
      fprintf(stderr,"%d %g %g Amplitude=%g\n",mode,acc,yc[100]-yi,zamp);//getchar();
      if(mode==3) break;
      else if(mode==0){
	zmav=0;zmsg=0;
	for(k=10;k<190;k++){
	  zm[3][k]=(zm[0][k]+zm[1][k]+zm[2][k])-(zm[0][k-1]+zm[1][k-1]+zm[2][k-1]+zm[0][k+1]+zm[1][k+1]+zm[2][k+1])/2;
	  zmav+=zm[3][k];
	  zmsg+=zm[3][k]*zm[3][k];
	}
	zmav/=180;zmsg=sqrt(zmsg/180-zmav*zmav);
	zamp=0;
	for(k=10;k<190;k++){
	  if(fabs(zm[3][k]-zmav)>zmsg*2){
	    fprintf(stderr,"Bright object #%d is found.(%g %g %g %g)\n",k+1,zm[3][k-1]-zmav,zm[3][k]-zmav,zm[3][k+1]-zmav,zmsg*2);
	    bo[k]=1;
	    zamp+=fabs(zm[0][k]-(zm[1][k]+zm[2][k])/2);
	    mode=1;
	  }
	}
	if(mode==0){
	  fprintf(stderr,"No bright object is found!\n");
	  break;
	}
	zampmax=zamp;
	fprintf(stderr,"%d %g %g Amplitude=%g\n",mode,acc,yc[100]-yi,zamp);
      }
      else if(zampmax>zamp){
	for(k=0;k<200;k++) yc[k]-=acc*inc2[mode];
	if(mode==2){
	  if(acc>accmin){acc/=2.;mode=0;}
	}
	mode++;
      }else zampmax=zamp;
    }
    
    if(loop==0) strcpy(fname,argv[2]);
    else{
      fname[strlen(argv[2])-5]='\0';
      sprintf(fname,"%s_%d.fits",fname,loop);
    }

    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);
    fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[1]);
    fpixel=1;
    nelements=naxes[0];
    for(i=0;i<naxes[1];i++){
      if (fits_write_img(fp, TFLOAT, fpixel, nelements, zs[i], &status))
	printerror(status);
      fpixel+=naxes[0];
    }
    if (fits_close_file(fp, &status))
      printerror(status);
    fprintf(stderr,"OK\n");

    if(bpm){
      fname[strlen(fname)-5]='\0';
      sprintf(fname,"%s_bpm.fits",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);
      fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[1]);
      fpixel=1;
      nelements=naxes[0];
      for(i=0;i<naxes[1];i++){
	if (fits_write_img(fp, TFLOAT, fpixel, nelements, bps[i], &status))
	  printerror(status);
	fpixel+=naxes[0];
      }
      if (fits_close_file(fp, &status))
	printerror(status);
      fprintf(stderr,"OK\n");
    }
    loop++;
  }
  fclose(fpd);

  if(argc>=4){
    sprintf(fname,"extract_%s",argv[3]);
    if((fpd=fopen(fname,"w"))==NULL){
      fprintf(stderr,"Can not open file: %s\n",fname);
      exit(1);
    }
    fprintf(fpd,"%g\n",yc[100]-yi);
    fclose(fpd);
  }

  if(1){
    for(k=0;k<9;k++){
      i=(int)((yc[k*20+19]+yc[k*20+20])/2-1);
      if(i<100) i=(int)((yc[k*20+39]+yc[k*20+40])/2-1);
      if(i>naxes1-100) i=(int)((yc[k*20-1]+yc[k*20])/2-1);
      fprintf(stderr,"%d %d\n",k,i);
      for(j=0;j<naxes[0];j++){
	for(di=-4;di<=4;di++){
	  bg[k*9+di+4][j]=z[i+di][j];
	}
      }
    }
    for(j=0;j<naxes[0];j++){
      zmav=0.;nd=0;
      for(i=0;i<81;i++){
	if(bg[i][j]>NOPIX*0.99){zmav+=bg[i][j];nd++;}
      }
      zmav/=nd;
      for(i=0;i<81;i++){
        if(bg[i][j]<NOPIX*0.99){bg[i][j]=zmav;}
      }
    }
    naxes[1]=9*9;
    fprintf(stderr,"Write bkg.fits  ... ");
    status=0;
    if (fits_create_file(&fp, "bkg.fits", &status))
      printerror(status);
    if (fits_create_img(fp, bitpix, naxis, naxes, &status))
      printerror(status);
    fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[1]);
    fpixel=1;
    nelements=naxes[0];
    for(i=0;i<naxes[1];i++){
      if (fits_write_img(fp, TFLOAT, fpixel, nelements, bg[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 */
}