/* 
   Find Prompt peaks:
*/

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

float *spc,*err,*tim;
int range=1024,
    tri=0,
    dev=50;

help()
{
  printf("finding prompt peaks in file (first approximation)\n");
  printf("fprompt file [options]\n");
  printf("print out prompt peak positions\n");
  printf("options:\n");
  printf("  -range  n  specifies the length of each spectrum (%d)\n",range);
  printf("  -dev n     specifies the maximum deviation of the prompt (%d)\n",dev);
  printf("  -tri [1/2] triangular spectrum beginning with minimum: 1\n");
  printf("                                                peak:    2\n");
  printf("  -get file  get old spectra definition file (for reading\n");
  printf("                              the spectrum phase and parity)\n");
  printf("  -def file  specifies another definition file for output\n\n");
  printf("the default output file is %s\n",TMPFILE);
  printf("This file is read by the RVT 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("The last two numbers are used to specify a range for\n");
  printf("the RVT routine and the exponential fit.\n");
  printf("\n(C) Rainer Kowallik\n");
  exit(0);
}

main(argc,argv)
int argc;
char *argv[];
{
int xx,yy,n,m,i,max,peakptr,peaks[30],phase[30],xend[30];
char z[80],defnam[80];
float x,y;
FILE *fp;

   spc= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   err= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   tim= (float *)calloc(_MAXSPCLEN+2,sizeof(float));

   strcpy(defnam,TMPFILE);

   guessphase(phase,tri);
   getphase(defnam,phase,&tri,&range);

   if(checkopt(argc,argv,"-range",z)) {
      range=atoi(z);
      if(tri>0) range = 2 * range;
   }
   if(checkopt(argc,argv,"-dev",z)) dev=atoi(z);
   if(checkopt(argc,argv,"-tri",z)) {
      tri=atoi(z);
      if(tri>2) {
         printf("don't fool me , or I will fool you too !\n");
         exit(0);
      }
      if(tri>0) range = 2 * range;
      guessphase(phase,tri);
   }
   if(checkopt(argc,argv,"-def",z)) strcpy(defnam,z);

   if(checkopt(argc,argv,"-get",z)) getphase(z,phase,&tri,&range);

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

   /* finding FIRST peak */
   
   yy=0; xx=0;
   for(n=0;n<range;n++) {
      if(spc[n]>yy) {
         yy=spc[n];
         xx=n;
      }
   }
   peakptr=0;
   peaks[peakptr++]=xx;
   if(tri>0) peaks[peakptr++]=xx;

   /* finding all other peaks */

   n=xx+range-dev;
   while(n<max) {
      xx=n; yy=0;
      for(i=0;i<(2*dev);i++) {
         if(spc[n]>yy) {
            yy=spc[n];
            xx=n;
         }
         n = n + 1;
      }
      peaks[peakptr++]=xx;
      if(tri>0) peaks[peakptr++]=xx;

      n=xx+range-dev;
   }

   /* --------------------------------------------------------------------
     finding arithmetic mean of peak heights. Used to figure out, if the
     peak is real or not
   --------------------------------------------------------------------- */
   
   xx=0;
   for(i=0;i<peakptr;i++) {
      xx=xx+spc[peaks[i]];
   }
   xx = xx/peakptr;
   xx = xx / 10;

   /* --------------------------------------------------------------------
      removing random "peaks"
      -------------------------------------------------------------------- */
  
   i=0;
   while(i<peakptr) {
      if(spc[peaks[i]] < xx) {
	 for(n=i+1;n<peakptr;n++) peaks[n-1] = peaks[n];
	 peakptr = peakptr - 1;
      } else {
	 i=i+1;
      }
   }

   /* ---------------------------------------
      finding end of usable data for spectra
      --------------------------------------- */

   findusable(peaks,range,tri,phase,xend,peakptr);   
   n=xend[0];
   for(i=0;i<peakptr;i++) {
      if(xend[i]<n) n=xend[i];
   }
   xend[peakptr-1] = n;                   /* minimum must be the last entry */


   /* ---------------------------------------------------
      writing data 
   ------------------------------------------------------ */

   fp=fopen(defnam,"w");
   for(i=0;i<peakptr;i++) {
      fprintf(fp,"%d   %2d   %d   %d\n",peaks[i],phase[i],5,xend[i]);
   }
   fclose(fp);
   printf("\ndata written to %s\n now start editor on this file !\n",TMPFILE);
   free(spc); free(err); free(tim);
   exit(0);
}

iabs(x)
int x;
{
if(x>=0) return(x);
return(-x);
}


/* -----------------------------------------------
   reading spectra definition file
   ----------------------------------------------- */
getphase(name,phase,tri,range)
char name[];
int phase[], *tri, *range;
{
int i,n,m,p[30];
FILE *fp;

   fp=fopen(name,"r");
   if(fp==NULL) {
      return(0);
   }
   i=0; *tri=0;
   while(!feof(fp)) {
      fscanf(fp,"%d %d %d %d\n",&p[i],&phase[i],&n,&n);
      i=i+1;
   }
   fclose(fp);
   for(n=0;n<i;n++) if(phase[i]<0) *tri=2;   /* triangular spectrum ? */
   if(phase[0]<0) *tri=1;        /* triangular, beginning with minimum */
   n=p[0];
   if(p[1]>(n+5)) {
      n = iabs(n - p[1]);
   } else {
      n = iabs(n - p[2]);
   }
   i=1;
   while(i<n) i= i << 1;
   m = i >> 1;
   if(iabs(n-m) < iabs(n-i)) {
      n = m;
   } else {
      n = i;
   }
   *range = n;
}

guessphase(phase,tri)
int phase[], tri;
{
int i,n,m;

   n=1;
   for(i=0;i<30;i++) {
      n= -1 * n;
      switch(tri) {
      case 0:
         if(n>0) phase[i]=1;
         if(n<0) phase[i]=2;
         break;
      case 1:
         if(n>0) {
            phase[i++] = -1;
            phase[i]=1;
         }
         if(n<0) {
            phase[i++]= -2;
            phase[i]=2;
         }
         break;
      case 2:
         if(n>0) {
            phase[i++]=1;
            phase[i]= -1;
         }
         if(n<0) {
            phase[i++]=2;
            phase[i]= -2;
         }
         break;
      }
   }
}

/* ------------------------------------------------
   finding usable data area
   ------------------------------------------------ */
findusable(peaks,range,tri,phase,xend,peakptr)
int peaks[],phase[],xend[],
    range,tri,peakptr;
{
int i,n,m,start,end;

   start = range / 2;
   end = range;
   if(tri>0) {
      start = range / 4;
      end = range / 2;
   }
   end = end - 3;
   for(n=0;n<peakptr;n++) {  /* scanning spectra */
      xend[n] = end;
      if(phase[n] > 0) {    /* positive time axis */
         m = peaks[n] + start;
         for(i=start;i<end;i++) {
            m = m + 1;
            if(!outrange(spc[m],err[m],spc[m+1],err[m+1])) continue;
            if(!outrange(spc[m],err[m],spc[m+2],err[m+2])) continue;
            if(!outrange(spc[m-1],err[m-1],spc[m+1],err[m+1])) continue;
            xend[n] = i - 2 ;
         }
      } else {
         m = peaks[n] - start;
         for(i=start;i<end;i++) {
            m = m - 1;
            if(!outrange(spc[m],err[m],spc[m-1],err[m-1])) continue;
            if(!outrange(spc[m],err[m],spc[m-2],err[m-2])) continue;
            if(!outrange(spc[m+1],err[m+1],spc[m-1],err[m-1])) continue;
            xend[n] = i - 2 ;
         }
      }
   }
}

outrange(x1,e1,x2,e2)
float x1,e1,x2,e2;
{
float a,b,c;

   a = (float) fabs((double)(x1-x2));
   b = 3 * (e1 + e2);
   if(a>b) return(TRUE);
   return(FALSE);
}