/*---------------------------------------------------------------------*/
/* ipar_arp.c -  functions to compute instantaneous PAR and ARP        */
/*                                                                     */
/* This is an implementation of MODIS ATBD 20, Version 7, by Carder,   */
/* Chen, and Hawes. Parts are based on the ipar-1.2.f code from modcol */
/* (the MODIS Level-2 code developed by Univ. of Miami) and parts of   */
/* that code are based on original work of W. Gregg.                   */
/*                                                                     */
/* References:                                                         */
/* Carder, K. L., Chen, R., Hawes, S., Instantaneous Photosynthetically*/
/* Available Radiation and Absorbed Radiation by Phytoplankton, MODIS  */
/* Ocean Science Team Algorithm Theoretical Basis Document (ATBD) 20,  */
/* Version 7, February 2003.                                           */
/*                                                                     */
/* Gregg, W.W., and Carder, K.L., A simple spectral soloar irradiance  */
/* model for cloudless maritime atmospheres, Limnol. Oceanogr., 35(8), */
/* 1657-1675, 1990.                                                    */
/*                                                                     */
/* Implementation for MSL12:  B. Franz, 21 Sep 2004.                   */
/*---------------------------------------------------------------------*/

#include "l12_proto.h"
#include "l2prod.h"

static float radeg = RADEG;
static float nw = 1.341;


/*---------------------------------------------------------------------*/
/* get_rhof - surface reflectance of foam from input windspeed         */
/*---------------------------------------------------------------------*/
float get_rhof(float ws)
{
    static float rhoair = 1.2E3;   /* density of air g/m3*/

    float Cd;
    float rhof;

    /*  Foam and diffuse reflectance */
    if (ws > 4.0) {
        if (ws <= 7.0) {
	    Cd  = 6.2E-4 + 1.56E-3/ws;
            rhof = rhoair*Cd*2.2E-5*ws*ws - 4.0E-4;
        } else {
	    Cd   = 0.49E-3 + 0.065E-3*ws;
            rhof = (rhoair*Cd*4.5E-5 - 4.0E-5)*ws*ws;
        }
    } else
        rhof   = 0.0;

    return(rhof);
}

/*---------------------------------------------------------------------*/
/* get_rhosps - surface reflectance for diffuse irrad. from windspeed  */
/*---------------------------------------------------------------------*/
float get_rhosps(float ws)
{
      if (ws > 4.0)
          return(0.057);
      else
          return(0.066);
}


/*---------------------------------------------------------------------*/
/* get_rhospd - surface reflectance for direct irrad. from input       */
/*              zenith angle (above surface) and windspeed.            */
/*---------------------------------------------------------------------*/
float get_rhospd(float theta, float ws)
{
    float rmin,rpls,sinp,tanp,a,b,rthetar,rtheta;
    float rhospd;

    /* Sub-surface angle via Snell's law */
    rtheta  = theta/radeg;
    rthetar = asin(sin(rtheta)/nw);

    /* Fresnel reflectance for theta < 40, ws < 2 m/s */
    if (theta < 40.0 || ws < 2.0) {
        if (theta == 0.0) {
	    rhospd = 0.0211;
        } else {
	    rmin = rtheta - rthetar;
	    rpls = rtheta + rthetar;
            sinp = (sin(rmin)*sin(rmin))/(sin(rpls)*sin(rpls));
            tanp = (tan(rmin)*tan(rmin))/(tan(rpls)*tan(rpls));
            rhospd = 0.5*(sinp + tanp);
        }
    } else {
        /* Empirical fit otherwise */
        a = 0.0253;
        b = -7.14E-4*ws + 0.0618;
        rhospd = a*exp(b*(theta-40.0));
    }

    return(rhospd);
}


/*---------------------------------------------------------------------*/
/* get_ipar - computes ipar for each pixel in L2 record                */
/*---------------------------------------------------------------------*/
void get_ipar(l2str *l2rec, float ipar[])
{
    static float badval = 0.0;
    static int   firstCall = 1;
    static int   iw400, iw700;
    static float wave[NBANDS];
    static float *wed;

    float *ws    = l2rec->ws;
    float *solz  = l2rec->solz;
    long  *bindx = l2rec->bindx;

    float rhof, rhospd, rho;
    float mu0, Ed0p, Ed0m;
    long  ip, ipb, iw, ib;

    if (firstCall) {
        firstCall = 0;
        for (iw=0; iw<l2rec->nbands; iw++)
	    wave[iw] = l2rec->fwave[l2rec->bindx[iw]];
        iw400 = windex(412.,wave,l2rec->nbands);
        iw700 = windex(670.,wave,l2rec->nbands);
        rdsensorinfo(l2rec->sensorID,"wed",(void **) &wed);
    }

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

        ipar[ip] = 0.;

        if (!l2rec->mask[ip]) {

            mu0 = cos(solz[ip]/radeg);

            /* diffuse reflectance for this geometry */
	    rhospd = get_rhospd(solz[ip],ws[ip]);
            rhof   = get_rhof(ws[ip]);
	    rho    = rhospd + rhof;

            for (iw=iw400; iw<=iw700; iw++) {

	        ib  = bindx[iw];
                ipb = ip*NBANDS+ib;

                Ed0p = l2rec->Fonom[ib]/100.    /* 100 mW/cm^2/um = 1 W/m^2/nm */
                     * l2rec->fsol
                     * l2rec->t_oz_sol[ipb]
                     * l2rec->t_sol   [ipb]
                     * mu0;

                Ed0m = Ed0p*(1-rho);

                /* 1 W/m2/nm = 119.3 uEin/m2/nm/sec ?????
                ipar[ip] += wave[iw]*Ed0m*wed[iw]/119.3;
		*/

                /* 1 W/m2 * 4.6 = 1 uEin/s/m^2     */
                /* wed is in nm, Ed is in W/m^2/nm */
                ipar[ip] += Ed0m*wed[iw]*4.6E-6;
	    }
	}
    }
}


/*---------------------------------------------------------------------*/
/* get_arp - computes arp for each pixel in L2 record                  */
/*---------------------------------------------------------------------*/
void get_arp(l2str *l2rec, float arp[])
{
    static float badval = 0.0;
    static int   firstCall = 1;
    static int   iw400, iw700;
    static float wave [NBANDS];
    static long  iwave[NBANDS];
    static float *wed;
    static float *waph;
    static float *aph675;
    static float *aph;
    static float *atot;
    static l2prodstr *p;

    float *ws      = l2rec->ws;
    float *solz    = l2rec->solz;
    float *senz    = l2rec->senz;
    long  *bindx   = l2rec->bindx;
    long  nband    = l2rec->nbands;
    long  npix     = l2rec->npix;
    long  nscan    = l2rec->nscans;
    long  sensorID = l2rec->sensorID;

    float rhof, rhosol, rhosen;
    float mu0, Ed0p, Ed0m;
    long  ip, ipb, iw, ib, ipw;
    char  name[32];
    float rrs[NBANDS];
    float solzp, senzp, deltaz;
    float aw685;
    float q, dterm1, uterm1, dterm2, uterm2;
    float avgcosd, avgcosu;
    float rkd, rku, rkexpd, rkexpu;

    if (firstCall) {
        firstCall = 0;
        aph675 = calloc(npix, sizeof(float));
        aph    = calloc(npix*nband, sizeof(float));
        atot   = calloc(npix*nband, sizeof(float));
        for (iw=0; iw<nband; iw++) {
	    wave [iw] = l2rec->fwave[bindx[iw]];
            iwave[iw] = (int) wave[iw];
            ib = l2rec->bindx[iw];
	}
        iw400 = windex(412.,wave,nband);
        iw700 = windex(670.,wave,nband);
        rdsensorinfo(l2rec->sensorID,"wed", (void **) &wed );
        rdsensorinfo(l2rec->sensorID,"waph",(void **) &waph);
    }

    /* get total and phytoplankton absorption using Carder model */
    p = get_l2prod_index("aph_carder",sensorID,nband,npix,nscan,bindx,iwave);
    get_carder(l2rec,p,aph675);
    for (iw=iw400; iw<=iw700; iw++) {
        sprintf(name,"aph_%3d_carder",iwave[iw]);
        p = get_l2prod_index(name,sensorID,nband,npix,nscan,bindx,iwave);
        get_carder(l2rec,p,&aph[iw*npix]);
        sprintf(name,"a_%3d_carder",iwave[iw]);
        p = get_l2prod_index(name,sensorID,nband,npix,nscan,bindx,iwave);
        get_carder(l2rec,p,&atot[iw*npix]);
    }

    for (ip=0; ip<npix; ip++) {

        arp[ip] = 0.;

        /* get Rrs, but skip pixel if Rrs is negative */
        for (iw=iw400; iw<=iw700; iw++) {
	    ib  = bindx[iw];
            ipb = ip*NBANDS+bindx[ib];
            rrs[iw] = l2rec->Rrs[ipb];
            if (rrs[iw] < 0.0) {
                l2rec->flags[ip] |= CHLFAIL;
                continue;
	    }
	}

        if (!l2rec->mask[ip] && aph675[ip] > 0.0) {

            mu0 = cos(solz[ip]/radeg);

            /* direct + diffuse reflectance for this geometry */
            rhof   = get_rhof(ws[ip]);
	    rhosol = get_rhospd(solz[ip],ws[ip]) + rhof;
	    rhosen = get_rhospd(senz[ip],ws[ip]) + rhof;

            /* first optical depth (deltaz) at 685 nm, satellite view */
            senzp  = asin(sin(senz[ip]/radeg)/nw);
            aw685  = 0.475;
            deltaz = cos(senzp)/(aw685 + aph675[ip]);

            /* avg. cosines for down/upwelling */
            solzp   = asin(sin(solz[ip]/radeg)/nw);
            avgcosd = 0.96*cos(solzp);
            avgcosu = 0.4;

	    /* down/upwelling terms for ARP integral */
            q = 4.0;
            dterm1 = 1.0/avgcosd;
            uterm1 = q*nw*nw/(1.0-rhosol)/(1.0-rhosen)/avgcosu;

            /* integrate over wavelengths */
            for (iw=iw400; iw<=iw700; iw++) {

	        ib  = bindx[iw];
                ipb = ip*NBANDS+ib;
                ipw = iw*npix+ip;

                Ed0p = l2rec->Fonom[ib]/100.
                     * l2rec->fsol
                     * l2rec->t_oz_sol[ipb]
                     * l2rec->t_sol   [ipb]
                     * mu0;

                Ed0m = Ed0p*(1-rhosol);

                rkd    = atot[ipw]/avgcosd;
                rku    = atot[ipw]/avgcosu;
 	        rkexpd = MIN(rkd*deltaz,60.);
  	        rkexpu = MIN(rku*deltaz,60.);
                dterm2 = dterm1*(1.0-exp(-rkexpd))/rkd;
                uterm2 = uterm1*(1.0-exp(-rkexpu))/rku*rrs[iw];

                arp[ip] += aph[ipw]*waph[iw]*Ed0m*wed[iw]*(dterm2 + uterm2)*4.6E-6;
	    }
	} else {
	    l2rec->flags[ip] |= CHLFAIL;
	}
    }
}