/******************************************************************************
** Copyright 2015 MBARI - Monterey Bay Aquarium Research Institute
*******************************************************************************
*******************************************************************************
** Summary  : Computes a least squares fit to a given linear or non-linear modelled absorbance spectra
** Filename : lsqft.c
** Author   : Luke Coletti
** Project  :
** Version  : 1.0
** Compiler : MetroWerks Code Warrior v8.3
** Created  : 12/15/14
** Archived :
*******************************************************************************
** Modification History:
** 01/02/15 : CF2 Port, LJC
*******************************************************************************/

#include        <cfxpico.h>             // Persistor PicoDOS Definitions

#include 		<stdio.h>
#include 		<math.h>
#include 		"lsqfit.h"
#include 		"nrutil.h"

#define NRANSI
#define TRUE  1
#define FALSE 0
#define ITMAX 100
#define EPS 3.0e-7
#define FPMIN 1.0e-30
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}


/******************************************************************************\
**	linear absorbance model
\******************************************************************************/

int LinFit(double **CoefVals, double *LambdaVals, double *AbsVals, double *AbsFitVals, double *ResVals, double *FitVals, ushort fit_pixels, ushort fit_concs, ushort bsl_model)
{


int    fitsal;
int    ret, i, MA, *ia;
static double *a;
double *sig, **covar;
double chisq, rr, sumrr, RMSerror;
double NO3=5.0, SAL=34.5, HS=1.0;

        //printf("\nLinFit()%lf,%lf,%.9lf,%.9lf,%.9lf\n", LambdaVals[1], AbsVals[1], CoefVals[1][0], CoefVals[1][1], CoefVals[1][2]);
        //printf("\nLinFit() FitData[0] = %lf\n", FitVals[SALINITY]);

        if( (bsl_model != 1) && (bsl_model != 2) )
          return(-1);

        if( (fit_concs != 2) && (fit_concs != 3) )
          return(-2);

        /* set number of fit parameters to minimize */
        if(bsl_model == 1)
          MA = fit_concs+2;
        else if(bsl_model == 2)
          MA = fit_concs+3;

        ia=ivector(1,MA);
        a=dvector(1,MA);
        sig=dvector(1,fit_pixels);
        covar=dmatrix(1,MA,1,MA);

        for(i = 1; i <= fit_pixels; i++) {
          sig[i]=0.001;
          }

        /* set which parameter to fit and */
        /* prime the fitted parameters to something reasonable */
        if(FitVals[SALINITY] != 0.0){
          SAL=FitVals[SALINITY];
          fitsal=FALSE;
          }
        else
          fitsal=TRUE;

        /* prime the pump */
        if(fit_concs == 2){
          if(bsl_model == 1){ // double conc w/2 term baseline fit
            ia[1]=fitsal; ia[2]=1; ia[3]=1; ia[4]=1;
            a[1]=SAL; a[2]=NO3; a[3]=0.0001; a[4]=0.05;
            }
          else if(bsl_model == 2){ // double conc w/3 term poly() baseline fit
            ia[1]=fitsal; ia[2]=1; ia[3]=1; ia[4]=1; ia[5]=1;
            a[1]=SAL; a[2]=NO3; a[3]=0.000001; a[4]=0.0001; a[5]=0.05;
            }
          }

        else if(fit_concs == 3){
          if(bsl_model == 1){ // triple conc w/2 term baseline fit
            ia[1]=fitsal; ia[2]=1; ia[3]=1; ia[4]=1; ia[5]=1;
            a[1]=SAL; a[2]=NO3; a[3]=HS; a[4]=0.0001; a[5]=0.05;
            }
          else if(bsl_model == 2){ // triple conc w/3 term poly() baseline fit
            ia[1]=fitsal; ia[2]=1; ia[3]=1; ia[4]=1; ia[5]=1; ia[6]=1;
            a[1]=SAL; a[2]=NO3; a[3]=HS; a[4]=0.000001; a[5]=0.0001; a[6]=0.05;
            }
          }

        /* run the fitting algorithim */
        ret = _lfit(LambdaVals,AbsVals,AbsFitVals,sig,CoefVals,a,ia,MA,covar,&chisq,fit_pixels,fit_concs,bsl_model,LinModelFuncs);

        /* was the fitting sucessful? */
        if( (ret == FALSE) || ((a[1] == SAL) && (a[2] == NO3)) ){
          if(ret == FALSE){
            for(i = 0; i < (MA+3); i++){
              FitVals[i] = -1.0;
              }
            }
          else{
            for(i = 0; i < (MA+3); i++){
              FitVals[i] = -2.0;
              }
            }
          free_dmatrix(covar,1,MA,1,MA);
          free_dvector(sig,1,fit_pixels);
          free_dvector(a,1,MA);
          free_ivector(ia,1,MA);
          return(FALSE);
          }

        /* compute the residuals */
        sumrr = 0;
        for(i=1; i<=fit_pixels; i++){
          ResVals[i] = AbsFitVals[i]-AbsVals[i];
          rr   = pow(ResVals[i], 2.0);
          sumrr= rr + sumrr;
          }
        RMSerror = sqrt(sumrr/fit_pixels);

        /* store the results */
        if(fit_concs == 2){
          if(bsl_model == 1){ // double conc w/2 term baseline fit
            FitVals[SALINITY]   = a[1];
            FitVals[NITRATE]    = a[2];
            FitVals[BISULFIDE]  = -1.0;
            FitVals[BASELINE_A] = a[3];
            FitVals[BASELINE_B] = a[4];
            FitVals[BASELINE_C] = -1.0;
            FitVals[FIT_ERROR]  = RMSerror;
            FitVals[FIT_CHISQ]  = chisq;
            FitVals[FIT_GAMMQ]  = _gammq(0.5*(fit_pixels-MA), 0.5*chisq);
            }
          else if(bsl_model == 2){ // double conc w/3 term poly() baseline fit
            FitVals[SALINITY]   = a[1];
            FitVals[NITRATE]    = a[2];
            FitVals[BISULFIDE]  = -1.0;
            FitVals[BASELINE_A] = a[3];
            FitVals[BASELINE_B] = a[4];
            FitVals[BASELINE_C] = a[5];
            FitVals[FIT_ERROR]  = RMSerror;
            FitVals[FIT_CHISQ]  = chisq;
            FitVals[FIT_GAMMQ]  = _gammq(0.5*(fit_pixels-MA), 0.5*chisq);
            }
          }

        else if(fit_concs == 3){
          if(bsl_model == 1){ // triple conc w/2 term baseline fit
            FitVals[SALINITY]   = a[1];
            FitVals[NITRATE]    = a[2];
            FitVals[BISULFIDE]  = a[3];
            FitVals[BASELINE_A] = a[4];
            FitVals[BASELINE_B] = a[5];
            FitVals[BASELINE_C] = -1.0;
            FitVals[FIT_ERROR]  = RMSerror;
            FitVals[FIT_CHISQ]  = chisq;
            FitVals[FIT_GAMMQ]  = _gammq(0.5*(fit_pixels-MA), 0.5*chisq);
            }
          else if(bsl_model == 2){ // triple conc w/3 term poly() baseline fit
            FitVals[SALINITY]   = a[1];
            FitVals[NITRATE]    = a[2];
            FitVals[BISULFIDE]  = a[3];
            FitVals[BASELINE_A] = a[4];
            FitVals[BASELINE_B] = a[5];
            FitVals[BASELINE_C] = a[6];
            FitVals[FIT_ERROR]  = RMSerror;
            FitVals[FIT_CHISQ]  = chisq;
            FitVals[FIT_GAMMQ]  = _gammq(0.5*(fit_pixels-MA), 0.5*chisq);
            }
          }

        free_dmatrix(covar,1,MA,1,MA);
        free_dvector(sig,1,fit_pixels);
        free_dvector(a,1,MA);
        free_ivector(ia,1,MA);

        return(TRUE);
}

int    _lfit(double x[], double y[], double yfit[], double sig[], double **Ecoefs, double a[],
        int ia[], int ma, double **covar, double *chisq, ushort fitpixels, ushort fitconcs, ushort bslmodel,
        void (*funcs)(double [], double [], double **, int, ushort, ushort) )
{

        int ret;
        int i,j,k,l,m,mfit=0;
        double ym,wt,sum,sig2i,**beta,*afunc;

        beta=dmatrix(1,ma,1,1);
        afunc=dvector(1,ma);
        for (j=1;j<=ma;j++)
                if (ia[j]) mfit++;
        if(mfit == 0){
          printf("lfit: no parameters to be fitted");
          free_dvector(afunc,1,ma);
          free_dmatrix(beta,1,ma,1,1);
          return(FALSE);
          }
        for (j=1;j<=mfit;j++) {
                for (k=1;k<=mfit;k++) covar[j][k]=0.0;
                beta[j][1]=0.0;
        }
        for (i=1;i<=fitpixels;i++) {
                (*funcs)(x,afunc,Ecoefs,i,fitconcs,bslmodel);
                ym=y[i];
                if (mfit < ma) {
                        for (j=1;j<=ma;j++)
                                if (!ia[j]) ym -= a[j]*afunc[j];
                }
                sig2i=1.0/DSQR(sig[i]);
                for (j=0,l=1;l<=ma;l++) {
                        if (ia[l]) {
                                wt=afunc[l]*sig2i;
                                for (j++,k=0,m=1;m<=l;m++)
                                        if (ia[m]) covar[j][++k] += wt*afunc[m];
                                beta[j][1] += ym*wt;
                        }
                }
        }
        for (j=2;j<=mfit;j++)
                for (k=1;k<j;k++)
                        covar[k][j]=covar[j][k];

        ret = _gaussj(covar,mfit,beta,1);
        if(ret == FALSE){   // test for singular matrix
          free_dvector(afunc,1,ma);
          free_dmatrix(beta,1,ma,1,1);
          return(FALSE);
          }

        for (j=0,l=1;l<=ma;l++)
                if (ia[l]) a[l]=beta[++j][1];
        *chisq=0.0;
        for (i=1;i<=fitpixels;i++) {
                (*funcs)(x,afunc,Ecoefs,i,fitconcs,bslmodel);
                for (sum=0.0,j=1;j<=ma;j++) sum += a[j]*afunc[j];
                yfit[i] = sum;
                *chisq += DSQR((y[i]-sum)/sig[i]);
        }
        _covsrt(covar,ma,ia,mfit);
        free_dvector(afunc,1,ma);
        free_dmatrix(beta,1,ma,1,1);
        return(TRUE);
}

#undef NRANSI


/******************************************************************************/
/** Summary  : Using General Least Squares, as published in NR in C 15.4,     */
/**            compute the basis functions, afunc (na of), for a selected     */
/**            model each of which varies in a linear manner to its parameter */
/**            (a) at x.                                                      */
/******************************************************************************/
void LinModelFuncs(double x[], double afunc[], double **Ecoefs, int i, ushort fitconcs, ushort bslmodel)
{
/*
   x          = lambda
   afunc[]    = basis functions
   Ecoefs[][] = molar extinction coefs
   i          = Ecoef index
   fitconcs   = number of concentrations in model
   bslmodel   = type of baseline in model (which defines # of parameter to fit)
*/

        if(fitconcs == 2){
          if(bslmodel == 1){ // double conc w/2 term linear baseline fit
            afunc[1] = Ecoefs[i][0];
            afunc[2] = Ecoefs[i][1];
            afunc[3] = x[i];
            afunc[4] = 1.0;
            //if(i == 1)
            //printf("%lf,%lf,%.9lf,%.9lf\n", x, afunc[1], Ecoefs[i][0], Ecoefs[i][1]);
            }
          else if(bslmodel == 2){ // double conc w/3 term quadratic baseline fit
            afunc[1] = Ecoefs[i][0];
            afunc[2] = Ecoefs[i][1];
            afunc[3] = x[i]*x[i];
            afunc[4] = x[i];
            afunc[5] = 1.0;
            //printf("%lf,%lf,%.9lf,%.9lf\n", x[i], afunc[1], Ecoefs[i][0], Ecoefs[i][1]);
            }
          }

        else if(fitconcs == 3){
          if(bslmodel == 1){ // triple conc w/2 term linear baseline fit
            afunc[1] = Ecoefs[i][0];
            afunc[2] = Ecoefs[i][1];
            afunc[3] = Ecoefs[i][2];
            afunc[4] = x[i];
            afunc[5] = 1.0;
            }
          else if(bslmodel == 2){ // triple conc w/3 term quadratic baseline fit
            afunc[1] = Ecoefs[i][0];
            afunc[2] = Ecoefs[i][1];
            afunc[3] = Ecoefs[i][2];
            afunc[4] = x[i]*x[i];
            afunc[5] = x[i];
            afunc[6] = 1.0;
            }
          }

}


/******************************************************************************\
**	non-linear absorbance model
\******************************************************************************/

#define NRANSI

int NLinFit(double **CoefVals, double *LambdaVals, double *AbsVals, double *AbsFitVals, double *ResVals, double *FitVals, ushort fit_pixels, ushort fit_concs, ushort bsl_model)
{

int    fitsal;
int    ret, i, k, MA, *ia;
static double *a;
double *sig, *yfitsav, **covar, **alpha;
double alamda, chisq, ochisq, rr, sumrr, RMSerror;
double NO3=5.0, SAL=34.5, HS=1.0;


        //printf("\nNLinFit()%lf,%lf,%.9lf,%.9lf,%.9lf\n", LambdaVals[1], AbsVals[1], CoefVals[1][0], CoefVals[1][1], CoefVals[1][2]);
        //printf("\nNLinFit() FitData[0] = %lf\n", FitVals[SALINITY]);

        if( (bsl_model != 3) && (bsl_model != 4) )
          return(-1);

        if( (fit_concs != 2) && (fit_concs != 3) )
          return(-2);

        /* set number of fit parameters to minimize */
        if(bsl_model == 3)
          MA = fit_concs+2;
        else if(bsl_model == 4)
          MA = fit_concs+3;

        ia=ivector(1,MA);
        a=dvector(1,MA);
        yfitsav=dvector(1,fit_pixels);
        sig=dvector(1,fit_pixels);
        covar=dmatrix(1,MA,1,MA);
        alpha=dmatrix(1,MA,1,MA);

        for(i = 1; i <= fit_pixels; i++) {
          sig[i]=0.001;
          }

        /* set which parameter to fit and */
        /* prime the fitted parameters to something reasonable */
        if(FitVals[SALINITY] != 0.0){
          SAL=FitVals[SALINITY];
          fitsal=FALSE;
          }
        else
          fitsal=TRUE;

        if(fit_concs == 2){
          if(bsl_model == 3){ // double conc w/2 term baseline fit
            ia[1]=fitsal; ia[2]=1; ia[3]=1; ia[4]=1;
            a[1]=SAL; a[2]=NO3; a[3]=-0.0025; a[4]=0.05;
            }
          else if(bsl_model == 4){ // double conc w/3 term poly() baseline fit
            ia[1]=fitsal; ia[2]=1; ia[3]=1; ia[4]=1; ia[5]=1;
            a[1]=SAL; a[2]=NO3; a[3]=0.000001; a[4]=-0.0025; a[5]=0.05;
            }
          }

        else if(fit_concs == 3){
          if(bsl_model == 3){ // triple conc w/2 term baseline fit
            ia[1]=fitsal; ia[2]=1; ia[3]=1; ia[4]=1; ia[5]=1;
            a[1]=SAL; a[2]=NO3; a[3]=HS; a[4]=-0.0025; a[5]=0.05;
            }
          else if(bsl_model == 4){ // triple conc w/3 term poly() baseline fit
            ia[1]=fitsal; ia[2]=1; ia[3]=1; ia[4]=1; ia[5]=1; ia[6]=1;
            a[1]=SAL; a[2]=NO3; a[3]=HS; a[4]=0.000001; a[5]=-0.0025; a[6]=0.05;
            }
          }

        /* initialize the fitting algorithim */
        alamda = -1.0;
        _mrqmin(LambdaVals,AbsVals,AbsFitVals,sig,CoefVals,a,ia,MA,covar,&chisq,&alamda,alpha,fit_pixels,fit_concs,bsl_model,NLinModelFuncs);

        /* run the fitting algorithim */
        k=1;
        for(;;) {
          k++;
          ochisq=chisq;
          for (i = 1; i <= fit_pixels; i++){
            yfitsav[i]=AbsFitVals[i]; /* save the fitted model function results */
            }
          ret=_mrqmin(LambdaVals,AbsVals,AbsFitVals,sig,CoefVals,a,ia,MA,covar,&chisq,&alamda,alpha,fit_pixels,fit_concs,bsl_model,NLinModelFuncs);
          if(ret == FALSE)
            break;
          if(chisq < ochisq)
            continue;
          /* generate stats from the fitting algorithim */
          alamda=0.0;
          _mrqmin(LambdaVals,AbsVals,AbsFitVals,sig,CoefVals,a,ia,MA,covar,&chisq,&alamda,alpha,fit_pixels,fit_concs,bsl_model,NLinModelFuncs);
          break;
          }

        /* was the fitting sucessful? */
        if( (ret == FALSE) || ((a[1] == SAL) && (a[2] == NO3)) ){
          if(ret == FALSE){
            for(i = 0; i < (MA+3); i++){
              FitVals[i] = -1.0;
              }
            }
          else{
            for(i = 0; i < (MA+3); i++){
              FitVals[i] = -2.0;
              }
            }
          free_dmatrix(covar,1,MA,1,MA);
          free_dvector(sig,1,fit_pixels);
          free_dvector(a,1,MA);
          free_ivector(ia,1,MA);
          return(FALSE);
          }

        /* compute the residuals */
        sumrr = 0;
        for(i=1; i<=fit_pixels; i++){
          ResVals[i] = AbsFitVals[i]-AbsVals[i];
          rr   = pow(ResVals[i], 2.0);
          sumrr= rr + sumrr;
          }
        RMSerror = sqrt(sumrr/fit_pixels);

        /* store the results */
        if(fit_concs == 2){
          if(bsl_model == 1){ // double conc w/2 term baseline fit
            FitVals[SALINITY]   = a[1];
            FitVals[NITRATE]    = a[2];
            FitVals[BISULFIDE]  = -1.0;
            FitVals[BASELINE_A] = a[3];
            FitVals[BASELINE_B] = a[4];
            FitVals[BASELINE_C] = -1.0;
            FitVals[FIT_ERROR]  = RMSerror;
            FitVals[FIT_CHISQ]  = chisq;
            FitVals[FIT_GAMMQ]  = _gammq(0.5*(fit_pixels-MA), 0.5*chisq);
            }
          else if(bsl_model == 2){ // double conc w/3 term poly() baseline fit
            FitVals[SALINITY]   = a[1];
            FitVals[NITRATE]    = a[2];
            FitVals[BISULFIDE]  = -1.0;
            FitVals[BASELINE_A] = a[3];
            FitVals[BASELINE_B] = a[4];
            FitVals[BASELINE_C] = a[5];
            FitVals[FIT_ERROR]  = RMSerror;
            FitVals[FIT_CHISQ]  = chisq;
            FitVals[FIT_GAMMQ]  = _gammq(0.5*(fit_pixels-MA), 0.5*chisq);
            }
          }

        else if(fit_concs == 3){
          if(bsl_model == 1){ // triple conc w/2 term baseline fit
            FitVals[SALINITY]   = a[1];
            FitVals[NITRATE]    = a[2];
            FitVals[BISULFIDE]  = a[3];
            FitVals[BASELINE_A] = a[4];
            FitVals[BASELINE_B] = a[5];
            FitVals[BASELINE_C] = -1.0;
            FitVals[FIT_ERROR]  = RMSerror;
            FitVals[FIT_CHISQ]  = chisq;
            FitVals[FIT_GAMMQ]  = _gammq(0.5*(fit_pixels-MA), 0.5*chisq);
            }
          else if(bsl_model == 2){ // triple conc w/3 term poly() baseline fit
            FitVals[SALINITY]   = a[1];
            FitVals[NITRATE]    = a[2];
            FitVals[BISULFIDE]  = a[3];
            FitVals[BASELINE_A] = a[4];
            FitVals[BASELINE_B] = a[5];
            FitVals[BASELINE_C] = a[6];
            FitVals[FIT_ERROR]  = RMSerror;
            FitVals[FIT_CHISQ]  = chisq;
            FitVals[FIT_GAMMQ]  = _gammq(0.5*(fit_pixels-MA), 0.5*chisq);
            }
          }

        free_dmatrix(covar,1,MA,1,MA);
        free_dvector(sig,1,fit_pixels);
        free_dvector(a,1,MA);
        free_ivector(ia,1,MA);

        //return(TRUE);
        return(k);
}


int    _mrqmin(double x[], double y[], double yfit[], double sig[], double **Ecoefs, double a[],
        int ia[], int ma, double **covar, double *chisq, double *alamda, double **alpha, ushort fitpixels, ushort fitconcs, ushort bslmodel,
        void (*funcs)(double [], double [], double [], double [], double **, int, ushort, ushort) )
{

        int ret;
        int j,k,l;
        static int mfit;
        static double ochisq,*atry,*beta,*da,**oneda;

        if (*alamda < 0.0) {  //initialization time
                atry=dvector(1,ma);
                beta=dvector(1,ma);
                da=dvector(1,ma);
                for (mfit=0,j=1;j<=ma;j++)
                        if (ia[j]) mfit++;
                oneda=dmatrix(1,mfit,1,1);
                *alamda=0.001;
                _mrqcof(x,y,yfit,sig,Ecoefs,a,ia,ma,alpha,beta,chisq,fitpixels,fitconcs,bslmodel,funcs);
                ochisq=(*chisq);
                for (j=1;j<=ma;j++) atry[j]=a[j];
        }

        for (j=1;j<=mfit;j++) {
                for (k=1;k<=mfit;k++) covar[j][k]=alpha[j][k];
                covar[j][j]=alpha[j][j]*(1.0+(*alamda));
                oneda[j][1]=beta[j];
        }
        ret = _gaussj(covar,mfit,oneda,1);
        if(ret == FALSE){   // test for singular matrix
          free_dmatrix(oneda,1,mfit,1,1);
          free_dvector(da,1,ma);
          free_dvector(beta,1,ma);
          free_dvector(atry,1,ma);
          return(FALSE);
        }

        for (j=1;j<=mfit;j++) da[j]=oneda[j][1];

        if (*alamda == 0.0) {  //process stats time
                _covsrt(covar,ma,ia,mfit);
                _covsrt(alpha,ma,ia,mfit);
                free_dmatrix(oneda,1,mfit,1,1);
                free_dvector(da,1,ma);
                free_dvector(beta,1,ma);
                free_dvector(atry,1,ma);
                return(TRUE);
        }
        for (j=0,l=1;l<=ma;l++)
                if (ia[l]) atry[l]=a[l]+da[++j];
        _mrqcof(x,y,yfit,sig,Ecoefs,atry,ia,ma,covar,da,chisq,fitpixels,fitconcs,bslmodel,funcs);
        if (*chisq < ochisq) {
                *alamda *= 0.1;
                ochisq=(*chisq);
                for (j=1;j<=mfit;j++) {
                        for (k=1;k<=mfit;k++) alpha[j][k]=covar[j][k];
                        beta[j]=da[j];
                }
                for (l=1;l<=ma;l++) a[l]=atry[l];
        } else {
                *alamda *= 10.0;
                *chisq=ochisq;
        }

        return(TRUE);

}

void _mrqcof(double x[], double y[], double yfit[], double sig[], double **Ecoefs, double a[],
        int ia[], int ma, double **alpha, double beta[], double *chisq, ushort fitpixels, ushort fitconcs, ushort bslmodel,
        void (*funcs)(double [], double [], double [], double [], double **, int, ushort, ushort) )
{
        int i,j,k,l,m,mfit=0;
        double wt,sig2i,dy,*dyda;

        dyda=dvector(1,ma);
        for (j=1;j<=ma;j++)
                if (ia[j]) mfit++;
        for (j=1;j<=mfit;j++) {
                for (k=1;k<=j;k++) alpha[j][k]=0.0;
                beta[j]=0.0;
        }
        *chisq=0.0;
        //printf("COF %lf,%lf,%lf,%lf\n", a[1],a[2],a[3],a[4]);
        for (i=1;i<=fitpixels;i++) {
                (*funcs)(x,a,yfit,dyda,Ecoefs,i,fitconcs,bslmodel);
                //printf("x = %lf y = %lf yfit = %lf\n", x[i], y[i], yfit[i]);
                sig2i=1.0/(sig[i]*sig[i]);
                dy=y[i]-yfit[i];
                for (j=0,l=1;l<=ma;l++) {
                        if (ia[l]) {
                                wt=dyda[l]*sig2i;
                                for (j++,k=0,m=1;m<=l;m++)
                                        if (ia[m]) alpha[j][++k] += wt*dyda[m];
                                beta[j] += dy*wt;
                        }
                }
                *chisq += dy*dy*sig2i;
                //printf("%lf\n", *chisq);
        }
        for (j=2;j<=mfit;j++)
                for (k=1;k<j;k++) alpha[k][j]=alpha[j][k];
        free_dvector(dyda,1,ma);
}

#undef NRANSI


/******************************************************************************/
/** Summary  : Using General Least Squares, as published in NR in C 15.5,     */
/**            compute the basis functions, afunc (na of), for a selected     */
/**            model each of which varies in a linear manner to its parameter */
/**            (a) at x.                                                      */
/******************************************************************************/

void NLinModelFuncs(double x[], double a[], double yfit[], double dyda[], double **Ecoefs, int i, ushort fitconcs, ushort bslmodel)
{
/*
   x[]        = lambda
   afunc[]    = parameters to be minimized
   yfit[]     = fitted absorbance, using current parameter vals
   dyda[]     = partial derivitive of the given parameter val a[n]
   Ecoefs[][] = molar extinction coefs
   i          = index to x, yfit and Ecoefs
   fitconcs   = number of concentrations in model
   bslmodel   = type of modeled baseline
*/

double temp;

        if(fitconcs == 2){
          if(bslmodel == 3){ // double conc w/2 term exp-linear baseline
            temp    = a[3]*x[i] + a[4];
            yfit[i] = Ecoefs[i][0]*a[1] + Ecoefs[i][1]*a[2] + exp(temp);
          //yfit[i] = Ecoefs[i][0]*a[1] + Ecoefs[i][1]*a[2] + temp;
            dyda[1] = Ecoefs[i][0];
            dyda[2] = Ecoefs[i][1];
            dyda[3] = x[i]*exp(temp);
            dyda[4] = exp(temp);
          //dyda[3] = x[i];
          //dyda[4] = 1.0;
          //if(i == 1)
            //printf("%lf,%lf,%lf,%.9lf,%.9lf\n", x[i], yfit[i], dyda[1], Ecoefs[i][0], Ecoefs[i][1]);
            }
          else if(bslmodel == 4){ // double conc w/3 term exp-quadratic baseline
            temp    = a[3]*x[i]*x[i] + a[4]*x[i] + a[5];
            yfit[i] = Ecoefs[i][0]*a[1] + Ecoefs[i][1]*a[2] + exp(temp);
          //yfit[i] = Ecoefs[i][0]*a[1] + Ecoefs[i][1]*a[2] + temp;
            dyda[1] = Ecoefs[i][0];
            dyda[2] = Ecoefs[i][1];
            dyda[3] = x[i]*x[i]*exp(temp);
            dyda[4] = x[i]*exp(temp);
            dyda[5] = exp(temp);
          //dyda[3] = x[i]*x[i];
          //dyda[4] = x[i];
          //dyda[5] = 1.0;
            }
          }

        else if(fitconcs == 3){
          if(bslmodel == 3){ // triple conc w/2 term exp-linear baseline
            temp = a[4]*x[i] + a[5];
            yfit[i] = Ecoefs[i][0]*a[1] + Ecoefs[i][1]*a[2] + Ecoefs[i][2]*a[3] + exp(temp);
            dyda[1] = Ecoefs[i][0];
            dyda[2] = Ecoefs[i][1];
            dyda[3] = Ecoefs[i][2];
            dyda[4] = x[i]*exp(temp);
            dyda[5] = exp(temp);
            }
          else if(bslmodel == 4){ // triple conc w/3 term exp-quadratic baseline
            temp    = a[4]*x[i]*x[i] + a[5]*x[i] + a[6];
            yfit[i] = Ecoefs[i][0]*a[1] + Ecoefs[i][1]*a[2] + Ecoefs[i][2]*a[3] + exp(temp);
            dyda[1] = Ecoefs[i][0];
            dyda[2] = Ecoefs[i][1];
            dyda[3] = Ecoefs[i][2];
            dyda[4] = x[i]*x[i]*exp(temp);
            dyda[5] = x[i]*exp(temp);
            dyda[6] = exp(temp);
            }
          }

}

/******************************************************************************\
**	NR in C math functons.
\******************************************************************************/

void _covsrt(double **covar, int ma, int ia[], int mfit)
{
        int i,j,k;
        double temp;

        for (i=mfit+1;i<=ma;i++)
                for (j=1;j<=i;j++) covar[i][j]=covar[j][i]=0.0;
        k=mfit;
        for (j=ma;j>=1;j--) {
                if (ia[j]) {
                        for (i=1;i<=ma;i++) SWAP(covar[i][k],covar[i][j])
                        for (i=1;i<=ma;i++) SWAP(covar[k][i],covar[j][i])
                        k--;
                }
        }
}


void _gcf(double *gammcf, double a, double x, double *gln)
{
        int i;
        double an,b,c,d,del,h;

        *gln=_gammln(a);
        b=x+1.0-a;
        c=1.0/FPMIN;
        d=1.0/b;
        h=d;
        for (i=1;i<=ITMAX;i++) {
                an = -i*(i-a);
                b += 2.0;
                d=an*d+b;
                if (fabs(d) < FPMIN) d=FPMIN;
                c=b+an/c;
                if (fabs(c) < FPMIN) c=FPMIN;
                d=1.0/d;
                del=d*c;
                h *= del;
                if (fabs(del-1.0) < EPS) break;
        }
//        if (i > ITMAX) nrerror("a too large, ITMAX too small in gcf");
        *gammcf=exp(-x+a*log(x)-(*gln))*h;
}


void _gser(double *gamser, double a, double x, double *gln)
{
        int n;
        double sum,del,ap;

        *gln=_gammln(a);
        if (x <= 0.0) {
                if (x < 0.0) nrerror("x less than 0 in routine gser");
                *gamser=0.0;
                return;
        } else {
                ap=a;
                del=sum=1.0/a;
                for (n=1;n<=ITMAX;n++) {
                        ++ap;
                        del *= x/ap;
                        sum += del;
                        if (fabs(del) < fabs(sum)*EPS) {
                                *gamser=sum*exp(-x+a*log(x)-(*gln));
                                return;
                        }
                }
//                nrerror("a too large, ITMAX too small in routine gser");
                return;
        }
}


int _gaussj(double **a, int n, double **b, int m)
{
        int *indxc,*indxr,*ipiv;
        int i,icol,irow,j,k,l,ll;
        double big,dum,pivinv,temp;

        indxc=ivector(1,n);
        indxr=ivector(1,n);
        ipiv=ivector(1,n);
        for (j=1;j<=n;j++) ipiv[j]=0;
        for (i=1;i<=n;i++) {
                big=0.0;
                for (j=1;j<=n;j++)
                        if (ipiv[j] != 1)
                                for (k=1;k<=n;k++) {
                                        if (ipiv[k] == 0) {
                                                if (fabs(a[j][k]) >= big) {
                                                        big=fabs(a[j][k]);
                                                        irow=j;
                                                        icol=k;
                                                }
                                        } else if (ipiv[k] > 1){ //nrerror("gaussj: Singular Matrix-1");
                                            free_ivector(ipiv,1,n);
                                            free_ivector(indxr,1,n);
                                            free_ivector(indxc,1,n);
                                            return(FALSE);
                                            }
                                }
                ++(ipiv[icol]);
                if (irow != icol) {
                        for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l])
                        for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
                }
                indxr[i]=irow;
                indxc[i]=icol;
                if (a[icol][icol] == 0.0){ //nrerror("gaussj: Singular Matrix-2");
                  free_ivector(ipiv,1,n);
                  free_ivector(indxr,1,n);
                  free_ivector(indxc,1,n);
                  return(FALSE);
                  }

                pivinv=1.0/a[icol][icol];
                a[icol][icol]=1.0;
                for (l=1;l<=n;l++) a[icol][l] *= pivinv;
                for (l=1;l<=m;l++) b[icol][l] *= pivinv;
                for (ll=1;ll<=n;ll++)
                        if (ll != icol) {
                                dum=a[ll][icol];
                                a[ll][icol]=0.0;
                                for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum;
                                for (l=1;l<=m;l++) b[ll][l] -= b[icol][l]*dum;
                        }
        }
        for (l=n;l>=1;l--) {
                if (indxr[l] != indxc[l])
                        for (k=1;k<=n;k++)
                                SWAP(a[k][indxr[l]],a[k][indxc[l]]);
        }
        free_ivector(ipiv,1,n);
        free_ivector(indxr,1,n);
        free_ivector(indxc,1,n);
        return(TRUE);

}


double _gammp(double a, double x)
{

        double gamser,gammcf,gln;

//        if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammp");
        if (x < (a+1.0)) {
                _gser(&gamser,a,x,&gln);
                return gamser;
        } else {
                _gcf(&gammcf,a,x,&gln);
                return 1.0-gammcf;
        }
}


double _gammq(double a, double x)
{

        double gamser,gammcf,gln;

//        if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammq");
        if (x < (a+1.0)) {
                _gser(&gamser,a,x,&gln);
                return 1.0-gamser;
        } else {
                _gcf(&gammcf,a,x,&gln);
                return gammcf;
        }
}


double _gammln(double xx)
{
        double x,y,tmp,ser;
        static double cof[6]={76.18009172947146,-86.50532032941677,
                24.01409824083091,-1.231739572450155,
                0.1208650973866179e-2,-0.5395239384953e-5};
        int j;

        y=x=xx;
        tmp=x+5.5;
        tmp -= (x+0.5)*log(tmp);
        ser=1.000000000190015;
        for (j=0;j<=5;j++) ser += cof[j]/++y;
        return -tmp+log(2.5066282746310005*ser/x);
}


#undef NRANSI
#undef ITMAX
#undef EPS
#undef FPMIN
#undef SWAP