--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/minfo/BedQualityCalculation.java	Fri Sep 14 14:14:46 2012 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/minfo/BedQualityCalculation.java	Fri Sep 14 14:20:42 2012 +0000
@@ -1,7 +1,10 @@
 package de.intevation.flys.artifacts.model.minfo;
 
+import gnu.trove.TDoubleArrayList;
+
 import java.util.LinkedList;
 import java.util.List;
+import java.util.Map;
 
 import org.apache.log4j.Logger;
 
@@ -20,6 +23,8 @@
     protected String river;
     protected double from;
     protected double to;
+    protected List<String> bedDiameter;
+    protected List<String> bedloadDiameter;
     protected List<DateRange> ranges;
 
     public BedQualityCalculation() {
@@ -31,6 +36,8 @@
         String river = access.getRiver();
         Double from = access.getFrom();
         Double to = access.getTo();
+        List<String> bedDiameter = access.getBedDiameter();
+        List<String> bedloadDiameter = access.getBedloadDiameter();
         List<DateRange> ranges = access.getDateRanges();
 
         if (river == null) {
@@ -58,6 +65,8 @@
             this.from = from;
             this.to = to;
             this.ranges = ranges;
+            this.bedDiameter = bedDiameter;
+            this.bedloadDiameter = bedloadDiameter;
 
             SedDBSessionHolder.acquire();
             try {
@@ -76,28 +85,224 @@
         List<BedQualityResult> results = new LinkedList<BedQualityResult>();
         // Calculate for all time periods.
         for (DateRange dr : ranges) {
-            QualityMeasurements bedMeasurements = QualityMeasurementFactory
-                .getBedMeasurements(river, from, to, dr.getFrom(), dr.getTo());
-            QualityMeasurements loadMeasurements = QualityMeasurementFactory
-                .getBedMeasurements(river, from, to, dr.getFrom(), dr.getTo());
+            QualityMeasurements loadMeasurements =
+                QualityMeasurementFactory.getBedMeasurements(
+                    river,
+                    from,
+                    to,
+                    dr.getFrom(),
+                    dr.getTo());
+            QualityMeasurements bedMeasurements =
+                QualityMeasurementFactory.getBedMeasurements(
+                    river,
+                    from,
+                    to,
+                    dr.getFrom(),
+                    dr.getTo());
+ 
+            for (String bd : bedDiameter) {
+                BedQualityResult bedResult =
+                    calculateBed(bedMeasurements, bd, dr);
 
-            BedQualityResult bedResult = calculateBed(bedMeasurements);
-            BedQualityResult loadResult = calculateBedload(loadMeasurements);
-            results.add(bedResult);
-            results.add(loadResult);
+                //Avoid adding empty result sets.
+                if (!bedResult.isEmpty()) {
+                    results.add(bedResult);
+                }
+            }
+            for (String bld : bedloadDiameter) {
+               BedQualityResult loadResult =
+                    calculateBedload(loadMeasurements, bld, dr);
+                results.add(loadResult);
+            }
         }
 
         return new CalculationResult(
             results.toArray(new BedQualityResult[results.size()]), this);
     }
 
-    protected BedQualityResult calculateBed(QualityMeasurements qm) {
-        // TODO
-        return new BedQualityResult();
+    protected BedQualityResult calculateBed(
+        QualityMeasurements qm,
+        String diameter,
+        DateRange range
+    ) {
+        List<Double> kms = qm.getKms();
+        TDoubleArrayList location = new TDoubleArrayList();
+        TDoubleArrayList avDiameterCap = new TDoubleArrayList();
+        TDoubleArrayList avDiameterSub = new TDoubleArrayList();
+        TDoubleArrayList porosityCap = new TDoubleArrayList();
+        TDoubleArrayList porositySub = new TDoubleArrayList();
+        TDoubleArrayList densityCap = new TDoubleArrayList();
+        TDoubleArrayList densitySub = new TDoubleArrayList();
+        for (double km : kms) {
+            //Filter cap and sub measurements.
+            QualityMeasurements capFiltered = filterCapMeasurements(qm);
+            QualityMeasurements subFiltered = filterSubMeasurements(qm);
+
+            List<QualityMeasurement> cm = capFiltered.getMeasurements(km);
+            List<QualityMeasurement> sm = subFiltered.getMeasurements(km);
+
+            double avCap = calculateAverage(cm, diameter);
+            double avSub = calculateAverage(sm, diameter);
+            double[] pCap = calculatePorosity(capFiltered, km, diameter);
+            double[] pSub = calculatePorosity(subFiltered, km, diameter);
+            double[] dCap = calculateDensity(capFiltered, pCap);
+            double[] dSub = calculateDensity(subFiltered, pSub);
+
+            double pCapRes = 0d;
+            double pSubRes = 0d;
+            double dCapRes = 0d;
+            double dSubRes = 0d;
+            for (int i = 0; i < pCap.length; i++) {
+                pCapRes += pCap[i];
+                dCapRes += dCap[i];
+            }
+            for (int i = 0; i < pSub.length; i++) {
+                pSubRes += pSub[i];
+                dSubRes += dSub[i];
+            }
+            location.add(km);
+            avDiameterCap.add(avCap);
+            avDiameterSub.add(avSub);
+            porosityCap.add((pCapRes / pCap.length) * 100 );
+            porositySub.add((pSubRes / pSub.length) * 100);
+            densityCap.add((dCapRes / dCap.length) / 1000);
+            densitySub.add((dSubRes / dSub.length) / 1000);
+        }
+        return new BedBedQualityResult(
+            diameter,
+            avDiameterCap,
+            avDiameterSub,
+            location,
+            range,
+            porosityCap,
+            porositySub,
+            densityCap,
+            densitySub);
     }
 
-    protected BedQualityResult calculateBedload(QualityMeasurements qm) {
-        // TODO
-        return new BedQualityResult();
+    private double[] calculateDensity(
+        QualityMeasurements capFiltered,
+        double[] porosity
+    ) {
+        double[] density = new double[porosity.length];
+        for (int i = 0; i < porosity.length; i++) {
+            density[i] = (1 - porosity[i]) * 2650;
+        }
+        return density;
+    }
+
+    private double[] calculatePorosity(
+        QualityMeasurements capFiltered,
+        double km,
+        String diameter
+    ) {
+        List<QualityMeasurement> list = capFiltered.getMeasurements(km);
+        double[] results = new double[list.size()];
+        int i = 0;
+        for (QualityMeasurement qm : list) {
+            double deviation = calculateDeviation(qm);
+            double p = calculateP(qm);
+            double porosity = 0.353 - 0.068 * deviation + 0.146 * p;
+            results[i] = porosity;
+            i++;
+        }
+
+        return results;
+    }
+
+    protected BedQualityResult calculateBedload(
+        QualityMeasurements qm,
+        String diameter,
+        DateRange range
+    ) {
+        List<Double> kms = qm.getKms();
+        TDoubleArrayList location = new TDoubleArrayList();
+        TDoubleArrayList avDiameter = new TDoubleArrayList();
+        for (double km : kms) {
+            List<QualityMeasurement> measurements = qm.getMeasurements(km);
+            double mid = calculateAverage(measurements, diameter);
+            location.add(km);
+            avDiameter.add(mid);
+        }
+        return new BedLoadBedQualityResult(
+            diameter,
+            avDiameter,
+            location,
+            range);
+    }
+
+    protected double calculateAverage(
+        List<QualityMeasurement> list,
+        String diameter
+    ) {
+        double av = 0;
+        for (QualityMeasurement qm : list) {
+            av += qm.getDiameter(diameter);
+        }
+        return av/list.size();
+    }
+
+    protected QualityMeasurements filterCapMeasurements(
+        QualityMeasurements qms
+    ) {
+        List<QualityMeasurement> result = new LinkedList<QualityMeasurement>();
+        for (QualityMeasurement qm : qms.getMeasurements()) {
+            if (qm.getDepth1() == 0d && qm.getDepth2() <= 0.3) {
+                result.add(qm);
+            }
+        }
+        return new QualityMeasurements(result);
+    }
+
+    protected QualityMeasurements filterSubMeasurements(
+        QualityMeasurements qms
+    ) {
+        List<QualityMeasurement> result = new LinkedList<QualityMeasurement>();
+        for (QualityMeasurement qm : qms.getMeasurements()) {
+            if (qm.getDepth1() > 0d && qm.getDepth2() <= 0.5) {
+                result.add(qm);
+            }
+        }
+        return new QualityMeasurements(result);
+    }
+
+    public double calculateDeviation(QualityMeasurement qm) {
+        Map<String, Double> dm = qm.getAllDiameter();
+        double phiM = 0;
+        double[] phis = new double[dm.size()];
+        double[] ps = new double[dm.size()];
+        int i = 0;
+        for (String key : dm.keySet()) {
+            double d = dm.get(key);
+            double phi = -Math.log(d)/Math.log(2);
+            phis[i] = phi;
+            double p = calculateWeight(qm, key);
+            ps[i] = p;
+            phiM += phi * p;
+            i++;
+        }
+
+        double sig = 0d;
+        for (i = 0; i < dm.size(); i++) {
+            sig += ps[i] * Math.exp(phis[i] - phiM);
+        }
+        double deviation = Math.sqrt(sig);
+        return deviation;
+    }
+
+    protected double calculateP(QualityMeasurement qm) {
+        return calculateWeight(qm, "dmin");
+    }
+
+    public double calculateWeight(QualityMeasurement qm, String diameter) {
+        Map<String, Double> dm = qm.getAllDiameter();
+        double value = qm.getDiameter(diameter);
+
+        double sum = 0d;
+        for (Double d : dm.values()) {
+            sum =+ d.doubleValue();
+        }
+        double weight = sum/100*value;
+        return weight;
     }
 }