/* =================================================================== */
/* module gsm01.c - Garver-Seigel-Maritorena 2001 Bio-optical model    */
/*                                                                     */
/* This module contains the functions to optimize and evaluate the     */
/* Garver-Seigel-Maritorena 2001 Bio-optical model.  The optimization  */
/* is performed for each set of Rrs values within a scanline, using    */
/* the Amoeba simplex minimization technique (Numerical Recipes).      */
/*                                                                     */
/* Reference:                                                          */
/* Maritorena S., D.A. Siegel & A. Peterson, Optimization of a Semi-   */
/* Analytical Ocean Color Model for Global Scale Applications, Applied */
/* Optics,  41(15): 2705-2714, 2002.                                   */
/*                                                                     */
/* Implementation:                                                     */
/* B. Franz, NASA/OCRPG/SAIC, September 2004 (rewrite of previous)     */
/*                                                                     */
/* =================================================================== */

#include <stdio.h>
#include <math.h>
#include "amoeba.h"
#include "l12_proto.h"

#define NPAR   3
#define MAXITR 500

static float lambda [NBANDS];
static float aw     [NBANDS];
static float bbw    [NBANDS];
static float aphstar[NBANDS];
static float admstar[NBANDS];
static float bbpstar[NBANDS];


/* ------------------------------------------------------------------- */
/* get_aphstar() - interpolate aph* to sensor wavelengths              */
/* ------------------------------------------------------------------- */
float get_aphstar(float wave) {

    static long  ntab = 6;
    static float twave[] = {412.0, 443.0, 490.0, 510.0, 555.0, 670.0};
    static float taphs[] = {0.00665,0.05582,0.02055,0.01910,0.01015,0.01424};

    return(linterp(twave,taphs,ntab,wave));
}


/* ------------------------------------------------------------------- */
/* ucsb_amb() GSM01 model, set-up to work with amoeba optimization     */
/*                                                                     */
/* Input:                                                              */
/* auxdata - amoeba interface structure                                */
/* par[]   - model parameters                                          */
/*                                                                     */
/* Output:                                                             */
/* function returns chi-squared, which is the quantity minimized.      */
/*                                                                     */
/* ------------------------------------------------------------------- */
double ucsb_amb(FITSTRUCT *auxdata, double par[])
{
  static int    first = 1;

  static float  adm_s = 0.02061;
  static float  bbp_s = 1.03373;
  static float  grd1  = 0.0949;
  static float  grd2  = 0.0794;

  int    iw;
  float  x, F, bb, ac;
  double merit = 0.0;
  int    nwave = auxdata->npnts;

  if (first) {
      for (iw=0; iw< nwave; iw++) {
          aphstar[iw] = get_aphstar(lambda[iw]);
          admstar[iw] = exp( - adm_s * (lambda[iw] - 443.0));
          bbpstar[iw] = pow((443.0/lambda[iw]), bbp_s);
      }
      first = 0;
  }

  for (iw=0; iw<nwave; iw++) {

      ac = aw [iw] + aphstar[iw]*par[0] + par[1]*admstar[iw];
      bb = bbw[iw] + par[2]*bbpstar[iw];
      x  = bb/(ac + bb);
      F  = grd1*x + grd2*pow(x,2);

      auxdata->yfit[iw] = F;

      merit += pow((auxdata->y[iw] - F), 2) * auxdata->wgt[iw];
  }

  auxdata->merit = merit;

  return(merit);
}


/* ------------------------------------------------------------------- */
/* gsm01_amb() - run amoeba optimization of GSM01 for one spectrum     */
/* ------------------------------------------------------------------- */
int gsm01_amb(double Rrs[], double wts[], long npts, double fitparms[], 
              double Rrs_fit[])
{
  static float tol = 1.e-6;           /* fractional change in chisqr   */
  static FITSTRUCT auxdata;           /* amoeba interface structure    */
  static double init[NPAR*(NPAR+1)];  /* initial simplex               */
  static double p0[] = {.1,.02,.0001};/* starting parameters (C,adg,bb)*/
  static double d0[] = {5.,1.0,0.01}; /* characteristic length         */

  int   i, j;
  short isml;
  int   status = 1;

  auxdata.niter = MAXITR;       /* max number of iterations allowed   */
  auxdata.nfunc = NPAR;         /* number of model parameters         */
  auxdata.npnts = npts;         /* number of wavelengths (Rrs values) */
  auxdata.y     = Rrs;          /* Input Rrs values                   */
  auxdata.wgt   = wts;          /* Input weights on Rrs values        */
  auxdata.yfit  = Rrs_fit;      /* Output model predicted Rrs values  */

  /* initialize simplex with first guess model parameters */
  for (j=0; j<NPAR+1; j++) for (i=0; i<NPAR; i++) init[j*NPAR+i] = p0[i];
  for (i=0; i<NPAR; i++) {
      init[NPAR+i*(NPAR+1)] += d0[i];
      fitparms[i] = 0.0;
  }

  /* run optimization */
  isml = amoeba(init, &auxdata, ucsb_amb, tol); 

  /* check convergence and record parameter results */
  if (auxdata.niter > MAXITR)
      /* failed to converge */
      status = 0;
  else
      for (i = 0; i < NPAR; i++)
          fitparms[i] = init[NPAR * isml + i];

  return(status);
}


/* ------------------------------------------------------------------- */
/* get_gsm01() - returns requested GSM01 product for full scanline     */ 
/* ------------------------------------------------------------------- */
void get_gsm01(l2str *l2rec, l2prodstr *p, float prod[])
{
    static int     firstCall = TRUE;
    static long    lastScan  = -1;
    static double *chl;
    static double *adg;
    static double *bbp;
    static long    nbands;
    static float   Fo[NBANDS];

    int band   = p->prod_ix;
    int prodID = p->cat_ix;

    long   ip, ib, iw, ipb;
    double model  [NPAR];
    double Rrs    [NBANDS];
    double wts    [NBANDS];
    double Rrs_fit[NBANDS];
    int    status;

    if (firstCall == TRUE) {

        firstCall = FALSE;

        if ((chl = calloc(l2rec->npix,sizeof(double))) == NULL) {
            printf("-E- %s line %d : error allocating memory for gsm01.\n",
                __FILE__,__LINE__);
            exit(1);
	}
        if ((adg = calloc(l2rec->npix,sizeof(double))) == NULL) {
            printf("-E- %s line %d : error allocating memory for gsm01.\n",
                __FILE__,__LINE__);
            exit(1);
	}
        if ((bbp = calloc(l2rec->npix,sizeof(double))) == NULL) {
            printf("-E- %s line %d : error allocating memory for gsm01.\n",
                __FILE__,__LINE__);
            exit(1);
	}

        /* build array of actual sensor wavelengths */
        for (iw=0; iw<l2rec->nbands; iw++) {
            ib = l2rec->bindx[iw];
            lambda[iw] = l2rec->fwave[ib]; 
	}

        /* limit to wavelengths less than ~600nm */
        nbands = MIN(windex(600.,lambda,l2rec->nbands)+1, l2rec->nbands);

        /* save associated irradiances and water absorption/backscatter     */
        /* use band-averaged vales if no nLw out-of-band correction applied */

        if (l2rec->input->outband_opt >= 2) {
            for (iw=0; iw<nbands; iw++) {
                ib = l2rec->bindx[iw];
                Fo [iw] = l2rec->Fonom[ib]; 
                aw [iw] = aw_spectra (lambda[iw],10);
                bbw[iw] = bbw_spectra(lambda[iw],10);
	    }
	} else {
            float *awptr, *bbwptr;
            rdsensorinfo(l2rec->sensorID,"aw" ,(void **) &awptr );
            rdsensorinfo(l2rec->sensorID,"bbw",(void **) &bbwptr);
            for (iw=0; iw<nbands; iw++) {
                ib = l2rec->bindx[iw];
                Fo [iw] = l2rec->Fobar[ib];
                aw [iw] = awptr [ib];
                bbw[iw] = bbwptr[ib];
	    }
	}
    }

    /* if this is a new scanline, fit model for every pixel of the scanline */

    if (l2rec->iscan != lastScan) {
        for (ip=0; ip<l2rec->npix; ip++) {

	    status   = 1;
            chl[ip] = -1.0;
            adg[ip] = -1.0;
            bbp[ip] = -1.0;

  	    if (l2rec->mask[ip] == 0) {

	        for (iw=0; iw<nbands; iw++) {
                    wts[iw] = 1.0;
                    Rrs[iw] = l2rec->nLw[ip*NBANDS+l2rec->bindx[iw]]/Fo[iw];
                    if (Rrs[iw] < 0.0) status=0;
		}
                
                /* if any negative reflectance, skip this pixel */
                if (status == 0) {
                    l2rec->flags[ip] |= CHLFAIL;
                    continue;
		}

                /* fit model for this pixel */
		status = gsm01_amb(Rrs,wts,nbands,model,Rrs_fit);

                /* if optimization was successful, save results */
                if (status == 1) {;
                    chl[ip] = model[0];
                    adg[ip] = model[1];
                    bbp[ip] = model[2];
		} else {
                    l2rec->flags[ip] |= CHLFAIL;
		}
   	    } 
	}            
        lastScan = l2rec->iscan;
    }


    /* return requested product, return -1 for abs/scatter in NIR */

    if (band >= 0) {
        if (lambda[band] > 700.0) {
            for (ip=0; ip<l2rec->npix; ip++)
                prod[ip] = -1.0;
            return;
	}
    }

    for (ip=0; ip<l2rec->npix; ip++) {

        switch (prodID) {

	  case CAT_chl_gsm01 : 
            prod[ip] = (float) chl[ip];
            break;

	  case CAT_adg_gsm01 : 
            if (band < 0)
              prod[ip] = (float) adg[ip];
            else 
              prod[ip] = (float) adg[ip]*admstar[band];
            break;

	  case CAT_bbp_gsm01 : 
            if (band < 0)
              prod[ip] = (float) bbp[ip];
            else 
              prod[ip] = (float) bbp[ip]*bbpstar[band];
            break;

  	  case CAT_aph_gsm01 : 
            prod[ip] = (float) aphstar[band]*chl[ip];
            break;

  	  case CAT_a_gsm01 : 
            prod[ip] = (float) aw[band] + aphstar[band]*chl[ip] + adg[ip]*admstar[band];
            break;

	  case CAT_bb_gsm01 : 
            prod[ip] = (float) bbw[band] + bbp[ip]*bbpstar[band];
            break;

          default:
            printf("-E- %s line %d : erroneous product ID %d passed to gsm01.\n",
                __FILE__,__LINE__,prodID);
            exit(1);
	}
    }

    return;
}