/*
   calibrate energy scale of auger spectrum
*/

#define ECALTMP "ecal.auger"

#include <stdio.h>
#ifdef AMIGA
#include <gfxamiga.h>
#else
#include <gfx.h>
#endif

float x,y,*spc,*err,*tim;
int   fdlockin, delay;
float stepv, startv, stopv;
float ecal;
int spcmax;
char *cname[110];
char *lname[110];
float *Eauger, *Ebind, *sens, *Ratom;
int offset;


int find_peak(Epos)
float Epos;
{
float    pmin, pmax, x, y, mean, diff;
int      k, ka, kmin, kmax, isl, islm, isr, isrm;

/* search minimum */
   ka = (int) (Epos / _tica) - offset;
   isl = (int) ((5.0 / _tica) / 2.0);
   islm = (int) (5.0 / _tica);
   isr = (int) ((5.0 / _tica) / 2.0);
   isrm = (int) (5.0 / _tica);

   while((isl < islm) && (isr < isrm)) {
      if(((ka - isl) < 1) || ((ka + isr) > spcmax)) return(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);

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

   kmax = kmin;
   pmax = pmin;
   while(isl < islm) {
      if((kmax - isl) < 1) return(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 = _tica * ((float) (kmin + offset));
   y = _tica * ((float) (kmax + offset));

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

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

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

   return(kmin);
}


float find_energy(z)
char *z;
{
int i,n,m;

   z[0] = z[0] & 223;               /* = toupper ! */
   if(z[1] > 0) z[1] = z[1] | 32;   /* = tolower ! */
   for(n = 0; n < 110; n++) {
      if(strcmp(cname[n],z) == 0) break;
   }
   if(n >= 109) {
      printf("sorry, element >%s< not found\n");
      exit(0);
   }
   return(Eauger[n]);
}


help()
{
  printf("augcal spectrumfile [-chan n] [-cur] -energy n.m\n");
  printf("       change energy calibration of specified file\n");
  printf("  -chan n           specify channel\n");
  printf("  -cur              use cursor to select the channel\n");
  printf("  -energy n.m       set energy of peak\n");
  printf("  -line Xx          use auger energy of element Xx as energy\n");
  printf("  -auto Xx          perform auto calibration for element Xx\n");
  printf("  -table filename   uses a custom table for elements\n");
  printf("  -v                verbouse: print out calibration factor\n");
}

float frline(stream)
FILE *stream;
{
int i,n;
char c,s[80];

   fgets(s,80,stream);
   return(atosf(s));
}

main(argc,argv)
int argc;
char *argv[];
{
int n,m,i;
char *z,*comment, *elements;
int channel;
int flag_v;
float energy;
float real_t,real_y;
float a,b;
FILE *fp;

   flag_v = 0;
   energy = 1.0;
   channel = 1;

   if(argc < 4) {help(); exit(0);}

   fdlockin = auxopen("auger");
   delay = atoi(auxparams[4]);
   stepv = atosf(auxparams[5]);
   startv = atosf(auxparams[6]);
   stopv = atosf(auxparams[7]);
   ecal = atosf(auxparams[8]);
   close(fdlockin);

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

   spcmax=readspec(argv[1],spc,err,tim,comment);
   a = tim[1] - tim[0]; b = tim[0];
   offset = (int) (b / a);

   if(checkopt(argc,argv,"-v",z)) flag_v = 1;
   if(checkopt(argc,argv,"-chan",z)) channel = atoi(z);
   if(checkopt(argc,argv,"-cur",z)) {
      getcur(&channel,&real_t,&real_y);
      if(flag_v == 1) printf("channel = %d\n",channel);
      channel = channel - offset;
   }
   if(checkopt(argc,argv,"-energy",z)) energy = atosf(z);
   if(checkopt(argc,argv,"-table",z)) strcpy(elements,z);
   if(checkopt(argc,argv,"-line",z)) {
      read_elements(elements);
      energy = find_energy(z);
   }
   if(checkopt(argc,argv,"-auto",z)) {
      read_elements(elements);
      energy = find_energy(z);
      channel = find_peak(energy);
      if(channel == 0) {
         printf("could not find peak !\nNo calibration performed !\n");
         exit(0);
      }
   }

   _tica = energy / (channel + offset);
   for(n = 0; n <= spcmax; n++) tim[n] = _tica * ((float) (n + offset));
   writespec(argv[1],spc,err,spcmax,2,comment);
   strcpy(z,argv[1]); strcat(z,".tim");
   writespec(z,tim,err,spcmax,2,comment);

   ecal = _tica / stepv;
   if(flag_v == 1) printf("calibration is %f eV/V\n",ecal);

   fp = fopen(ECALTMP,"w");
   if(fp != NULL) {
      fprintf(fp,"%f\n",ecal);
   }
   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);
}

getcur(channel,time,value)
int *channel;
float *time, *value;
{
float tmin,
      tmax,
      ymin,
      ymax,
      real_t,
      real_y,
      yfak,
      tfak,
      xoffset,
      yoffset;
int   leftmargin,
      bottommargin,
      rightmargin,
      topmargin,
      _win_flg;
int x,xx,y,yy;
FILE *fp;

      tekopen();
      getxy(&x,&y);
      close(_tek4014);
      /* calculating channel and time equivalents */
      fp=fopen(ASHBDRY,"r");
      leftmargin = frline(fp);
      bottommargin = frline(fp);
      rightmargin = frline(fp);
      topmargin = frline(fp);
      _win_flg = frline(fp);
      tmin = frline(fp);
      tmax = frline(fp);
      ymin = frline(fp);
      ymax = frline(fp);
      _tica = frline(fp);
      fclose(fp);
      yfak=(topmargin-bottommargin)/(ymax-ymin);
      tfak=(rightmargin-leftmargin)/(tmax-tmin);
      xoffset= ((- tmin) * tfak) + leftmargin;
      yoffset= ((- ymin) * yfak) + bottommargin;

      real_t= (x - xoffset)/tfak;
      real_y= (y - yoffset)/yfak;
      xx = real_t / _tica;
      *channel = xx;
      *time = real_t;
      *value = real_y;
}

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