/*
  split one spectra file, containing 4 or 8 spectra, into
  single spectra, possibly reverse spectra, subtract Background
  Calculate PAC ratio function:

                                        1/4
          2   (product(Si(180) - BACKi))
   R(t) = - * --------------------------
          3   (product(Sj(90) - BACKj))^(1/4)

  Input is read from the file specified.
  The peakpositions, phase and parity are read from *.def or global.def
*/

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

char deffile[80];
float *spc, *err, *tim, *usp1, *usp2, *uer1, *uer2;
int range=1024,
    tri=0,
    flg_v=FALSE,
    flg_ba=0,
    lastn=10,
    dev=50;

int *back;

help()
{
  printf("Splitting spectrum file into single spectra, possibly reversing\n");
  printf("them, subtracting Background\n");
  printf("Calculating PAC Ratio function:\n");
  printf("                                        1/4       \n");
  printf("        2   (product(Si(180,t) - BACKi))          \n");
  printf(" R(t) = - * --------------------------------- - 1 \n");
  printf("        3   (product(Sj(90,t) - BACKj))^(1/4)     \n");
  printf("\n");
  printf("               1     2   [ 1                               \n"); 
  printf(" Error(t) = ------ * - * [-- * sum(1/(Si(180,t)-BACKi)) +  \n");
  printf("            R(t)+1   3   [16                               \n");
  printf("                                                          1/2 \n");
  printf("                            1                            ]    \n"); 
  printf("                         + -- * sum(1/(Sj(180,t)-BACKj)) ]    \n");
  printf("                           16                            ]    \n");
  printf("\n");
  printf("RVT2 file [-o file] [-back n] [-last n]\n");
  printf("options:\n");
  printf("   -o  file     output spectrum instead of stdout\n");
  printf("   -back n1,n2,n3,..  background for all spectra\n");
  printf("   -last n      takes arithmetic mean of last n channels as\n");
  printf("                background. Default is, to find the background\n");
  printf("                by fitting an exponential curve to the spectrum\n");
  printf("   -v           verbose mode. Prints fitted background\n");
  printf("   -def file    use spectra definition file different from %s\n",TMPFILE);
  printf("The definition of the spectrum is read from %s\n",deffile);
  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);
}

float sd(x,y)
float x,y;
{
if(y!=0.0) return((float)(x / y));
return(0.0);
}

float root(n,x)
float x;
int n;
{
float y;

y = fabs(exp((double)(log(fabs((double)x))/((double)n))));
return(y);
}

main(argc,argv)
int argc;
char *argv[];
{
int  xbeg,n,m,i,max,peakptr,pacmax,peaks[30],phase[30];
int  wurzel1,wurzel2,merk;
char *z,*s,*comment,*defnam;
FILE *fp;

   if(argc == 1) help();

   z = (char *) malloc(256);
   s = (char *) malloc(256);
   comment = (char *) malloc(256);
   defnam = (char *) malloc(256);
   back = (int *) calloc(32,sizeof(int));

   strcpy(deffile,TMPFILE);
   flg_ba=0;            /* defaults to exponential fit */
   for(n = 0; n < 32; n++) back[n] = -1;
   flg_v=FALSE;

   if(checkopt(argc,argv,"-def",z)) strcpy(deffile,z);
   strcpy(defnam,deffile);
   if(checkopt(argc,argv,"-back",z)) {
      i = 0; n = 0; m = 0; flg_ba=1;
      while(z[n] != 0) {
         s[m++] = z[n++];
         if(z[n] == ',') {
            s[m] = 0; m = 0; n = n + 1;
            back[i++] = atoi(s);
         }
      }
      while(i < 32) back[i++] = atoi(s);
   }
   if(checkopt(argc,argv,"-last",z)) {
      flg_ba=2; lastn = atoi(z);
   }
   if(checkopt(argc,argv,"-v",z)) flg_v=TRUE;


   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);

/* -----------------------------------------------
   generating error spectrum 
   ----------------------------------------------- */

   for(n=0;n<max;n++) err[n] = sqrt(spc[n]);

/* -----------------------------------------------
   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,deffile);
      fp=fopen(defnam,"r");
   }
   if(fp==NULL) {
      home = (char *) getenv("LISEPRG");
      if(home != NULL) {strcpy(defnam,home); strcat(defnam,"/"); }
      strcat(defnam,deffile);
      fp=fopen(defnam,"r");
   }
#endif
#ifdef AMIGA
   if(fp==NULL) {
      strcpy(defnam,"LISE:"); strcat(defnam,deffile);
      fp=fopen(defnam,"r");
   }
#endif
   if(fp == NULL) {
      printf("rvt2: could not open >%s<\n",deffile);
      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);

/* -------------------------------------------------
   initializing usp1, usp2 ,uer1 and uer2
   ------------------------------------------------- */

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

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

   wurzel1 = 0; wurzel2 = 0;

   for(n=0;n<peakptr;n++) {
      xbeg = peaks[n];
      m = phase[n];
/* if negative, go and mirror spectrum */
      if(m < 0) {
         for(i=pacmax/2 ; i<pacmax ; i++) {
            merk = spc[xbeg - i];
            spc[xbeg - i] = spc[xbeg + i - pacmax];
            spc[xbeg + i - pacmax] = merk;
            err[xbeg - i] = sqrt(spc[xbeg - i]);
            err[xbeg + i - pacmax] = sqrt(spc[xbeg + i - pacmax]);
         }
         xbeg = xbeg - pacmax; m = 0 - m;
      }

      SubBackground(spc,err,xbeg,xbeg+pacmax,n);

      switch(m) {   /* which spectrum to add ? */
      case 1:       /* 180 degree spectrum */
         for(i=0;i<pacmax;i++) {
         usp1[i] = usp1[i] * spc[xbeg+i];
         uer1[i] = uer1[i] + sd(1.0,spc[xbeg+i]);
         }
         wurzel1 = wurzel1 + 1;
         break;
      case 2:       /* 90 degree spectrum */
         for(i=0;i<pacmax;i++) {
         usp2[i] = usp2[i] * spc[xbeg+i];
         uer2[i] = uer2[i] + sd(1.0,spc[xbeg+i]);
         }
         wurzel2 = wurzel2 + 1;
         break;
      }
   }

/* --------------------------------------------------------------
         calculate ratio function and error spectrum 
   -------------------------------------------------------------- */

   for(i=0;i<pacmax;i++) {
      spc[i] = sd(2.0 * root(wurzel1,usp1[i]) , 3.0 * root(wurzel2,usp2[i]));
      spc[i] = spc[i] - 0.666667;
      err[i] = sd(1.0,1.0 + spc[i]) * (2.0 / 3.0);
      err[i] = err[i] *
         sqrt(sd(1.0,(float)(wurzel1 * wurzel1)) * uer1[i] + 
              sd(1.0,(float)(wurzel2 * wurzel2)) * uer2[i]);
   }


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

   strcpy(z,comment);
   n=instr("|",comment);
   if(n>0) midstr(z,comment,0,n-1);
   strcat(z," | RVT(PAC)");
   writespec("",spc,err,pacmax,2,z);
   free(uer1); free(uer2); free(usp1); free(usp2);
   free(spc); free(err); free(tim);
   free(back); free(defnam); free(comment);
   free(s); free(z);
   exit(0);
}
/* --------------------------------------------
   finding and subtracting correct background
   -------------------------------------------- */
SubBackground(spc,err,xbeg,xend,j)
float *spc, *err;
int xbeg, xend, j;
{
int n,ba;

   switch(flg_ba) {
   case 0:
      ba = expbafit(spc,err,xbeg,xend);
      break;
   case 1: /* allready done when parsing arguments */
      ba = back[j];
      break;
   case 2:
      ba = sumspc(spc,xend-lastn,xend) / lastn;
      break;
   }
   if(flg_v) printf("background=%d\n",ba);

/*    Subtracting background  */
   for(n=xbeg ; n<xend ; n++) spc[n] = spc[n] - ba;
   return(0);
}

sumspc(spc,a,z)
float spc[];
int a,z;
{
int i,erg;

   erg=0;
   for(i=a;i<z;i++) erg=erg+spc[i];
   return(erg);
}

/* ---------------------------------------------------------
   Linear regression after logarithm of spectrum
   This is only to calculate the chi**2, all
   other values are lost.
   UP TO NOW, THE ERRORBARS ARE NOT INCLUDED !
   I HAVE FORGOTTEN MY OLD SUPERBASIC SOURCE AT HOME, SORRY.
   --------------------------------------------------------- */

float linreg(spc,err,n,a,z)
float spc[],err[];
int n,a,z;
{
int i;
float ysoll,x,y,sumx,sumxx,sumy,sumyy,sumxy,nges,
      steigung,offset,erg,
      *ylog;

ylog=(float *)calloc(_MAXSPCLEN+2,sizeof(float));

nges = (float) (z-a);    /* this we know for shure */
sumx = 0.0;
sumxx = 0.0;
sumy = 0.0;
sumyy = 0.0;
sumxy = 0.0;

for(i=a;i<z;i++) {
   x = (float) i;
   y = spc[i];
   y = y - ((float) n);
   if(y<0.1) y=1.0;                 /* trap overflow on log */
   y = (float) log((double) y);     /* linearization of exp function */
   ylog[i] = y;                     /* store for later use */
   sumx = sumx + x;
   sumxx = sumxx + (x*x);
   sumy = sumy + y;
   sumyy = sumyy + (y*y);
   sumxy = sumxy + (x*y);
}
sumx = sumx / nges;
sumy = sumy / nges;
sumxx = sumxx - (nges * sumx * sumx);
sumxy = sumxy - (nges * sumx * sumy);
sumyy = sumyy - (nges * sumy * sumy);
steigung = sumxy / sumxx;
offset = sumy - (steigung * sumx);

/* calculating chi**2 */

erg=0.0;
for(i=a;i<z;i++) {
   x = (float) i;
   y = spc[i];
   y = ylog[i];
   ysoll = (steigung * x) + offset;
   y = y - ysoll;
   erg = erg + (y * y); 
}
free(ylog);
return(erg);
}

   
/* ---------------------------------------------------------
   Fitting background to exponential function.
   This is done, by calculating the chi**2 with
   linear regression for the logarithm of the spectrum.
   We than look for a minimum of chi**2 for different
   backgrounds.
   UP TO NOW, THE ERRORBARS ARE NOT INCLUDED !
   I HAVE FORGOTTEN MY OLD SUPERBASIC SOURCE AT HOME, SORRY.
   --------------------------------------------------------- */
expbafit(spc,err,a,z)
float spc[],err[];
int a,z;
{
int i,n,bamin,bamax,babest,inc,xraster,xdepth;
float chisq,chibest;

xraster = 10;
xdepth = 5;
chibest = 1E10;
bamax = sumspc(spc,z-10,z) / 10 ; /* get maximum background */
bamin = bamax / 2;                /* get minimum background */
for(i=1;i<xdepth;i++) {
   inc = (bamax - bamin) / xraster;
   if(inc<1) break;
   n = bamin; while(n<=bamax) {
      chisq = linreg(spc,err,n,a,z); /* just calculate chi**2 */
      if(chisq < chibest) {
	 chibest = chisq;
	 babest = n;
	 if(flg_v) {
	 printf("new best: background = %d   chi**2 = %f\n",babest,chibest);
	 }
      }
      n = n + inc;
   }
   bamin = babest - inc;
   bamax = babest + inc;
}
return(babest);
}