#include <string.h>
#include <stdio.h>
#include <tblsys.h>
#include <tblerr.h>
#include <tbldef.h>
#include <macrogen.h>
#include <atype.h>
#include <midas_def.h>
#include <math.h>
/* #include <sunmath.h> */
#include "nrutil.h"

#define TRUE 1
#define FALSE 0
#define MLEN  1024
#define IMAX  300

int   Mlev;
char  Mes[MLEN];


int nonlfit(x,y,sig,ndata,param,npar,listp,mfit,funcs)
float *x, *y, *sig, *param;
int   *listp;
int   ndata, npar, mfit;
void (*funcs)();
{
  int    i, err;
  float  **covar;
  float  **alpha;
  float  ochisq, chisq, alamda;

  covar = matrix(1,npar,1,npar);
  alpha = matrix(1,npar,1,npar);


  alamda = -1.0;   /* initialization */
  err = mrqmin(x,y,sig,ndata,param,npar,listp,mfit,
	       covar,alpha,&chisq,funcs,&alamda);
  if (err < 0) return(err);

  i = 0;
  do
    {
      ochisq = chisq;
      /*
      sprintf(Mes,"Iteration: %d\tChi^2  = %f",i+1,chisq);
      sprintf(&Mes[strlen(Mes)],"\nParam  : %f %f %f %f",
	      param[1],param[2],param[3],param[4]);
      messout(Mlev,4,"nonlfit",Mes);
      SCTPUT(Mes);
      */

      err=mrqmin(x,y,sig,ndata,param,npar,listp,mfit,
		 covar,alpha,&chisq,funcs,&alamda);
      i++;
      /*      printf("%f %f %g %d\n",ochisq,chisq-ochisq,alamda,err); */
      if (err < 0) return(err);
    }while ((alamda > 1.E-6) && (i<IMAX) && (fabs(chisq-ochisq)>1.E-6));

  if (i>1){
   /*
   sprintf(Mes,"Iteration: %d\tChi^2  = %f",i,chisq);
   sprintf(&Mes[strlen(Mes)],"\nParam  : %f %f %f %f",
	   param[1],param[2],param[3],param[4]);
   messout(Mlev,4,"nonlfit",Mes);
   SCTPUT(Mes);
   */
  }

  /* Calculate covariance and curvature matrix. */
  /*
  sprintf(Mes,"Save fit.");
  messout(Mlev,4,"nonlfit",Mes);
  SCTPUT(Mes);
  */
  alamda=0.0;
  err=mrqmin(x,y,sig,ndata,param,npar,listp,mfit,
	     covar,alpha,&chisq,funcs,&alamda);

  return(err);
}


/* the following macro has to be done for FORTRAN codes
   as the linker tries to find nonlft instead of non_lfit!
   This is due to the esoext.exe pre-compiler!!
 */
int nonlft(y,sig,ndata,param,npar,listp,mfit,funcs)
float *y, *sig, *param;
int   *listp;
int   *ndata, *npar, *mfit;
void (*funcs)();
{
  non_lfit(y,sig,*ndata,param,*npar,listp,*mfit,funcs);
}

int non_lfit(y,sig,ndata,param,npar,listp,mfit,funcs)
float *y, *sig, *param;
int   *listp;
int   ndata, npar, mfit;
void (*funcs)();
{
  int    i, err;
  float  **covar;
  float  **alpha;
  float  ochisq, chisq, alamda;

  covar = matrix(1,npar,1,npar);
  alpha = matrix(1,npar,1,npar);
  /*
  printf("non_lfit: %f %f %f %f, %f %f, %d %d\n",
	 param[1],param[2],param[3],param[4],y[1],sig[1],ndata,npar);
  printf("non_lfit: %d %d %d %d\n",
	 listp[1],listp[2],listp[3],listp[4]);
	 */
  alamda = -1.0;   /* initialization */
  err = mrq_min(y,sig,ndata,param,npar,listp,mfit,
	       covar,alpha,&chisq,funcs,&alamda);
  if (err < 0) return(err);
  /*
  printf("result: %f %f %f %f, %f\n",
	 param[1],param[2],param[3],param[4],chisq);
  exit(1);
  */

  i = 0;
  do
    {
      ochisq = chisq;
      /*
      sprintf(Mes,"Iteration: %d\tChi^2  = %f",i+1,chisq);
      sprintf(&Mes[strlen(Mes)],"\nParam  : %f %f %f %f",
	      param[1],param[2],param[3],param[4]);
      messout(Mlev,4,"nonlfit",Mes);
      SCTPUT(Mes);
      */

      err=mrq_min(y,sig,ndata,param,npar,listp,mfit,
		 covar,alpha,&chisq,funcs,&alamda);
      i++;
      /*      printf("%f %f %f %g %d\n",chisq,ochisq,chisq-ochisq,alamda,i); */
      if (err < 0) return(err);
    }while ((alamda > 1.E-6) && (i<IMAX));
  /*
    }while ((alamda > 1.E-6) && (alamda < 1.E38) && (i<IMAX));
    }while ((alamda > 1.E-6) && (i<IMAX) && (fabs(chisq-ochisq)>1.E-8));
    */  
  /*
  if (i>1){
   sprintf(Mes,"Iteration: %d\talamda = %f \tChi^2  = %f",i,alamda,chisq);
   sprintf(&Mes[strlen(Mes)],"\nParam  : %f %f %f %f",
	   param[1],param[2],param[3],param[4]);
   messout(Mlev,4,"nonlfit",Mes);
   SCTPUT(Mes);
  }
  */

  /* Calculate covariance and curvature matrix. */
  /*
  sprintf(Mes,"Save fit.");
  messout(Mlev,4,"nonlfit",Mes);
  SCTPUT(Mes);
  */
  alamda=0.0;
  err=mrq_min(y,sig,ndata,param,npar,listp,mfit,
	     covar,alpha,&chisq,funcs,&alamda);

  return(err);
}


int calculfit(x, y, fit, res, funcs, param, num, npar)
float *x, *y, *fit, *res, *param;
int   num, npar;
void (*funcs)();
{
  int i;
  float *dy;

  dy = vector(1,npar);
  for (i=1;i<=num;i++){
    (*funcs)(x[i], param, &fit[i], dy);
    /*    printf("%d %f %f %f\n",i,x[i],y[i],fit[i]); */
    res[i] = y[i] - fit[i];
  }
}



int calcul_fit(y, fit, res, funcs, param, num, npar)
float *y, *fit, *res, *param;
int   num, npar;
void (*funcs)();
{
  int i;
  float x, *dy;

  dy = vector(1,npar);
  for (i=1;i<=num;i++){
    x = i;
    (*funcs)(x, param, &fit[i], dy);
    /*    printf("%d %f %f %f\n",i,x[i],y[i],fit[i]); */
    res[i] = y[i] - fit[i];
  }
}


int show_param(param,listp,num,mfit)
float *param;
int   *listp, num, mfit;
{
  int i;
  
  sprintf(Mes,"Initial parameter set:\n\t\t\tVALUE\t\tUSED-FLAG");
  for (i=1;i<=num;i++){
    sprintf(&Mes[strlen(Mes)],"\nparameter_%d:\t%f\t%d",i,param[i],listp[i]);
  }
    
  sprintf(&Mes[strlen(Mes)],"\n(USED-FLAG: 0 indicates a fixed parameter)");
  /*  messout(Mlev,4,"nonlfit",Mes); */
  SCTPUT(Mes);
}