/*
 * apf_math.c
 *
 *  Created on: Sep 1, 2024
 *      Author: tm
 */


/******************************************************************************
** 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"

#include <stdint.h>
#include <stdbool.h>
#include <time.h>
#include <sys/types.h>
//#include "stat.h"       // for struct stat down below

#include "system.h"
#include "string.h"
#include "uartstdio.h"
#include "user_io.h"

#include "apf.h"
#include "apf_util.h"
#include "apf_math.h"
#include "apf_numbr.h"
#include "board_util.h"
#include "serial_cmd_fxn.h"

#define NR_MATH  1   // turn math ON (11 Mar 2025), was off=0, did not fit 123g flash or ram

#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;}

extern int mmsDbgFlg; //defd in mms_spec.c

#define MATRIX_BUF_ROW_MAX 3
#define MATRIX_BUF_COL_MAX 512
#define MAX_MATRIX_BUF_SIZE  (MATRIX_BUF_ROW_MAX * MATRIX_BUF_COL_MAX)

//double Matrixbuffer[MATRIX_BUF_ROW_MAX][MATRIX_BUF_COL_MAX];   //2k per row = 6k

//APFTODO - move to apf_math.c
//void InitLambdaArray(double *Lambda, double *ZeissCoefs)
void InitLambdaArray(void)      //globals stored in FRAM, computed once, then just copied at init time
{
    uint32_t indx;

    double p;
    double ZC0 = apfArray.ZeissCoefs[0];
    double ZC1 = apfArray.ZeissCoefs[1];
    double ZC2 = apfArray.ZeissCoefs[2];
    double ZC3 = apfArray.ZeissCoefs[3];
    double ZC4 = apfArray.ZeissCoefs[4];

    //fram_apfArray_retrieve();       // DEBUG
    //uprintf("DBG: InitLambdaArray"); prn_nl();
    //apf_zcoef_print();


// both of these are 10 times slower than below. The pow() seems to be the issue.
//    for( i = 1; i <= 256; i++ ){
//      Lambda[i] = CellNumToLambda(i, ZeissCoefs);
//      Lambda[i] = ZC0 + ZC1*i + ZC2*pow(i, 2.0) + ZC3*pow(i, 3.0) + ZC4*pow(i, 4.0);
//      }

    indx = 0;
    for( p = 1.0; p <= 256.0; p++ )
    {
        apfArray.LambdaData[indx] = ZC0 + ZC1*p + ZC2*p*p + ZC3*p*p*p + ZC4*p*p*p*p;
        //uprintf("%.0lf, %lf\n", p, Lambda[(int)p]);
        indx++;
        wdog_pet();
    }

}


// function to update the values of apf.pixel_beg and _end by feeding in wavelength, aka, lambda values
//int MatchPixelsToLambda(double *Lambda, int lambda_beg, int lambda_end, unsigned short *pixel_beg, unsigned short *pixel_end)  // 1 Apr 2025, mod to uint32_t
int MatchPixelsToLambda(double *Lambda, uint32_t lambda_beg, uint32_t lambda_end, uint32_t *pixel_beg, uint32_t *pixel_end)
{
    unsigned short i;
    double min1, min2;
    double delta1, delta2;
    // JLD - performance enhancement. Stop looking at the entire table.
    double prevDelta1 = 1000.0;
    double prevDelta2 = 1000.0;

    min1=min2=delta1=delta2=100.0;

    for(i = 1; i <= 256; i++) {
        delta1 = fabs((double)(Lambda[i]-lambda_beg));
        wdog_pet();
        delta2 = fabs((double)(Lambda[i]-lambda_end));
        if(delta1 < min1) {
            min1 = delta1;
            prevDelta1 = delta1;
            *pixel_beg = i;               // The index used not off by one for direct array access.  The calculation function was off by one
        }
        if(delta2 < min2) {
            min2 = delta2;
            prevDelta2 = delta2;
            *pixel_end = i;
        }
        // performance enhancement. The table is in order so as soon as the deltas begin
        // to get larger we are done checking.
        if ((delta2 > prevDelta2) && (delta1 > prevDelta1))
            break;

    } //end for
    return( TRUE );

}

//double CellNumToLambda(int cell, double *ZeissCoefs)
//double CellNumToLambda(int cell)
double CellNumToLambda(uint32_t cell)
{
    double lambda = 0.0;
    double dCell;

    dCell = (double)cell;

    // pow(double x,double y)
    lambda = apfArray.ZeissCoefs[0] + apfArray.ZeissCoefs[1]*cell
             + apfArray.ZeissCoefs[2]*pow(dCell, 2.0)
             + apfArray.ZeissCoefs[3]*pow(dCell, 3.0)
             + apfArray.ZeissCoefs[4]*pow(dCell, 4.0);

    return( lambda );

}


void FloatVectorMinMax(unsigned short end, float data[], float *min, float *max)
{
    unsigned short i;

    *max = data[0];
    *min = data[0];

    for(i = 0; i < end; i++) {
        if(data[i] > *max)
            *max = data[i];
        else if(data[i]< *min)
            *min = data[i];
    }

}


void FloatVectorMeanStdDev(unsigned short end, float data[], float *mean, float *stdev)
{
    unsigned short i;
    float sum, stdvar;

    sum = stdvar = 0.0;
    *mean = *stdev = 0.0;

    for(i = 0; i < end; i++)
        sum += data[i];
    *mean = (sum/end);

    for(i = 0; i < end; i++)
        stdvar += pow((data[i] - *mean), 2.0);
    stdvar /= (end - 1);
    *stdev = sqrt(stdvar);

}


void ShortVectorMinMax(unsigned short end, unsigned short data[], unsigned short *min, unsigned short *max)
{
    unsigned short i;

    *max = data[0];
    *min = data[0];

    for(i = 0; i < end; i++) {
        if(data[i] > *max)
            *max = data[i];
        else if(data[i]< *min)
            *min = data[i];
    }

}

void ShortVectorMeanStdDev(unsigned short end, unsigned short data[], float *mean, float *stdev)
{
    unsigned short i;
    float sum, stdvar;

    sum = stdvar = 0.0;
    *mean = *stdev = 0.0;

    for(i = 0; i < end; i++)
        sum += (float)data[i];
    *mean = (sum/end);

    for(i = 0; i < end; i++)
        stdvar += pow((data[i] - *mean), 2.0);
    stdvar /= (end - 1);
    *stdev = sqrt(stdvar);

}

float CtoF(float degC)
{
    return( 9.0/5.0 * degC + 32.0 );
}


/**
 * Evaluates p(x) for a polynom p.
 * Calculates the value of polynom p at x accordings to
 * Horners schema.
 * @param p polynom.
 * @param x value to be inserted into the polynom.
 * @return calculated polynom value.
 */
double poly(double x, int degree, double p[])
{
    double retval = 0.0;
    int i;

    for (i = degree; i >= 0; i--)
        retval = retval * x + p[i];
    return retval;
}

/* Conversion from resistance to temperature.
 * Calculates and returns temperature for given resistance.
 * @param t resistance (in Ohm).
 * @return corresponding temperature.
 */
double rtot(double r, double a[])
{
    double ti;

    ti = poly(log(r), 3, a);
    ti = 1.0 / ti + TABS;
    return ti;
}

/**
   Conversion from temperature to resistance.
 * Calculates and returns resistance for given temperature.
 * @param t temperature (in degree Celsius).
 * @return corresponding resistance.
 */
double ttor(double t, double a[])
{
    double r;
    double u, v, p, q, b, c, d;

    t = t - TABS;
    d = (a[0] - 1.0 / t) / a[3];
    c = a[1] / a[3];
    b = a[2] / a[3];
    q = 2.0 / 27.0 * b * b * b - 1.0 / 3.0 * b * c + d;
    p = c - 1.0 / 3.0 * b * b;
    v = -pow(q / 2.0 + sqrt(q * q / 4.0 + p * p * p / 27.0), 1.0 / 3.0);
    u =   pow(-q / 2.0 + sqrt(q * q / 4.0 + p * p * p / 27.0), 1.0 / 3.0);
    r  = exp(u + v - b / 3.0);
    return r;
}

/*
 * updcrc macro derived from article Copyright (C) 1986 Stephen Satchell.
 *  NOTE: First srgument must be in range 0 to 255.
 *        Second argument is referenced twice.
 *
 * Programmers may incorporate any or all code into their programs,
 * giving proper credit within the source. Publication of the
 * source routines is permitted so long as proper credit is given
 * to Stephen Satchell, Satchell Evaluations and Chuck Forsberg,
 * Omen Technology.
 */
/* crctab calculated by Mark G. Mendel, Network Systems Corporation */
static unsigned short icrctab[256] = {
    0x0000,  0x1021,  0x2042,  0x3063,  0x4084,  0x50a5,  0x60c6,  0x70e7,
    0x8108,  0x9129,  0xa14a,  0xb16b,  0xc18c,  0xd1ad,  0xe1ce,  0xf1ef,
    0x1231,  0x0210,  0x3273,  0x2252,  0x52b5,  0x4294,  0x72f7,  0x62d6,
    0x9339,  0x8318,  0xb37b,  0xa35a,  0xd3bd,  0xc39c,  0xf3ff,  0xe3de,
    0x2462,  0x3443,  0x0420,  0x1401,  0x64e6,  0x74c7,  0x44a4,  0x5485,
    0xa56a,  0xb54b,  0x8528,  0x9509,  0xe5ee,  0xf5cf,  0xc5ac,  0xd58d,
    0x3653,  0x2672,  0x1611,  0x0630,  0x76d7,  0x66f6,  0x5695,  0x46b4,
    0xb75b,  0xa77a,  0x9719,  0x8738,  0xf7df,  0xe7fe,  0xd79d,  0xc7bc,
    0x48c4,  0x58e5,  0x6886,  0x78a7,  0x0840,  0x1861,  0x2802,  0x3823,
    0xc9cc,  0xd9ed,  0xe98e,  0xf9af,  0x8948,  0x9969,  0xa90a,  0xb92b,
    0x5af5,  0x4ad4,  0x7ab7,  0x6a96,  0x1a71,  0x0a50,  0x3a33,  0x2a12,
    0xdbfd,  0xcbdc,  0xfbbf,  0xeb9e,  0x9b79,  0x8b58,  0xbb3b,  0xab1a,
    0x6ca6,  0x7c87,  0x4ce4,  0x5cc5,  0x2c22,  0x3c03,  0x0c60,  0x1c41,
    0xedae,  0xfd8f,  0xcdec,  0xddcd,  0xad2a,  0xbd0b,  0x8d68,  0x9d49,
    0x7e97,  0x6eb6,  0x5ed5,  0x4ef4,  0x3e13,  0x2e32,  0x1e51,  0x0e70,
    0xff9f,  0xefbe,  0xdfdd,  0xcffc,  0xbf1b,  0xaf3a,  0x9f59,  0x8f78,
    0x9188,  0x81a9,  0xb1ca,  0xa1eb,  0xd10c,  0xc12d,  0xf14e,  0xe16f,
    0x1080,  0x00a1,  0x30c2,  0x20e3,  0x5004,  0x4025,  0x7046,  0x6067,
    0x83b9,  0x9398,  0xa3fb,  0xb3da,  0xc33d,  0xd31c,  0xe37f,  0xf35e,
    0x02b1,  0x1290,  0x22f3,  0x32d2,  0x4235,  0x5214,  0x6277,  0x7256,
    0xb5ea,  0xa5cb,  0x95a8,  0x8589,  0xf56e,  0xe54f,  0xd52c,  0xc50d,
    0x34e2,  0x24c3,  0x14a0,  0x0481,  0x7466,  0x6447,  0x5424,  0x4405,
    0xa7db,  0xb7fa,  0x8799,  0x97b8,  0xe75f,  0xf77e,  0xc71d,  0xd73c,
    0x26d3,  0x36f2,  0x0691,  0x16b0,  0x6657,  0x7676,  0x4615,  0x5634,
    0xd94c,  0xc96d,  0xf90e,  0xe92f,  0x99c8,  0x89e9,  0xb98a,  0xa9ab,
    0x5844,  0x4865,  0x7806,  0x6827,  0x18c0,  0x08e1,  0x3882,  0x28a3,
    0xcb7d,  0xdb5c,  0xeb3f,  0xfb1e,  0x8bf9,  0x9bd8,  0xabbb,  0xbb9a,
    0x4a75,  0x5a54,  0x6a37,  0x7a16,  0x0af1,  0x1ad0,  0x2ab3,  0x3a92,
    0xfd2e,  0xed0f,  0xdd6c,  0xcd4d,  0xbdaa,  0xad8b,  0x9de8,  0x8dc9,
    0x7c26,  0x6c07,  0x5c64,  0x4c45,  0x3ca2,  0x2c83,  0x1ce0,  0x0cc1,
    0xef1f,  0xff3e,  0xcf5d,  0xdf7c,  0xaf9b,  0xbfba,  0x8fd9,  0x9ff8,
    0x6e17,  0x7e36,  0x4e55,  0x5e74,  0x2e93,  0x3eb2,  0x0ed1,  0x1ef0
};

#define updcrc(cp, crc) ( icrctab[((crc >> 8) & 255)] ^ (crc << 8) ^ cp)

unsigned short calcCRC(unsigned char val, unsigned short crc)
{
    crc = updcrc((0xff & val), crc);

    return(crc);
}

unsigned short fCrc16Bit(const unsigned char *msg)
{
    unsigned short crc = 0xffff;

    if(mmsDbgFlg == 1) { uprintf("fCrc16Bit needs imp, apf_util.c"); prn_nl(); }
//#ifdef APFTODO
    // dangerous - need to insure *msg is null term
    while (*msg)
    {
        crc = updcrc((0xff & *msg), crc);
        msg++;
    }

    crc = updcrc(0,updcrc(0,crc));  // macro icrtab.c  see zmodem
//#endif
    return( (crc & 0xffff) );

}


//1L,256L,1L,3L
double **statmatrix(unsigned char *arrayRowPtr)         //, uint32_t arrayByteCnt)
/* build a table of row pointers to the static arrays in apfArray for use with LinFit */
{

    long nrl;
    long nrh;
    long ncl;
    long nch;

    long i, ncol;
    double **m;


    nrl = 1;
    nrh = 256;
    ncl = 1;
    nch = 3;

    //nrow = nrh - nrl + 1;       //nrow = 256 - 1 + 1 = 256
    ncol = nch - ncl + 1;       //ncol = 3 - 1 + 1 = 3

#define NR_END 1
    /* allocate pointers to rows */

    m = (double **)arrayRowPtr;  //*EcoefFileRowPtr[260]; or *EcoefDataRowPtr[260];

    for(i=nrl+1;i<=nrh;i++)
    {
        m[i]=m[i-1]+ncol;
    }


    /* return pointer to array of pointers to rows */
    return m;
}

//remove args: &apfArray.EcoefData[apf_eeprom.isus_pixel_beg-1], &apfArray.LambdaData[apf_eeprom.isus_pixel_beg-1],
//int LinFit(double *AbsVals, double *AbsFitVals, double *ResVals, double *FitVals, ushort fit_pixels, ushort fit_concs, ushort bsl_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;


    //uprintf("\nLinFit()%lf,%lf,%.9lf,%.9lf,%.9lf\n", LambdaVals[1], AbsVals[1], CoefVals[1][0], CoefVals[1][1], CoefVals[1][2]);
    //uprintf("\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);
}



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


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

    }

}

#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);
        }
    }

}



