--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoadCalculation.java	Thu Apr 25 15:23:37 2013 +0200
@@ -0,0 +1,343 @@
+package org.dive4elements.river.artifacts.model.minfo;
+
+import gnu.trove.TDoubleArrayList;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.apache.log4j.Logger;
+
+import org.dive4elements.river.artifacts.access.SedimentLoadAccess;
+import org.dive4elements.river.artifacts.model.Calculation;
+import org.dive4elements.river.artifacts.model.CalculationResult;
+
+
+/** Calculate sediment load. */
+public class SedimentLoadCalculation
+extends Calculation
+{
+
+    /** Private logger. */
+    private static final Logger logger = Logger
+        .getLogger(SedimentLoadCalculation.class);
+
+    protected String river;
+    protected String yearEpoch;
+    protected double kmUp;
+    protected double kmLow;
+    protected int[] period;
+    protected int[][] epoch;
+    protected String unit;
+
+    public SedimentLoadCalculation() {
+    }
+
+    public CalculationResult calculate(SedimentLoadAccess access) {
+        logger.info("SedimentLoadCalculation.calculate");
+
+        String river = access.getRiver();
+        String yearEpoch = access.getYearEpoch();
+        String unit = access.getUnit();
+        int[] period = null;
+        int[][] epoch = null;
+        double kmUp = access.getUpperKM();
+        double kmLow = access.getLowerKM();
+        if (yearEpoch.equals("year")) {
+            period = access.getPeriod();
+            epoch = null;
+        }
+        else if (yearEpoch.equals("epoch") || yearEpoch.equals("off_epoch")) {
+            epoch = access.getEpochs();
+            period = null;
+        }
+        else {
+            addProblem("minfo.missing.year_epoch");
+        }
+
+        if (river == null) {
+            // TODO: i18n
+            addProblem("minfo.missing.river");
+        }
+
+        if (period == null && epoch == null) {
+            addProblem("minfo.missing.time");
+        }
+
+        if (!hasProblems()) {
+            this.river = river;
+            this.yearEpoch = yearEpoch;
+            this.unit = unit;
+            this.period = period;
+            this.epoch = epoch;
+            this.kmUp = kmUp;
+            this.kmLow = kmLow;
+            return internalCalculate();
+        }
+
+        return new CalculationResult();
+    }
+
+    private CalculationResult internalCalculate() {
+        logger.debug("internalCalulate; mode:" + yearEpoch);
+        if (yearEpoch.equals("year")) {
+            List<SedimentLoadResult> results =
+                new ArrayList<SedimentLoadResult>();
+            for (int i = period[0]; i <= period[1]; i++) {
+                SedimentLoadResult res = calculateYear(i);
+                results.add(res);
+            }
+            return new CalculationResult(
+                results.toArray(new SedimentLoadResult[results.size()]), this);
+        }
+        else if (yearEpoch.equals("epoch")) {
+            List<SedimentLoadResult> results =
+                new ArrayList<SedimentLoadResult>();
+            for (int i = 0; i < epoch.length; i++) {
+                SedimentLoadResult res = calculateEpoch(i);
+                results.add(res);
+            }
+            return new CalculationResult(
+                results.toArray(new SedimentLoadResult[results.size()]), this);
+        }
+        else if (yearEpoch.equals("off_epoch")) {
+            List<SedimentLoadResult> results =
+                new ArrayList<SedimentLoadResult>();
+            for (int i = 0; i < epoch.length; i++) {
+                SedimentLoadResult res = calculateOffEpoch(i);
+                results.add(res);
+            }
+            return new CalculationResult(
+                results.toArray(new SedimentLoadResult[results.size()]), this);
+        }
+        return null;
+    }
+
+    private SedimentLoadResult calculateEpoch(int i) {
+        List<SedimentLoad> epochLoads = new ArrayList<SedimentLoad>();
+        for (int j = epoch[i][0]; j < epoch[i][1]; j++) {
+            epochLoads.add(SedimentLoadFactory.getLoadWithData(
+                this.river,
+                this.yearEpoch,
+                this.kmLow,
+                this.kmUp,
+                j,
+                j));
+        }
+
+        SedimentLoad resLoad = new SedimentLoad();
+        TDoubleArrayList kms = new TDoubleArrayList();
+
+        for (SedimentLoad load : epochLoads) {
+            for (double km : load.getKms()) {
+                if (!kms.contains(km)) {
+                    kms.add(km);
+                }
+            }
+        }
+
+        for (int j = 0; j < kms.size(); j++) {
+            int cSum = 0;
+            int fmSum = 0;
+            int sSum = 0;
+            int ssSum = 0;
+            int ssbSum = 0;
+            int sseSum = 0;
+            double km = kms.get(j);
+            for (SedimentLoad load : epochLoads) {
+                SedimentLoadFraction f = load.getFraction(km);
+                if (f.getCoarse() > 0d) {
+                    double c = resLoad.getFraction(km).getCoarse();
+                    resLoad.setCoarse(km, c + f.getCoarse());
+                    cSum++;
+                }
+                if (f.getFine_middle() > 0d) {
+                    double fm = resLoad.getFraction(km).getFine_middle();
+                    resLoad.setFineMiddle(km, fm + f.getFine_middle());
+                    fmSum++;
+                }
+                if (f.getSand() > 0d) {
+                    double s = resLoad.getFraction(km).getSand();
+                    resLoad.setSand(km, s + f.getSand());
+                    sSum++;
+                }
+                if (f.getSusp_sand() > 0d) {
+                    double s = resLoad.getFraction(km).getSusp_sand();
+                    resLoad.setSuspSand(km, s + f.getSusp_sand());
+                    ssSum++;
+                }
+                if (f.getSusp_sand_bed() > 0d) {
+                    double s = resLoad.getFraction(km).getSusp_sand_bed();
+                    resLoad.setSuspSandBed(km, s + f.getSusp_sand_bed());
+                    ssbSum++;
+                }
+                if (f.getSusp_sediment() > 0d) {
+                    double s = resLoad.getFraction(km).getSusp_sediment();
+                    resLoad.setSuspSediment(km, s + f.getSusp_sediment());
+                    sseSum++;
+                }
+            }
+            SedimentLoadFraction fr = resLoad.getFraction(km);
+            resLoad.setCoarse(km, fr.getCoarse()/cSum);
+            resLoad.setFineMiddle(km, fr.getFine_middle()/fmSum);
+            resLoad.setSand(km, fr.getSand()/sSum);
+            resLoad.setSuspSand(km, fr.getSusp_sand()/ssSum);
+            resLoad.setSuspSandBed(km, fr.getSusp_sand_bed()/ssbSum);
+            resLoad.setSuspSediment(km, fr.getSusp_sediment()/sseSum);
+        }
+        resLoad.setDescription("");
+        resLoad.setEpoch(true);
+
+        SedimentLoadResult result;
+        SedimentLoad sl = calculateTotalLoad(resLoad, this.epoch[i][0]);
+        if (this.unit.equals("m3_per_a")) {
+            SedimentLoad slu = calculateUnit(sl, this.epoch[i][0]);
+            result = new SedimentLoadResult(
+                this.epoch[i][0],
+                this.epoch[i][1],
+                slu);
+        }
+        else {
+            result = new SedimentLoadResult(
+                this.epoch[i][0],
+                this.epoch[i][1],
+                sl);
+        }
+
+        return result;
+    }
+
+    private SedimentLoadResult calculateOffEpoch(int i) {
+        SedimentLoad load = SedimentLoadFactory.getLoadWithData(
+            this.river,
+            this.yearEpoch,
+            kmLow,
+            kmUp,
+            this.epoch[i][0],
+            this.epoch[i][1]);
+        SedimentLoadResult result;
+        SedimentLoad sl = calculateTotalLoad(load, this.epoch[i][0]);
+        if (unit.equals("m3_per_a")) {
+            SedimentLoad slu = calculateUnit(sl, epoch[i][0]);
+            result = new SedimentLoadResult(
+                this.epoch[i][0],
+                this.epoch[i][1],
+                slu);
+        }
+        else {
+            result = new SedimentLoadResult(
+                this.epoch[i][0],
+                this.epoch[i][1],
+                sl);
+        }
+
+        return result;
+    }
+
+    /** Fetch loads for a single year, calculate total and
+     * return the result containing both. */
+    private SedimentLoadResult calculateYear(int y) {
+        SedimentLoad load = SedimentLoadFactory.getLoadWithData(
+            this.river,
+            this.yearEpoch,
+            this.kmLow,
+            this.kmUp,
+            y,
+            y);
+
+        SedimentLoadResult result;
+        SedimentLoad sl = calculateTotalLoad(load, y);
+        if (unit.equals("m3_per_a")) {
+            SedimentLoad slu = calculateUnit(sl, y);
+            result = new SedimentLoadResult(y, 0, slu);
+        }
+        else {
+            result = new SedimentLoadResult(y, 0, sl);
+        }
+        return result;
+    }
+
+    private SedimentLoad calculateTotalLoad(SedimentLoad load, int year) {
+        logger.debug("calculateTotalLoad");
+        boolean problemThisYear = false;
+        if (!load.hasCoarse()) {
+            addProblem("missing.fraction.coarse", Integer.toString(year));
+            problemThisYear = true;
+        }
+        if (!load.hasFineMiddle()) {
+            addProblem("missing.fraction.fine_middle", Integer.toString(year));
+            problemThisYear = true;
+        }
+        if (!load.hasSand()) {
+            addProblem("missing.fraction.sand", Integer.toString(year));
+            problemThisYear = true;
+        }
+        if (!load.hasSuspSand()) {
+            addProblem("missing.fraction.susp_sand", Integer.toString(year));
+            problemThisYear = true;
+        }
+        if (!load.hasSuspSediment()) {
+            addProblem("missing.fraction.susp_sediment", Integer.toString(year));
+            problemThisYear = true;
+        }
+        if (problemThisYear) {
+            logger.warn("Some problem, not calculating total load.");
+            return load;
+        }
+        for(double km : load.getKms()) {
+            SedimentLoadFraction fraction = load.getFraction(km);
+            double total = 0d;
+            if ((fraction.getCoarse() <= 0d && load.hasCoarse())){
+                addProblem(km, "missing.data.coarse");
+                continue;
+            }
+            if (fraction.getFine_middle() <= 0d && load.hasFineMiddle()) {
+                addProblem(km, "missing.data.fine_middle");
+                continue;
+            }
+            if (fraction.getSand() <= 0d && load.hasSand()) {
+                addProblem(km, "missing data.sand");
+                continue;
+            }
+            if (fraction.getSusp_sand() <= 0d && load.hasSuspSand()) {
+                addProblem(km, "missing.data.susp_sand");
+                continue;
+            }
+            if (fraction.getSusp_sediment() <= 0d && load.hasSuspSediment()) {
+                addProblem(km, "missing.data.susp_sediment");
+                continue;
+            }
+            total += fraction.getCoarse() +
+                fraction.getFine_middle() +
+                fraction.getSand() +
+                fraction.getSusp_sand() +
+                fraction.getSusp_sediment();
+            load.setTotal(km, total);
+        }
+        return load;
+    }
+
+    private SedimentLoad calculateUnit(SedimentLoad load, int year) {
+        SedimentDensity density =
+            SedimentDensityFactory.getSedimentDensity(river, kmLow, kmUp, year);
+        for (double km: load.getKms()) {
+            double dens = density.getDensity(km, year);
+            SedimentLoadFraction fraction = load.getFraction(km);
+            double coarse = fraction.getCoarse();
+            double fineMiddle = fraction.getFine_middle();
+            double sand = fraction.getSand();
+            double suspSand = fraction.getSusp_sand();
+            double bedSand = fraction.getSusp_sand_bed();
+            double sediment = fraction.getSusp_sediment();
+            double total = fraction.getTotal();
+            load.setCoarse(km, (coarse * dens));
+            load.setFineMiddle(km, (fineMiddle * dens));
+            load.setSand(km, (sand * dens));
+            load.setSuspSand(km, (suspSand * dens));
+            load.setSuspSandBed(km, (bedSand * dens));
+            load.setSuspSediment(km, (sediment * dens));
+            load.setTotal(km, (total * dens));
+        }
+        return load;
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :