/*
   ratio function:
   sp = (usp1 - usp2) / (usp1 + usp2)
*/

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

float *spc, *err, *tim, *usp1, *usp2, *uer1, *uer2;
int   flg_v;

help()
{
  printf("ratio function: sp = (a*usp1 - b*usp2) / (a*usp1 + b*usp2)\n");
  printf("rvt [-s1 file1] [-s2 file2] [-o file] [-ba1 n] [-ba2 n] [-last n]\n");
  printf("options:\n");
  printf("   -s1 file1    0 degree spctrum instead of USP1\n");
  printf("   -s2 file2    90 degree spectrum instead of USP2\n");
  printf("   -o  file     output spectrum instead of stdout\n");
  printf("   -ba1 n       background for 0 degree spectrum\n");
  printf("   -ba2 n       background for 90 degree spectrum\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    set spectra definition file other then global.def\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);
}

main(argc,argv)
int argc;
char *argv[];
{
int  n,i,max,flg_ba,lastn,ba1,ba2,xbeg,xend;
char z[80],comment[80],nam1[80],nam2[80],defnam[80];
float x,alpha,sum1,sum2;
FILE *fp;

   strcpy(nam1,"USP1");
   strcpy(nam2,"USP2");
   strcpy(defnam,TMPFILE);
   flg_ba=0;            /* defaults to exponential fit */
   ba1 = -1; ba2 = -1;
   flg_v=FALSE;

   if(checkopt(argc,argv,"-s1",z)) strcpy(nam1,z); 
   if(checkopt(argc,argv,"-s2",z)) strcpy(nam2,z); 
   if(checkopt(argc,argv,"-def",z)) strcpy(defnam,z);
   if(checkopt(argc,argv,"-ba1",z)) {
      ba1=atoi(z); flg_ba=1;
      ba2=ba1;
   }
   if(checkopt(argc,argv,"-ba2",z)) {
      ba2=atoi(z); flg_ba=1;
      if(ba1==-1) ba1=ba2;
   }
   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(nam1,usp1,uer1,tim,comment);
   max=readspec(nam2,usp2,uer2,tim,comment);

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

/* -----------------------------------------------
   reading spectra definition file
   ----------------------------------------------- */

   fp=fopen(defnam,"r");
   while(!feof(fp)) {
      fscanf(fp,"%d %d %d %d\n",&i,&i,&xbeg,&xend);
   }
   if(flg_v) printf("xbeg = %d, xend = %d\n",xbeg,xend);
   fclose(fp);


/* --------------------------------------------
   finding correct background
   -------------------------------------------- */
 
   switch(flg_ba) {
   case 0:
      ba1=expbafit(usp1,uer1,xbeg,xend);
      ba2=expbafit(usp2,uer2,xbeg,xend);
      break;
   case 1: /* allready done when parsing arguments */
      break;
   case 2:
      ba1 = sumspc(usp1,max-lastn,max) / lastn;
      ba2 = sumspc(usp2,max-lastn,max) / lastn;
      break;
   }
   if(flg_v) printf("ba1=%d   ba2=%d\n",ba1,ba2);

/* --------------------------------------------
   Subtracting background
   -------------------------------------------- */
   for(n=0;n<max;n++) {
      usp1[n] = usp1[n] - ba1;
      usp2[n] = usp2[n] - ba2;
   }

/* --------------------------------------------
   calculating normalization
   -------------------------------------------- */
   sum1 = (float) sumspc(usp1,xbeg,max);
   sum2 = (float) sumspc(usp2,xbeg,max);
   alpha = sd(sum1 , sum2);
   if(flg_v) {
      printf("sum usp1 = %f\nsum usp2 = %f\n",sum1,sum2);
      printf("alpha    = %f\n",alpha);
   }
/* --------------------------------------------
   Ratio function
   -------------------------------------------- */
   for(n=0;n<max;n++) {

      usp2[n] = alpha * usp2[n];
      x = usp1[n] + usp2[n];
      spc[n] = sd(usp1[n] - usp2[n] , x);

      sum1 = sd(2*usp2[n] , x*x);
      sum2 = sd(2*alpha*usp1[n] , x*x);
      sum1 = sum1 * sum1;
      sum2 = sum2 * sum2;

      err[n] = sqrt((sum1 * usp1[n]) + (sum2 * usp2[n]));
   }
/* -------------------------------------------
   now we can go and save the new spectra 
   ------------------------------------------- */

   strcpy(z,comment);
   n=instr("|",comment);
   if(n>0) midstr(z,comment,0,n-1);
   strcat(z," | RVT");
   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);
}

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