/*
  split one spectra file, containing 4 or 8 spectra, into
  single spectra, possibly reverse spectra, and sum them together, to
  for two spectra files, one for each phase.
  Input is read from the file specified.
  The peakpositions, phase and parity are read from *.def or global.def
  Output goes to USP1.spc/.err and USP2.spc/.err , if not specified
  with -o1 -o2.
*/

#include <stdio.h>
#include <spec.h>
#define TMPFILE "global.def"

float *spc, *err, *tim, *usp1, *usp2, *uer1, *uer2;
int range=1024,
    tri=0,
    dev=50;

help()
{
  printf("Splitting spectrum file into sinngle spectra, possibly reversing\n");
  printf("them, and adding them to 2 spectra (one for each phase)\n");
  printf("pacsplit file [options]\n");
  printf("options:\n");
  printf("  -o1 file    specifies output file (0). Default: USP1\n");
  printf("  -o2 file    specifies output file (90). Default: USP2\n");
  printf("  -def file   specifies another definition file than %s\n",TMPFILE);
  printf("The definition of the spectrum is read from %s\n",TMPFILE);
  printf("This file should be in a format as generated by the FPROMPT \n");
  printf("routine and is organized as follows:\n");
  printf("Each line consists of 4 integers.\n");
  printf("The first specifies the Time 0 position.\n");
  printf("The second can be +1 -1 +2 or -2 and specifies, if it is\n");
  printf("     to be mirrored (-) and if it is a 0 degree (1)\n");
  printf("     or 90 degree (2) spectrum\n");
  printf("     If any other number is given, the spectrum is ignored\n");
  printf("The last two numbers are used to specify the range\n");
  printf("of the usable data area. Up to now, these parameters\n");
  printf("are only significant in the last line !\n");
  printf("\n(C) Rainer Kowallik\n");
  exit(0);
}

main(argc,argv)
int argc;
char *argv[];
{
int  xbeg,xend,n,m,i,max,peakptr,pacmax,peaks[30],phase[30];
char z[80],comment[80],nam1[80],nam2[80],defnam[80];
char *home;
float x;
FILE *fp;

   strcpy(nam1,"USP1");
   strcpy(nam2,"USP2");
   strcpy(defnam,TMPFILE);

   if(checkopt(argc,argv,"-o1",z)) strcpy(nam1,z);
   if(checkopt(argc,argv,"-o2",z)) strcpy(nam2,z);
   if(checkopt(argc,argv,"-def",z)) strcpy(defnam,z);

   spc=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
   err=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
   tim=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
   usp1=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
   usp2=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
   uer1=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
   uer2=(float *)calloc(_MAXSPCLEN+2,sizeof(float));

/* -----------------------------------------------
   reading spectrum 
   ----------------------------------------------- */

   max=readspec(argv[1],spc,err,tim,comment);

/* -----------------------------------------------
   reading spectra definition file
   ----------------------------------------------- */

   fp=fopen(defnam,"r");
#ifdef UNIX
   if(fp==NULL) {
      home = (char *) getenv("HOME");
      if(home != NULL) {strcpy(defnam,home); strcat(defnam,"/"); }
      strcat(defnam,TMPFILE);
      fp=fopen(defnam,"r");
   }
   if(fp==NULL) {
      home = (char *) getenv("LISEPRG");
      if(home != NULL) {strcpy(defnam,home); strcat(defnam,"/"); }
      strcat(defnam,TMPFILE);
      fp=fopen(defnam,"r");
   }
#endif
#ifdef AMIGA
   if(fp==NULL) {
      strcpy(defnam,"LISE:"); strcat(defnam,TMPFILE);
      fp=fopen(defnam,"r");
   }
#endif
   if(fp == NULL) {
      printf("rvt2: could not open >%s<\n",TMPFILE);
      exit(0);
   }
   peakptr=0;
   while(!feof(fp)) {
      fscanf(fp,"%d %d %d %d\n",&peaks[peakptr],&phase[peakptr],&i,&pacmax);
      peakptr=peakptr+1;
   }
   fclose(fp);

/* -------------------------------------------------
   setting usp1, usp2 ,uer1 and uer2 to 0
   ------------------------------------------------- */

   for(n=0;n<pacmax;n++) {
      usp1[n]=0.0;
      usp2[n]=0.0;
      uer1[n]=0.0;
      uer2[n]=0.0;
   }

/* ----------------------------------------------
   here we start calculating the new sum spectra 
  -----------------------------------------------*/

   for(n=0;n<peakptr;n++) {
      xbeg = peaks[n];
      m = phase[n];
      switch(m) {     /* which spectrum to add ? reversing ? */
      case -2:        /* mirror, 90 degree spectrum */
	 for(i=0;i<pacmax;i++) {
	    usp2[i] = usp2[i] + spc[xbeg-i];
	    uer2[i] = uer2[i] + err[xbeg-i];
         }
	 break;
      case -1:       /* mirror 0 degree spectrum */
	 for(i=0;i<pacmax;i++) {
	    usp1[i] = usp1[i] + spc[xbeg-i];
	    uer1[i] = uer1[i] + err[xbeg-i];
         }
	 break;
      case 1:       /* normal 0 degree spectrum */
	 for(i=0;i<pacmax;i++) {
	    usp1[i] = usp1[i] + spc[xbeg+i];
	    uer1[i] = uer1[i] + err[xbeg+i];
         }
	 break;
      case 2:       /* normal 90 degree spectrum */
	 for(i=0;i<pacmax;i++) {
	    usp2[i] = usp2[i] + spc[xbeg+i];
	    uer2[i] = uer2[i] + err[xbeg+i];
         }
	 break;
      }
   }
   for(i=0;i<pacmax;i++) { /* generate new error spectrum */
      uer1[i] = (float) sqrt((double) usp1[i]);
      uer2[i] = (float) sqrt((double) usp2[i]);
   }

/* -------------------------------------------
   now we can go and save the new spectra 
   ------------------------------------------- */

   strcpy(z,comment);
   strcat(z," | 0");
   writespec(nam1,usp1,uer1,pacmax,1,z);
   strcpy(z,comment);
   strcat(z," | 90");
   writespec(nam2,usp2,uer2,pacmax,1,z);
   free(uer1); free(uer2); free(usp1); free(usp2);
   free(spc); free(err); free(tim);
   exit(0);
}