/*
 * JGrass - Free Open Source Java GIS http://www.jgrass.org 
 * (C) HydroloGIS - www.hydrologis.com 
 * 
 * This library is free software; you can redistribute it and/or modify it under
 * the terms of the GNU Library General Public License as published by the Free
 * Software Foundation; either version 2 of the License, or (at your option) any
 * later version.
 * 
 * This library is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more
 * details.
 * 
 * You should have received a copy of the GNU Library General Public License
 * along with this library; if not, write to the Free Foundation, Inc., 59
 * Temple Place, Suite 330, Boston, MA 02111-1307 USA
 */
package org.jgrasstools.hortonmachine.modules.hydrogeomorphology.energyindexcalculator;

import static java.lang.Math.PI;
import static java.lang.Math.acos;
import static java.lang.Math.asin;
import static java.lang.Math.cos;
import static java.lang.Math.sin;
import static java.lang.Math.tan;
import static org.jgrasstools.gears.libs.modules.JGTConstants.intNovalue;
import static org.jgrasstools.gears.libs.modules.JGTConstants.isNovalue;
import static org.jgrasstools.gears.libs.modules.JGTConstants.omega;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_AUTHORCONTACTS;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_AUTHORNAMES;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_DESCRIPTION;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_KEYWORDS;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_LABEL;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_LICENSE;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_NAME;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_STATUS;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_inAspect_DESCRIPTION;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_inBasins_DESCRIPTION;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_inCurvatures_DESCRIPTION;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_inElev_DESCRIPTION;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_inSlope_DESCRIPTION;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_outAltimetry_DESCRIPTION;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_outArea_DESCRIPTION;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_outEnergy_DESCRIPTION;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_pDt_DESCRIPTION;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_pEi_DESCRIPTION;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSENERGYINDEXCALCULATOR_pEs_DESCRIPTION;

import java.awt.image.RenderedImage;
import java.awt.image.WritableRaster;
import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;

import javax.media.jai.iterator.RandomIter;
import javax.media.jai.iterator.RandomIterFactory;

import oms3.annotations.Author;
import oms3.annotations.Description;
import oms3.annotations.Execute;
import oms3.annotations.In;
import oms3.annotations.Keywords;
import oms3.annotations.Label;
import oms3.annotations.License;
import oms3.annotations.Name;
import oms3.annotations.Out;
import oms3.annotations.Status;

import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.geometry.jts.JTS;
import org.geotools.referencing.CRS;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.jgrasstools.gears.io.eicalculator.EIAltimetry;
import org.jgrasstools.gears.io.eicalculator.EIAreas;
import org.jgrasstools.gears.io.eicalculator.EIEnergy;
import org.jgrasstools.gears.libs.exceptions.ModelsIllegalargumentException;
import org.jgrasstools.gears.libs.modules.JGTModel;
import org.jgrasstools.gears.utils.coverage.CoverageUtilities;
import org.jgrasstools.hortonmachine.i18n.HortonMessageHandler;
import org.opengis.referencing.operation.MathTransform;

import com.vividsolutions.jts.geom.Coordinate;

@Description(OMSENERGYINDEXCALCULATOR_DESCRIPTION)
@Author(name = OMSENERGYINDEXCALCULATOR_AUTHORNAMES, contact = OMSENERGYINDEXCALCULATOR_AUTHORCONTACTS)
@Keywords(OMSENERGYINDEXCALCULATOR_KEYWORDS)
@Label(OMSENERGYINDEXCALCULATOR_LABEL)
@Name(OMSENERGYINDEXCALCULATOR_NAME)
@Status(OMSENERGYINDEXCALCULATOR_STATUS)
@License(OMSENERGYINDEXCALCULATOR_LICENSE)
public class OmsEnergyIndexCalculator extends JGTModel {

    @Description(OMSENERGYINDEXCALCULATOR_inElev_DESCRIPTION)
    @In
    public GridCoverage2D inElev = null;

    @Description(OMSENERGYINDEXCALCULATOR_inBasins_DESCRIPTION)
    @In
    public GridCoverage2D inBasins = null;

    @Description(OMSENERGYINDEXCALCULATOR_inCurvatures_DESCRIPTION)
    @In
    public GridCoverage2D inCurvatures = null;

    @Description(OMSENERGYINDEXCALCULATOR_inAspect_DESCRIPTION)
    @In
    public GridCoverage2D inAspect = null;

    @Description(OMSENERGYINDEXCALCULATOR_inSlope_DESCRIPTION)
    @In
    public GridCoverage2D inSlope = null;

    @Description(OMSENERGYINDEXCALCULATOR_pEs_DESCRIPTION)
    @In
    public int pEs = -1;

    @Description(OMSENERGYINDEXCALCULATOR_pEi_DESCRIPTION)
    @In
    public int pEi = -1;

    @Description(OMSENERGYINDEXCALCULATOR_pDt_DESCRIPTION)
    @In
    public double pDt = -1;

    @Description(OMSENERGYINDEXCALCULATOR_outAltimetry_DESCRIPTION)
    @Out
    public List<EIAltimetry> outAltimetry;

    @Description(OMSENERGYINDEXCALCULATOR_outEnergy_DESCRIPTION)
    @Out
    public List<EIEnergy> outEnergy;

    @Description(OMSENERGYINDEXCALCULATOR_outArea_DESCRIPTION)
    @Out
    public List<EIAreas> outArea;

    private static final int NOVALUE = intNovalue;
    private HashMap<Integer, Integer> id2indexMap = null;
    private HashMap<Integer, Integer> index2idMap = null;

    private double avgLatitude = -1;
    private double dx;
    private double dy;
    private int[][] eibasinID;
    private int[][] outputShadow;
    private double[][][] eibasinEmonth;
    private double[][] eibasinESrange;
    private double[][] eibasinES;
    private double[][] eibasinEI_mean;
    private double[][][] eibasinEI;
    private double[][] eibasinE;
    private double[][][] eibasinA;
    private final GeomorphUtilities geomorphUtilities = new GeomorphUtilities();
    private int eibasinNum;

    private RandomIter idbasinImageIterator;

    private RandomIter elevImageIterator;

    private WritableRaster curvatureImage;

    private RandomIter aspectImageIterator;

    private RandomIter slopeImageIterator;

    private int rows;

    private int cols;

    private HortonMessageHandler msg = HortonMessageHandler.getInstance();

    @Execute
    public void process() throws Exception {
        HashMap<String, Double> regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inBasins);
        cols = regionMap.get(CoverageUtilities.COLS).intValue();
        rows = regionMap.get(CoverageUtilities.ROWS).intValue();
        dx = regionMap.get(CoverageUtilities.XRES);
        dy = regionMap.get(CoverageUtilities.YRES);
        double n = regionMap.get(CoverageUtilities.NORTH);
        double s = regionMap.get(CoverageUtilities.SOUTH);
        double w = regionMap.get(CoverageUtilities.WEST);
        double e = regionMap.get(CoverageUtilities.EAST);

        double meanX = w + (e - w) / 2.0;
        double meanY = s + (n - s) / 2.0;

        Coordinate tmp = new Coordinate(meanX, meanY);
        MathTransform mathTransform = CRS.findMathTransform(inAspect.getCoordinateReferenceSystem(), DefaultGeographicCRS.WGS84);
        Coordinate newC = JTS.transform(tmp, null, mathTransform);
        avgLatitude = newC.y;

        RenderedImage idbasinImage = inBasins.getRenderedImage();
        idbasinImageIterator = RandomIterFactory.create(idbasinImage, null);

        RenderedImage elevImage = inElev.getRenderedImage();
        elevImageIterator = RandomIterFactory.create(elevImage, null);

        RenderedImage tmpImage = inCurvatures.getRenderedImage();
        curvatureImage = CoverageUtilities
                .createDoubleWritableRaster(tmpImage.getWidth(), tmpImage.getHeight(), null, null, null);
        RandomIter tmpIterator = RandomIterFactory.create(tmpImage, null);
        // TODO check what this is for?!?!?
        for( int i = 0; i < tmpImage.getHeight(); i++ ) {
            for( int j = 0; j < tmpImage.getWidth(); j++ ) {
                double value = tmpIterator.getSampleDouble(j, i, 0);
                curvatureImage.setSample(j, i, 0, value);
            }
        }

        RenderedImage aspectImage = inAspect.getRenderedImage();
        aspectImageIterator = RandomIterFactory.create(aspectImage, null);

        RenderedImage slopeImage = inSlope.getRenderedImage();
        slopeImageIterator = RandomIterFactory.create(slopeImage, null);

        avgLatitude *= (PI / 180.0);

        pm.message(msg.message("eicalculator.preparing_inputs")); //$NON-NLS-1$
        eibasinNum = prepareInputsOutputs();

        pm.beginTask(msg.message("eicalculator.computing"), 6); //$NON-NLS-1$
        for( int m = 0; m < 6; m++ ) {
            pm.worked(1);
            compute_EI(m + 1);
        }
        pm.done();

        average_EI(10, 6);

        pm.beginTask(msg.message("eicalculator.calc_areas"), eibasinNum); //$NON-NLS-1$
        for( int i = 0; i < eibasinNum; i++ ) {
            pm.worked(1);
            area(i);
        }
        pm.done();

        /*
         * putting the results together
         */
        outAltimetry = new ArrayList<EIAltimetry>();
        outEnergy = new ArrayList<EIEnergy>();
        outArea = new ArrayList<EIAreas>();
        for( int i = 0; i < eibasinNum; i++ ) {
            int realBasinId = index2idMap.get(i + 1);
            /*
             * ENERGY BANDS
             * 
             * Cycle over the virtual months:
             * 0: 22 DICEMBRE - 20 GENNAIO
             * 1: 21 GENNAIO - 20 FEBBRAIO
             * 2: 21 FEBBRAIO - 22 MARZO
             * 3: 23 MARZO - 22 APRILE
             * 4: 23 APRILE - 22 MAGGIO
             * 5: 23 MAGGIO - 22 GIUGNO
             */
            for( int j = 0; j < 6; j++ ) {
                for( int k = 0; k < pEi; k++ ) {
                    EIEnergy tmpEi = new EIEnergy();
                    // the basin id
                    tmpEi.basinId = realBasinId;
                    tmpEi.energeticBandId = k;
                    tmpEi.virtualMonth = j;
                    tmpEi.energyValue = eibasinEI[0][k][i];
                    outEnergy.add(tmpEi);
                }
            }

            /*
             * ALTIMETRIC BANDS
             */
            for( int k = 0; k < pEs; k++ ) {
                EIAltimetry tmpAl = new EIAltimetry();
                tmpAl.basinId = realBasinId;
                tmpAl.altimetricBandId = k;
                tmpAl.elevationValue = eibasinES[k][i];
                tmpAl.bandRange = eibasinESrange[k][i];
                outAltimetry.add(tmpAl);
            }

            /*
             * AREAS
             */
            for( int j = 0; j < pEs; j++ ) {
                for( int k = 0; k < pEi; k++ ) {
                    EIAreas tmpAr = new EIAreas();
                    tmpAr.basinId = realBasinId;
                    tmpAr.altimetricBandId = j;
                    tmpAr.energyBandId = k;
                    tmpAr.areaValue = eibasinA[j][k][i];
                    outArea.add(tmpAr);
                }
            }

        }
    }

    private int prepareInputsOutputs() {

        List<Integer> idList = new ArrayList<Integer>();
        eibasinID = new int[rows][cols];

        for( int r = 0; r < rows; r++ ) {
            for( int c = 0; c < cols; c++ ) {
                // get the value
                if (isNovalue(idbasinImageIterator.getSampleDouble(c, r, 0))) {
                    eibasinID[r][c] = NOVALUE;
                } else {
                    eibasinID[r][c] = (int) idbasinImageIterator.getSampleDouble(c, r, 0);
                    // put the value in the id list, if it isn't already there
                    if (!idList.contains(eibasinID[r][c])) {
                        idList.add(eibasinID[r][c]);
                    }
                }

            }
        }
        // sort the id list
        Collections.sort(idList);
        /*
         * now the number of involved subbasins is known
         */
        int eibasinNum = idList.size();

        /*
         * now substitute the numbers in the ID matrix with a sequential index without wholes
         */
        // first create the mapping
        id2indexMap = new HashMap<Integer, Integer>();
        index2idMap = new HashMap<Integer, Integer>();
        for( int i = 1; i <= idList.size(); i++ ) {
            id2indexMap.put(idList.get(i - 1), i);
            index2idMap.put(i, idList.get(i - 1));
        }
        for( int r = 0; r < eibasinID.length; r++ ) {
            for( int c = 0; c < eibasinID[0].length; c++ ) {
                if (eibasinID[r][c] != NOVALUE) {
                    eibasinID[r][c] = id2indexMap.get(eibasinID[r][c]);
                }
            }
        }

        pm.message(msg.message("eicalculator.subbasinsnum") + eibasinNum); //$NON-NLS-1$

        /*
         * prepare outputs
         */
        outputShadow = new int[rows][cols];
        for( int r = 0; r < outputShadow.length; r++ ) {
            for( int c = 0; c < outputShadow[0].length; c++ ) {
                outputShadow[r][c] = NOVALUE;
            }
        }

        eibasinE = new double[rows][cols];
        for( int r = 0; r < eibasinE.length; r++ ) {
            for( int c = 0; c < eibasinE[0].length; c++ ) {
                eibasinE[r][c] = NOVALUE;
            }
        }

        eibasinEmonth = new double[6][rows][cols];
        for( int r = 0; r < eibasinEmonth.length; r++ ) {
            for( int c = 0; c < eibasinEmonth[0].length; c++ ) {
                for( int t = 0; t < eibasinEmonth[0][0].length; t++ ) {
                    eibasinEmonth[r][c][t] = NOVALUE;
                }
            }
        }

        eibasinES = new double[pEs][eibasinNum];
        eibasinESrange = new double[pEs][eibasinNum];
        eibasinEI_mean = new double[pEi][eibasinNum];

        eibasinEI = new double[6][pEi][eibasinNum];

        eibasinA = new double[pEs][pEi][eibasinNum];

        return eibasinNum;
    }

    private void compute_EI( int month ) {

        int[] day_beg = new int[1], day_end = new int[1], daymonth = new int[1], monthyear = new int[1];
        int day;
        double hour;
        double[] Rad_morpho = new double[1], Rad_flat = new double[1];
        double[] E0 = new double[1], alpha = new double[1], direction = new double[1];
        double[][] Rad_morpho_cum, Rad_flat_cum;

        find_days(month, day_beg, day_end);

        hour = 0.5 * pDt;
        day = day_beg[0];

        Rad_morpho_cum = new double[rows][cols];
        Rad_flat_cum = new double[rows][cols];

        get_date(day, monthyear, daymonth);

        printReport(daymonth, monthyear, hour);

        do {
            sun(hour, day, E0, alpha, direction);

            for( int r = 0; r < eibasinID.length; r++ ) {
                for( int c = 0; c < eibasinID[0].length; c++ ) {
                    if (eibasinID[r][c] != NOVALUE) {
                        radiation(Rad_morpho, Rad_flat, E0[0], alpha[0], direction[0],
                                aspectImageIterator.getSampleDouble(c, r, 0), slopeImageIterator.getSampleDouble(c, r, 0),
                                outputShadow[r][c]);
                        Rad_morpho_cum[r][c] += Rad_morpho[0];
                        Rad_flat_cum[r][c] += Rad_flat[0];
                    }
                }
            }

            hour += pDt;
            if (hour >= 24) {
                hour -= 24.0;
                day += 1;
            }

            get_date(day, monthyear, daymonth);

            printReport(daymonth, monthyear, hour);

        } while( day <= day_end[0] );

        for( int r = 0; r < eibasinID.length; r++ ) {
            for( int c = 0; c < eibasinID[0].length; c++ ) {
                if (eibasinID[r][c] != NOVALUE) {
                    pm.message("Bacino: " + index2idMap.get(eibasinID[r][c]));
                    if (Rad_flat_cum[r][c] == 0) {
                        pm.message("Rad flat nulla");
                        Rad_morpho_cum[r][c] = 1;
                        Rad_flat_cum[r][c] = 1;
                    }
                    eibasinEmonth[month - 1][r][c] = Rad_morpho_cum[r][c] / Rad_flat_cum[r][c];
                    pm.message("Rad morfo: " + Rad_morpho_cum[r][c]);
                    pm.message("Rad flat: " + Rad_flat_cum[r][c]);
                } else {
                    eibasinEmonth[month - 1][r][c] = NOVALUE;
                }
            }
        }

    }

    private void printReport( int[] daymonth, int[] monthyear, double hour ) {
        StringBuilder sb = new StringBuilder();
        sb.append(msg.message("hm.day")); //$NON-NLS-1$
        sb.append(": "); //$NON-NLS-1$
        sb.append(daymonth[0]);
        sb.append("/"); //$NON-NLS-1$
        sb.append(monthyear[0]);
        sb.append(msg.message("hm.hour")); //$NON-NLS-1$
        sb.append(": "); //$NON-NLS-1$
        sb.append(hour);
        if ((hour - (long) hour) * 60 < 10) {
            sb.append(":0"); //$NON-NLS-1$
        } else {
            sb.append(":"); //$NON-NLS-1$
        }
        sb.append(((hour - (long) hour) * 60));
        pm.message(sb.toString());
    }

    private void find_days( int month, int[] day_begin, int[] day_end ) {

        if (month == 1) {
            day_begin[0] = -9;
            day_end[0] = 20;
        } else if (month == 2) {
            day_begin[0] = 21;
            day_end[0] = 51;
        } else if (month == 3) {
            day_begin[0] = 52;
            day_end[0] = 81;
        } else if (month == 4) {
            day_begin[0] = 82;
            day_end[0] = 112;
        } else if (month == 5) {
            day_begin[0] = 113;
            day_end[0] = 142;
        } else if (month == 6) {
            day_begin[0] = 143;
            day_end[0] = 173;
        } else {
            throw new ModelsIllegalargumentException("Incorrect in find_days", "OmsEnergyIndexCalculator", pm);
        }
    }

    private void get_date( int julianday, int[] month, int[] daymonth ) {

        if (julianday <= 0) {
            month[0] = 12;
            daymonth[0] = 31 + julianday;
        } else if (julianday >= 1 && julianday <= 31) {
            month[0] = 1;
            daymonth[0] = julianday;
        } else if (julianday >= 32 && julianday <= 59) {
            month[0] = 2;
            daymonth[0] = julianday - 31;
        } else if (julianday >= 60 && julianday <= 90) {
            month[0] = 3;
            daymonth[0] = julianday - 59;
        } else if (julianday >= 91 && julianday <= 120) {
            month[0] = 4;
            daymonth[0] = julianday - 90;
        } else if (julianday >= 121 && julianday <= 151) {
            month[0] = 5;
            daymonth[0] = julianday - 120;
        } else if (julianday >= 152 && julianday <= 181) {
            month[0] = 6;
            daymonth[0] = julianday - 151;
        } else if (julianday >= 182 && julianday <= 212) {
            month[0] = 7;
            daymonth[0] = julianday - 181;
        } else if (julianday >= 213 && julianday <= 243) {
            month[0] = 8;
            daymonth[0] = julianday - 212;
        } else if (julianday >= 244 && julianday <= 273) {
            month[0] = 9;
            daymonth[0] = julianday - 243;
        } else if (julianday >= 274 && julianday <= 304) {
            month[0] = 10;
            daymonth[0] = julianday - 273;
        } else if (julianday >= 305 && julianday <= 334) {
            month[0] = 11;
            daymonth[0] = julianday - 304;
        } else if (julianday >= 335 && julianday <= 365) {
            month[0] = 12;
            daymonth[0] = julianday - 334;
        }

    }

    private void sun( double hour, int day, double[] E0, double[] alpha, double[] direction ) {

        // latitudine, longitudine in [rad]

        double G, Et, local_hour, D, Thr, beta;

        // correction sideral time
        G = 2.0 * PI * (day - 1) / 365.0;
        Et = 0.000075 + 0.001868 * cos(G) - 0.032077 * sin(G) - 0.014615 * cos(2 * G) - 0.04089 * sin(2 * G);

        // local time
        local_hour = hour + Et / omega; // Iqbal: formula 1.4.2

        // earth-sun distance correction
        E0[0] = 1.00011 + 0.034221 * cos(G) + 0.00128 * sin(G) + 0.000719 * cos(2 * G) + 0.000077 * sin(2 * G);

        // solar declination
        D = 0.006918 - 0.399912 * cos(G) + 0.070257 * sin(G) - 0.006758 * cos(2 * G) + 0.000907 * sin(2 * G) - 0.002697
                * cos(3 * G) + 0.00148 * sin(3 * G);

        // Sunrise and sunset with respect to midday [hour]
        Thr = (acos(-tan(D) * tan(avgLatitude))) / omega;

        if (local_hour >= 12.0 - Thr && local_hour <= 12.0 + Thr) {

            // alpha: solar height (complementar to zenith angle), [rad]
            alpha[0] = asin(sin(avgLatitude) * sin(D) + cos(avgLatitude) * cos(D) * cos(omega * (12.0 - local_hour)));

            // direction: azimuth angle (0 Nord, clockwise) [rad]
            if (local_hour <= 12) {
                if (alpha[0] == PI / 2.0) { /* sole allo zenit */
                    direction[0] = PI / 2.0;
                } else {
                    direction[0] = PI - acos((sin(alpha[0]) * sin(avgLatitude) - sin(D)) / (cos(alpha[0]) * cos(avgLatitude)));
                }
            } else {
                if (alpha[0] == PI / 2.0) { /* sole allo zenit */
                    direction[0] = 3 * PI / 2.0;
                } else {
                    direction[0] = PI + acos((sin(alpha[0]) * sin(avgLatitude) - sin(D)) / (cos(alpha[0]) * cos(avgLatitude)));
                }
            }

            // CALCOLO OMBRE
            /*
             * Chiama Orizzonte# Inputs: dx: dim. pixel (funziona solo per pixel quadrati)
             * 2(basin.Z.nch + basin.Z.nrh): dimensione matrice alpha: altezza solare Z: matrice
             * elevazioni curv: matrice curvature beta: azimuth +#PI/4 NOVALUE: novalue per Z0
             * Outputs: shadow: matrice ombre (1 ombra 0 sole)
             */

            if (direction[0] >= 0. && direction[0] <= PI / 4.) {
                beta = direction[0];
                geomorphUtilities.orizzonte1(dx, 2 * (cols + rows), beta, alpha[0], elevImageIterator, curvatureImage,
                        outputShadow);

            } else if (direction[0] > PI / 4. && direction[0] <= PI / 2.) {
                beta = (PI / 2. - direction[0]);
                geomorphUtilities.orizzonte2(dx, 2 * (cols + rows), beta, alpha[0], elevImageIterator, curvatureImage,
                        outputShadow);

            } else if (direction[0] > PI / 2. && direction[0] <= PI * 3. / 4.) {
                beta = (direction[0] - PI / 2.);

                geomorphUtilities.orizzonte3(dx, 2 * (cols + rows), beta, alpha[0], elevImageIterator, curvatureImage,
                        outputShadow);

            } else if (direction[0] > PI * 3. / 4. && direction[0] <= PI) {
                beta = (PI - direction[0]);
                geomorphUtilities.orizzonte4(dx, 2 * (cols + rows), beta, alpha[0], elevImageIterator, curvatureImage,
                        outputShadow);

            } else if (direction[0] > PI && direction[0] <= PI * 5. / 4.) {
                beta = (direction[0] - PI);
                geomorphUtilities.orizzonte5(dx, 2 * (cols + rows), beta, alpha[0], elevImageIterator, curvatureImage,
                        outputShadow);

            } else if (direction[0] > PI * 5. / 4. && direction[0] <= PI * 3. / 2.) {
                beta = (PI * 3. / 2. - direction[0]);
                geomorphUtilities.orizzonte6(dx, 2 * (cols + rows), beta, alpha[0], elevImageIterator, curvatureImage,
                        outputShadow);

            } else if (direction[0] > PI * 3. / 2. && direction[0] <= PI * 7. / 4.) {
                beta = (direction[0] - PI * 3. / 2.);
                geomorphUtilities.orizzonte7(dx, 2 * (cols + rows), beta, alpha[0], elevImageIterator, curvatureImage,
                        outputShadow);

            } else if (direction[0] > PI * 7. / 4. && direction[0] < 2. * PI) {
                beta = (2. * PI - direction[0]);
                /*
                 * here we should have Orizzonte8, but the routine has an error. So the 1 is called.
                 * Explanation: quello è un errore dovuto al fatto che la routine orizzonte8 è
                 * sbagliata e dà errore, allora ci ho messo una pezza e ho richiamato la
                 * orizzonte1, invece che la orizzonte 8. tuttavia le orizzonte 1 e 8 vengono
                 * chiamate solo quando il sole è a nord e da noi questo non capita mai.. quindi
                 * puoi lasciare così com'è
                 */
                geomorphUtilities.orizzonte1(dx, 2 * (cols + rows), beta, alpha[0], elevImageIterator, curvatureImage,
                        outputShadow);
                // error!!!
            }

        } else {

            for( int r = 0; r < eibasinID.length; r++ ) {
                for( int c = 0; c < eibasinID[0].length; c++ ) {
                    if (eibasinID[r][c] != (int) NOVALUE)
                        outputShadow[r][c] = 1;
                }
            }
            alpha[0] = 0.0;
            direction[0] = 0.0;

        }

    }

    private void radiation( double[] Rad_morpho, double[] Rad_flat, double E0, double alpha, double direction, double aspect,
            double slope, int shadow ) {

        Rad_flat[0] = E0 * sin(alpha);
        if (shadow == 1 || alpha == 0.0) { // in ombra o di notte
            Rad_morpho[0] = 0.0;
        } else {
            Rad_morpho[0] = E0 * (cos(slope) * sin(alpha) + sin(slope) * cos(alpha) * cos(-aspect + direction));
        }

        if (Rad_morpho[0] < 0)
            Rad_morpho[0] = 0.0;
    }

    private void average_EI( int month_begin, int month_end ) {

        int m, month;

        if (month_end < month_begin)
            month_end += 12;

        for( int r = 0; r < eibasinID.length; r++ ) {
            for( int c = 0; c < eibasinID[0].length; c++ ) {
                if (eibasinID[r][c] != NOVALUE) {
                    eibasinE[r][c] = 0.0;
                    for( m = month_begin - 1; m < month_end; m++ ) {
                        if (m > 11) {
                            month = m - 12;
                        } else {
                            month = m;
                        }
                        if (month < 6) {
                            eibasinE[r][c] += eibasinEmonth[month][r][c];
                        } else {
                            eibasinE[r][c] += eibasinEmonth[12 - month][r][c];
                        }
                    }
                    eibasinE[r][c] /= (double) (month_end - month_begin + 1);
                }
            }
        }

    }

    private void area( int i ) {

        double minES, maxES, minEI, maxEI;

        maxES = Double.NEGATIVE_INFINITY;
        minES = Double.POSITIVE_INFINITY;
        maxEI = Double.NEGATIVE_INFINITY;
        minEI = Double.POSITIVE_INFINITY;

        for( int r = 0; r < eibasinID.length; r++ ) {
            for( int c = 0; c < eibasinID[0].length; c++ ) {
                if (eibasinID[r][c] != NOVALUE) {
                    // System.out.println("Bacino: " + eibasinID[r][c]);
                    if (eibasinID[r][c] == i + 1) {
                        double value = elevImageIterator.getSampleDouble(c, r, 0);
                        if (value < minES)
                            minES = value;
                        if (value > maxES)
                            maxES = value;
                        if (eibasinE[r][c] < minEI)
                            minEI = eibasinE[r][c];
                        if (eibasinE[r][c] > maxEI)
                            maxEI = eibasinE[r][c];
                    }
                }
            }
            // System.out.println("minEi: " + minEI);
            // System.out.println("maxEi: " + maxEI);
        }

        for( int r = 0; r < eibasinID.length; r++ ) {
            for( int c = 0; c < eibasinID[0].length; c++ ) {
                if (eibasinID[r][c] == (i + 1)) {
                    for( int j = 0; j < pEs; j++ ) {
                        double minCurrentAltimetricBand = minES + (j) * (maxES - minES) / (double) pEs;
                        double maxCurrentAltimetricBand = minES + (j + 1) * (maxES - minES) / (double) pEs;
                        double value = elevImageIterator.getSampleDouble(c, r, 0);
                        if ((value > minCurrentAltimetricBand && value <= maxCurrentAltimetricBand) || (j == 0 && value == minES)) {
                            for( int k = 0; k < pEi; k++ ) {
                                double minCurrentEnergeticBand = minEI + (k) * (maxEI - minEI) / (double) pEi;
                                double maxCurrentEnergeticBand = minEI + (k + 1) * (maxEI - minEI) / (double) pEi;
                                if ((eibasinE[r][c] > minCurrentEnergeticBand && eibasinE[r][c] <= maxCurrentEnergeticBand)
                                        || (k == 0 && eibasinE[r][c] == minEI)) {
                                    eibasinA[j][k][i] += 1.0;
                                    break;
                                }
                            }
                            break;
                        }
                    }
                }
            }
        }

        for( int j = 0; j < pEs; j++ ) {
            for( int k = 0; k < pEi; k++ ) {
                eibasinA[j][k][i] *= (dx * dy * 1.0E-6); // in [km2]
            }
        }

        for( int j = 0; j < pEs; j++ ) {
            eibasinESrange[j][i] = (maxES - minES) / (double) pEs;
            eibasinES[j][i] = minES + (j + 1 - 0.5) * (maxES - minES) / (double) pEs;
        }

        int cont = 0;
        for( int m = 0; m < 6; m++ ) {
            for( int k = 0; k < pEi; k++ ) {
                cont = 0;
                for( int r = 0; r < eibasinID.length; r++ ) {
                    for( int c = 0; c < eibasinID[0].length; c++ ) {
                        if (eibasinID[r][c] != NOVALUE) {
                            double test1 = minEI + (k) * (maxEI - minEI) / (double) pEi;
                            double test2 = minEI + (k + 1) * (maxEI - minEI) / (double) pEi;
                            if ((eibasinE[r][c] > test1 && eibasinE[r][c] <= test2) || (k == 0 && eibasinE[r][c] == minEI)) {
                                cont += 1;
                                eibasinEI[m][k][i] += eibasinEmonth[m][r][c];
                            }
                        }
                    }
                }
                if (cont == 0) {
                    eibasinEI[m][k][i] = 0.00001;
                } else {
                    eibasinEI[m][k][i] /= (double) cont;
                }
            }
        }

    }

    // private void output( int i ) {
    //
    // String filename = outputPath + formatter.format(reverseidMappings.get(i + 1)) + ".txt";
    //
    // out.println(MessageFormat.format("Writing output {0} to file {1}", i, filename));
    // if (new File(filename).exists()) {
    // // copy file to backup
    // try {
    // // Create channel on the source
    // FileChannel srcChannel = new FileInputStream(filename).getChannel();
    //
    // // Create channel on the destination
    // FileChannel dstChannel = new FileOutputStream(filename + ".old").getChannel();
    //
    // // Copy file contents from source to destination
    // dstChannel.transferFrom(srcChannel, 0, srcChannel.size());
    //
    // // Close the channels
    // srcChannel.close();
    // dstChannel.close();
    //
    // // remove the file
    // boolean success = (new File(filename)).delete();
    // if (!success) {
    // out.println("Cannot remove file: " + filename + "!");
    // return;
    // }
    //
    // } catch (IOException e) {
    // e.printStackTrace();
    // return;
    // }
    // }
    //
    // BufferedWriter bw;
    // try {
    // bw = new BufferedWriter(new FileWriter(filename, true));
    //
    // bw.write("# generated by EIcalculator*/");
    // bw.write("\n@15");
    // bw.write("\n");
    // bw
    // .write("\n# 1 block - INDICE ENERGETICO PER LE BANDE ENERGETICHE PER OGNI MESE DELL'ANNO (-)")
    // ;
    // bw.write("\n");
    //
    // bw.write("\n# 22 DICEMBRE - 20 GENNAIO ");
    // bw.write("\n% " + numEi + "\n");
    // for( int k = 0; k < numEi - 1; k++ ) {
    // bw.write(eibasinEI[0][k][i] + " ");
    // }
    // bw.write("" + eibasinEI[0][numEi - 1][i]);
    // bw.write("\n");
    //
    // bw.write("\n# 21 GENNAIO - 20 FEBBRAIO");
    // bw.write("\n% " + numEi + "\n");
    // for( int k = 0; k < numEi - 1; k++ ) {
    // bw.write(eibasinEI[1][k][i] + " ");
    // }
    // bw.write("" + eibasinEI[1][numEi - 1][i]);
    // bw.write("\n");
    //
    // bw.write("\n# 21 FEBBRAIO - 22 MARZO");
    // bw.write("\n% " + numEi + "\n");
    // for( int k = 0; k < numEi - 1; k++ ) {
    // bw.write(eibasinEI[2][k][i] + " ");
    // }
    // bw.write("" + eibasinEI[2][numEi - 1][i]);
    // bw.write("\n");
    //
    // bw.write("\n# 23 MARZO - 22 APRILE");
    // bw.write("\n% " + numEi + "\n");
    // for( int k = 0; k < numEi - 1; k++ ) {
    // bw.write(eibasinEI[3][k][i] + " ");
    // }
    // bw.write("" + eibasinEI[3][numEi - 1][i]);
    // bw.write("\n");
    //
    // bw.write("\n# 23 APRILE - 22 MAGGIO");
    // bw.write("\n% " + numEi + "\n");
    // for( int k = 0; k < numEi - 1; k++ ) {
    // bw.write(eibasinEI[4][k][i] + " ");
    // }
    // bw.write("" + eibasinEI[4][numEi - 1][i]);
    // bw.write("\n");
    //
    // bw.write("\n#23 MAGGIO - 22 GIUGNO");
    // bw.write("\n% " + numEi + "\n");
    // for( int k = 0; k < numEi - 1; k++ ) {
    // bw.write(eibasinEI[5][k][i] + " ");
    // }
    // bw.write("" + eibasinEI[5][numEi - 1][i]);
    // bw.write("\n");
    //
    // bw.write("\n");
    // bw.write("\n");
    // bw.write("\n");
    // bw
    // .write(
    // "\n# 2 block - QUOTA DEL BARICENTRO DELLE FASCIE ALTIMETRICHE e RANGE DI QUOTA PER OGNI FASCIA (m)"
    // );
    // bw.write("\n");
    //
    // bw.write("\n% " + numEs + "\n");
    // for( int j = 0; j < numEs - 1; j++ ) {
    // bw.write(eibasinES[j][i] + " ");
    // }
    // bw.write("" + eibasinES[numEs - 1][i]);
    //
    // bw.write("\n% " + numEs + "\n");
    // for( int j = 0; j < numEs - 1; j++ ) {
    // bw.write(eibasinESrange[j][i] + " ");
    // }
    // bw.write("" + eibasinESrange[numEs - 1][i]);
    //
    // bw.write("\n");
    // bw.write("\n");
    // bw.write("\n");
    //
    // bw
    // .write("\n# 3 block - AREE PER FASCIA ALTIMETRICA (riga) E BANDA ENERGETICA (colonna) (km2)");
    // bw.write("\n% " + numEs + " " + numEi + "\n");
    // for( int j = 0; j < numEs; j++ ) {
    // for( int k = 0; k < numEi; k++ ) {
    // bw.write(eibasinA[j][k][i] + " ");
    // }
    // bw.write("\n");
    // }
    //
    // bw.close();
    //
    // } catch (IOException e) {
    // e.printStackTrace();
    // }
    //
    // }
}
