/* analyse auger spectrum
*/

#include <stdio.h>
#include <spec.h>

float x,y,*spc,*err,*tim;
float ecal;
int spcmax;
char *cname[110];
char *lname[110];
float *Eauger, *Ebind, *sens, *Ratom;
float intensity[110];
float deviation = 0.1;
int   elementnr[110];
int   iptr = 0;
int   nel = 4;

/* ---------------------------------------
      function returns peakheight or 0
   --------------------------------------- */
float augpeak(Epos,Emeas)
float Epos;
float *Emeas;
{
float    pmin, pmax, x, y, mean, diff;
float    a,b;
int      offset;
int      k, ka, kmin, kmax, isl, islm, isr, isrm;

   a = tim[1] - tim[0]; b = tim[0];
   offset = (int) (b / a);

/* search minimum */
   ka = ((int) (Epos / a)) - offset;
   if(ka <= 0) return(0.0);
   isl = (int) ((5.0 / a) / 2.0);
   islm = (int) (5.0 / a);
   isr = (int) ((5.0 / a) / 2.0);
   isrm = (int) (5.0 / a);

   while((isl < islm) && (isr < isrm)) {
      if(((ka - isl) < 1) || ((ka + isr) > spcmax)) return(0.0);
      pmin = spc[ka];
      kmin = ka; mean = 0.0;
      for(k = (ka - isl); k <= (ka + isr); k++) {
         mean = mean + spc[k];
         if(spc[k] < pmin) {
            pmin = spc[k];
            kmin = k;
         }
      }
      mean = mean / ((float) (isr + isl + 1));
      if(spc[ka - isl - 1] <= pmin) {
         isl = isl + 1;
         continue;
      }
      if(spc[ka + isr + 1] <= pmin) {
         isr = isr + 1;
         continue;
      }
      diff = ((mean - pmin) / mean); if(diff < 0.0) diff = -1.0 * diff;
      if(diff < 0.01) {
         isr = isr + 1;
         isl = isl + 1;
         continue;
      }
      break;
   }

   diff = ((mean - pmin) / mean); if(diff < 0.0) diff = -1.0 * diff;
   if(diff < 0.01) return(0.0);

/* search maximum */
   isl = (int) (8.0 / a);
   islm = (int) (12.0 / a);

   kmax = kmin;
   pmax = pmin;
   while(isl < islm) {
      if((kmax - isl) < 1) return(0.0);
      for(k = (kmax - isl); k <= kmax; k++) {
         if(spc[k] > pmax) {
            pmax = spc[k];
            kmax = k;
         }
      }
      if(spc[kmin - isl - 1] >= pmax) {
         isl = isl + 1;
         continue;
      }
      break;
   }
   x = tim[kmin];
   y = tim[kmax];

   *Emeas = x;

/* Difference between maximum and minimum should not exceed 20 eV*/
   diff = x - y ; if(diff < 0.0) diff = -1.0 * diff;
   if(diff > 20.0) return(0.0);

/* Difference between measured peak and value of table should not exceed 10% */
   diff = x - Epos; if(diff < 0.0) diff = -1.0 * diff;
   if((diff / Epos) > deviation) return(0.0);

   if(kmax == kmin) return(0.0);

   y = pmax - pmin;
   return(y);
}


help()
{
  printf("anaaug spectrumfile [-cal n.m]\n");
  printf("analyse auger spectrum using peakpositions and\n");
  printf("peak intensities\n");
  printf("   -cal  n.m        specifies a alternate energy calibration\n");
  printf("   -table filename  uses a custom table for elements\n");
  printf("   -find list       give a list of elements to search for\n");
  printf("   -nel n           sets number of elements which contribute\n");
  printf("   -dev n.m         sets maximum deviation between energy from\n");
  printf("                    table and energy found by peak finder.\n");
  printf("                    default is %f\n\n",deviation);
}

main(argc,argv)
int argc;
char *argv[];
{
int n,m,i;
char *z,*comment,*elements;
FILE *fp;

   Eauger = (float *) calloc(110,sizeof(float));
   Ebind = (float *) calloc(110,sizeof(float));
   sens = (float *) calloc(110,sizeof(float));
   Ratom = (float *) calloc(110,sizeof(float));

   z = (char *) malloc(80);
   comment = (char *) malloc(80);
   elements = (char *) malloc(80); strcpy(elements,"elements");
   spc= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   err= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   tim= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   if(tim==NULL) {
      printf("sorry, not enough memory\n");
      exit(-1);
   }

   for(n = 0; n < 110; n++) {
      cname[n] = (char *) malloc(4); cname[n][0] = 0;
      lname[n] = (char *) malloc(16); lname[n][0] = 0;
      sens[n] = 0.0;
      Eauger[n] = 0.0;
      Ebind[n] = 0.0;
      Ratom[n] = 0.0;
   }

   if(checkopt(argc,argv,"-cal",z)) ecal = atoi(z);
   if(checkopt(argc,argv,"-table",z)) strcpy(elements,z);
   nel = 4;
   if(checkopt(argc,argv,"-nel",z)) nel = atoi(z);
   deviation = 0.1;
   if(checkopt(argc,argv,"-dev",z)) deviation = atosf(z);

   read_elements(elements);

   if(checkopt(argc,argv,"-find",z)) restrictto(z);

   spcmax = readspec(argv[1],spc,err,tim,comment);
   analyse();
   free(err); free(spc); free(tim);
   free(elements); free(comment); free(z);
   for(n = 0; n < 110; n++) {
      free(lname[n]);
      free(cname[n]);
   }
   free(Ratom); free(sens); free(Ebind); free(Eauger);
   exit(0);
}

analyse()
{
int i,j,n,m;
float Peakh, Epos, Emeas;
float z, sumintens;

   iptr = 0;
   sumintens = 0.0;
   for(n = 0; n < 89; n++) {
      if(Eauger[n] > 0.0) {
         Epos = Eauger[n];
         Peakh = augpeak(Epos,&Emeas);
         if(Peakh > 0.0) {
            z = (Peakh / sens[n]);
            intensity[iptr] = z;
            elementnr[iptr] = n;
            iptr = iptr + 1;
            sumintens += z;
/*
            printf("\n%f %s at %f (%f)\n",Peakh,cname[n],Emeas,Epos);
*/
         }
      }
   }

/* sort peaks for intensity */
   for(n = 0; n < (iptr - 1); n++) {
      i = 1;
      for(m = 0; m < ((iptr - n) - 1); m++) {
         if(intensity[m] < intensity[m+1]) {
            z = intensity[m];
            intensity[m] = intensity[m+1];
            intensity[m+1] = z;
            j = elementnr[m];
            elementnr[m] = elementnr[m+1];
            elementnr[m+1] = j;
            i = 0;
         }
      }
      if(i == 1) break;
   }
/* print out result */
   n = iptr;
   if(n > nel) n = nel;
   sumintens = 0; for(i = 0; i < n; i++) sumintens += intensity[i];
   for(i = 0; i < n; i++) {
      z = 100.0 * intensity[i] / sumintens;
      j = elementnr[i];
      augpeak(Eauger[j],&Emeas);
      printf("%3.1f%s  %14s  Eauger = %4.2f  Epeak = %4.2f\n",
               z,"%",lname[j],Eauger[j],Emeas);
   }
}

read_elements(elements)
char *elements;
{
int i,n,m;
char z[80];
FILE *fp;

   fp = fopen(elements,"r");
   if(fp == NULL) {
      fprintf(stderr,"file >%s< not found !\n",elements);
      exit(0);
   }
   while(!feof(fp)) {
      fgets2(z,80,fp);
      if(feof(fp)) break;
      if(z[0] == ';') continue;
      if(z[0] < 33) continue;
      n = atoi(z);
      if((n < 1) || (n > 109)) {
         fprintf(stderr,"error in >%s< unexpected >%s<\n",elements,z);
         exit(0);
      }
      fgets2(z,80,fp); strcpy(cname[n],z);
      fgets2(z,80,fp); strcpy(lname[n],z);
      fgets2(z,80,fp); sens[n] = atosf(z);
      fgets2(z,80,fp); Eauger[n] = atosf(z);
      fgets2(z,80,fp); Ebind[n] = atosf(z);
      fgets2(z,80,fp); Ratom[n] = atosf(z);
   }
   fclose(fp);
}


fgets2(s,n,fp)
char *s;
int n;
FILE *fp;
{
int i,l;

   fgets(s,n,fp);
   l = strlen(s);
   for(i = 0; i <= l; i++) {
      if(s[i] < 33) s[i] = 0;
      if(s[i] == ';') s[i] = 0;
   }
}

/* disable all elements, which are not contained in the list */
restrictto(s)
char *s;
{
char elist[10][4];
char c;
int l,i,n,m;

   for(n = 0; n < 10; n++) for(i = 0; i < 4; i++) elist[n][i] = 0;
   l = strlen(s);
   n = 0; i = 0; m = 0;
   if(l < 1) error1(s,i);
   while(c = s[i++] , c != 0) {
      elist[n][m++] = c;
      if(m > 3) error1(s,i);
      if((c == ',') || (c == '.') || (c == ' ') || (c == '-')) {
         elist[n][m - 1] = 0;
         m = 0; n = n + 1;
      }
   }
   nel = n + 1;   /* set number of elements to find */

   for(n = 0; n < 110; n++) {
      m = 0;
      for(i = 0; i < nel; i++) {
         if(strcmp(elist[i],cname[n]) == 0) m = 1;
      }
      if(m == 0) {
         Eauger[n] = 0.0;
         sens[n] = 999999.9;
      }
   }
}

error1(s,i)
char *s;
int i;
{
printf("error in element list >%s< at %d\n",s,i);
printf("String must consist of chemical names separated by colon\n");
exit(0);
}