class SmacAlgorithm {

    public final float[] run(float[] sza, float[] saa, float[] vza, float[] vaa,
                             float[] taup550, float[] uh2o, float[] uo3,
                             float[] airPressure, boolean[] process, float invalid, float[] r_toa,
                             float[] r_surfRecycle) {
        // array to be returned
        float[] r_return;
        double us, invUs, us2, uv, invUv, usTimesuv, invUsTimesUv;
        double dphi, Peq, m, s, cksi, ksiD;
        double taup, tautot, Res_6s;
        double uo2, uco2, uch4, uno2, uco;
        double to3, th2o, to2, tco2, tch4, tno2, tco;
        double ttetas, ttetav;
        double ray_phase, ray_ref, taurz, Res_ray;
        double aer_phase, aer_ref, Res_aer;
        double atm_ref, tg;
        double d, del, dp, e, f, ss;
        double q1, q2, q3;
        double c1, c2, cp1, cp2;
        double x, y, z;
        double aa1, aa2, aa3;
        double temp;

        // try to reuse the recyle array to prevent memory waste. We can reuse if
        // a) it's present and
        // b) has the same size as the input vector
        if ((r_surfRecycle == null) || (r_surfRecycle.length != r_toa.length)) {
            r_return = new float[r_toa.length];
        } else {
            r_return = r_surfRecycle;
        }

        // loop over vectors
        // -----------------
        int n;
        for (n = 0; n < r_toa.length; n++) {
            // check for process flag. If set to false we must set the default value for
            // invalid pixels and process the next pixel
            if (!process[n]) {
                r_return[n] = invalid;
                continue;
            }
            // parameter setup
            us = Math.cos(sza[n] * _cdr);
            invUs = 1.0 / us;
            us2 = us * us;

            uv = Math.cos(vza[n] * _cdr);
            invUv = 1.0 / uv;
            usTimesuv = us * uv;
            invUsTimesUv = 1.0 / usTimesuv;

            dphi = (saa[n] - vaa[n]) * _cdr;
            Peq = airPressure[n] * _invMaxPressure;

            /*------ 1) air mass */
            m = invUs + invUv;

            /*------ 2) aerosol optical depth in the spectral band, taup  */
            taup = _A0taup + _A1taup * taup550[n];

            /*------ 3) gaseous transmissions (downward and upward paths)*/
            uo2 = 1.0;
            if (_calcUo2) {
                uo2 = Math.pow(Peq, _po2);
            }
            uco2 = 1.0;
            if (_calcUco2) {
                uco2 = Math.pow(Peq, _pco2);
            }
            uch4 = 1.0;
            if (_calcUch4) {
                uch4 = Math.pow(Peq, _pch4);
            }
            uno2 = 1.0;
            if (_calcUno2) {
                uno2 = Math.pow(Peq, _pno2);
            }
            uco = 1.0;
            if (_calcUco) {
                uco = Math.pow(Peq, _pco);
            }

            /*------ 4) if uh2o <= 0 and uo3 <= 0 no gaseous absorption is computed*/
            to3 = 1.0;
            th2o = 1.0;
            to2 = 1.0;
            tco2 = 1.0;
            tch4 = 1.0;
            tno2 = 1.0;
            tco = 1.0;
            if ((uh2o[n] > 0.) || (uo3[n] > 0.)) {
                if (_calcTo3) {
                    to3 = Math.exp(_ao3 * Math.pow((uo3[n] * m), _no3));
                }
                if (_calcTh2o) {
                    th2o = Math.exp(_ah2o * Math.pow((uh2o[n] * m), _nh2o));
                }
                if (_calcTo2) {
                    to2 = Math.exp(_ao2 * Math.pow((uo2 * m), _no2));
                }
                if (_calcTco2) {
                    tco2 = Math.exp(_aco2 * Math.pow((uco2 * m), _nco2));
                }
                if (_calcTch4) {
                    tch4 = Math.exp(_ach4 * Math.pow((uch4 * m), _nch4));
                }
                if (_calcTno2) {
                    tno2 = Math.exp(_ano2 * Math.pow((uno2 * m), _nno2));
                }
                if (_calcTco) {
                    tco = Math.exp(_aco * Math.pow((uco * m), _nco));
                }
            }

            /*------  5) Total scattering transmission */
            temp = _a2T * Peq + _a3T;
            /* downward */
            ttetas = _a0T + _a1T * taup550[n] * invUs + temp / (1.0 + us);
            /* upward   */
            ttetav = _a0T + _a1T * taup550[n] * invUv + temp / (1.0 + uv);

            /*------ 6) spherical albedo of the atmosphere */
            s = _a0s * Peq + _a3s + _a1s * taup550[n] + _a2s * taup550[n] * taup550[n];

            /*------ 7) scattering angle cosine */
            cksi = -(usTimesuv + (Math.sqrt(1.0 - us2) * Math.sqrt(1.0 - uv * uv) * Math.cos(dphi)));
            if (cksi < -1) {
                cksi = -1.0;
            }

            /*------ 8) scattering angle in degree */
            ksiD = _crd * Math.acos(cksi);

            /*------ 9) rayleigh atmospheric reflectance */
            /* pour 6s on a delta = 0.0279 */
            ray_phase = 0.7190443 * (1.0 + (cksi * cksi)) + 0.0412742;
            taurz = _taur * Peq;
            ray_ref = (taurz * ray_phase) * _oneQuarter * invUsTimesUv;

            /*-----------------Residu Rayleigh ---------*/
            temp = taurz * ray_phase * invUsTimesUv;
            Res_ray = _resr1 + _resr2 * temp + _resr3 * temp * temp;

            /*------ 10) aerosol atmospheric reflectance */
            temp = ksiD * ksiD;
            aer_phase = _a0P + _a1P * ksiD + _a2P * temp + _a3P * temp * ksiD + _a4P * temp * temp;

            // now the uncommented block :-)
            // -----------------------------
            temp = 1.0 / (4.0 * (1.0 - _ak2 * us2));
            e = -3.0 * us2 * _wo * temp;
            f = -_onemwo * 3.0 * _gc * us2 * _wo * temp;
            dp = e / (3.0 * us) + us * f;
            d = e + f;
            del = Math.exp(_ak * taup) * _onepb2 - Math.exp(-_ak * taup) * _onemb2;
            ss = us / (1.0 - _ak2 * us2);
            temp = 3.0 * us;
            q1 = 2.0 + temp + _onemwo * temp * _gc * (1.0 + 2.0 * us);
            q2 = 2.0 - temp - _onemwo * temp * _gc * (1.0 - 2.0 * us);
            q3 = q2 * Math.exp(-taup * invUs);
            temp = (_ww * ss) / del;
            c1 = temp * (q1 * Math.exp(_ak * taup) * _onepb + q3 * _onemb);
            c2 = -temp * (q1 * Math.exp(-_ak * taup) * _onemb + q3 * _onepb);
            cp1 = c1 * _pfac;
            cp2 = -c2 * _pfac;
            temp = _wo * 3.0 * _gc * uv;
            z = d - temp * dp + _wo * aer_phase * _oneQuarter;
            x = c1 - temp * cp1;
            y = c2 - temp * cp2;
            temp = _ak * uv;
            aa1 = uv / (1.0 + temp);
            aa2 = uv / (1.0 - temp);
            aa3 = usTimesuv / (us + uv);

            aer_ref = x * aa1 * (1.0 - Math.exp(-taup / aa1));
            aer_ref += y * aa2 * (1.0 - Math.exp(-taup / aa2));
            aer_ref += z * aa3 * (1.0 - Math.exp(-taup / aa3));
            aer_ref *= invUsTimesUv;

            /*--------Residu Aerosol --------*/
            temp = taup * m * cksi;
            Res_aer = (_resa1 + _resa2 * temp + _resa3 * temp * temp)
                    + _resa4 * temp * temp * temp;

            /*---------Residu 6s-----------*/
            tautot = taup + taurz;
            temp = tautot * m * cksi;
            Res_6s = (_rest1 + _rest2 * temp + _rest3 * temp * temp)
                    + _rest4 * temp * temp * temp;

            /*------ 11) total atmospheric reflectance */
            atm_ref = ray_ref - Res_ray + aer_ref - Res_aer + Res_6s;

            /*-------- reflectance at toa*/
            tg = th2o * to3 * to2 * tco2 * tch4 * tco * tno2;

            /* reflectance at surface */
            /*------------------------ */
            temp = r_toa[n] - (atm_ref * tg);
            temp = temp / ((tg * ttetas * ttetav) + (temp * s));
            r_return[n] = (float) temp;
        }

        return r_return;
    }

}