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

#define MAXC 10
#define CEXP -1
#define NEXP -2

float *x, *y, *y_calc, *resid;
float coef[MAXC], sig[MAXC];
int   nrow, ncol;
float correl_coef;
float *spc, *err, *tim;

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

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

square(x, y, a, g, nrow, ncol)
   float  *x[MAXC];
   float  *y;
   float  a[MAXC][MAXC];
   float  g[MAXC];
   int   nrow, ncol;
{
int  i,l,k;

   for (k = 0; k <= ncol; k++) {
      for (l = 0; l <= k; l++) {
         a[k][l] = 0.0;
         for (i = 0 ; i <= nrow ; i++) {
            a[k][l] = a[k][l] + x[i][l] * x[i][k];
            if (k != l)   a[l][k] = a[k][l];
         }
      }
      g[k] = 0.0;
      for (i = 0; i <= nrow ; i++)  g[k] = g[k] + y[i] * x[i][k];
   }
}

swap(a,b)
   float   *a, *b;
{
   float   hold;

   hold = *a; *a = *b;
   *b = hold;
}

gaussj(b, y, coef, ncol, error)
   float  b[MAXC][MAXC];
   float  y[MAXC];
   float  coef[MAXC];
   int   ncol;
   int   *error;
{
float w[MAXC][MAXC];
int   index[MAXC][MAXC];
int   irow, icol, n;
int   l1, l, k, j, i;
float determ, pivot, hold, sum, t, ab, big;

   *error = FALSE;
   n = ncol;
   for (i = 0; i <= n ; i++) {
      w[i][0] = y[i];
      index[i][2] = 0;
   }
   determ = 1.0;
   for (i = 0; i <= n ; i++) {
      big = 0.0;
      for (j = 0; j <= n ; j++) {
         if (index[j][2] != 1) {
            for (k = 0; k <= n ; k++) {
               if (index[k][2] > 1) {
                  printf("fehler: Matrix singulaer\n");
                  *error = TRUE;
                  return(0);
               }
               if (index[k][2] < 1) {
                  if (((float)fabs((double)b[j][k])) > ((float)big)) {
                     irow = j; icol = k;
                     big = (float)fabs((double)b[j][k]);
                  }
               }
            }
         }
      }
      index[icol][2] = index[icol][2] + 1;
      index[i][0] = irow;
      index[i][1] = icol;
      if (irow != icol) {
         determ = -determ;
         swap(b[irow][0], b[icol][0]);
         swap(w[irow][0], w[icol][0]);
      }
      pivot = b[icol][icol];
      determ = determ * pivot;
      b[icol][icol] = 1.0;
      for (l = 0; l <= n; l++) b[icol][l] = b[icol][l] / pivot;
      w[icol][0] = w[icol][0] / pivot;
      for (l1 = 0; l1 <= n ; l1++) {
         if (l1 != icol) {
            t = b[l1][icol];
            b[l1][icol] = 0.0;
            for (l = 0; l <= n ; l++) b[l1][l] = b[l1][l] - b[icol][l] * t;
            w[l1][0] = w[l1][0] - w[icol][0] * t;
         }
      }
   }
   for (i = 0 ; i <= n ; i++) {
      l = n - i;
      if (index[l][0] != index[l][1]) {
         irow = index[l][0];
         icol = index[l][1];
         for (k = 0; k <= n ; k++) swap(b[k][irow], b[k][icol]);
      }
   }
   for (k = 0; k <= n ; k++) {
      if (index[k][2] != 1) {
         printf("Fehler: Matrix singulaer\n");
         *error = TRUE;
         return(0);
      }
   }
   for (i = 0 ; i <= n ; i++) coef[i] = w[i][0];
}

get_data(xarr, yarr)
   float xarr[], yarr[];
{
int   nrow,i,n,x,y,xx,yy,
      leftmargin,
      bottommargin,
      rightmargin,
      topmargin,
      _win_flg;
float tmin,
      tmax,
      ymin,
      ymax,
      real_t,
      real_y,
      yfak,
      tfak,
      xoffset,
      yoffset;
float a,b;
char  s[80];
FILE *fp;

   nrow = -1;
   while(1 == 1) {
      printf("%d: X [fp-number | Q<uit> | C<ursor>] x: ",nrow+1); fflush(stdout);
      fgets(s,80,stdin);
      if(toupper(s[0]) == 'Q') break;
      if(toupper(s[0]) == 'C') {

         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;
         a= (x - xoffset)/tfak;
         b= (y - yoffset)/yfak;
         xx= a / _tica;
      }
      if(s[0] <'A') {
         a = atosf(s);
         printf("%d: Y [fp-number]                  y: "); fflush(stdout);
         fgets(s,80,stdin);
         b = atosf(s);
      }
      nrow = nrow +1;
      xarr[nrow] = a;
      yarr[nrow] = b;
printf("x = %f    y = %f\n",a,b);
   }
   return(nrow);
}

write_data()
{
    int   i;

/* ------------------
   printf("nrow = %d\n\n", nrow);
   printf(" I     X         Y      Y Ber.   Resid\n");
   for (i = 0; i <= nrow ; i++) {
      printf("%d %f %f %f %f\n", i, x[i], y[i], y_calc[i], resid[i]);
   }
------------------- */
   printf("Koeffizienten  Fehler\n");
   printf("%f     %f    konstanter term\n", coef[0], sig[0]);
   for (i = 1; i <= ncol; i++) {
      printf("%f     %f\n", coef[i], sig[i]);
   }
   printf("korrelationskoeffizient %f\n", correl_coef);
}


polyfit(x, y, y_calc, resid, coef, sig, nrow, ncol)
   float  *x, *y, *y_calc, *resid, *coef, *sig;
   int   nrow, ncol;
{
float *xmatr[MAXC];
float a[MAXC][MAXC];
float g[MAXC];
int   error;
int   nm;
int   j;
int   i;
float   xi, yi, yc, srs, see, sum_y, sum_y2;

   for(i = 0; i < MAXC; i++) {
      xmatr[i] = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
      if(xmatr[i] == NULL) {
         fprintf(stderr,"not enough memory, exiting\n");
         exit(-1);
      }
   }

   for (i = 0; i <= nrow; i++) {
      xi = x[i];
      xmatr[i][0] = 1.0;
      for (j = 1; j <= ncol; j++) xmatr[i][j] = xmatr[i][j - 1] * xi;
   }
   square(xmatr, y, a, g, nrow, ncol);
   gaussj(a, g, coef, ncol, &error);
   sum_y = 0.0;
   sum_y2 = 0.0;
   srs = 0.0;
   for (i = 0; i <= nrow; i++) {
      yi = y[i];
      yc = 0.0;
      for (j = 0; j <= ncol; j++) yc = yc + coef[j] * xmatr[i][j];
      y_calc[i] = yc;
      resid[i] = yc - yi;
      srs = srs + (resid[i] * resid[i]);
      sum_y = sum_y + yi;
      sum_y2 = sum_y2 + yi * yi;
   }
   correl_coef = (float)sqrt((double)(1.0 - srs / (sum_y2 - (sum_y * sum_y) / nrow)));
   if (nrow == ncol) {
      nm = 1;
   } else {
      nm = nrow - ncol;
   }
   see = (float) sqrt((double)(srs / (float)nm));
   for (i = 0; i <= ncol; i++) {
       sig[i] = see * (float)sqrt((double)a[i][i]);
   }
   for(i = 0; i < MAXC; i++) free(xmatr[i]);
}

help()
{
printf("Linear or Polynomial Fit to Data.\n");
printf("                        Highest possible order: %d\n",MAXC);
printf("LinFIT [-in file] [-o file] [-dif file] [-v] [-poly n]\n");
printf("options:\n");
printf("    -in filename   reads data from spectrum file\n");
printf("    -o  filename   output calculated fitcurve to spectrum file\n");
printf("    -dif filename  diferentiate fitcurve and write data to spectrum\n");
printf("    -v             verbouse: prints out resulting coeffitients\n");
printf("    -poly n        calculate polynom of order n (default is 3)\n");
printf("    -cexp n.m      calculate critical exponent for Tc = n.m\n");
printf("                    Formula: y = a * (Tc - x) ^b\n");
printf("                    if Tc is set to 0, x is assumed to be (Tc - T)\n");
printf("    -nexp n.m      calculate normal exponential curve for Background n.m\n");
printf("                    Formula: y = a * x ^b + n.m\n\n");
}

gen_curve(spc,tim,x,y,coef,numchan,ncol,nrow)
float *spc,*tim,*x,*y,*coef;
int numchan,ncol,nrow;
{
float a,b,xx;
int n,m,i;

/* --------------- finding maximum and minimum x ------------- */
   a = x[0]; b = a;
   for(i = 0; i <= nrow; i++) {
      if(x[i] < a) a = x[i];
      if(x[i] > b) b = x[i];
   }
/* generating curve */
   for(i = 0 ; i <= numchan ; i++) {
      xx = 1.0; spc[i] = 0.0;
      tim[i] = a + (b - a) * ((float)i)/((float)numchan);
      for(n = 0; n <= ncol; n++) {
         spc[i] = spc[i] + (coef[n] * xx);
         xx = xx * tim[i];
      }
   }
}

fit_cexp(x,y,spc,tim,coef,nrow,ncol,tc,numchan)
float *x, *y, *spc, *tim, *coef;
int nrow, ncol;
float tc;
int numchan;
{
int i,n;
float x1,x2,a,b;
   for(i = 0; i <= nrow; i++) {
      x[i] = tim[i]; y[i] = spc[i];
      if(tc != 0.0) x[i] = tc - tim[i];
      if(x[i] <= 0.0) x[i] = 0.00001;
      if(y[i] <= 0.0) y[i] = 0.00001;
      x[i] = (float) log((double) x[i]);
      y[i] = (float) log((double) y[i]);
   }
   polyfit(x, y, y_calc, resid, coef, sig, nrow, 1);
   a = coef[0] ; b = coef[1] ;
   if(b == 0) {
      fprintf(stderr,"sorry, b=0 . exiting.\n");
      exit(0);
   }
   a = (float) exp((double) (a / b));
   x1 = tim[0]; x2 = x1;
   for(i = 0; i <= nrow; i++) {
      if(tim[i] < x1) x1 = tim[i];
      if(tim[i] > x2) x2 = tim[i];
   }
/* generating curve */
   for(i = 0 ; i <= numchan ; i++) {
      tim[i] = x1 + (x2 - x1) * ((float)i)/((float)numchan);
      if(tc != 0.0) {
         spc[i] = a * pow((tc - tim[i]), b);
      } else {
         spc[i] = a * pow(tim[i],b);
      }
   }
}

main(argc,argv)
int argc;
char *argv[];
{
int  i,n,m,numchan;
char s[80],z[80];
float xx,fparam,a,b;

   checkopt(argc,argv,"-dummy",z);

   x      = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   y      = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   y_calc = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   resid  = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   if(resid == NULL) {
      fprintf(stderr,"not enough memory, exiting...\n");
      exit(-1);
   }

   spc    = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   err    = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   tim    = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   if(tim == NULL) {
      fprintf(stderr,"not enough memory, exiting...\n");
      exit(-1);
   }

   ncol = 3; numchan = 2048;
/*  ------------------ read order of polynom ------------------ */
   if(checkopt(argc,argv,"-poly",s)) {
      ncol = atoi(s);
      if(ncol > MAXC) {
         printf("you should not try to use more than %d terms\n",MAXC);
         exit(0);
      }
   }

   if(checkopt(argc,argv,"-cexp",s)) {
      ncol = CEXP;
      fparam = atosf(s);
   }

   if(checkopt(argc,argv,"-nexp",s)) {
      ncol = NEXP;
      fparam = atosf(s);
   }

/*  ------------------ read data points ------------------ */
   if(checkopt(argc,argv,"-in",s)) {
      nrow = readspec(s,spc,err,tim,z);
      nrow = nrow - 1;
      for(i = 0; i <= nrow; i++) {
         x[i] = tim[i]; y[i] = spc[i];
      }
   } else {
      nrow = get_data(tim,spc);
      for(i = 0; i <= nrow; i++) {
         x[i] = tim[i]; y[i] = spc[i];
      }
      strcpy(_spc_onam,"");    /* signal "no output file names */
      writespec("userdata",spc,err,nrow+1,'2',"user data");
      writespec("userdata.tim",tim,err,nrow+1,'2',"user data");
      printf("Your data were stored in userdata.spc and userdata.tim\n");
   }

/* ------------------ going to start fit ...  ------------------ */
/* --- 1. simple polynomial fit --- */
   if(ncol > 0) {
      polyfit(x, y, y_calc, resid, coef, sig, nrow, ncol);
      gen_curve(spc,tim,x,y,coef,numchan,ncol,nrow);
      if(checkopt(argc,argv,"-v",s)) write_data();
   }
/* --- 2. critical exponent fit --- */
   if(ncol == CEXP) {
      ncol = 1;
      fit_cexp(x,y,spc,tim,coef,nrow,ncol,fparam,numchan);
      a = coef[0] ; b = coef[1] ;
      a = (float) exp((double) (a / b));
      if(checkopt(argc,argv,"-v",s)) printf("Beta = %f  Alpha = %f\n",b,a);
   }
/* --- 3. exponential fit --- */
   if(ncol == NEXP) {
printf("NOT IMPLEMENTED YET !!!!!\n");
      for(i = 0; i <= nrow; i++) {
         x[i] = (float) log((double)(tim[i] - fparam));
         y[i] = (float) log((double) y[i]);
      }
      ncol = 1;
      polyfit(x, y, y_calc, resid, coef, sig, nrow, ncol);
      gen_curve(spc,tim,x,y,coef,numchan,ncol,nrow);
      for(i = 0 ; i <= numchan ; i++) spc[i] = (float) exp((double)spc[i]);
   }

/* ------------------ save curve ------------------ */
   if(checkopt(argc,argv,"-o",s)) {
      writespec(s,spc,err,numchan,2,"poly fit curve");
      strcat(s,".tim");
      writespec(s,tim,err,numchan,2,"poly fit curve");
   }
/* ------------------ generate differentiated curve ------------------ */
   if(checkopt(argc,argv,"-dif",s)) {
      for(i = 0 ; i <= (numchan - 1) ; i++) {
         spc[i] = (spc[i] - spc[i+1]) / (tim[i] - tim[i+1]);
      }
      spc[numchan] = spc[numchan-1];
      writespec(s,spc,err,numchan,2,"poly fit differentiated curve");
      strcat(s,".tim");
      writespec(s,tim,err,numchan,2,"poly fit differentiated curve");
   }
   free(x); free(y); free(y_calc); free(resid); free(spc); free(err); free(tim);
   exit(0);
}