#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "fitsio.h"
#include <pthread.h>
#include <sched.h>
#include <unistd.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

typedef struct _thread_arg{
  int id, i, ihdu, ihdumax;
  int *zs[LINEPIX], *xzs[LINEPIX], *cr[LINEPIX];
  double *zzs[LINEPIX];
  unsigned short *buffer, *buffers;
  int sx,sx2,sigma;
  long naxes[2];
  int cpuid, cpu_num, lock;
} thread_arg_t;

pthread_mutex_t mutex;
pthread_cond_t cond;

// multi core thread ------------------------------------
void thread_func(thread_arg_t *arg) {
  thread_arg_t* targ = arg;
  int j;
  int c_start, c_end, range;
  long naxes[2] = {targ->naxes[0], targ->naxes[1]};
  int zex, zre;
  double a, b, stddev;

#ifdef __CPU_ZERO
  cpu_set_t mask;
  __CPU_ZERO(&mask);
  __CPU_SET(targ->cpuid, &mask);
  if(sched_setaffinity(0, sizeof(mask), &mask) == -1)
    fprintf(stderr, "WARNING: failed to set CPU affinity...(cpuid=%d)\n");
#endif

  // determine ranges of calculations for this thread
  range = naxes[0] / targ->cpu_num;
  c_start = targ->id * range;
  c_end = (targ->id+1) * range;
  if(c_end > naxes[0]) c_end = naxes[0];

  pthread_mutex_lock(&mutex);
  while(targ->lock < 0) pthread_cond_wait(&cond, &mutex);
  pthread_mutex_unlock(&mutex);

  while(targ->lock < 2){
    
    for(j=c_start;j<c_end;j++){
      if(targ->xzs[targ->i][j]>MARK_SATULATE){
	a=((double)targ->xzs[targ->i][j]*targ->ihdu-(double)targ->zs[targ->i][j]*targ->sx)/targ->sigma;
	b=(targ->zs[targ->i][j]-a*targ->sx)/targ->ihdu;
	if(targ->ihdu==0){
	}
	if(targ->buffer[j]<SATULATE && targ->buffers[j]<SATULATE){
	  zex=(int)(a*targ->ihdu+b);
	  zre=(int)targ->buffer[j]-(int)targ->buffers[j];
	  stddev=sqrt((targ->zzs[targ->i][j]-2*a*targ->xzs[targ->i][j]-2*b*targ->zs[targ->i][j]+a*a*targ->sx2+2*a*b*targ->sx)/targ->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(targ->ihdu>3 && targ->cr[targ->i][j]<=ZSHIFT/2){
	    if(targ->ihdu>4 && abs(zre-zex)>COSMIC_RAY && abs((zre-zex)-targ->cr[targ->i][j])>COSMIC_RAY && abs(zre-zex)>stddev*10 && abs((zre-zex)-targ->cr[targ->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);
	      targ->cr[targ->i][j]=zre-zex+ZSHIFT;
	    }
	    else
	      targ->cr[targ->i][j]=zre-zex;
	  }
	  //if(i==0&&j==0)getchar();
	  if(targ->cr[targ->i][j]>ZSHIFT/2) zre=(int)targ->buffer[j]-(int)targ->buffers[j]-targ->cr[targ->i][j]+ZSHIFT;
	  else zre=(int)targ->buffer[j]-(int)targ->buffers[j];
	  targ->zs[targ->i][j]+=zre;
	  targ->xzs[targ->i][j]+=targ->ihdu*zre;
	  targ->zzs[targ->i][j]+=zre*zre;
	}
	else{
	  if(targ->sigma==0){
	    targ->zs[targ->i][j]=BAD_PIXEL;
	    targ->xzs[targ->i][j]=MARK_BAD;
	  }
	  else{
	    targ->zs[targ->i][j]=(int)(a*(targ->ihdumax-1));
	    targ->xzs[targ->i][j]=MARK_SATULATE;
	  }
	}
      }
    }

    pthread_mutex_lock(&mutex);
    targ->lock = -1;
    pthread_cond_broadcast(&cond);
    while(targ->lock < 0) pthread_cond_wait(&cond, &mutex);
    pthread_mutex_unlock(&mutex);
    
  }

  return;
}

// end of the multi core thread -------------------------

void printerror(int status);

int main(int argc,char *argv[])
{
  fitsfile *fp,*fps;
  int i,j,k,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,irs,star=0,scr=0,ncr;
  char buf[256],fname[3][128],*ret,mod[8],mod2[4];
  //int zex,zre;
  //double a,b,stddev;
  FILE *fpt;

  int cpu_num = sysconf(_SC_NPROCESSORS_CONF);
  int core;

  pthread_t handle[cpu_num];
  thread_arg_t targ[cpu_num];
  int masterkey;

  fprintf(stderr,"CPU_num=%i\n",cpu_num);

  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;

  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;

  // Creating thread ------------------------------
  if(ihdumax != 1){
    for(core = 0; core < cpu_num; core++){

      targ[core].id = core;
      targ[core].cpuid = core%cpu_num;
      targ[core].cpu_num = cpu_num;
      
      targ[core].naxes[0] = naxes[0];
      targ[core].naxes[1] = naxes[1];
      
      targ[core].ihdumax = ihdumax;

      targ[core].buffer = buffer;
      targ[core].buffers = buffers;
      
      for(k=0; k<LINEPIX; k++){
	targ[core].zs[k] = zs[k];
	targ[core].xzs[k] = xzs[k];
	targ[core].cr[k] = cr[k];
	targ[core].zzs[k] = zzs[k];
      }
      
      targ[core].lock = -1;
      
      pthread_create(&handle[core], NULL, (void *)thread_func, (void *)&targ[core]);
    }
  }
  // ----------------------------------------------

  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);
	
	// program for multicore -----------------
	for(core = 0; core < cpu_num; core++){
	  targ[core].ihdu = ihdu;
	  targ[core].i = i;
	  targ[core].sx = sx;
	  targ[core].sx2 = sx2;
	  targ[core].sigma = sigma;
	  
	  targ[core].lock = 1;
	}
	// ------------------------------------
	
	/*
	  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(ihdu==0){
	    }
	    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;
	      }
	    }
	  }
	}
	*/
	pthread_mutex_lock(&mutex);
	pthread_cond_broadcast(&cond);
	masterkey = 0;
	while(masterkey > -1*cpu_num) {
	  pthread_cond_wait(&cond, &mutex);
	  masterkey = 0;
	  for(core = 0; core < cpu_num; core++){
	    masterkey += targ[core].lock;
	  }
	}
	pthread_mutex_unlock(&mutex);
      }
      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){
    for(core = 0; core < cpu_num; core++){
	  targ[core].lock = 2;
    }
    pthread_cond_broadcast(&cond);

    //fprintf(stderr,"waiting for thread end");
    for(i = 0; i < cpu_num; i++)
      pthread_join(handle[i], NULL);
    //fprintf(stderr," --- all threads joined!\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 */
}