/* =================================================================== */
/* Algorithms for computing diffuse attenuation coefficient for MSl12. */
/*                                                                     */
/* B. Franz, NASA Ocean Biology Processing Group, SAIC, March 2005.    */
/* =================================================================== */

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

#define KD_BAD -1.0
#define KD_MAX  6.4

/* ------------------------------------------------------------------- */
/* Kd490_mueller() -  diffuse attenuation at 490nm (J. Mueller).       */
/*                                                                     */
/* Inputs:                                                             */
/*     l2rec - level-2 structure containing one complete scan after    */
/*             atmospheric correction.                                 */
/* Outputs:                                                            */
/*     k490 - diffuse attenuation coefficient, 1 value per pixel.      */
/*                                                                     */
/* Algorithm provided by: J. Mueller                                   */
/* OCTS/POLDER coefficents: S. Bailey, 16 July 2001                    */
/* ------------------------------------------------------------------- */
void Kd490_mueller(l2str *l2rec, float k490[])
{
    long   ip, ipb;

    static float a  =  0.15645;
    static float b  = -1.5401;
    static float Kw =  0.016;

    static float badval = 0.0;
    static float maxval = 6.4;

    if (l2rec->sensorID == OCTS || l2rec->sensorID == POLDER) {
        /* Coefs for 490/565 band combination */ 
        a =  0.2166;
        b = -1.6355;
    }

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

        ipb = NBANDS*ip;

        if (l2rec->mask[ip]) {                   /* pixel was masked          */
            k490[ip] = badval;
        } else if (l2rec->nLw[ipb+2] <= 0.0 ) {  /* unknown, high attenuation */
            k490[ip] = maxval;
        } else if (l2rec->nLw[ipb+4] <= 0.0 ) {  /* unknown, low  attenuation */
            k490[ip] = Kw;
        } else {
            k490[ip] = a*pow(l2rec->nLw[ipb+2]/l2rec->nLw[ipb+4],b) + Kw;
            if (k490[ip] > maxval) {
                k490[ip] = maxval;
	    }
	}
    }
}


/* ------------------------------------------------------------------- */
/* Kd490_obpg -  diffuse attenuation at 490nm (Mueller & Werdell).     */
/*                                                                     */
/* Inputs:                                                             */
/*     l2rec - level-2 structure containing one complete scan after    */
/*             atmospheric correction.                                 */
/* Outputs:                                                            */
/*     k490 - diffuse attenuation coefficient, 1 value per pixel.      */
/*                                                                     */
/* Algorithm provided by: J. Mueller                                   */
/* New coefficents and updated form: P.J.Werdell, February 2005.       */
/* ------------------------------------------------------------------- */
void Kd490_obpg(l2str *l2rec, float k490[])
{
    long   ip, ipb;

    static float a;
    static float b;
    static long ib490;
    static long ib555;

    static firstCall = 1;

    if (firstCall) {

        long *bindx = l2rec->bindx;
        long  nwave = l2rec->nbands;
        float wave[NBANDS];
        long  iw555, iw565, iw;

        /* select 490/555 or 490/565 fit, depending on proximity of */
        /* sensor bands to fitted bands                             */

        if ((ib555 = bindex_get(555)) > 0) {
            a =  0.1853;
            b = -1.349;
        } else if ((ib555 = bindex_get(565)) > 0) {
            a =  0.1787;
            b = -1.122;
	} else {
	  printf("Kd_obpg: incompatible sensor wavelengths (no 555 or 565).\n");
          exit(1);
	}

        if ((ib490 = bindex_get(490)) < 0) {
	  printf("Kd_obpg: incompatible sensor wavelengths (no 490).\n");
          exit(1);
	}

        firstCall = 0;
    }

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

        ipb = NBANDS*ip;

        if (l2rec->mask[ip] || 
            l2rec->nLw[ipb+ib490] <= 0.0 || l2rec->nLw[ipb+ib555] <= 0.0) {
            k490[ip] = KD_BAD;
            l2rec->flags[ip] |= CHLFAIL;
        } else {
            k490[ip] = a*pow(l2rec->nLw[ipb+ib490]/l2rec->nLw[ipb+ib555],b);
            if (k490[ip] > KD_MAX) {
                k490[ip] = KD_MAX;
                l2rec->flags[ip] |= CHLRANGE;
	    }
	}

    }
}


/* ------------------------------------------------------------------- */
/* Kd_lee() - spectral diffuse attenuation using Lee, et. (2005)       */
/*                                                                     */
/* Inputs:                                                             */
/*     l2rec - level-2 structure containing one complete scan after    */
/*             atmospheric correction.                                 */
/*     band  - waveband number (0 - nbands-1) at which Kd computed.    */
/*                                                                     */
/* Outputs:                                                            */
/*     Kd - diffuse attenuation at specified band, 1 value per pixel.  */
/*                                                                     */
/* Description:                                                        */
/*  This produces the estimate of diffuse attenation at the given band */
/*  using the satellite derived total absorption and backscattering.   */
/*                                                                     */
/* Reference:                                                          */
/*  Lee, Z.P., M. Darecki, K.L. Carder, C.O.Davis,                     */
/*  D. Stramski, W.J. Rhea, "Diffuse Attenuation coefficient of        */
/*  downwelling irradiance:  An evalution of remote sensing methods".  */
/*                                                                     */
/* Original Implementation: P. Martinolich, NRL/Stennis, March 2005    */
/*---------------------------------------------------------------------*/
void Kd_lee(l2str *l2rec, int band, float *Kd)
{

    const float m1 =   4.18;
    const float m2 =   0.52;
    const float m3 = -10.80;

    float m0;
    int ip, ipb;

    if (l2rec->input->iop_opt == IOPNONE) {
        printf("IOP-based Kd_lee product requires iop model selection (iop_opt).  ");
        printf("Using default model.\n");
        l2rec->input->iop_opt = IOPDEFAULT;
        get_iops(l2rec,l2rec->input->iop_opt);
    }

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

        ipb = ip*NBANDS + band;

        if (l2rec->mask[ip] || 
            l2rec->a[ipb] <= 0.0 || l2rec->bb[ipb] <= 0.0 ) {
	    Kd[ip] = KD_BAD;
            l2rec->flags[ip] |= CHLFAIL;
	} else {
	    m0 = 1.0 + 0.005 * l2rec->solz[ip];
	    Kd[ip] = m0 * l2rec->a[ipb]
		       + m1 * (1.0 - m2 * exp( m3 * l2rec->a[ipb])) * l2rec->bb[ipb];
        }
    }

    return;
}


/* ------------------------------------------------------------------- */
/* get_Kd() - l2_hdf_generic interface for Kd                          */
/* ------------------------------------------------------------------- */
void get_Kd(l2str *l2rec, l2prodstr *p, float prod[])
{
    switch (p->cat_ix) {
	case CAT_Kd_mueller:
	    Kd490_mueller(l2rec,prod);
            break;
	case CAT_Kd_obpg:
	    Kd490_obpg(l2rec,prod);
            break;
	case CAT_Kd_lee:
	    Kd_lee(l2rec,p->prod_ix,prod);
            break;
    }

    return;
}