/*
   Exponential Fit
*/

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

float *spc, *err, *tim, *usp1, *usp2, *uer1, *uer2;
int   flg_v;
int   bgrnd,xbeg,xend;
float chisq,fit_a, fit_b;

help()
{
  printf("\n\nExponential Fit y = c + b * exp(a * x)\n");
  printf("ExpFit file [-o file] [-ba n] [-xbeg n] [-xend n] [-last n]\n");
  printf("options:\n");
  printf("   file      spectrum to be fitted\n");
  printf("   -o file   output spectrum instead of stdout\n");
  printf("   -ba  n    background for spectrum (to avoid Background fit)\n");
  printf("   -xbeg n   first channel to fit\n");
  printf("   -xend n   last channel to fit\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 parameters\n");
  printf("   -vv       very verbose mode. Also Prints fitting progress\n");
  printf("\n(C) Rainer Kowallik\n");
  exit(0);
}

/* ---------------------------------------------------------
   Liner 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,tim,err,n,a,z)
float spc[],tim[],err[];
int n,a,z;
{
int i,m,k;
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);    /* Number of Elements */
   sumx = 0.0;
   sumxx = 0.0;
   sumy = 0.0;
   sumyy = 0.0;
   sumxy = 0.0;

   for(i=a;i<z;i++) {
      x = tim[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 = tim[i];
      y = ylog[i];
      ysoll = (steigung * x) + offset;
      y = y - ysoll;
      erg = erg + (y * y); 
   }

   fit_a = steigung;
   fit_b = offset;

   free(ylog);
   return(erg);
}

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

main(argc,argv)
int argc;
char *argv[];
{
int  n,m,i,max,flg_ba,lastn;
char z[80],comment[80],nam[80];
float x,alpha,sum1,sum2;
FILE *fp;

   flg_ba=0;            /* defaults to exponential fit */
   bgrnd = -1;
   flg_v=FALSE;
   xbeg = 0;
   xend = -1;

   if(argc < 2) {help(); exit(0);}
   strcpy(nam,argv[1]); 
   if(checkopt(argc,argv,"-ba",z)) {
      bgrnd=atoi(z); flg_ba=1;
   }
   if(checkopt(argc,argv,"-xbeg",z)) {
      xbeg = atoi(z);
   }
   if(checkopt(argc,argv,"-xend",z)) {
      xend = atoi(z);
   }
   if(checkopt(argc,argv,"-last",z)) {
      flg_ba=2; lastn = atoi(z);
   }
   if(checkopt(argc,argv,"-v",z)) flg_v=1;
   if(checkopt(argc,argv,"-vv",z)) flg_v=2;

   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(nam,usp1,uer1,tim,comment);
   if(xend < 0) xend = max - 1;

/* -----------------------------------------------
   generating error spectrum 
   ----------------------------------------------- */
   for(n=0;n<max;n++) {
      uer1[n] = sqrt(usp1[n]);
   }

/* --------------------------------------------
   finding correct background
   -------------------------------------------- */
 
   switch(flg_ba) {
   case 0:
      bgrnd = expbafit(usp1,tim,uer1,xbeg,xend);
      chisq = linreg(usp1,tim,uer1,bgrnd,xbeg,xend);
      break;
   case 1: /* allready done when parsing arguments */
      chisq = linreg(usp1,tim,uer1,bgrnd,xbeg,xend);
      break;
   case 2:
      bgrnd = sumspc(usp1,max-lastn,max) / lastn;
      chisq = linreg(usp1,tim,uer1,bgrnd,xbeg,xend);
      break;
   }

   fit_b = exp(fit_b);

   if(flg_v) {
      printf("CHI**2 = %f\n",chisq);
      printf("background = %d\n",bgrnd);
      printf("b (t0)     = %f\n",fit_b);
      printf("a (Tau)    = %f\n",-1.0/fit_a);
   }

/* --------------------------------------------
   calculating Fit curve
   -------------------------------------------- */
   
   for(n = 0; n < max; n++) {
      spc[n] = bgrnd + (fit_b * exp(fit_a * tim[n]));
   }


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

   strcpy(z,comment);
   n=instr("|",comment);
   if(n>0) midstr(z,comment,0,n-1);
   strcat(z," | ExpFit");
   writespec("",spc,err,max,2,z);
   free(uer1); free(uer2); free(usp1); free(usp2);
   free(spc); free(err); free(tim);
   exit(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);
}

   
/* ---------------------------------------------------------
   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,tim,err,a,z)
float spc[],tim[],err[];
int a,z;
{
int i,n,erg,bamin,bamax,babest,inc,xraster,xdepth;
float 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,tim,err,n,a,z); /* just calculate chi**2 */
      if(chisq < chibest) {
	 chibest = chisq;
	 babest = n;
	 if(flg_v > 1) {
	 printf("new best: background = %d   chi**2 = %f\n",babest,chibest);
	 }
      }
      n = n + inc;
   }
   bamin = babest - inc;
   bamax = babest + inc;
}
return(babest);
}