/disk01/web/ocssw/build/src/l2gen/qaa.c (r8106/r7411)

Go to the documentation of this file.

#include <stdlib.h>
#include <math.h>
#include <stdarg.h>
#include <assert.h>
#include "qaa.h"

static int idx410 = -1;
static int idx440 = -1;
static int idx490 = -1;
static int idx555 = -1;
static int idx670 = -1;
static int initialized  = 0;
static int aph_check = 1;
static double S = 0.015;
static double chl = -1;
static double acoefs[3];

int
qaa_is_initialized( void )
{
    return initialized;
}

int
qaa_init( int i410, int i440, int i490, int i555, int i670 )
{
    idx410 = i410;
    idx440 = i440;
    idx490 = i490;
    idx555 = i555;
    idx670 = i670;

//  SeaWiFS coefficients

    acoefs[0] = -1.146;
    acoefs[1] = -1.366;
    acoefs[2] = -0.469;

    initialized = 1;

    return 0;
}

int
qaa_set_param( int param, ... )
{
    va_list  ap;
    va_start( ap, param );
    switch (param) {
    case QAA_S_PARAM:
        S = va_arg(ap,double);
        break;
    case QAA_CHL_PARAM:
        chl = va_arg(ap,double);
        break;
    case QAA_COEFS_PARAM:
        acoefs[0] = va_arg(ap,double);
        acoefs[1] = va_arg(ap,double);
        acoefs[2] = va_arg(ap,double);
        break;
    case QAA_APH_CHECK:
        aph_check = va_arg(ap,int);
        break;
    }
    va_end(ap);
    return 0;
}

int
qaa_v6( int nbands, double *wavel, double *Rrs, double *aw, double *bbw,
         double *rrs, double *u, double *a, double *bb, unsigned char *flags  )
{

    const double g0 = 0.08945;
    const double g1 = 0.1247;

    int     i, idxref;

    double  rat, aref;
    double  bbpref;
    double  Y;
    double  rho, numer, denom;

    assert( idx440 >= 0 );
    assert( idx490 >= 0 );
    assert( idx555 >= 0 );
    assert( nbands >= 0 );

    /* pre-test 670 */
    if ( (Rrs[idx670] > 20.0*pow(Rrs[idx555],1.5)) ||
         (Rrs[idx670] < 0.9*pow(Rrs[idx555],1.7)) ) {
        Rrs[idx670] = 1.27*pow(Rrs[idx555],1.47) + 0.00018*pow(Rrs[idx490]/Rrs[idx555],-3.19);
        if ( flags ) *flags |= 0x02;
    }

    /* Step 0 */
    for ( i = 0; i < nbands; i++ ) {
        rrs[i] = Rrs[i] / (0.52 + 1.7 * Rrs[i] );
        if ( Rrs[i] < 0.0 )
            if ( flags ) *flags |= 0x01;
    }

    /* Step 1 */
    for ( i = 0; i < nbands; i++ )
        u[i]   = (sqrt(g0*g0 + 4.0*g1*rrs[i]) - g0) / (2.0 * g1);

    /* Step 2 */
    if ( Rrs[idx670] >= 0.0015 ) {
        aref    = aw[idx670] + 0.07*pow(Rrs[idx670]/Rrs[idx440],1.1);
        idxref  = idx670;
    } else {
        numer   = Rrs[idx440] + Rrs[idx490];
        denom   = Rrs[idx555] + 5*Rrs[idx670]*(Rrs[idx670]/Rrs[idx490]);
        rho     = log10( numer / denom );
        rho     = acoefs[0] + acoefs[1]*rho + acoefs[2]*rho*rho;
        aref    = aw[idx555] + pow(10.0,rho);
        idxref  = idx555;
    }

    /* Step 3 */
    bbpref = ((u[idxref] * aref) / (1.0 - u[idxref])) - bbw[idxref];

    /* Step 4 */
    rat    = rrs[idx440] / rrs[idx555];
    Y = 2.0 * (1.0 - 1.2 * exp( -0.9*rat) );

    /* Step 5 */
    for ( i = 0; i < nbands; i++ ) {
        bb[i] = bbpref * pow((wavel[idxref]/wavel[i]),Y) + bbw[i];
        if ( bb[i] < 0.0 )
            if ( flags ) *flags |= 0x04;
    }

    /* Step 6 */
    for ( i = 0; i < nbands; i++ ) {
        a[i] = ((1.0 - u[i]) * bb[i]) / u[i];
        if ( a[i] < 0.0 )
            if ( flags ) *flags |= 0x08;
    }

    return 0;
}

int
qaaf_v6( int nbands, float *wavel, float *Rrs, float *aw, float *bbw,
         float *rrs, float *u, float *a, float *bb,
         unsigned char *flags  )
{

    const float g0 = 0.08945;
    const float g1 = 0.1247;

    float   rho, numer, denom;

    float   rat, aref;
    float   bbpref;
    float   Y;

    int     i, idxref;

    assert( idx440 >= 0 );
    assert( idx490 >= 0 );
    assert( idx555 >= 0 );
    assert( idx670 >= 0 );
    assert( nbands >= 0 );

    /* Test for bad Rrs at idx555 */
    if ( Rrs[idx555] <= 0.0 )
        Rrs[idx555] = 0.001;

    /* pre-test 670 */
    if ( (Rrs[idx670] > 20.0*powf(Rrs[idx555],1.5)) ||
         (Rrs[idx670] < 0.9*powf(Rrs[idx555],1.7)) ) {
        Rrs[idx670] = 1.27*powf(Rrs[idx555],1.47) + 0.00018*powf(Rrs[idx490]/Rrs[idx555],-3.19);
        if ( flags ) *flags |= 0x02;
    }

    /* Step 0 */
    for ( i = 0; i < nbands; i++ ) {
        rrs[i] = Rrs[i] / (0.52 + 1.7 * Rrs[i] );
        if ( Rrs[i] < 0.0 )
            if ( flags ) *flags |= 0x01;
    }

    /* Step 1 */
    for ( i = 0; i < nbands; i++ )
        u[i]   = (sqrt(g0*g0 + 4.0*g1*rrs[i]) - g0) / (2.0 * g1);

    /* Step 2 */
    if ( Rrs[idx670] >= 0.0015 ) {
        aref   = aw[idx670] + 0.07*pow(Rrs[idx670]/Rrs[idx440],1.1);
        idxref = idx670;
    } else {
        numer  = Rrs[idx440] + Rrs[idx490];
        denom  = Rrs[idx555] + 5*Rrs[idx670]*(Rrs[idx670]/Rrs[idx490]);
        rho    = log10f( numer / denom );
        rho    = acoefs[0] + acoefs[1]*rho + acoefs[2]*rho*rho;
        aref   = aw[idx555] + powf(10.0,rho);
        idxref = idx555;
    }

    /* Step 3 */
    bbpref = ((u[idxref] * aref) / (1.0 - u[idxref])) - bbw[idxref];

    /* Step 4 */
    rat   = rrs[idx440] / rrs[idx555];
    Y = 2.0 * (1.0 - 1.2 * expf( -0.9*rat) );

    /* Step 5 */
    for ( i = 0; i < nbands; i++ ) {
        bb[i] = bbpref * pow((wavel[idxref]/wavel[i]),Y) + bbw[i];
        if ( bb[i] < 0.0 )
            if ( flags ) *flags |= 0x04;
    }

    /* Step 6 */
    for ( i = 0; i < nbands; i++ ) {
        a[i] = ((1.0 - u[i]) * bb[i]) / u[i];
        if ( a[i] < 0.0 )
            if ( flags ) *flags |= 0x08;
    }

    return 0;
}

int
qaa_decomp( int nbands, double *wavel, double *rrs, double *a, double *aw,
            double *adg, double *aph, unsigned char *flags )
{
    int     i;
    double  symbol, x1, x2;
    double  zeta, denom, dif1, dif2;
    double  rat, ag440;
    double  Sr;

    assert( idx410 >= 0 );
    assert( idx440 >= 0 );
    assert( idx555 >= 0 );
    assert( nbands >= 0 );

    /* step 7 */
    rat    = rrs[idx440] / rrs[idx555];
    symbol = 0.74 + ( 0.2 / ( 0.8 + rat ) );

    /* step 8 */
    Sr   = S + 0.002 / ( 0.6 + rat );
    zeta = exp( Sr * (wavel[idx440] - wavel[idx410]) );

    /* step 9 */
    denom = zeta - symbol;
    dif1  = a[idx410]  - symbol * a[idx440];
    dif2  = aw[idx410] - symbol * aw[idx440];
    ag440 = (dif1 - dif2) / denom;

    for ( i = 0; i < nbands; i++ ) {
        adg[i] = ag440 * exp( Sr * (wavel[idx440] - wavel[i]));
    aph[i] = a[i] - adg[i] - aw[i];
    }

    /* check aph443 range */

    if ( aph_check ) {

    x1 = aph[idx440] / a[idx440];
    if ( x1 < 0.15 || x1 > 0.6 ) {

        if ( flags ) *flags |= 0x10;

        x2 = -0.8 + 1.4 * (a[idx440] - aw[idx440])/(a[idx410] - aw[idx410]);
        if ( x2 < 0.15 ) {
        x2 = 0.15;
        if ( flags ) *flags |= 0x20;
        }
        if ( x2 > 0.6 ) {
        x2 = 0.6;
        if ( flags ) *flags |= 0x40;
        }

        aph[idx440] = a[idx440] * x2;
        ag440 = a[idx440] - aph[idx440] - aw[idx440];

        for ( i = 0; i < nbands; i++ ) {
            adg[i] = ag440 * exp( Sr * (wavel[idx440] - wavel[i]));
            aph[i] = a[i] - adg[i] - aw[i];
        }

        }
    }
    return 0;
}

int
qaaf_decomp( int nbands, float *wavel, float *rrs, float *a, float *aw,
             float *adg, float *aph, unsigned char *flags )
{
    int    i;
    float  symbol, x1, x2;
    float  zeta, denom, dif1, dif2;
    float  rat, ag440;
    float  Sr;

    assert( idx410 >= 0 );
    assert( idx440 >= 0 );
    assert( idx555 >= 0 );
    assert( nbands >= 0 );

    /* step 7 */
    rat    = rrs[idx440] / rrs[idx555];
    symbol = 0.74 + ( 0.2 / ( 0.8 + rat ) );

    /* step 8 */
    Sr   = S + 0.002 / ( 0.6 + rat );
    zeta = expf( Sr * (wavel[idx440] - wavel[idx410]) );

    /* step 9 */
    denom  = zeta - symbol;
    dif1   = a[idx410]  - symbol * a[idx440];
    dif2   = aw[idx410] - symbol * aw[idx440];
    ag440  = (dif1 - dif2) / denom;

    for ( i = 0; i < nbands; i++ ) {
        adg[i] = ag440 * expf( Sr * (wavel[idx440] - wavel[i]));
    aph[i] = a[i] - adg[i] - aw[i];
    }

    /* check aph443 range */

    if ( aph_check ) {

        x1 = aph[idx440] / a[idx440];
        if ( x1 < 0.15 || x1 > 0.6 ) {

            if ( flags ) *flags |= 0x10;

            x2 = -0.8 + 1.4 * (a[idx440] - aw[idx440])/(a[idx410] - aw[idx410]);
            if ( x2 < 0.15 ) {
                x2 = 0.15;
                if ( flags ) *flags |= 0x20;
            }
            if ( x2 > 0.6 ) {
                x2 = 0.6;
                if ( flags ) *flags |= 0x40;
            }

            aph[idx440] = a[idx440] * x2;
            ag440 = a[idx440] - aph[idx440] - aw[idx440];

            for ( i = 0; i < nbands; i++ ) {
                adg[i] = ag440 * expf( Sr * (wavel[idx440] - wavel[i]));
                aph[i] = a[i] - adg[i] - aw[i];
            }

        }

    }

    return 0;
}

#ifdef TEST_QAA

#include <stdio.h>

/*
 * to compile:
 * cc -o qaa -g -DTEST_QAA -I. qaa.c -lm
 * ./qaa
 */



static void print_out( int n, float *fwl, float *Rrs, float *rrs, float *u,
           float *a, float *aph, float *adg, float *aw, float *bb, float *bbw )
{

    int i;

    printf("lamda ");
    for ( i = 0; i < n; i++ )
    printf("%9.0f ", fwl[i]);
    printf("\n");

    printf("Rrs : ");
    for ( i = 0; i < n; i++ )
    printf("%9.6f ", Rrs[i]);
    printf("\n");

    printf("rrs : ");
    for ( i = 0; i < n; i++ )
    printf("%9.6f ", rrs[i]);
    printf("\n");

    printf("u   : ");
    for ( i = 0; i < n; i++ )
    printf("%9.6f ", u[i]);
    printf("\n");

    printf("a   : ");
    for ( i = 0; i < n; i++ )
    printf("%9.6f ", a[i]);
    printf("\n");

    printf("aph : ");
    for ( i = 0; i < n; i++ )
    printf("%9.6f ", aph[i]);
    printf("\n");

    printf("adg : ");
    for ( i = 0; i < n; i++ )
    printf("%9.6f ", adg[i]);
    printf("\n");

    printf("aw  : ");
    for ( i = 0; i < n; i++ )
    printf("%9.6f ", aw[i]);
    printf("\n");

    printf("bb  : ");
    for ( i = 0; i < n; i++ )
    printf("%9.6f ", bb[i]);
    printf("\n");

    printf("bbw : ");
    for ( i = 0; i < n; i++ )
    printf("%9.6f ", bbw[i]);
    printf("\n");

}

int
main( int argc, char *argv[] )
{

    /* the 5th position here will be changed for 640nm later down in code */

#define NBANDS      6
#define NUM_SPECTRA 24
    float Rrs_insitu[NUM_SPECTRA][NBANDS] = {
    /*    412       443       490       510       555       670 */
        { 0.001919, 0.002297, 0.004420, 0.005547, 0.009138, 0.004110},
        { 0.001753, 0.002657, 0.004348, 0.005212, 0.007641, 0.003950},
        { 0.003968, 0.003374, 0.002932, 0.002142, 0.001214, 0.000136},
        { 0.001332, 0.001414, 0.002116, 0.002464, 0.003891, 0.001669},
        { 0.001572, 0.001537, 0.002132, 0.002504, 0.004157, 0.001909},
        { 0.004830, 0.004084, 0.003540, 0.002648, 0.001584, 0.000194},
        { 0.001145, 0.001457, 0.002501, 0.002945, 0.004180, 0.001524},
        { 0.000921, 0.001004, 0.001465, 0.001456, 0.001363, 0.000252},
        { 0.003120, 0.003226, 0.003530, 0.002857, 0.001882, 0.000196},
        { 0.005170, 0.004682, 0.003794, 0.002535, 0.001343, 0.000131},
        { 0.004718, 0.004378, 0.003763, 0.002627, 0.001472, 0.000151},
        { 0.003503, 0.003394, 0.003358, 0.002510, 0.001507, 0.000181},
        { 0.001005, 0.001131, 0.001879, 0.002219, 0.003316, 0.001100},
        { 0.007704, 0.006917, 0.005540, 0.003721, 0.002129, 0.000389},
        { 0.003311, 0.003071, 0.002920, 0.002315, 0.001508, 0.000285},
        { 0.003476, 0.003285, 0.003073, 0.002347, 0.001436, 0.000200},
        { 0.004661, 0.005739, 0.008028, 0.007809, 0.006632, 0.001035},
        { 0.004212, 0.004859, 0.006455, 0.006117, 0.004895, 0.000676},
        { 0.002090, 0.002280, 0.003091, 0.002780, 0.002177, 0.000444},
        { 0.003237, 0.003042, 0.002947, 0.002317, 0.001440, 0.000166},
        { 0.003125, 0.003444, 0.004119, 0.003777, 0.002875, 0.000413},
        { 0.002637, 0.002896, 0.003477, 0.002987, 0.002085, 0.000251},
        { 0.003472, 0.003596, 0.003834, 0.003067, 0.001969, 0.000331},
        { 0.004554, 0.003772, 0.002074, 0.002082, 0.001338, 0.000120}
    };

    float fwl[NBANDS] ={ 412,      443,      490,      510,      555,       670 };
    float aw[NBANDS]  ={ 0.004994, 0.007512, 0.025010, 0.040020, 0.077080, 0.445600 };
    float bbw[NBANDS] ={ 0.003271, 0.002421, 0.001568, 0.001339, 0.000939, 0.000426 };

    int   i, j;
    float Rrs[NBANDS];
    float rrs[NBANDS];
    float u[NBANDS];
    float a[NBANDS];
    float bb[NBANDS];
    float aph[NBANDS];
    float adg[NBANDS];
    unsigned char flags;

    int   idx410 = 0;
    int   idx440 = 1;
    int   idx490 = 2;
    int   idx555 = 4;
    int   idx670 = 5;
    int   nbands;

    // Ping's bbw using 0.0038 * pow((400.0/lambda),4.32);
//    for ( i = 0; i< NBANDS; i++ )
//        bbw[i] = 0.0038 * pow(400.0/fwl[i],4.32);

    qaa_init( idx410, idx440, idx490, idx555, idx670 );
    qaa_set_param( QAA_APH_CHECK, 0 );

    printf("QAA v6\n" );
    for ( j = 0; j < NUM_SPECTRA; j++ ) {

    flags  = 0;
        nbands = NBANDS;
    /* 412 to 670 */
    for ( i = 0; i < nbands; i++ )
        Rrs[i] = Rrs_insitu[j][i];

    qaaf_v6( nbands, fwl, Rrs, aw, bbw, rrs, u, a, bb, NULL );
    qaaf_decomp( nbands, fwl, rrs, a, aw, adg, aph, NULL );

    for ( i = 0; i < 6; i++ )
        if ( a[i] < aw[i] )
        flags |= 0x08;

    for ( i = 0; i < 6; i++ )
        if ( bb[i] < bbw[i] )
        flags |= 0x80;

    print_out( nbands, fwl, Rrs, rrs, u, a, aph, adg, aw, bb, bbw );

    if ( flags & 0x10 )
        printf("original aph/a ratio was out of range (0.15 to 0.6)\n");
    if ( flags & 0x20 )
        printf("    and was forced to minimum (0.15)\n");
    if ( flags & 0x40 )
        printf("    and was forced to maximum (0.6)\n");
    printf("\n");
    }

    return 0;
}
#endif