class RayleighCorrAlgorithm {

    double[] getRhoBrr(RayleighAux rayleighAux, double[] rayleighOpticalThickness, double[] corrOzoneRefl) {
        final double[] airMasses = rayleighAux.getAirMass();
        final double[] aziDiffs = rayleighAux.getAziDifferent();
        final double[] cosSZARads = rayleighAux.getCosSZARads();
        final double[] cosOZARads = rayleighAux.getCosOZARads();
        final double[] sARay = rayleighAux.getInterpolateRayleighThickness(rayleighOpticalThickness);
        final double[] tau_ray = rayleighAux.getTaur();

        final Map<Integer, double[]> fourier = rayleighAux.getFourier();
        final Map<Integer, List<double[]>> interpolation = rayleighAux.getInterpolation();
        final int length = cosOZARads.length;

        final double[] rho_BRR = new double[length];

        for (int index = 0; index < length; index++) {

            double corrOzone = corrOzoneRefl[index];
            if (Double.isNaN(corrOzone)) {
                rho_BRR[index] = RayleighConstants.INVALID_VALUE;
                continue;
            }
            if (corrOzone <= 0) {
                rho_BRR[index] = 0.0;
                continue;
            }
            double taurVal = rayleighOpticalThickness[index];
            if (Double.isNaN(taurVal)) {
                rho_BRR[index] = taurVal;
                continue;
            }
            double aziDiff = aziDiffs[index];
            double massAir = airMasses[index];

            double cosOZARad = cosOZARads[index];
            double cosSZARad = cosSZARads[index];

            List<double[]> interpolateValues = interpolation.get(index);
            if (interpolateValues == null || interpolateValues.size() == 0) {
                // this might happen if we are out of range of RayleighAux data (e.g. SZA > 80deg)
                // --> set to NaN (CB, 20171026)
                rho_BRR[index] = RayleighConstants.INVALID_VALUE;
                continue;
            }
            double[] fourierSeries = fourier.get(index);

            double[] rho_Rm = getFourierSeries(taurVal, massAir, cosOZARad, cosSZARad, interpolateValues, fourierSeries);

            double rho_R = rho_Rm[0] + 2.0 * rho_Rm[1] * Math.cos(aziDiff) + 2.0 * rho_Rm[2] * Math.cos(2.0 * aziDiff);

            // polynomial coefficients tpoly0, tpoly1 and tpoly2 from MERIS LUT

            double tRs = ((2.0 / 3.0 + cosSZARad) + (2.0 / 3.0 - cosSZARad) * Math.exp(-taurVal / cosSZARad)) / (4.0 / 3.0 + taurVal);

            double tR_thetaS = tau_ray[0] + tau_ray[1] * tRs + tau_ray[2] * Math.pow(tRs, 2);
            //#Rayleigh Transmittance sun - surface
            double tRv = ((2.0 / 3.0 + cosOZARad) + (2.0 / 3.0 - cosOZARad) * Math.exp(-taurVal / cosOZARad)) / (4.0 / 3.0 + taurVal);
            //#Rayleigh Transmittance surface - sensor
            double tR_thetaV = tau_ray[0] + tau_ray[1] * tRv + tau_ray[2] * Math.pow(tRv, 2);

            double saRay = sARay[index];


            double rho_toaR = (corrOzone - rho_R) / (tR_thetaS * tR_thetaV); //toa corrOzoneRefl corrected for Rayleigh scattering
            double sphericalFactor = 1.0 / (1.0 + saRay * rho_toaR); //#factor used in the next equation to account for the spherical albedo
            //#top of aerosol reflectance, which is equal to bottom of Rayleigh reflectance
            rho_BRR[index] = (rho_toaR * sphericalFactor);
        }
        return rho_BRR;
    }

    public  RayleighCorrAlgorithm(String bandPattern, int numBands);
    public  RayleighCorrAlgorithm(Sensor sensor);

}