/*---------------------------------------------------------------------*/
/* calcite.c -  get calcium carbonate concentration.                   */
/*                                                                     */
/* Inputs:                                                             */
/*     l2rec - level-2 structure containing one complete scan after    */
/*             atmospheric correction.                                 */
/* Outputs:                                                            */
/*     caco3 - calcium carbonate concentration, per pixel        .     */
/*                                                                     */
/* Written by: W. Robinson, GSC, 7 Jun 2000.                           */
/*             S. Bailey, OCDPG, July 2004, conversion to C.           */
/*             B. Franz, OCDPG, Sep 2004, sensor generalization and    */
/*                 implementation of 2-Band algorithm.                 */
/*                                                                     */
/*---------------------------------------------------------------------*/

#include <stdlib.h>
#include <math.h>
#include "l12_proto.h"
#include "l2_flags.h"

static float pi        = PI;
static float radeg     = RADEG;
static long  caco3_msk = ATMFAIL | LAND | HISUNGLINT | CLOUD;
static float caco3min  = 0.003;

/* --------------------------------------------------------------------- */
/* calcite_3b() - calcium carbonate concentration from 3-Band algorith.. */
/*                                                                       */
/* Gordon, H.R. Boynton, G.C., Balch, W.M., Groom, S.B., Harbour, D.S.,  */
/* Smyth, T.J., Retrieval of Coccolithophore Calcite Concentration from  */
/* SeaWiFS Imagery, GRL, 28, 8, 1587-1590.                               */
/*                                                                       */
/* --------------------------------------------------------------------- */

void calcite_3b(l2str *l2rec, float caco3[])
{
    static int   firstCall = 1;
    static int   maxiter   = 10;
    static float bbw546    = 0.000939;
    static float wave[3]   = {670.,760.,860.};    /* approx. wavelengths */
    static int   bx  [3];
    static float aw  [3];
    static float bbw [3];
    static float bbc [3];
    static float t   [3];
    static float b68diff;
    static float b78diff;
    static float fw1,fw2;
    
    static float oobswf[3][8] = { 
    {0.000313529, 0.000770558, 0.00152194,  0.000155573, 
     0.00116455,  0.0,         0.000445433, 0.000124172},
    {0.000201709, 6.96143e-05, 7.00147e-06, 2.28957e-07, 
     4.17788e-05, 0.00159814,  0.0,         0.00536827},
    {0.000463807, 8.54003e-05, 2.47401e-05, 0.000755890, 
     0.00587073,  0.00021686,  0.0111331,   0.0} };

    long   ip, ipb, ib, iw, i;
    float  rhoaw[NBANDS];
    float  rho[3];
    float  bbc_cclth, r8_cclth, aeps_cclth;
    float  mu, mu0, airmass;
    float  bbcinit, bbctol;
    int    numiter;
    long   status = 0;
    float *awptr, *bbwptr;

    if (firstCall) {

        /* get band-averaged absorption and backscatter for pure water */
        rdsensorinfo(l2rec->sensorID,"aw" ,(void **) &awptr );
        rdsensorinfo(l2rec->sensorID,"bbw",(void **) &bbwptr);

        /* save coeffs for the three bands, resolve actual sensor wave */
        for (i=0; i<3; i++) {
            bx  [i] = windex(wave[i],l2rec->fwave,NBANDS);
            wave[i] = l2rec->fwave[bx[i]];
            aw  [i] = awptr [bx[i]];
	    bbw [i] = bbwptr[bx[i]];
            bbc [i] = 0.0;
	}

        b68diff = wave[2]-wave[0];
        b78diff = wave[1]-wave[0];

        /* spectral dependence of bbc */
        fw1 = pow(wave[0]/wave[1],1.35);
        fw2 = pow(wave[0]/wave[2],1.35);

        firstCall = 0;
    }

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

        status    = 0;
        numiter   = 0;
        bbctol    = 100.;
        bbcinit   = 0.00085;
        caco3[ip] = 0.0;
        bbc[0]    = 0.000;

        /* skip pixel if already masked */
        if ( (l2rec->flags[ip] & caco3_msk) != 0 ) {
	    l2rec->flags[ip] |= CHLFAIL;
            continue;
	}

        mu0 = cos( l2rec->solz[ip]/radeg );
        mu  = cos( l2rec->senz[ip]/radeg );
        airmass = 1.0/mu0 + 1.0/mu;

        /* compute the aerosol/water reflectance (include out-of-band correction for SeaWiFS) */
        if (l2rec->sensorID == SEAWIFS) {
            for( iw = 0; iw < l2rec->nbands; iw++ ) {
	        ib  = l2rec->bindx[iw];
                ipb = NBANDS*ip + ib;
                rhoaw[ib] = ( ( (l2rec->Lt[ipb] - l2rec->tLf[ipb])/l2rec->t_oz_sol[ipb]/l2rec->t_oz_sen[ipb]
                          - l2rec->Lr[ipb] ) / l2rec->t_o2[ipb] - l2rec->TLg[ipb]) * pi/l2rec->Fo[ib]/mu0;
            }
            for (i = 0; i < 3; i++) {
	        rho[i] = rhoaw[bx[i]];
                for( iw = 0; iw < l2rec->nbands; iw++ ) {
		    ib = l2rec->bindx[iw];
                    rho[i] -= rhoaw[ib]*oobswf[i][iw]; 
                    if (rho[i] <= 0.0) status = 1;
		}
            }
        } else {
            for( i = 0; i < 3; i++ ) {
	        ib  = bx[i];
                ipb = NBANDS*ip + ib;
                rho[i] = ( ( (l2rec->Lt[ipb] - l2rec->tLf[ipb])/l2rec->t_oz_sol[ipb]/l2rec->t_oz_sen[ipb]
                       - l2rec->Lr[ipb] ) / l2rec->t_o2[ipb] - l2rec->TLg[ipb]) * pi/l2rec->Fo[ib]/mu0;
                if (rho[i] <= 0.0) status = 1;
            }
	}

        /* skip pixel on negative surface reflectance */
        if (status != 0) {
 	    l2rec->flags[ip] |= CHLFAIL;
 	    continue;
	}


	/* compute total transmittance */
        for (i=0; i<3; i++) {        
            ipb  = NBANDS*ip + bx[i];
            t[i] = l2rec->t_oz_sol[ipb]*l2rec->t_oz_sen[ipb]*l2rec->t_sol[ipb]*l2rec->t_sen[ipb];
	}

        /* compute backscatter at 546 nm */
        while (bbctol > 5. && numiter < maxiter) {

            numiter++;

            bbc[1] = bbc[0]*fw1;
            bbc[2] = bbc[0]*fw2;

            /* reflectance at longest wavelength */        
            r8_cclth = rho[2] - (bbw[2]+bbc[2])/(aw[2]+bbw[2]+bbc[2])/6.0*t[2];

            if ((r8_cclth > 0.09) || (r8_cclth < 0)){
                status = 1;
                bbc[0] = 0;
                break;
            }

            /* atmospheric epsilon at two longest wavelengths */        
            aeps_cclth = log((rho[1] - (bbw[1]+bbc[1])/(aw[1]+bbw[1]+bbc[1])/6.0*t[1])/r8_cclth)/b78diff;

            if (aeps_cclth > 0.4){
                status = 1;
                bbc[0] = 0;
                break;
            }

            /* --------------- */
            bbc[0] = (rho[0] - r8_cclth * exp(aeps_cclth*b68diff))/t[0] * (aw[0]+bbw[0]+bbc[0])*6.0 - bbw[0];

            if ((bbc[0] <= 0) || isnan(bbc[0])) {
                status = 1;
                bbc[0] = 0;
                break;
            }
       
            bbctol  = fabs((bbcinit - bbc[0])/bbcinit)*100.;
            bbcinit = bbc[0];                        
        }


        if (status == 0) {
            bbc_cclth = bbc[0]/pow((546./wave[0]),1.35) + bbw546;    
            caco3[ip] = bbc_cclth * 0.625 + 0.00225; 
            if (caco3[ip] < caco3min ) {
	        l2rec->flags[ip] |= CHLRANGE;
	    }
        } else {
	    l2rec->flags[ip] |= CHLFAIL;
	}

    } /* pixel loop */

    return;
}


/* --------------------------------------------------------------------- */
/* calcite_2b() - calcium carbonate concentration from 2-Band algorith.. */
/*                                                                       */
/* Gordon, H.R. and Balch, W.M., MODIS Detached Coccolith Concentration  */
/* Algorithm Theoretical Basis Document,  April 30, 1999                 */
/*                                                                       */
/* --------------------------------------------------------------------- */

#define N443 490
#define N550 456

void calcite_2b(l2str *l2rec, float caco3[])
{
    static int   firstCall = 1;
    static int   n443 = N443;
    static int   n550 = N550;
    static float t443[N443]; 
    static float t550[N550]; 
    static float tbb [N443][N550];
    static float tpig[N443][N550];
    static long  ib443;
    static long  ib550;

    int   i443, i550;
    float x,y,a,b;
    long  i,nc, ip;
    float x443, x550;
    float bb1, bb2, bb;
    float pig1, pig2, pig;

    if (firstCall) {

        FILE *fp;
        char *filedir;
        char filename[FILENAME_MAX];
        char line [80];

        ib443 = windex(443.,l2rec->fwave,NBANDS);
        ib550 = windex(550.,l2rec->fwave,NBANDS);

        /* Locate coccolith table using environment variable */
        if ((filedir = getenv("MSL12_DATA")) == NULL) {
            printf("-E- %s: MSL12_DATA env variable undefined.\n", __FILE__);
            exit(1);
        }
        strcpy(filename, filedir); strcat(filename, "/common/coccolith.dat");

        printf("Loading coccolithophore table %s\n",filename);

        if ( (fp = fopen(filename,"r")) == NULL ) {
            fprintf(stderr,"-E- %s line %d: unable to open %s for reading\n",
                    __FILE__,__LINE__,filename);
            exit(1);
        }


        /* Skip comment lines */
        nc = 0;
        while ( fgets(line,80,fp ) ) {
	    if (line[0] != '#' && line[0] != '\n' ) {
                break;
	    }
	    nc++;
	}
        rewind(fp); for (i=0; i<nc; i++) fgets(line,80,fp);

        /* Load table */
        for (i443=0; i443<n443; i443++) for (i550=0; i550<n550; i550++) {
            fscanf(fp,"%f %f %f %f\n",&x,&y,&a,&b);
            t443[i443] = x;
            t550[i550] = y;
            tbb [i443][i550] = a;
            tpig[i443][i550] = b;
	}

        firstCall = 0;
    }

    for (ip=0; ip<l2rec->npix; ip++) {
    
        caco3[ip] = 0.0;

        if ( (l2rec->flags[ip] & caco3_msk) != 0 || l2rec->mask[ip]) {
	    l2rec->flags[ip] |= CHLFAIL;
            continue;
	}

        x443 = l2rec->nLw[ip*NBANDS+ib443];
        x550 = l2rec->nLw[ip*NBANDS+ib550];

        if (x443 <= 0.0 || x550 <= 0.0) {
	    l2rec->flags[ip] |= CHLFAIL;
            continue;
	}

        /* locate bounding table indices */
        for (i=0; i<n443; i++) {
	    if (x443 < t443[i]) {
	        i443 = i;
		break;
	    }
	}
        for (i=0; i<n550; i++) {
	    if (x550 < t550[i]) {
	        i550 = i;
		break;
	    }
	}

        if (i443 <=0 || i443 >= n443 || i550 <=0 || i550 >= n550) {
	    l2rec->flags[ip] |= CHLFAIL;
            continue;
	}

        if (tbb[i443-1][i550-1] > 998.9 || tbb[i443  ][i550-1] > 998.9 || 
            tbb[i443-1][i550  ] > 998.9 || tbb[i443  ][i550  ] > 998.9) {
	    l2rec->flags[ip] |= CHLFAIL;
            continue;
	}

        /* interpolate to get bb(546) */
        bb1 = tbb[i443-1][i550-1] + (x443 - t443[i443-1])*
          (tbb[i443][i550-1] - tbb[i443-1][i550-1])/(t443[i443] - t443[i443-1]);

        bb2 = tbb[i443-1][i550  ] + (x443 - t443[i443-1])*
          (tbb[i443][i550  ] - tbb[i443-1][i550  ])/(t443[i443] - t443[i443-1]);

        bb = bb1 + (x550 - t550[i550-1]) * (bb2 - bb1)/(t550[i550] - t550[i550-1]);

	/* coccolith pigment, currently disabled
        pig1 = tpig[i443-1][i550-1] + (x443 - t443[i443-1])* 
          (tpig[i443][i550-1] - tpig[i443-1][i550-1])/ (t443[i443] - t443[i443-1]);

        pig2 = tpig[i443-1][i550] + (x443 - t443[i443-1])* 
          (tpig[i443][i550  ] - tpig[i443-1][i550  ])/ (t443[i443] - t443[i443-1]);

        pig = pig1 + (x550 - t550[i550-1]) * (pig2 - pig1)/(t550[i550] - t550[i550-1]);
	*/

        /* convert to calcite in moles/m^3 (ref. HRG 1999) */
        caco3[ip] = bb * 0.625 + 0.00225;
        if (caco3[ip] < caco3min ) {
            l2rec->flags[ip] |= CHLRANGE;
	}
    }
}