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

static float badval =    -1;
static float minval =   0.0;
static float maxval = 100.0;

float get_tsm_clark_seawifs(float nLw[])
{
    float tsm = badval;
    float x;

    static float a[6] = {0.51897,-2.24106,1.20113,-4.35315,9.07162,-5.10552};

    if (nLw[1] <= 0.0 || nLw[2] <= 0.0 || nLw[4] <= 0.0)
            return(tsm);

    x   = log10((nLw[1]+nLw[2])/nLw[4]);
    tsm = pow(10.0,(a[0] + x*(a[1] + x*(a[2] + x*(a[3] + x*(a[4] + x*a[5]))))));

    return (tsm);
}    


float get_tsm_clark_modis(float nLw[])
{
    float tsm = badval;
    float x;

    static float a[6] = {0.490330,-2.712882,3.412666,-8.336478,12.111023,-5.961926};

    if (nLw[1] <= 0.0 || nLw[2] <= 0.0 || nLw[4] <= 0.0)
            return(tsm);

    x   = log10((nLw[1]+nLw[2])/nLw[4]);
    tsm = pow(10.0,(a[0] + x*(a[1] + x*(a[2] + x*(a[3] + x*(a[4] + x*a[5]))))));

    return (tsm);
}    


void get_tsm(l2str *l2rec, int prodnum, float prod[])
{
    long  ip;

    /*                                                      */
    /* Compute desired products at each pixel               */
    /*                                                      */
    for (ip=0; ip<l2rec->npix; ip++) {

        switch (prodnum) {

	case CAT_tsm_clark:
            switch (l2rec->sensorID) {
                case SEAWIFS:
                    prod[ip] = get_tsm_clark_seawifs(&l2rec->nLw[ip*NBANDS]);
                    break;
                case MODISA:
                case MODIST:
                    prod[ip] = get_tsm_clark_modis(&l2rec->nLw[ip*NBANDS]);
                    break;
                default:
		  prod[ip] = badval;
                    break;
            }
            if (prod[ip] == badval)
	        l2rec->flags[ip] |= CHLFAIL;
            else if (prod[ip] < minval || prod[ip] > maxval) 
	        l2rec->flags[ip] |= CHLRANGE;
            break;

           default:
            printf("Error: %s : Unknown product specifier: %d\n",__FILE__,prodnum);
            exit(FATAL_ERROR);
            break;
        }

    }

}